% SUpDEq - Spatial Upsampling by Directional Equalization
function [eqHRTFdataset, HRTF_equalized_L, HRTF_equalized_R] = supdeq_eq(sparseHRTFdataset, eqDataset, N, samplingGrid, tikhEps, phaseOnly)
This function performs the equalization of a sparse HRTF dataset (in
SH-domain / stored as SH-coefficients) with the pre-defined eqDataset
(dataset for equalization in SH-domain / stored as SH-coefficients). The
result, the eqHRTFdataset, is the equalized counterpart of the input HRTF
dataset, stored as SH-coefficients. Transform order N of the
eqHRTFdataset is set to the maximum order with respect to the sparse
sampling grid. For example, a lebedev grid with 86 nodes has Nmax = 7.
Output:
eqHRTFdataset - Struct with the equalized HRTF dataset as
SH-coefficients for the left (Hl_nm) and
right (Hr_nm) channel/ear, equalized HRTFs
for the left (HRTF_L) and right (HRTF_R) ear,
absolute frequency scale f, transform order N,
and FFToversize
HRTF_equalized_L/R - Equalized complex HRTFs (single-sided spectrum)
Input:
sparseHRTFdataset - Struct with sparse HRTF dataset in frequency
domain ('normal' HRTFs...). Struct needs to
provide HRTF_L/HRTF_R, the samplingGrid, and Nmax
eqDataset - Struct with equalization dataset as
SH-coefficients. Can be the output of
supdeq_getEqDataset.
N - Transform order N
samplingGrid - Spatial sampling grid (Q x 2 or Q x 3 matrix),
where the first column holds the azimuth,
the second the elevation (both in degree), and
optionally the third the sampling weights.
Azimuth in degree
(0=front, 90=left, 180=back, 270=right)
(0 points to positive x-axis, 90 to positive y-axis)
Elevations in degree
(0=North Pole, 90=front, 180=South Pole)
(0 points to positive z-axis, 180 to negative z-axis)
tikhEps - Define epsilon of Tikhonov regularization if
regularization should be applied
Applying the Tikhonov regularization will always result
in a least-square fit solution for the SH transform.
Variable 'transformCore' is neglected when 'tikhEps' is
defined as the regularized least-square spherical Fourier
transform is applied directly without any third party
toolbox.
Default: 0 (no Tikhonov regularization)
phaseOnly - Set to 1 if only phase response of eqDataset
should be applied for equalization and not
the magnitude response too
Default: 0 (equalize with magnitude and phase)
Dependencies: SOFiA toolbox, AKTools
(C) 2018/2019 by JMA, Johannes M. Arend
CP, Christoph Pörschmann
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 [eqHRTFdataset, HRTF_equalized_L, HRTF_equalized_R] = supdeq_eq(sparseHRTFdataset, eqDataset, N, samplingGrid, tikhEps, phaseOnly) 0004 % 0005 % This function performs the equalization of a sparse HRTF dataset (in 0006 % SH-domain / stored as SH-coefficients) with the pre-defined eqDataset 0007 % (dataset for equalization in SH-domain / stored as SH-coefficients). The 0008 % result, the eqHRTFdataset, is the equalized counterpart of the input HRTF 0009 % dataset, stored as SH-coefficients. Transform order N of the 0010 % eqHRTFdataset is set to the maximum order with respect to the sparse 0011 % sampling grid. For example, a lebedev grid with 86 nodes has Nmax = 7. 0012 % 0013 % Output: 0014 % eqHRTFdataset - Struct with the equalized HRTF dataset as 0015 % SH-coefficients for the left (Hl_nm) and 0016 % right (Hr_nm) channel/ear, equalized HRTFs 0017 % for the left (HRTF_L) and right (HRTF_R) ear, 0018 % absolute frequency scale f, transform order N, 0019 % and FFToversize 0020 % HRTF_equalized_L/R - Equalized complex HRTFs (single-sided spectrum) 0021 % 0022 % Input: 0023 % sparseHRTFdataset - Struct with sparse HRTF dataset in frequency 0024 % domain ('normal' HRTFs...). Struct needs to 0025 % provide HRTF_L/HRTF_R, the samplingGrid, and Nmax 0026 % eqDataset - Struct with equalization dataset as 0027 % SH-coefficients. Can be the output of 0028 % supdeq_getEqDataset. 0029 % N - Transform order N 0030 % samplingGrid - Spatial sampling grid (Q x 2 or Q x 3 matrix), 0031 % where the first column holds the azimuth, 0032 % the second the elevation (both in degree), and 0033 % optionally the third the sampling weights. 0034 % Azimuth in degree 0035 % (0=front, 90=left, 180=back, 270=right) 0036 % (0 points to positive x-axis, 90 to positive y-axis) 0037 % Elevations in degree 0038 % (0=North Pole, 90=front, 180=South Pole) 0039 % (0 points to positive z-axis, 180 to negative z-axis) 0040 % tikhEps - Define epsilon of Tikhonov regularization if 0041 % regularization should be applied 0042 % Applying the Tikhonov regularization will always result 0043 % in a least-square fit solution for the SH transform. 0044 % Variable 'transformCore' is neglected when 'tikhEps' is 0045 % defined as the regularized least-square spherical Fourier 0046 % transform is applied directly without any third party 0047 % toolbox. 0048 % Default: 0 (no Tikhonov regularization) 0049 % phaseOnly - Set to 1 if only phase response of eqDataset 0050 % should be applied for equalization and not 0051 % the magnitude response too 0052 % Default: 0 (equalize with magnitude and phase) 0053 % 0054 % Dependencies: SOFiA toolbox, AKTools 0055 % 0056 % (C) 2018/2019 by JMA, Johannes M. Arend 0057 % CP, Christoph Pörschmann 0058 % TH Köln - University of Applied Sciences 0059 % Institute of Communications Engineering 0060 % Department of Acoustics and Audio Signal Processing 0061 0062 function [eqHRTFdataset, HRTF_equalized_L, HRTF_equalized_R] = supdeq_eq(sparseHRTFdataset, eqDataset, N, samplingGrid, tikhEps, phaseOnly) 0063 0064 if nargin < 5 || isempty(tikhEps) 0065 tikhEps = 0; 0066 end 0067 0068 if nargin < 6 || isempty(phaseOnly) 0069 phaseOnly = 0; 0070 end 0071 0072 %Get fs 0073 fs = sparseHRTFdataset.f(end)*2; 0074 0075 %Transform eqDataset to Fourier domain at sparse sampling grid points 0076 %(inverse spherical Fourier transform) 0077 fprintf('Extracting %d eq transfer functions. This may take some time...\n',size(samplingGrid,1)); 0078 [eqTF_L,eqTF_R] = supdeq_getArbHRTF(eqDataset,samplingGrid,[],[],'ak'); %Use AKtools... 0079 fprintf('%d eq transfer functions extracted...\n',size(samplingGrid,1)) 0080 0081 if phaseOnly 0082 fprintf('Phase only equalization...\n') 0083 %Get only phase response of eqTF 0084 eqTF_L_phase = angle(eqTF_L); 0085 eqTF_R_phase = angle(eqTF_R); 0086 eqTF_L_mag = ones(size(eqTF_L,1),size(eqTF_L,2)); 0087 eqTF_R_mag = ones(size(eqTF_R,1),size(eqTF_R,2)); 0088 eqTF_L = eqTF_L_mag.*exp(1i*eqTF_L_phase); 0089 eqTF_R = eqTF_R_mag.*exp(1i*eqTF_R_phase); 0090 end 0091 0092 %Perform equalization (division in frequency domain) 0093 HRTF_equalized_L = zeros(size(samplingGrid,1),length(eqDataset.f)); 0094 HRTF_equalized_R = zeros(size(samplingGrid,1),length(eqDataset.f)); 0095 for kk = 1:size(samplingGrid,1) 0096 HRTF_equalized_L(kk,:)=sparseHRTFdataset.HRTF_L(kk,:)./eqTF_L(kk,:); 0097 HRTF_equalized_R(kk,:)=sparseHRTFdataset.HRTF_R(kk,:)./eqTF_R(kk,:); 0098 0099 if ~mod(kk,10) 0100 fprintf('%d of %d HRTFs equalized...\n',kk,size(samplingGrid,1)) 0101 end 0102 end 0103 0104 % Transform equalized HRTFs to SH-domain 0105 fprintf('Performing spherical Fourier transform with N = %d. This may take some time...\n',N); 0106 eqHRTFdataset = supdeq_hrtf2sfd(HRTF_equalized_L,HRTF_equalized_R,N,samplingGrid,fs,'ak',tikhEps); %Use AKtools... 0107 eqHRTFdataset.HRTF_L = HRTF_equalized_L; 0108 eqHRTFdataset.HRTF_R = HRTF_equalized_R; 0109 eqHRTFdataset.FFToversize = sparseHRTFdataset.FFToversize; 0110 eqHRTFdataset.eqEarDistance = eqDataset.earDistance; 0111 eqHRTFdataset.eqWaveType = eqDataset.waveType; 0112 if eqDataset.waveType == 1 0113 eqHRTFdataset.sourceDistance = eqDataset.sourceDistance; 0114 end 0115 if phaseOnly 0116 eqHRTFdataset.phaseOnly = 1; 0117 end 0118 if isfield(eqDataset,'limited') 0119 if eqDataset.limited 0120 eqHRTFdataset.limited = true; 0121 eqHRTFdataset.limitInfo = eqDataset.limitInfo; 0122 end 0123 end 0124 fprintf('Done with equalization...\n'); 0125 0126 end 0127