% SUpDEq - Spatial Upsampling by Directional Equalization
function results = supdeq_calcSpectralDifference_HRIR( testHRIRdataset, referenceHRIRdataset, fs, FFToversize, filterMode, plotResults)
This function calculates the spectral differences between a test and a
reference HRTF dataset at given spatial sampling points
Output:
results - Results struct with spectral difference per
frequency and sampling point in dB (20*log10)
(specDifferencePerPoint), spectral
difference averaged over all spatial sampling
points in dB (specDifference) as well
as mean of specDifference in dB
--> Spectral difference averaged over all
spatial sampling points and frequencies
(meanSpecDifference). Additionally,
meanSpecDifference_freqRange provides the mean
spectral difference in a frequency range between
100Hz and 12kHz. "f" is a frequency vector
Optional: filterMode = true
Results struct with spectral difference averaged
over all spatial sampling points in dB
(20*log10) and for every of the 40 filter bands
("specDiffPerBand"). Additionally,
"specDiffL_meanPerBand" returns the mean spectral
difference per band (mean per band of
specDiffPerBand), fc returns the center
frequencies of the auditory filterbank, and f
is a frequency vector.
Input:
testHRIRdataset - [N x M] array with HRIRs of the test HRIR
dataset, with N samples and M channels
referenceHRIRdataset - [N x M] array with HRIRs of the reference HRIR
dataset, with N samples and M channels
fs - Sampling Rate
Default: 48000
FFToversize - FFToversize rises the FFT Blocksize [default = 1]
A FFT of the blocksize (FFToversize*NFFT) is applied
to the time domain data, where NFFT is determinded
as the next power of two of the signalSize which is
signalSize = (lastSample-firstSample).
filterMode - Filter mode true/false. Use auditory filterbank
with 40 bands (fLow = 50 Hz, fHigh = 20kHz)
Default: false
plotResults - plotResults true/false. Simple plot...
Default: false
Dependencies: Auditory Modeling Toolbox (AMT), if filterMode = true
(C) 2018 by JMA, Johannes M. Arend
TH Köln - University of Applied Sciences
Institute of Communications Engineering
Department of Acoustics and Audio Signal Processing0001 %% SUpDEq - Spatial Upsampling by Directional Equalization 0002 % 0003 % function results = supdeq_calcSpectralDifference_HRIR( testHRIRdataset, referenceHRIRdataset, fs, FFToversize, filterMode, plotResults) 0004 % 0005 % This function calculates the spectral differences between a test and a 0006 % reference HRTF dataset at given spatial sampling points 0007 % 0008 % Output: 0009 % results - Results struct with spectral difference per 0010 % frequency and sampling point in dB (20*log10) 0011 % (specDifferencePerPoint), spectral 0012 % difference averaged over all spatial sampling 0013 % points in dB (specDifference) as well 0014 % as mean of specDifference in dB 0015 % --> Spectral difference averaged over all 0016 % spatial sampling points and frequencies 0017 % (meanSpecDifference). Additionally, 0018 % meanSpecDifference_freqRange provides the mean 0019 % spectral difference in a frequency range between 0020 % 100Hz and 12kHz. "f" is a frequency vector 0021 % 0022 % Optional: filterMode = true 0023 % Results struct with spectral difference averaged 0024 % over all spatial sampling points in dB 0025 % (20*log10) and for every of the 40 filter bands 0026 % ("specDiffPerBand"). Additionally, 0027 % "specDiffL_meanPerBand" returns the mean spectral 0028 % difference per band (mean per band of 0029 % specDiffPerBand), fc returns the center 0030 % frequencies of the auditory filterbank, and f 0031 % is a frequency vector. 0032 % 0033 % Input: 0034 % testHRIRdataset - [N x M] array with HRIRs of the test HRIR 0035 % dataset, with N samples and M channels 0036 % referenceHRIRdataset - [N x M] array with HRIRs of the reference HRIR 0037 % dataset, with N samples and M channels 0038 % fs - Sampling Rate 0039 % Default: 48000 0040 % FFToversize - FFToversize rises the FFT Blocksize [default = 1] 0041 % A FFT of the blocksize (FFToversize*NFFT) is applied 0042 % to the time domain data, where NFFT is determinded 0043 % as the next power of two of the signalSize which is 0044 % signalSize = (lastSample-firstSample). 0045 % filterMode - Filter mode true/false. Use auditory filterbank 0046 % with 40 bands (fLow = 50 Hz, fHigh = 20kHz) 0047 % Default: false 0048 % plotResults - plotResults true/false. Simple plot... 0049 % Default: false 0050 % 0051 % Dependencies: Auditory Modeling Toolbox (AMT), if filterMode = true 0052 % 0053 % (C) 2018 by JMA, Johannes M. Arend 0054 % TH Köln - University of Applied Sciences 0055 % Institute of Communications Engineering 0056 % Department of Acoustics and Audio Signal Processing 0057 0058 function results = supdeq_calcSpectralDifference_HRIR( testHRIRdataset, referenceHRIRdataset, fs, FFToversize, filterMode, plotResults) 0059 0060 %Default sampling grid (Lebedev grid with 2702 nodes) 0061 if nargin < 3 || isempty(fs) 0062 fs = 48000; 0063 end 0064 0065 if nargin < 4 || isempty(FFToversize) 0066 FFToversize = 1; 0067 end 0068 0069 if nargin < 5 || isempty(filterMode) 0070 filterMode = false; 0071 end 0072 0073 if nargin < 6 0074 plotResults = false; 0075 end 0076 0077 %% Get NFFT based on number of samples in HRIR arrays and FFToversize 0078 NFFT = size(testHRIRdataset,1); 0079 if size(referenceHRIRdataset,1) > NFFT 0080 NFFT = size(referenceHRIRdataset,1); 0081 end 0082 0083 NFFT = (2^nextpow2(NFFT)); 0084 if FFToversize > 1 0085 NFFT = NFFT*FFToversize; 0086 end 0087 %Get frequency vector 0088 f = linspace(0,fs/2,NFFT/2+1); 0089 0090 if size(testHRIRdataset,2) ~= size(referenceHRIRdataset,2) 0091 error('Number of sampling points in both HRIR datasets must be equal!'); 0092 end 0093 0094 numPoints = size(testHRIRdataset,2); 0095 0096 %% Filter mode off 0097 if ~filterMode 0098 0099 %Transform HRIRs to frequency domain 0100 fprintf('Transforming 2 x %d HRIRs to frequency domain. This may take some time...\n',size(testHRIRdataset,2)); 0101 HRTFs_Test = abs(fft(testHRIRdataset,NFFT)); 0102 HRTFs_Test = HRTFs_Test(1:end/2+1,:); 0103 HRTFs_Ref = abs(fft(referenceHRIRdataset,NFFT)); 0104 HRTFs_Ref = HRTFs_Ref(1:end/2+1,:); 0105 0106 %Calculate spectral difference 0107 disp('Calculating spectral differences...'); 0108 specDiff = abs( 20*log10(abs(HRTFs_Ref)) - 20*log10(abs(HRTFs_Test))); 0109 specDiffPerPoint = specDiff; 0110 specDiff = sum(specDiff,2); 0111 specDiff = specDiff/numPoints; 0112 0113 %Get values from 100Hz - 12 kHz 0114 lowBin = round(size(f,2) / f(end) * 100) + 1; 0115 highBin = round(size(f,2) / f(end) * 12000) + 1; 0116 0117 %Write results struct 0118 results.specDifferencePerPoint = specDiffPerPoint; 0119 results.specDifference = specDiff; 0120 results.meanSpecDifference = mean(results.specDifference); 0121 results.meanSpecDifference_freqRange = mean(results.specDifference(lowBin:highBin,:)); 0122 results.f = f; 0123 0124 disp('Done with calculation of spectral differences...'); 0125 0126 if plotResults 0127 figure; 0128 semilogx(results.f,results.specDifference(:,1),'k','LineWidth',1.5); 0129 grid on; 0130 xlim([20,20000]); 0131 xlabel('Frequency [Hz]'); 0132 ylabel('Spectral difference [dB]'); 0133 end 0134 0135 end 0136 0137 %% Filter mode on 0138 if filterMode 0139 0140 %Transform HRIRs to 3D array 0141 HRIRs_Test = zeros(size(testHRIRdataset,1), 1, numPoints); 0142 HRIRs_Ref = zeros(size(referenceHRIRdataset,1), 1, numPoints); 0143 0144 %Kind of a workaround, but 4-D arrays needed later.... 0145 HRIRs_Test(:,1,:) = reshape(testHRIRdataset,[size(HRIRs_Test,1),1,size(HRIRs_Test,3)]); 0146 HRIRs_Ref(:,1,:) = reshape(referenceHRIRdataset,[size(HRIRs_Ref,1),1,size(HRIRs_Ref,3)]); 0147 0148 %Zero-pad to 512 samples (if length is smaller) to get higher frequency 0149 %resolution. This only works if HRIRs_Test and HRIRs_Ref have equal 0150 %length! 0151 if NFFT < 512 0152 newLength = 512; 0153 else 0154 newLength = NFFT; 0155 end 0156 zeroPad = zeros(newLength-size(HRIRs_Test,1),1,numPoints); 0157 HRIRs_Test = [HRIRs_Test;zeroPad]; 0158 HRIRs_Ref = [HRIRs_Ref;zeroPad]; 0159 %Calc f needed later for results struct 0160 f = linspace(0,fs/2,newLength/2+1); 0161 0162 %Filter HRIRs with auditory filter bank (settings lead to 40 bands) 0163 %Output is a 4D array with [samples X filterBand x channel x 0164 %samplingPoint] 0165 disp('Applying auditory filter bank. This may take some time...'); 0166 HRIRs_Test_filt = auditoryfilterbank(HRIRs_Test,fs,'flow',50,'fhigh',20000); 0167 [HRIRs_Ref_filt,fc] = auditoryfilterbank(HRIRs_Ref,fs,'flow',50,'fhigh',20000); 0168 disp('Done with filtering...'); 0169 0170 %Perform FFT 0171 disp('Calculating spectral differences...'); 0172 HRTFs_Test_filt = fft(HRIRs_Test_filt); 0173 HRTFs_Test_filt = HRTFs_Test_filt(1:size(HRTFs_Test_filt,1)/2+1,:,:,:); 0174 HRTFs_Ref_filt = fft(HRIRs_Ref_filt); 0175 HRTFs_Ref_filt = HRTFs_Ref_filt(1:size(HRTFs_Ref_filt,1)/2+1,:,:,:); 0176 0177 %Calculate spectral difference for each band 0178 specDiff = squeeze(abs(20*log10(abs(HRTFs_Ref_filt(:,:,1,:))) - 20*log10(abs(HRTFs_Test_filt(:,:,1,:))))); 0179 specDiffPerPoint = specDiff; 0180 %Average across spatial sampling points 0181 specDiff = sum(specDiff,3); 0182 for i = 1:size(specDiff,2) 0183 specDiff(:,i) = specDiff(:,i)/numPoints; 0184 end 0185 %Get mean per band 0186 for i = 1:size(specDiff,2) 0187 specDiff_meanPerBand(i,:) = mean(specDiff(:,i)); 0188 end 0189 0190 %Write results struct 0191 results.specDiffPerPoint = specDiffPerPoint; 0192 results.specDiffPerBand = specDiff; 0193 results.specDiff_meanPerBand = specDiff_meanPerBand; 0194 results.fc = fc; 0195 results.f = f; 0196 0197 disp('Done with calculation of spectral differences...'); 0198 0199 if plotResults 0200 figure; 0201 semilogx(results.fc,results.specDiff_meanPerBand,'k','LineWidth',1.5); 0202 grid on; 0203 xlim([20,20000]); 0204 xlabel('Frequency [Hz]'); 0205 ylabel('Spectral Difference [dB]'); 0206 end 0207 0208 end 0209 0210 0211 end 0212