% SUpDEq - Spatial Upsampling by Directional Equalization
function results = supdeq_calcSpectralDifference( testHRTFdataset, referenceHRTFdataset, samplingGrid, fs, 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) for left and right channel
spectral difference averaged over all spatial
sampling points in dB (20*log10) for left and
right channel (specDifference) as well as mean of
specDifference in dB for left and right channel
--> 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 per
frequency and sampling point in dB (20*log10)
(specDifferencePerPoint) for left and right
channel, spectral difference averaged over all spatial
sampling points in dB (20*log10) for left and
right channel and for every of the 40 filter bands
("specDiffPerBand_L/R"). Additionally,
"specDiffL_meanPerBand" returns the mean spectral
difference per band (mean per band of
specDiffPerBand_L/R), fc returns the center
frequencies of the auditory filterbank, and f
is a frequency vector.
Input:
testHRTFdataset - Struct with (SH-Coefficients) of the test
HRTF dataset (for example the de-equalized HRTF set
or any other HRTF set based on the referenceHRTFdataset)
for the left (Hl_nm) and right (Hr_nm), channel/ear,
absolute frequency scale f, transform order N, and FFToversize
referenceHRTFdataset - Struct with (SH-Coefficients) of the reference
HRTF dataset for the left (Hl_nm) and right (Hr_nm), channel/ear,
absolute frequency scale f, transform order N, and FFToversize
NOTE: testHRTFdataset and referenceHRTFdataset
need to have the same FFT size and FFToversize!!
samplingGrid - Spatial sampling grid. Must be a Qx2 matrix where
the first column holds the azimuth and the second
the elevation.
If no samplingGrid (or []) is passed, the default
sampling grid is passed (Lebedev grid with 2702
nodes)
fs - Sampling Rate
Default: 48000
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( testHRTFdataset, referenceHRTFdataset, samplingGrid, fs, 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) for left and right channel 0012 % spectral difference averaged over all spatial 0013 % sampling points in dB (20*log10) for left and 0014 % right channel (specDifference) as well as mean of 0015 % specDifference in dB for left and right channel 0016 % --> Spectral difference averaged over all spatial 0017 % sampling points and frequencies 0018 % (meanSpecDifference). Additionally, 0019 % meanSpecDifference_freqRange provides the mean 0020 % spectral difference in a frequency range between 0021 % 100Hz and 12kHz. "f" is a frequency vector 0022 % 0023 % Optional: filterMode = true 0024 % Results struct with spectral difference per 0025 % frequency and sampling point in dB (20*log10) 0026 % (specDifferencePerPoint) for left and right 0027 % channel, spectral difference averaged over all spatial 0028 % sampling points in dB (20*log10) for left and 0029 % right channel and for every of the 40 filter bands 0030 % ("specDiffPerBand_L/R"). Additionally, 0031 % "specDiffL_meanPerBand" returns the mean spectral 0032 % difference per band (mean per band of 0033 % specDiffPerBand_L/R), fc returns the center 0034 % frequencies of the auditory filterbank, and f 0035 % is a frequency vector. 0036 % 0037 % Input: 0038 % testHRTFdataset - Struct with (SH-Coefficients) of the test 0039 % HRTF dataset (for example the de-equalized HRTF set 0040 % or any other HRTF set based on the referenceHRTFdataset) 0041 % for the left (Hl_nm) and right (Hr_nm), channel/ear, 0042 % absolute frequency scale f, transform order N, and FFToversize 0043 % referenceHRTFdataset - Struct with (SH-Coefficients) of the reference 0044 % HRTF dataset for the left (Hl_nm) and right (Hr_nm), channel/ear, 0045 % absolute frequency scale f, transform order N, and FFToversize 0046 % NOTE: testHRTFdataset and referenceHRTFdataset 0047 % need to have the same FFT size and FFToversize!! 0048 % samplingGrid - Spatial sampling grid. Must be a Qx2 matrix where 0049 % the first column holds the azimuth and the second 0050 % the elevation. 0051 % If no samplingGrid (or []) is passed, the default 0052 % sampling grid is passed (Lebedev grid with 2702 0053 % nodes) 0054 % fs - Sampling Rate 0055 % Default: 48000 0056 % filterMode - Filter mode true/false. Use auditory filterbank 0057 % with 40 bands (fLow = 50 Hz, fHigh = 20kHz) 0058 % Default: false 0059 % plotResults - plotResults true/false. Simple plot... 0060 % Default: false 0061 % 0062 % Dependencies: Auditory Modeling Toolbox (AMT), if filterMode = true 0063 % 0064 % (C) 2018 by JMA, Johannes M. Arend 0065 % TH Köln - University of Applied Sciences 0066 % Institute of Communications Engineering 0067 % Department of Acoustics and Audio Signal Processing 0068 0069 function results = supdeq_calcSpectralDifference( testHRTFdataset, referenceHRTFdataset, samplingGrid, fs, filterMode, plotResults) 0070 0071 %Default sampling grid (Lebedev grid with 2702 nodes) 0072 if nargin < 3 || isempty(samplingGrid) 0073 samplingGrid = supdeq_lebedev(2702); 0074 end 0075 0076 if nargin < 4 || isempty(fs) 0077 fs = 48000; 0078 end 0079 0080 if nargin < 5 || isempty(filterMode) 0081 filterMode = false; 0082 end 0083 0084 if nargin < 6 0085 plotResults = false; 0086 end 0087 0088 0089 %% Filter mode off 0090 if ~filterMode 0091 0092 %Get HRTFs according to sampling grid 0093 fprintf('Extracting 2 x %d HRTFs. This may take some time...\n',size(samplingGrid,1)); 0094 [HRTFs_Test_L,HRTFs_Test_R] = supdeq_getArbHRTF(testHRTFdataset,samplingGrid,[],[],'ak'); 0095 [HRTFs_Ref_L,HRTFs_Ref_R] = supdeq_getArbHRTF(referenceHRTFdataset,samplingGrid,[],[],'ak'); 0096 fprintf('2 x %d HRTFs extracted...\n',size(samplingGrid,1)) 0097 0098 %Calculate spectral difference 0099 disp('Calculating spectral differences...'); 0100 specDiffL = abs( 20*log10(abs(HRTFs_Ref_L.')) - 20*log10(abs(HRTFs_Test_L.'))); 0101 specDiffR = abs( 20*log10(abs(HRTFs_Ref_R.')) - 20*log10(abs(HRTFs_Test_R.'))); 0102 specDiffLperPoint = specDiffL; 0103 specDiffRperPoint = specDiffR; 0104 specDiffL = sum(specDiffL,2); 0105 specDiffR = sum(specDiffR,2); 0106 specDiffL = specDiffL/size(samplingGrid,1); 0107 specDiffR = specDiffR/size(samplingGrid,1); 0108 0109 %Get values from 100Hz - 12 kHz 0110 lowBin = round(size(referenceHRTFdataset.f,2) / referenceHRTFdataset.f(end) * 100) + 1; 0111 highBin = round(size(referenceHRTFdataset.f,2) / referenceHRTFdataset.f(end) * 12000) + 1; 0112 0113 %Write results struct 0114 results.specDifferencePerPoint = {specDiffLperPoint, specDiffRperPoint}; 0115 results.specDifference = [specDiffL,specDiffR]; 0116 results.meanSpecDifference = mean(results.specDifference); 0117 results.meanSpecDifference_freqRange = mean(results.specDifference(lowBin:highBin,:)); 0118 results.f = referenceHRTFdataset.f; 0119 0120 disp('Done with calculation of spectral differences...'); 0121 0122 if plotResults 0123 figure; 0124 semilogx(results.f,results.specDifference(:,1),'k','LineWidth',1.5); 0125 hold on; 0126 semilogx(results.f,results.specDifference(:,2),'r','LineWidth',1.5); 0127 grid on; 0128 xlim([20,20000]); 0129 xlabel('Frequency [Hz]'); 0130 ylabel('Spectral difference [dB]'); 0131 legend('Left','Right','Location','NorthWest'); 0132 end 0133 0134 end 0135 0136 %% Filter mode on 0137 if filterMode 0138 0139 %Get HRIRs according to sampling grid 0140 HRIRs_Test = zeros((size(testHRTFdataset.f,2)*2-2) / testHRTFdataset.FFToversize,2,size(samplingGrid,1)); 0141 HRIRs_Ref = zeros((size(referenceHRTFdataset.f,2)*2-2) / referenceHRTFdataset.FFToversize,2,size(samplingGrid,1)); 0142 0143 fprintf('Extracting 2 x %d HRIRs. This may take some time...\n',size(samplingGrid,1)); 0144 [HRIRs_Test_L,HRIRs_Test_R] = supdeq_getArbHRIR(testHRTFdataset,samplingGrid,[],[],'ak'); 0145 [HRIRs_Ref_L,HRIRs_Ref_R] = supdeq_getArbHRIR(referenceHRTFdataset,samplingGrid,[],[],'ak'); 0146 %Kind of a workaround, but 4-D arrays needed later.... 0147 HRIRs_Test(:,1,:) = reshape(HRIRs_Test_L,[size(HRIRs_Test,1),1,size(HRIRs_Test,3)]); 0148 HRIRs_Test(:,2,:) = reshape(HRIRs_Test_R,[size(HRIRs_Test,1),1,size(HRIRs_Test,3)]); 0149 HRIRs_Ref(:,1,:) = reshape(HRIRs_Ref_L,[size(HRIRs_Ref,1),1,size(HRIRs_Ref,3)]); 0150 HRIRs_Ref(:,2,:) = reshape(HRIRs_Ref_R,[size(HRIRs_Ref,1),1,size(HRIRs_Ref,3)]); 0151 fprintf('2 x %d HRIRs extracted...\n',size(samplingGrid,1)) 0152 0153 %Zero-pad to 512 samples (if length is smaller) to get higher frequency 0154 %resolution. This only works if HRIRs_Test and HRIRs_Ref have equal 0155 %length! 0156 newLength = 512; 0157 if size(HRIRs_Test,1) < newLength 0158 zeroPad = zeros(newLength-size(HRIRs_Test,1),2,size(samplingGrid,1)); 0159 HRIRs_Test = [HRIRs_Test;zeroPad]; 0160 HRIRs_Ref = [HRIRs_Ref;zeroPad]; 0161 %Calc f needed later for results struct 0162 f = linspace(0,fs/2,newLength/2+1); 0163 else 0164 f = linspace(0,fs/2,size(HRIRs_Test,1)/2+1); 0165 end 0166 0167 %Filter HRIRs with auditory filter bank (settings lead to 40 bands) 0168 %Output is a 4D array with [samples X filterBand x channel x 0169 %samplingPoint] 0170 disp('Applying auditory filter bank. This may take some time...'); 0171 HRIRs_Test_filt = auditoryfilterbank(HRIRs_Test,fs,'flow',50,'fhigh',20000); 0172 [HRIRs_Ref_filt,fc] = auditoryfilterbank(HRIRs_Ref,fs,'flow',50,'fhigh',20000); 0173 disp('Done with filtering...'); 0174 0175 %Perform FFT 0176 disp('Calculating spectral differences...'); 0177 HRTFs_Test_filt = fft(HRIRs_Test_filt); 0178 HRTFs_Test_filt = HRTFs_Test_filt(1:size(HRTFs_Test_filt,1)/2+1,:,:,:); 0179 HRTFs_Ref_filt = fft(HRIRs_Ref_filt); 0180 HRTFs_Ref_filt = HRTFs_Ref_filt(1:size(HRTFs_Ref_filt,1)/2+1,:,:,:); 0181 0182 %Calculate spectral difference for each band 0183 specDiffL = squeeze(abs(20*log10(abs(HRTFs_Ref_filt(:,:,1,:))) - 20*log10(abs(HRTFs_Test_filt(:,:,1,:))))); 0184 specDiffR = squeeze(abs(20*log10(abs(HRTFs_Ref_filt(:,:,2,:))) - 20*log10(abs(HRTFs_Test_filt(:,:,2,:))))); 0185 specDiffLperPoint = specDiffL; 0186 specDiffRperPoint = specDiffR; 0187 %Average across spatial sampling points 0188 specDiffL = sum(specDiffL,3); 0189 specDiffR = sum(specDiffR,3); 0190 for i = 1:size(specDiffL,2) 0191 specDiffL(:,i) = specDiffL(:,i)/size(samplingGrid,1); 0192 specDiffR(:,i) = specDiffR(:,i)/size(samplingGrid,1); 0193 end 0194 %Get mean per band 0195 for i = 1:size(specDiffL,2) 0196 specDiffL_meanPerBand(i,:) = mean(specDiffL(:,i)); 0197 specDiffR_meanPerBand(i,:) = mean(specDiffR(:,i)); 0198 end 0199 0200 %Write results struct 0201 results.specDiffPerPointAndBand_L = specDiffLperPoint; 0202 results.specDiffPerPointAndBand_R = specDiffRperPoint; 0203 results.specDiffPerBand_L = specDiffL; 0204 results.specDiffPerBand_R = specDiffR; 0205 results.specDiffL_meanPerBand = specDiffL_meanPerBand; 0206 results.specDiffR_meanPerBand = specDiffR_meanPerBand; 0207 results.fc = fc; 0208 results.f = f; 0209 0210 disp('Done with calculation of spectral differences...'); 0211 0212 if plotResults 0213 figure; 0214 semilogx(results.fc,results.specDiffL_meanPerBand,'k','LineWidth',1.5); 0215 hold on; 0216 semilogx(results.fc,results.specDiffR_meanPerBand,'r','LineWidth',1.5); 0217 grid on; 0218 xlim([20,20000]); 0219 xlabel('Frequency [Hz]'); 0220 ylabel('Spectral Difference [dB]'); 0221 legend('Left','Right','Location','NorthWest'); 0222 end 0223 0224 end 0225 0226 0227 end 0228