-
Notifications
You must be signed in to change notification settings - Fork 7
/
Copy pathsupdeq_demo_MCA.m
executable file
·166 lines (132 loc) · 6.77 KB
/
supdeq_demo_MCA.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
%% 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/2023 by JMA, Johannes M. Arend
% TU Berlin
% Audio Communication Group
%% (1) - Define sparse and dense grid
%Sparse Lebedev Grid
%Azimuth, Colatitude, Weight
Ns = 3;
sgS = supdeq_lebedev([],Ns);
%Dense Lebedev Grid (Target grid for upsampling)
%Azimuth, Colatitude, Weight
Nd = 44;
sgD = supdeq_lebedev([],Nd);
%% (2) - Get sparse HRTF set
%Get sparse KU100 HRTF set based on "HRIR_L2702.sofa" dataset in SH domain
sparseHRTF = supdeq_getSparseDataset(sgS,Ns,44,'ku100');
fs = sparseHRTF.fs;
%% (3) - Perform spatial upsampling to dense grid
headRadius = 0.0875; %Also default value in function
%Interpolation / upsampling with SH interpolation but without any
%pre/post-processing ('none'), called SH here
interpHRTF_sh = supdeq_interpHRTF(sparseHRTF,sgD,'None','SH',nan,headRadius);
%Interpolation / upsampling with SUpDEq time alignment and SH
%interpolation, called 'conventional'
interpHRTF_con = supdeq_interpHRTF(sparseHRTF,sgD,'SUpDEq','SH',nan,headRadius);
%Interpolation / upsampling with MCA interpolation, i.e., in this example
%SUpDEq time alignment and SH interpolation plus post-interpolation
%magnitude correction, as described in the paper
%Maximum boost of magnitude correction set to inf (unlimited gain),
%resulting in no soft-limiting (as presented in the paper)
%The other MCA settings are by default as described in the paper, i.e.:
% (a) Magnitude-correction filters are set to 0dB below spatial aliasing
% frequency fA
% (b) Soft-limiting knee set to 0dB (no knee)
% (c) Magnitude-correction filters are designed as minimum phase filters
interpHRTF_mca = supdeq_interpHRTF(sparseHRTF,sgD,'SUpDEq','SH',inf,headRadius);
%The resulting datasets can be saved as SOFA files using the function
%supdeq_writeSOFAobj
%% (4) - Optional: Get reference HRTF set
%Get reference dataset (using the "supdeq_getSparseDataset" function for
%convenience here)
refHRTF = supdeq_getSparseDataset(sgD,Nd,44,'ku100');
%Also get HRIRs
refHRTF.HRIR_L = ifft(AKsingle2bothSidedSpectrum(refHRTF.HRTF_L.'));
refHRTF.HRIR_L = refHRTF.HRIR_L(1:end/refHRTF.FFToversize,:);
refHRTF.HRIR_R = ifft(AKsingle2bothSidedSpectrum(refHRTF.HRTF_R.'));
refHRTF.HRIR_R = refHRTF.HRIR_R(1:end/refHRTF.FFToversize,:);
%% (5) - Optional: Plot HRIRs
%Plot frontal left-ear HRIR (Az = 0, Col = 90) of conventional and MCA
%interpolated dataset. Differences are small.
idFL = 16; %ID in sgD
supdeq_plotIR(interpHRTF_con.HRIR_L(:,idFL),interpHRTF_mca.HRIR_L(:,idFL),[],[],8);
%Plot contralateral left-ear HRIR (Az = 270, Col = 90) of SH-only and MCA
%interpolated dataset. Differences are strong
idCL = 2042; %ID in sgD
supdeq_plotIR(interpHRTF_sh.HRIR_L(:,idCL),interpHRTF_mca.HRIR_L(:,idCL),[],[],8);
%Plot contralateral left-ear HRIR (Az = 270, Col = 90) of conventional and MCA
%interpolated dataset. Differences are still significant.
supdeq_plotIR(interpHRTF_con.HRIR_L(:,idCL),interpHRTF_mca.HRIR_L(:,idCL),[],[],8);
%Plot contralateral left-ear HRIR (Az = 270, Col = 90) of SH-only
%interpolated dataset and reference.
supdeq_plotIR(interpHRTF_sh.HRIR_L(:,idCL),refHRTF.HRIR_L(:,idCL),[],[],8);
%Plot contralateral left-ear HRIR (Az = 270, Col = 90) of conventional
%interpolated dataset and reference.
supdeq_plotIR(interpHRTF_con.HRIR_L(:,idCL),refHRTF.HRIR_L(:,idCL),[],[],8);
%Plot contralateral left-ear HRIR (Az = 270, Col = 90) of MCA
%interpolated dataset and reference.
%MCA interpolated HRTF is much closer to the reference than HRTF from
%conventional interpolation. The "bump" above 10 kHz is corrected through
%the magnitude correction. Interpolation errors (in this case spatial
%aliasing errors at the contralateral ear) are reduced
supdeq_plotIR(interpHRTF_mca.HRIR_L(:,idCL),refHRTF.HRIR_L(:,idCL),[],[],8);
%% (6) - Optional: Calculate log-spectral difference
%Calculate left-ear log-spectral differences in dB.
%Log-spectral difference is often used as a measure for spectral distance
%between a reference and an interpolated HRTF set. Not the same as
%the magnitude error in auditory filters presented in the paper!
lsd_sh = supdeq_calcLSD_HRIR(interpHRTF_sh.HRIR_L,refHRTF.HRIR_L,fs,16);
lsd_con = supdeq_calcLSD_HRIR(interpHRTF_con.HRIR_L,refHRTF.HRIR_L,fs,16);
lsd_mca = supdeq_calcLSD_HRIR(interpHRTF_mca.HRIR_L,refHRTF.HRIR_L,fs,16);
%Quick plot of log-spectral difference averaged over all directions of the
%2702-point Lebedev grid. As shown in the paper, MCA provides the most
%significant benefit in the critical contralateral region. Thus, when
%averaging over all positions, the benefit of MCA seems smaller. The
%analysis in the paper as well as the provided audio examples reveal
%in much more detail the considerable improvements that the proposed
%magnitude correction provides.
AKf(18,9);
semilogx(lsd_sh.f,lsd_sh.lsd_freq,'LineWidth',1.5);
hold on;
semilogx(lsd_con.f,lsd_con.lsd_freq,'LineWidth',1.5);
semilogx(lsd_mca.f,lsd_mca.lsd_freq,'LineWidth',1.5);
xlim([500 20000])
legend('SH W/O MC','SH SUpDEq W/O MC','SH SUpDEq W/ MC','Location','NorthWest');
xlabel('Frequency in Hz');
ylabel('Log-Spectral Difference in dB');
grid on;
%% (7) - Optional: Calculate magnitude error in auditory filters (similar to paper)
[erb_sh,fc_erb] = AKerbError(interpHRTF_sh.HRIR_L,refHRTF.HRIR_L, [50 fs/2], fs);
erb_con = AKerbError(interpHRTF_con.HRIR_L,refHRTF.HRIR_L, [50 fs/2], fs);
erb_mca = AKerbError(interpHRTF_mca.HRIR_L,refHRTF.HRIR_L, [50 fs/2], fs);
%Quick plot of absolute magnitude error in auditory filters averaged over all
%directions of the 2702-point Lebedev grid.
AKf(18,9);
semilogx(fc_erb,mean(abs(erb_sh),2),'LineWidth',1.5);
hold on;
semilogx(fc_erb,mean(abs(erb_con),2),'LineWidth',1.5);
semilogx(fc_erb,mean(abs(erb_mca),2),'LineWidth',1.5);
xlim([500 20000])
legend('SH W/O MC','SH SUpDEq W/O MC','SH SUpDEq W/ MC','Location','NorthWest');
xlabel('Frequency in Hz');
ylabel('\DeltaG(f_c) in dB');
grid on;