Description of supdeq

0001 %% SUpDEq - Spatial Upsampling by Directional Equalization
0002 %
0003 % function [denseHRTFdataset, denseHRIRdataset, denseHRTFdataset_sfd] = supdeq_deq(eqHRTFdataset, deqDataset, N, samplingGrid, windowLength, normalization, parallax)
0004 %
0005 % This function performs the de-equalization of a equalized sparse HRTF
0006 % dataset, here called eqHRTFdataset, (in SH-domain / stored as SH-coefficients)
0007 % with the pre-defined deqDataset (dataset for de-equalization in SH-domain /
0008 % stored as SH-coefficients). The result is a dense hrtf dataset (according
0009 % to the (dense) spatial sampling grid), saved as in structs as HRTFs
0010 % (denseHRTFdataset), HRIRs (denseHRIRdataset) or as SH-coefficients
0011 % (denseHRTFdataset_sfd)
0012 %
0013 % Output:
0014 % denseHRTFdataset      - Dense HRTF dataset as a struct
0015 % denseHRIRdataset      - Dense HRIR dataset as a struct
0016 % denseHRTFdataset_sfd  - Dense HRIR dataset as SH-coefficients,
0017 %                         transformed with given order N
0018 %
0019 % Input:
0020 % eqHRTFdataset         - Struct with (SH-Coefficients) of the equalized
0021 %                         HRTF dataset for the left (Hl_nm) and right (Hr_nm)
0022 %                         channel/ear, absolute frequency scale f,
0023 %                         transform order N, and FFToversize
0024 % deqDataset            - Struct with de-equalization dataset in SH-domain / as
0025 %                         SH-coefficients. Can be the output of supdeq_getEqDataset.
0026 % N                     - Transform order N
0027 % samplingGrid          - Spatial sampling grid (Q x 2 or Q x 3 matrix),
0028 %                         where the first column holds the azimuth,
0029 %                         the second the elevation (both in degree), and
0030 %                         optionally the third the sampling weights.
0031 %                         Azimuth in degree
0032 %                         (0=front, 90=left, 180=back, 270=right)
0033 %                         (0 points to positive x-axis, 90 to positive y-axis)
0034 %                         Elevations in degree
0035 %                         (0=North Pole, 90=front, 180=South Pole)
0036 %                         (0 points to positive z-axis, 180 to negative z-axis)
0037 % windowLength          - Values of head and tail window in samples [a,b].
0038 %                         Head and tail window will be a half sided hann
0039 %                         window with the defined length. Windowing will be
0040 %                         applied to the HRIRs. The HRTF dataset in frequency-
0041 %                         and SH-domain will be calculated based on these
0042 %                         windowed HRIRs. If the matrix is empty ([]), no
0043 %                         windowing will be applied!
0044 %                         Default: [], no windowing
0045 % normalization         - Linear value (for example 0.99) which can be
0046 %                         applied in order to normalize the de-equalized
0047 %                         HRIRs. Can be useful if de-equalization with
0048 %                         waveType == 1 (spherical wave) is applied.
0049 %                         Default: [], no normalization.
0050 % parallax              - Boolean with true or false. Only applies if
0051 %                         waveType in deq is 1 (spherical wave) and a
0052 %                         distance shift is applied (DVF).
0053 %                         false: acoustic parallax not considered
0054 %                         true (Default: acoustic parallax considered
0055 %
0056 % Dependencies: SOFiA toolbox
0057 %
0058 % (C) 2018 by JMA, Johannes M. Arend
0059 %             CP,  Christoph Pörschmann
0060 %             TH Köln - University of Applied Sciences
0061 %             Institute of Communications Engineering
0062 %             Department of Acoustics and Audio Signal Processing
0063 
0064 function [denseHRTFdataset, denseHRIRdataset, denseHRTFdataset_sfd] = supdeq_deq(eqHRTFdataset, deqDataset, N, samplingGrid, windowLength, normalization, parallax)
0065 
0066 if nargin < 5 || isempty(windowLength)
0067     windowLength = [];
0068 end
0069 
0070 if nargin <6 || isempty(normalization)
0071     normalization = [];
0072 end
0073 
0074 if nargin < 7 || isempty(parallax)
0075     parallax = true;
0076 end
0077 
0078 %Check if eqDataset was with limiting but deqDataset is without limiting
0079 if isfield(eqHRTFdataset,'limited')
0080     if eqHRTFdataset.limited
0081         if isfield(deqDataset,'limited')
0082             if ~deqDataset.limited
0083                 error('Limited eq transfer functions applied for equalization, but not applied for de-equalization!');
0084             end
0085         else
0086             error('Limited eq transfer functions applied for equalization, but not applied for de-equalization!');
0087         end
0088     end
0089 end
0090 
0091 %Check if applied radius and fA is the same
0092 if isfield(eqHRTFdataset,'limited') && isfield(deqDataset,'limited')
0093     if eqHRTFdataset.limited && deqDataset.limited
0094         if eqHRTFdataset.limitInfo.appliedRadius ~= deqDataset.limitInfo.appliedRadius
0095             error('Applied radius of limited equalization and de-equalization dataset must be the same!');
0096         end
0097         if eqHRTFdataset.limitInfo.fA ~= deqDataset.limitInfo.fA
0098             error('fA of limited equalization and de-equalization dataset must be the same!');
0099         end
0100     end
0101 end
0102 
0103 %Check if phaseOnly equalization was applied
0104 if isfield(eqHRTFdataset,'phaseOnly')
0105     phaseOnly = 1;
0106 else
0107     phaseOnly = 0;
0108 end
0109 
0110 %Get fs
0111 fs = eqHRTFdataset.f(end)*2;
0112 
0113 %% PROCESSING
0114 %Processing for de-equalization with plane wave deqDataset
0115 if deqDataset.waveType == 0
0116 
0117     %Transform deqDataset to Fourier domain at dense sampling grid points
0118     %(inverse spherical Fourier transform)
0119     fprintf('Extracting %d deq transfer functions. This may take some time...\n',size(samplingGrid,1));
0120     [deqTF_L,deqTF_R] = supdeq_getArbHRTF(deqDataset,samplingGrid,[],[],'ak'); %Use AKtools...
0121     fprintf('%d deq transfer functions extracted...\n',size(samplingGrid,1))     
0122     
0123     if phaseOnly    
0124         fprintf('Phase only de-equalization...\n')
0125         %Get only phase response of eqTF
0126         deqTF_L_phase = angle(deqTF_L);
0127         deqTF_R_phase = angle(deqTF_R);
0128         deqTF_L_mag = ones(size(deqTF_L,1),size(deqTF_L,2));
0129         deqTF_R_mag = ones(size(deqTF_R,1),size(deqTF_R,2));
0130         deqTF_L = deqTF_L_mag.*exp(1i*deqTF_L_phase);
0131         deqTF_R = deqTF_R_mag.*exp(1i*deqTF_R_phase);
0132     end
0133     
0134     %Get HRTF from equalized HRTF dataset by inverse spherical Fourier transform at
0135     %sampling points of dense sampling grid
0136     fprintf('Extracting %d HRTFs. This may take some time...\n',size(samplingGrid,1));
0137     [eqHRTF_L,eqHRTF_R] = supdeq_getArbHRTF(eqHRTFdataset,samplingGrid,[],[],'ak'); %Use AKtools...
0138     fprintf('%d HRTFs extracted...\n',size(samplingGrid,1))
0139 
0140     %Perform de-equalization (multiplication in frequency domain)
0141     HRTF_deequalized_L = zeros(size(samplingGrid,1),length(deqDataset.f));
0142     HRTF_deequalized_R = zeros(size(samplingGrid,1),length(deqDataset.f));
0143     for kk = 1:size(samplingGrid,1)
0144         HRTF_deequalized_L(kk,:)=eqHRTF_L(kk,:).*deqTF_L(kk,:);
0145         HRTF_deequalized_R(kk,:)=eqHRTF_R(kk,:).*deqTF_R(kk,:); 
0146 
0147         if ~mod(kk,100)
0148             fprintf('%d of %d HRTFs de-equalized...\n',kk,size(samplingGrid,1))
0149         end
0150 
0151     end
0152 
0153 end
0154 
0155 
0156 %Processing for de-equalization with spherical wave deqDataset (DVF - Distance Variation Function)
0157 if deqDataset.waveType == 1
0158     
0159     fprintf('Applying DVF...\n');
0160 
0161     %Transform deqDataset to Fourier domain at dense sampling grid points
0162     %(inverse spherical Fourier transform)
0163     fprintf('Extracting %d deq transfer functions. This may take some time...\n',size(samplingGrid,1));
0164     [deqTF_L,deqTF_R] = supdeq_getArbHRTF(deqDataset,samplingGrid,[],[],'ak'); %Use AKtools...
0165     fprintf('%d deq transfer functions extracted...\n',size(samplingGrid,1))
0166     
0167     if parallax %Consider acoustic parallax
0168         
0169         fprintf('Considering acoustic parallax...\n');
0170     
0171         %Get parallax shifted dense sampling grid for left and right ear
0172         [samplingGridParL,samplingGridParR] = supdeq_parallax(samplingGrid,deqDataset.sourceDistance,deqDataset.radius,deqDataset.earPosition);
0173 
0174         %Get HRTF from equalized HRTF dataset by inverse spherical Fourier transform at
0175         %sampling points of dense sampling grid
0176         fprintf('Extracting %d HRTFs. This may take some time...\n',size(samplingGridParL,1));
0177         %Get left and right HRTFs separately
0178         [eqHRTF_L,~] = supdeq_getArbHRTF(eqHRTFdataset,samplingGridParL,'DEG',0,'ak'); %Use AKtools...
0179         [~,eqHRTF_R] = supdeq_getArbHRTF(eqHRTFdataset,samplingGridParR,'DEG',1,'ak'); %Use AKtools...
0180         fprintf('%d HRTFs extracted...\n',size(samplingGrid,1))
0181     
0182     else %parallax = false, do not consider acoustic parallax
0183         
0184         fprintf('Not considering acoustic parallax...\n');
0185         
0186         %Get HRTF from equalized HRTF dataset by inverse spherical Fourier transform at
0187         %sampling points of dense sampling grid
0188         fprintf('Extracting %d HRTFs. This may take some time...\n',size(samplingGrid,1));
0189         [eqHRTF_L,eqHRTF_R] = supdeq_getArbHRTF(eqHRTFdataset,samplingGrid,[],[],'ak'); %Use AKtools...
0190         fprintf('%d HRTFs extracted...\n',size(samplingGrid,1))
0191         
0192     end
0193 
0194     %Perform de-equalization (multiplication in frequency domain)
0195     HRTF_deequalized_L = zeros(size(samplingGrid,1),length(deqDataset.f));
0196     HRTF_deequalized_R = zeros(size(samplingGrid,1),length(deqDataset.f));
0197     for kk = 1:size(samplingGrid,1)
0198         HRTF_deequalized_L(kk,:)=eqHRTF_L(kk,:).*deqTF_L(kk,:);
0199         HRTF_deequalized_R(kk,:)=eqHRTF_R(kk,:).*deqTF_R(kk,:); 
0200 
0201         if ~mod(kk,100)
0202             fprintf('%d of %d HRTFs de-equalized...\n',kk,size(samplingGrid,1))
0203         end
0204 
0205     end
0206 
0207 end
0208 
0209 %% Store in structs
0210 
0211 % Transform HRTFs to time domain --> get HRIRs
0212 %Get both sided spectrum
0213 HRTF_deequalized_L_double = conj([HRTF_deequalized_L(:,:), conj(fliplr(HRTF_deequalized_L(:,2:end-1)))])';
0214 HRTF_deequalized_R_double = conj([HRTF_deequalized_R(:,:), conj(fliplr(HRTF_deequalized_R(:,2:end-1)))])';
0215 %Transform back to time domain
0216 HRIR_deequalized_L = real(ifft(HRTF_deequalized_L_double));
0217 HRIR_deequalized_R = real(ifft(HRTF_deequalized_R_double));
0218 %Cut if FFToversize of eqHRTFdataset > 1
0219 if eqHRTFdataset.FFToversize > 1
0220     newLength = size(HRIR_deequalized_L,1)/eqHRTFdataset.FFToversize;
0221     HRIR_deequalized_L = HRIR_deequalized_L(1:newLength,:);
0222     HRIR_deequalized_R = HRIR_deequalized_R(1:newLength,:);
0223 end
0224 
0225 %Window if windowLength is not empty
0226 if ~isempty(windowLength) %windowing
0227     disp('Applying windows...');
0228     
0229     %Apply normalization if desired
0230     if ~isempty(normalization)
0231         disp('Applying normalization...');
0232         
0233         %Get abs max of left/right channel
0234         maxL = max(max(abs(HRIR_deequalized_L)));
0235         maxR = max(max(abs(HRIR_deequalized_R)));
0236         
0237         %Calculate normalization factor
0238         if maxL >= maxR
0239             normFactor = normalization/maxL;
0240         elseif maxR > maxL
0241             normFactor = normalization/maxR;
0242         end
0243         
0244         %Apply to HRIRs
0245         HRIR_deequalized_L = HRIR_deequalized_L*normFactor;
0246         HRIR_deequalized_R = HRIR_deequalized_R*normFactor; 
0247     end
0248     
0249     %Windowing
0250     %Get window function in time domain
0251     win = supdeq_win(size(HRIR_deequalized_L,1),windowLength);
0252     %Apply to HRIRs
0253     for kk = 1:size(HRIR_deequalized_L,2)
0254         HRIR_deequalized_L(:,kk) = HRIR_deequalized_L(:,kk) .* win;
0255         HRIR_deequalized_R(:,kk) = HRIR_deequalized_R(:,kk) .* win;
0256     end
0257     
0258     %Store in struct denseHRIRdataset (window, normalization optionally)
0259     denseHRIRdataset.HRIR_L = HRIR_deequalized_L'; %Flip again to have same array settings as in SOFiA
0260     denseHRIRdataset.HRIR_R = HRIR_deequalized_R';
0261     denseHRIRdataset.fs = eqHRTFdataset.f(end)*2;
0262     denseHRIRdataset.samplingGrid = samplingGrid;
0263     denseHRIRdataset.deqEarDistance = deqDataset.earDistance;
0264     denseHRIRdataset.deqWaveType  = deqDataset.waveType;
0265     if deqDataset.waveType == 1
0266         denseHRIRdataset.sourceDistance = deqDataset.sourceDistance;
0267     end
0268     if ~isempty(normalization)
0269         denseHRIRdataset.normFactor = normFactor;
0270     end
0271     denseHRIRdataset.windowLength = windowLength;
0272     
0273     %Transform de-equalized windowed denseHRIRdataset to Fourier domain to get
0274     %denseHRTFdataset (with same FFToversize)
0275     HRTF_deequalized_win_L = fft(HRIR_deequalized_L',size(HRIR_deequalized_L,1)*eqHRTFdataset.FFToversize,2);
0276     HRTF_deequalized_win_L = HRTF_deequalized_win_L(:,1:end/2+1);
0277     HRTF_deequalized_win_R = fft(HRIR_deequalized_R',size(HRIR_deequalized_R,1)*eqHRTFdataset.FFToversize,2);
0278     HRTF_deequalized_win_R = HRTF_deequalized_win_R(:,1:end/2+1);
0279     
0280     %Store in struct denseHRTFdataset (window, normalization optionally)
0281     denseHRTFdataset.HRTF_L = HRTF_deequalized_win_L;
0282     denseHRTFdataset.HRTF_R = HRTF_deequalized_win_R;
0283     denseHRTFdataset.f      = eqHRTFdataset.f;
0284     denseHRTFdataset.fs     = eqHRTFdataset.f(end)*2;
0285     denseHRTFdataset.FFToversize = eqHRTFdataset.FFToversize;
0286     denseHRTFdataset.samplingGrid = samplingGrid;
0287     denseHRTFdataset.deqEarDistance = deqDataset.earDistance;
0288     denseHRTFdataset.deqWaveType  = deqDataset.waveType;
0289     if deqDataset.waveType == 1
0290         denseHRTFdataset.sourceDistance = deqDataset.sourceDistance;
0291     end
0292     if ~isempty(normalization)
0293         denseHRTFdataset.normFactor = normFactor;
0294     end
0295     denseHRTFdataset.windowLength = windowLength;
0296     
0297     % Transform de-equalized HRTFs to SH-domain --> get
0298     % denseHRTFdataset_sfd (window, normalization optionally)
0299     fprintf('Performing spherical Fourier transform with N = %d. This may take some time...\n',N);
0300     denseHRTFdataset_sfd = supdeq_hrtf2sfd(HRTF_deequalized_win_L,HRTF_deequalized_win_R,N,samplingGrid,fs,'ak'); %Use AKtools...
0301     denseHRTFdataset_sfd.FFToversize = eqHRTFdataset.FFToversize;
0302     denseHRTFdataset_sfd.deqEarDistance = deqDataset.earDistance;
0303     denseHRTFdataset_sfd.deqWaveType  = deqDataset.waveType;
0304     if deqDataset.waveType == 1
0305         denseHRTFdataset_sfd.sourceDistance = deqDataset.sourceDistance;
0306     end
0307     if ~isempty(normalization)
0308         denseHRTFdataset_sfd.normFactor = normFactor;
0309     end
0310     denseHRTFdataset_sfd.windowLength = windowLength;
0311     
0312     fprintf('Done with de-equalization...\n');
0313       
0314 else %no windowing
0315  
0316 %Apply normalization if desired
0317     if ~isempty(normalization)
0318         disp('Applying normalization...');
0319         
0320         %Get abs max of left/right channel
0321         maxL = max(max(abs(HRIR_deequalized_L)));
0322         maxR = max(max(abs(HRIR_deequalized_R)));
0323         
0324         %Calculate normalization factor
0325         if maxL >= maxR
0326             normFactor = normalization/maxL;
0327         elseif maxR > maxL
0328             normFactor = normalization/maxR;
0329         end
0330         
0331         %Apply to HRIRs
0332         HRIR_deequalized_L = HRIR_deequalized_L*normFactor;
0333         HRIR_deequalized_R = HRIR_deequalized_R*normFactor;   
0334         
0335         %Store in struct denseHRIRdataset (no window, normalization)
0336         denseHRIRdataset.HRIR_L = HRIR_deequalized_L'; %Flip again to have same array settings as in SOFiA
0337         denseHRIRdataset.HRIR_R = HRIR_deequalized_R';
0338         denseHRIRdataset.fs = eqHRTFdataset.f(end)*2;
0339         denseHRIRdataset.samplingGrid = samplingGrid;
0340         denseHRIRdataset.deqEarDistance = deqDataset.earDistance;
0341         denseHRIRdataset.deqWaveType  = deqDataset.waveType;
0342         if deqDataset.waveType == 1
0343             denseHRIRdataset.sourceDistance = deqDataset.sourceDistance;
0344         end
0345         denseHRIRdataset.normFactor = normFactor;
0346         
0347         %Transform de-equalized normalized denseHRIRdataset to Fourier domain to get
0348         %denseHRTFdataset (with same FFToversize)
0349         HRTF_deequalized_L = fft(HRIR_deequalized_L',size(HRIR_deequalized_L,1)*eqHRTFdataset.FFToversize,2);
0350         HRTF_deequalized_L = HRTF_deequalized_L(:,1:end/2+1);
0351         HRTF_deequalized_R = fft(HRIR_deequalized_R',size(HRIR_deequalized_R,1)*eqHRTFdataset.FFToversize,2);
0352         HRTF_deequalized_R = HRTF_deequalized_R(:,1:end/2+1);
0353         
0354         %Store in struct denseHRTFdataset (no window, normalization)
0355         denseHRTFdataset.HRTF_L = HRTF_deequalized_L;
0356         denseHRTFdataset.HRTF_R = HRTF_deequalized_R;
0357         denseHRTFdataset.f      = eqHRTFdataset.f;
0358         denseHRTFdataset.fs     = eqHRTFdataset.f(end)*2;
0359         denseHRTFdataset.FFToversize = eqHRTFdataset.FFToversize;
0360         denseHRTFdataset.samplingGrid = samplingGrid;
0361         denseHRTFdataset.deqEarDistance = deqDataset.earDistance;
0362         denseHRTFdataset.deqWaveType  = deqDataset.waveType;
0363         if deqDataset.waveType == 1
0364             denseHRTFdataset.sourceDistance = deqDataset.sourceDistance;
0365         end
0366         denseHRTFdataset.normFactor = normFactor;
0367         
0368         % Transform de-equalized HRTFs to SH-domain --> get
0369         % denseHRTFdataset_sfd (no window, normalization)
0370         fprintf('Performing spherical Fourier transform with N = %d. This may take some time...\n',N);
0371         denseHRTFdataset_sfd = supdeq_hrtf2sfd(HRTF_deequalized_L,HRTF_deequalized_R,N,samplingGrid,fs,'ak'); %Use AKtools...
0372         denseHRTFdataset_sfd.FFToversize = eqHRTFdataset.FFToversize;
0373         denseHRTFdataset_sfd.deqEarDistance = deqDataset.earDistance;
0374         denseHRTFdataset_sfd.deqWaveType  = deqDataset.waveType;
0375         if deqDataset.waveType == 1
0376             denseHRTFdataset_sfd.sourceDistance = deqDataset.sourceDistance;
0377         end
0378         denseHRTFdataset_sfd.normFactor = normFactor;
0379         fprintf('Done with de-equalization...\n');
0380     
0381         
0382     else % No normalization
0383         %Store in struct denseHRIRdataset (no window, no normalization)
0384         denseHRIRdataset.HRIR_L = HRIR_deequalized_L'; %Flip again to have same array settings as in SOFiA
0385         denseHRIRdataset.HRIR_R = HRIR_deequalized_R';
0386         denseHRIRdataset.fs = eqHRTFdataset.f(end)*2;
0387         denseHRIRdataset.samplingGrid = samplingGrid;
0388         denseHRIRdataset.deqEarDistance = deqDataset.earDistance;
0389         denseHRIRdataset.deqWaveType  = deqDataset.waveType;
0390         if deqDataset.waveType == 1
0391             denseHRIRdataset.sourceDistance = deqDataset.sourceDistance;
0392         end
0393 
0394         %Store in struct denseHRTFdataset (no window, no normalization)
0395         denseHRTFdataset.HRTF_L = HRTF_deequalized_L;
0396         denseHRTFdataset.HRTF_R = HRTF_deequalized_R;
0397         denseHRTFdataset.f      = eqHRTFdataset.f;
0398         denseHRTFdataset.fs     = eqHRTFdataset.f(end)*2;
0399         denseHRTFdataset.FFToversize = eqHRTFdataset.FFToversize;
0400         denseHRTFdataset.samplingGrid = samplingGrid;
0401         denseHRTFdataset.deqEarDistance = deqDataset.earDistance;
0402         denseHRTFdataset.deqWaveType  = deqDataset.waveType;
0403         if deqDataset.waveType == 1
0404             denseHRTFdataset.sourceDistance = deqDataset.sourceDistance;
0405         end
0406 
0407         % Transform de-equalized HRTFs to SH-domain --> get
0408         % denseHRTFdataset_sfd (no window, no normalization)
0409         fprintf('Performing spherical Fourier transform with N = %d. This may take some time...\n',N);
0410         denseHRTFdataset_sfd = supdeq_hrtf2sfd(HRTF_deequalized_L,HRTF_deequalized_R,N,samplingGrid,fs,'ak'); %Use AKtools...
0411         denseHRTFdataset_sfd.FFToversize = eqHRTFdataset.FFToversize;
0412         denseHRTFdataset_sfd.deqEarDistance = deqDataset.earDistance;
0413         denseHRTFdataset_sfd.deqWaveType  = deqDataset.waveType;
0414         if deqDataset.waveType == 1
0415             denseHRTFdataset_sfd.sourceDistance = deqDataset.sourceDistance;
0416         end
0417         fprintf('Done with de-equalization...\n');
0418 
0419     end
0420 end
0421 
0422 end
0423
supdeq_deq

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
