% SUpDEq - Spatial Upsampling by Directional Equalization
script supdeq_demo_MCA
This script presents magnitude-corrected and time-aligned interpolation,
hereafter abbreviated as MCA interpolation. The method combines a
perceptually motivated post-interpolation magnitude correction with
time-aligned interpolation. This demo script uses SUpDEq processing for
time alignment and spherical harmonics (SH) interpolation for spatial
upsampling of a sparse HRTF set to a dense HRTF set. The script compares
the results of MCA interpolation and conventional time-aligned
interpolation. The proposed magnitude correction as well as various
time-alignment and interpolation approaches are implemented in the
function supdeq_interpHRTF, used throughout this demo script.
Reference:
J. M. Arend, C. Pörschmann, S. Weinzierl, F. Brinkmann,
"Magnitude-Corrected and Time-Aligned Interpolation of
Head-Related Transfer Functions," (Manuscript submitted for publication).
(C) 2022 by JMA, Johannes M. Arend
TU Berlin
Audio Communication Group0001 %% SUpDEq - Spatial Upsampling by Directional Equalization 0002 % 0003 % script supdeq_demo_MCA 0004 % 0005 % This script presents magnitude-corrected and time-aligned interpolation, 0006 % hereafter abbreviated as MCA interpolation. The method combines a 0007 % perceptually motivated post-interpolation magnitude correction with 0008 % time-aligned interpolation. This demo script uses SUpDEq processing for 0009 % time alignment and spherical harmonics (SH) interpolation for spatial 0010 % upsampling of a sparse HRTF set to a dense HRTF set. The script compares 0011 % the results of MCA interpolation and conventional time-aligned 0012 % interpolation. The proposed magnitude correction as well as various 0013 % time-alignment and interpolation approaches are implemented in the 0014 % function supdeq_interpHRTF, used throughout this demo script. 0015 % 0016 % Reference: 0017 % J. M. Arend, C. Pörschmann, S. Weinzierl, F. Brinkmann, 0018 % "Magnitude-Corrected and Time-Aligned Interpolation of 0019 % Head-Related Transfer Functions," (Manuscript submitted for publication). 0020 % 0021 % (C) 2022 by JMA, Johannes M. Arend 0022 % TU Berlin 0023 % Audio Communication Group 0024 0025 %% (1) - Define sparse and dense grid 0026 0027 %Sparse Lebedev Grid 0028 %Azimuth, Colatitude, Weight 0029 Ns = 3; 0030 sgS = supdeq_lebedev([],Ns); 0031 0032 %Dense Lebedev Grid (Target grid for upsampling) 0033 %Azimuth, Colatitude, Weight 0034 Nd = 44; 0035 sgD = supdeq_lebedev([],Nd); 0036 0037 %% (2) - Get sparse HRTF set 0038 0039 %Get sparse KU100 HRTF set based on "HRIR_L2702.sofa" dataset in SH domain 0040 0041 sparseHRTF = supdeq_getSparseDataset(sgS,Ns,44,'ku100'); 0042 fs = sparseHRTF.f(end)*2; 0043 0044 %% (3) - Perform spatial upsampling to dense grid 0045 0046 headRadius = 0.0875; %Also default value in function 0047 0048 %Interpolation / upsampling with SH interpolation but without any 0049 %pre/post-processing ('none'), called SH here 0050 interpHRTF_sh = supdeq_interpHRTF(sparseHRTF,sgD,'None','SH',nan,headRadius); 0051 0052 %Interpolation / upsampling with SUpDEq time alignment and SH 0053 %interpolation, called 'conventional' 0054 interpHRTF_con = supdeq_interpHRTF(sparseHRTF,sgD,'SUpDEq','SH',nan,headRadius); 0055 0056 %Interpolation / upsampling with MCA interpolation, i.e., in this example 0057 %SUpDEq time alignment and SH interpolation plus post-interpolation 0058 %magnitude correction, as described in the paper 0059 %Maximum boost of magnitude correction set to inf (unlimited gain), 0060 %resulting in no soft-limiting (as presented in the paper) 0061 %The other MCA settings are by default as described in the paper, i.e.: 0062 % (a) Magnitude-correction filters are set to 0dB below spatial aliasing 0063 % frequency fA 0064 % (b) Soft-limiting knee set to 0dB (no knee) 0065 % (c) Magnitude-correction filters are designed as minimum phase filters 0066 interpHRTF_mca = supdeq_interpHRTF(sparseHRTF,sgD,'SUpDEq','SH',inf,headRadius); 0067 0068 %The resulting datasets can be saved as SOFA files using the function 0069 %supdeq_writeSOFAobj 0070 0071 %% (4) - Optional: Get reference HRTF set 0072 0073 %Get reference dataset (using the "supdeq_getSparseDataset" function for 0074 %convenience here) 0075 refHRTF = supdeq_getSparseDataset(sgD,Nd,44,'ku100'); 0076 0077 %Also get HRIRs 0078 refHRTF.HRIR_L = ifft(AKsingle2bothSidedSpectrum(refHRTF.HRTF_L.')); 0079 refHRTF.HRIR_L = refHRTF.HRIR_L(1:end/refHRTF.FFToversize,:); 0080 refHRTF.HRIR_R = ifft(AKsingle2bothSidedSpectrum(refHRTF.HRTF_R.')); 0081 refHRTF.HRIR_R = refHRTF.HRIR_R(1:end/refHRTF.FFToversize,:); 0082 0083 %% (5) - Optional: Plot HRIRs 0084 0085 %Plot frontal left-ear HRIR (Az = 0, Col = 90) of conventional and MCA 0086 %interpolated dataset. Differences are small. 0087 idFL = 16; %ID in sgD 0088 supdeq_plotIR(interpHRTF_con.HRIR_L(:,idFL),interpHRTF_mca.HRIR_L(:,idFL),[],[],8); 0089 0090 %Plot contralateral left-ear HRIR (Az = 270, Col = 90) of SH-only and MCA 0091 %interpolated dataset. Differences are strong 0092 idCL = 2042; %ID in sgD 0093 supdeq_plotIR(interpHRTF_sh.HRIR_L(:,idCL),interpHRTF_mca.HRIR_L(:,idCL),[],[],8); 0094 0095 %Plot contralateral left-ear HRIR (Az = 270, Col = 90) of conventional and MCA 0096 %interpolated dataset. Differences are still significant. 0097 supdeq_plotIR(interpHRTF_con.HRIR_L(:,idCL),interpHRTF_mca.HRIR_L(:,idCL),[],[],8); 0098 0099 %Plot contralateral left-ear HRIR (Az = 270, Col = 90) of SH-only 0100 %interpolated dataset and reference. 0101 supdeq_plotIR(interpHRTF_sh.HRIR_L(:,idCL),refHRTF.HRIR_L(:,idCL),[],[],8); 0102 0103 %Plot contralateral left-ear HRIR (Az = 270, Col = 90) of conventional 0104 %interpolated dataset and reference. 0105 supdeq_plotIR(interpHRTF_con.HRIR_L(:,idCL),refHRTF.HRIR_L(:,idCL),[],[],8); 0106 0107 %Plot contralateral left-ear HRIR (Az = 270, Col = 90) of MCA 0108 %interpolated dataset and reference. 0109 %MCA interpolated HRTF is much closer to the reference than HRTF from 0110 %conventional interpolation. The "bump" above 10 kHz is corrected through 0111 %the magnitude correction. Interpolation errors (in this case spatial 0112 %aliasing errors at the contralateral ear) are reduced 0113 supdeq_plotIR(interpHRTF_mca.HRIR_L(:,idCL),refHRTF.HRIR_L(:,idCL),[],[],8); 0114 0115 %% (6) - Optional: Calculate spectral difference 0116 0117 %Calculate left-ear spectral differences in dB by spectral division of upsampled 0118 %HRTF set with reference HRTF set. Not the same as the magnitude error in 0119 %auditory bands presented in the paper! 0120 0121 specDiff_sh = supdeq_calcSpectralDifference_HRIR(interpHRTF_sh.HRIR_L,refHRTF.HRIR_L,fs,16); 0122 specDiff_con = supdeq_calcSpectralDifference_HRIR(interpHRTF_con.HRIR_L,refHRTF.HRIR_L,fs,16); 0123 specDiff_mca = supdeq_calcSpectralDifference_HRIR(interpHRTF_mca.HRIR_L,refHRTF.HRIR_L,fs,16); 0124 0125 %Quick plot of spectral difference averaged over all directions of the 0126 %2702-point Lebedev grid. As shown in the paper, MCA provides the most 0127 %significant benefit in the critical contralateral region. Thus, when 0128 %averaging over all positions, the benefit of MCA seems smaller. The 0129 %analysis in the paper as well as the provided audio examples reveal 0130 %in much more detail the considerable improvements that the proposed 0131 %magnitude correction provides. 0132 AKf(18,9); 0133 semilogx(specDiff_sh.f,specDiff_sh.specDifference,'LineWidth',1.5); 0134 hold on; 0135 semilogx(specDiff_con.f,specDiff_con.specDifference,'LineWidth',1.5); 0136 semilogx(specDiff_mca.f,specDiff_mca.specDifference,'LineWidth',1.5); 0137 xlim([500 20000]) 0138 legend('SH','SH SUpDEq W/O MC','SH SUpDEq W/ MC','Location','NorthWest'); 0139 xlabel('Frequency in Hz'); 0140 ylabel('Magnitude Error in dB'); 0141 grid on;