supdeq_alfe

% SUpDEq - Spatial Upsampling by Directional Equalization

 function [eqHRTFDataset_alfe] = supdeq_alfe(eqHRTFDataset, samplingGrid, xover, fs)

 This function applies the ALFE (Adaptive Low Frequency Extention)
 algorithm in frequency domain (ALFE_fd) to HRTFs. Typically, it is 
 applied on the equalized HRTF dataset as a processing step of SUpDEq.

 Input:
 eqHRTFDataset         - Equalized HRTF dataset in the SH domain 
 xover                 - Crossover frequency in Hz
 fs                    - Sampling rate
                         Default: 48000
 Output:       
 eqHRTFDataset_ALFE    - Equalized and ALFE-processeed HRTF dataset in the SH domain 
 
 Dependencies: ALFE_fd

 References:
 Xie, Bosun: 'On the low frequency characteristics of head-related
 transfer function'. In: Chinese J. Acoust. 28(2):1-13 (2009).

 (C) 2019 by CP,  Christoph Pörschmann
             JMA, Johannes M. Arend
             Technische Hochschule Köln
             University of Applied Sciences
             Institute of Communications Engineering
             Department of Acoustics and Audio Signal Processing

0001 %% SUpDEq - Spatial Upsampling by Directional Equalization
0002 %
0003 % function [eqHRTFDataset_alfe] = supdeq_alfe(eqHRTFDataset, samplingGrid, xover, fs)
0004 %
0005 % This function applies the ALFE (Adaptive Low Frequency Extention)
0006 % algorithm in frequency domain (ALFE_fd) to HRTFs. Typically, it is
0007 % applied on the equalized HRTF dataset as a processing step of SUpDEq.
0008 %
0009 % Input:
0010 % eqHRTFDataset         - Equalized HRTF dataset in the SH domain
0011 % xover                 - Crossover frequency in Hz
0012 % fs                    - Sampling rate
0013 %                         Default: 48000
0014 % Output:
0015 % eqHRTFDataset_ALFE    - Equalized and ALFE-processeed HRTF dataset in the SH domain
0016 %
0017 % Dependencies: ALFE_fd
0018 %
0019 % References:
0020 % Xie, Bosun: 'On the low frequency characteristics of head-related
0021 % transfer function'. In: Chinese J. Acoust. 28(2):1-13 (2009).
0022 %
0023 % (C) 2019 by CP,  Christoph Pörschmann
0024 %             JMA, Johannes M. Arend
0025 %             Technische Hochschule Köln
0026 %             University of Applied Sciences
0027 %             Institute of Communications Engineering
0028 %             Department of Acoustics and Audio Signal Processing
0029 
0030 function [eqHRTFDataset_alfe] = supdeq_alfe(eqHRTFDataset, samplingGrid, xover, fs)
0031 
0032 if nargin < 4 || isempty(fs)
0033     fs = 48000;
0034 end
0035 
0036 if length(xover)~=2
0037     warning('ALFE Crossover frequency wrong defined');
0038 end
0039 
0040 if or(min(xover)<300,max(xover)>1300) 
0041     warning('ALFE Crossover frequencies in not in typical range between 300 and 1300 Hz');
0042 end
0043 
0044 %Get equalized HRTFs from SH coefficients of equalized HRTF dataset
0045 %[HRTF_equalized_L, HRTF_equalized_R] = supdeq_getArbHRTF(eqHRTFDataset,samplingGrid,'DEG',2,'ak');
0046 %Get equalized HRTFs
0047 HRTF_equalized_L = eqHRTFDataset.HRTF_L;
0048 HRTF_equalized_R = eqHRTFDataset.HRTF_R;
0049 %Get mirror spectrum
0050 HlTwoSided = AKsingle2bothSidedSpectrum(HRTF_equalized_L.');
0051 HrTwoSided = AKsingle2bothSidedSpectrum(HRTF_equalized_R.');
0052 %Transform back to time domain
0053 hl = real(ifft(HlTwoSided));
0054 hr = real(ifft(HrTwoSided));
0055 
0056 %% Determine average value in crossover range
0057 
0058 NFFT = size(hl,1);
0059 n_link = round(xover / (fs/NFFT)) + 1;
0060 ALFE_abs = mean(mean(abs(HlTwoSided(n_link(1):n_link(2),:)))+mean(abs(HrTwoSided(n_link(1):n_link(2),:))))/2;
0061 
0062 %% Apply ALFE to all HRIRs
0063 
0064 fprintf('Applying ALFE\n');
0065 
0066 for kk = 1:length(hl(1,:))
0067     
0068     %Get ALFE processed HRIR
0069     hl_ALFE  = ALFE_fd(hl(:,kk),'f_link', xover,'fs',fs,'L_target',20*log10(ALFE_abs));
0070     hr_ALFE  = ALFE_fd(hr(:,kk),'f_link', xover,'fs',fs,'L_target',20*log10(ALFE_abs));
0071 
0072     %Interchange original HRIRs with ALFE processed HRIRs
0073     hl(:,kk)=hl_ALFE;
0074     hr(:,kk)=hr_ALFE;
0075     
0076     %Transform back to Fourier Domain
0077     Hl=fft(hl(:,kk));
0078     Hr=fft(hr(:,kk));
0079     
0080     %Get rid of mirror spectrum
0081     HRTF_equalized_ALFE_L(kk,:)=Hl(1:length(HRTF_equalized_L(1,:)));
0082     HRTF_equalized_ALFE_R(kk,:)=Hr(1:length(HRTF_equalized_R(1,:)));
0083     
0084     %Status prints
0085     if ~mod(kk,10)
0086         fprintf('|');
0087     end
0088     if ~mod(kk,500)
0089         fprintf(' %d\n',kk); 
0090     end
0091 end
0092 
0093 %Transform ALFE processed HRTFs to SH domain
0094 fprintf('\nPerforming spherical Fourier transform with N = %d. This may take some time...\n',eqHRTFDataset.N);
0095 eqHRTFDataset_temp = supdeq_hrtf2sfd(HRTF_equalized_ALFE_L,HRTF_equalized_ALFE_R,eqHRTFDataset.N,samplingGrid,fs,'ak');
0096 
0097 %Write new output struct
0098 eqHRTFDataset_alfe = eqHRTFDataset;
0099 eqHRTFDataset_alfe.Hl_nm = eqHRTFDataset_temp.Hl_nm;
0100 eqHRTFDataset_alfe.Hr_nm = eqHRTFDataset_temp.Hr_nm;
0101 eqHRTFDataset_alfe.HRTF_L = HRTF_equalized_ALFE_L;
0102 eqHRTFDataset_alfe.HRTF_R = HRTF_equalized_ALFE_R;
0103 eqHRTFDataset_alfe.alfe = 1;
0104 eqHRTFDataset_alfe.alfeXover = xover;
0105 clear eqHRTFDataset_temp
0106 
0107 fprintf('\nDone with ALFE...\n');