Tidying up

core/applyCovConst.m

@@ -111,14 +111,6 @@
 end
 fprintf('\n')
 
-% Hnm_match2(:,:,1) = Hnm(:,:,1).*G1(:,1) + Hnm(:,:,2).*G2(:,1);
-% Hnm_match2(:,:,2) = Hnm(:,:,2).*G1(:,2) + Hnm(:,:,1).*G2(:,2);
-% 
-% Hnm_match3 = Hnm.*permute(G1,[1,3,2]);
-% Hnm_match3(:,:,1) = Hnm_match3(:,:,1) + Hnm(:,:,2).*G2(:,1);
-% Hnm_match3(:,:,2) = Hnm_match3(:,:,2) + Hnm(:,:,1).*G2(:,2);
-
-
 Hnm = Hnm_match;
 
 %% Plots

core/fromSH.m

@@ -39,7 +39,7 @@
 kr = k*r;
 
 %% Process
-Y = AKsh(N, [], az*180/pi, el*180/pi, 'real').'; % SH coefficients
+Y = getRealSHmatrix(N,az,el); % SH coefficients
 Hnm = ffth(hnm,nfft,1); % to frequency domain
 H = mult3(Hnm,Y); % interpolate
 if isaligned % re-align if specified

core/getYinv.m

@@ -0,0 +1,63 @@
+function Y_inv = getYinv(N,az,el,w,reg_eps)
+% Get transformation matrix to obtain the SH coefficients of an HRTF.
+% Available options:
+%   1. Use provided integration weights (see w)
+%   2. Use pseudoinverse method (least squares solution)
+%   3. Use Tikhonov regularisation which is useful if there are missing
+%       directions, e.g. at the bottom of the sphere [1]
+%
+% SIMPLE USAGE EXAMPLE:
+%   Y_inv = getYinv(15,az,el);
+%
+% INPUT:
+%   N = target SH order
+%   az = HRIR azimuth (ndirs x 1) in rad
+%   el = HRIR elevation (ndirs x 1) in rad (0=top, pi/2=front)
+%   w = quadrature weights (ndirs x 1); if empty, use pseudoinverse
+%   reg_eps = epsilon used for Tikhonov regularisation according to [1]; if
+%       empty, set to 0 (don't regularise)
+%
+% OUTPUT:
+%   Y_inv = transformation matrix (ndirs x (N+1)^2)
+%
+% REFERENCES:
+%   [1] Duraiswaini, R., Dmitry N. Zotkin, and Nail A. Gumerov.
+%       "Interpolation and range extrapolation of HRTFs." 2004 IEEE 
+%       International Conference on Acoustics, Speech, and Signal
+%       Processing. Vol. 4. IEEE, 2004.
+%
+% AUTHOR: Isaac Engel - isaac.engel(at)imperial.ac.uk
+% August 2021
+
+if ~exist('w','var')
+    w = [];
+end
+
+if ~exist('reg_eps','var')
+    reg_eps = 0; % if epsilon not provided, don't regularise
+end
+
+Y = getRealSHmatrix(N,az,el);
+
+if ~isempty(w) && reg_eps == 0
+
+    % If integration weights are provided, use them
+%     Y_inv = 4*pi*w.*Y'; % not compatible with old Matlab versions
+    Y_inv = mult2(4*pi*w,Y');
+
+else
+
+    if reg_eps == 0
+
+        % If no regularisation requested, just apply the pseudoinverse
+        Y_inv = pinv(Y); 
+
+    else
+
+        % If regualarisation requested, apply it according to [12]
+        D = (1 + N*(N+1)) * eye((N+1)^2); % reg. matrix
+        Y_inv = Y'*(Y*Y'+reg_eps*D)^(-1); % regularised weights
+
+    end
+end
+

core/toSH.m

@@ -65,6 +65,8 @@
 %   dualBand = whether to use dual-band tapering (def=false)
 %   Hnm_ref = precomputed high-order Hnm used to calculate diffuse field
 %       EQ. To speed things up if processing many HRTFs
+%   reg_eps = epsilon used for Tikhonov regularisation according to [12]
+%       (def=0=don't regularise)
 %
 % OUTPUTS:
 %   hnm = HRIRs' SH coefficients (nfftOut x (N+1)^2 x 2 ears)
@@ -119,16 +121,14 @@
 %   [11] Archontis Politis, Microphone array processing for parametric
 %       spatial audio techniques, 2016 Doctoral Dissertation, Department of
 %       Signal Processing and Acoustics, Aalto University, Finland.
+%   [12] Duraiswaini, R., Dmitry N. Zotkin, and Nail A. Gumerov.
+%       "Interpolation and range extrapolation of HRTFs." 2004 IEEE 
+%       International Conference on Acoustics, Speech, and Signal
+%       Processing. Vol. 4. IEEE, 2004.
 %
 % AUTHOR: Isaac Engel - isaac.engel(at)imperial.ac.uk
 % February 2021
 
-% TODO: implement Tikhonov regularization for irregular grids
-% According to Duraiswami 2004, it goes like this:
-% reg_eps = 1e-6; % epsilon
-% D = (1 + N*(N+1)) * eye((N+1)^2); % Tikhonov regularization matrix
-% Y_inv = Y' * (Y * Y' + reg_eps*D)^(-1); % this instead of pseudoinverse
-
 %% Parse inputs
 p = inputParser;
 addParameter(p,'mode','Trunc',@ischar) 
@@ -152,6 +152,7 @@
 addParameter(p,'dualBand',false,@isscalar)
 addParameter(p,'tapering',0,@isscalar)
 addParameter(p,'Hnm_ref',[])
+addParameter(p,'reg_eps',0,@isscalar)
 
 parse(p,varargin{:})
 mode = p.Results.mode;
@@ -175,6 +176,7 @@
 dualBand = p.Results.dualBand;
 tapering = p.Results.tapering;
 Hnm_ref = p.Results.Hnm_ref;
+reg_eps = p.Results.reg_eps;
 
 if isempty(fc) && ~contains(mode,'BiMagLS','IgnoreCase',1)
     fc = N*c/(2*pi*r); % if fc not provided, use aliasing frequency
@@ -226,15 +228,9 @@
 alignOption = 1; % 1=phase delay, 2=onset detection, 0=nothing
 if alignOption == 1 % Option 1: phase delay 
     H = ffth(h,nfft,1); % to frequency domain
-    Y0 = AKsh(0, [], az*180/pi, el*180/pi, 'real').';
-    if ~isempty(w)
-%         Y0_inv = 4*pi*w.*Y0'; % if integrations weights are provided, use them
-        Y0_inv = mult2(4*pi*w,Y0'); % use integrations weights if provided
-    else
-        Y0_inv = pinv(Y0); % if not, the pseudoinverse will do just fine
-    end
+    Y0_inv = getYinv(0,az,el,w,reg_eps);
     H0 = mult3(H,Y0_inv); % 0th order SH signal 
-%     pd = -unwrap(angle(H0(2:end,1,:)))./(2*pi*f(2:end)); % phase delay in s
+%     pd = -unwrap(angle(H0(2:end,1,:)))./(2*pi*f(2:end)); % phase delay
 %     pd_mean = mean(pd,'all'); % avg delay
 %     H = H.*exp(1i*2*pi*f*pd_mean); % subtract delay
     pd = div2(-unwrap(angle(H0(2:end,1,:))),(2*pi*f(2:end))); % phase delay
@@ -244,7 +240,6 @@
     onsL = AKonsetDetect(h(:,:,1)); % left ear onsets
     onsR = AKonsetDetect(h(:,:,2)); % right rear onsets
     ons = min([onsL;onsR]); % ipsilateral ear onsets
-    %     nshift = -floor(min(ons))); % detect earliest onset
     nshift = -round(median(ons)); % median onset
     h = circshift(h,nshift);
     H = ffth(h,nfft,1);
@@ -255,30 +250,30 @@
 %% Preprocess using the indicated method
 
 if strcmpi(mode,'Trunc') % just order truncation
-    % Order-truncated SH-HRTF (faster than order Nmax + truncation)
-    Y = AKsh(N,[],az*180/pi,el*180/pi,'real').';
-    if ~isempty(w)
-%         Y_inv = 4*pi*w.*Y'; % use integrations weights if provided
-        Y_inv = mult2(4*pi*w,Y'); % use integrations weights if provided
-    else
-        Y_inv = pinv(Y); % if not, pseudoinverse will do just fine
-    end
+    % Order-truncated SH-HRTF (faster than order Nmax + truncation)  
+    Y_inv = getYinv(N,az,el,w,reg_eps);
     Hnm = mult3(H,Y_inv);
 
 elseif strcmpi(mode,'SpSub') % subsampling aka virtual loudspeaker
-    Hnm = toSH_SpSub(H,N,az,el,w,Nmax);
+    Hnm = toSH_SpSub(H,N,az,el,w,Nmax,reg_eps);
 
 elseif strcmpi(mode,'FDTA') % frequency-dependent time-alignment
     if fc>=fs/2
         warning('fc (%0.2f Hz) >= fs/2 (%0.2f Hz). FDTA preprocessing has no effect...',fc,fs/2)
     end
-    [Hnm,varOut.fc] = toSH_FDTA(H,N,az,el,fs,w,fc,r,earAz,earEl);
+    [Hnm,varOut.fc] = toSH_FDTA(H,N,az,el,fs,w,fc,r,earAz,earEl,reg_eps);
 
 elseif strcmpi(mode,'MagLS') % magnitude least squares
     if fc>=fs/2
         warning('fc (%0.2f Hz) >= fs/2 (%0.2f Hz). MagLS preprocessing has no effect...',fc,fs/2)
     end
-    [Hnm,varOut.fc] = toSH_MagLS(H,N,az,el,fs,w,fc,k,r);
+    [Hnm,varOut.fc] = toSH_MagLS(H,N,az,el,fs,w,fc,k,r,reg_eps);
+
+elseif strcmpi(mode,'MagSun') % magnitude optimisation by Sun
+    if fc>=fs/2
+        warning('fc (%0.2f Hz) >= fs/2 (%0.2f Hz). MagSun preprocessing has no effect...',fc,fs/2)
+    end
+    [Hnm,varOut.fc] = toSH_MagSun(H,N,az,el,fs,w,fc,k,r,reg_eps);
 
 elseif strcmpi(mode,'MagLSold') % old smoothing, kept to reproduce old data
     if fc>=fs/2
@@ -287,13 +282,13 @@
     [Hnm,varOut.fc] = toSH_MagLS_old(H,N,az,el,fs,w,fc,k,r);
 
 elseif strcmpi(mode,'TA') % time alignment (ear alignment)
-    Hnm = toSH_TA(H,N,az,el,fs,w,r,earAz,earEl);
+    Hnm = toSH_TA(H,N,az,el,fs,w,r,earAz,earEl,reg_eps);
 
 elseif strcmpi(mode,'BiMagLS') % ear alignment + MagLS
     if fc>=fs/2
         warning('fc (%0.2f Hz) >= fs/2 (%0.2f Hz). BiMagLS preprocessing has no effect...',fc,fs/2)
     end
-    [Hnm,varOut.fc] = toSH_BiMagLS(H,N,az,el,fs,w,fc,k,r,earAz,earEl);
+    [Hnm,varOut.fc] = toSH_BiMagLS(H,N,az,el,fs,w,fc,k,r,earAz,earEl,reg_eps);
 
 elseif strcmpi(mode,'BiMagLSold') % old smoothing, kept to reproduce old data
     if fc>=fs/2
@@ -305,7 +300,7 @@
     if fc>=fs/2
         warning('fc (%0.2f Hz) >= fs/2 (%0.2f Hz). SpSubMod preprocessing is the same as SpSub...',fc,fs/2)
     end
-    [Hnm,varOut.fc] = toSH_SpSubMod(H,N,az,el,fs,w,fc,r,earAz,earEl,Nmax);
+    [Hnm,varOut.fc] = toSH_SpSubMod(H,N,az,el,fs,w,fc,r,earAz,earEl,Nmax,reg_eps);
 
 else
     error('Unknown method %s. Use one of these: ''Trunc'', ''SpSub'', ''FDTA'', ''MagLS'', ''TA'', ''BiMagLS'', ''SpSubMod''',mode)
@@ -364,7 +359,7 @@
 
     if EQ==1 || EQ==2 % these modes require a high-order reference
         if isempty(Hnm_ref)
-            Y = AKsh(Nmax,[],az*180/pi,el*180/pi,'real').';
+            Y = getRealSHmatrix(Nmax,az,el);
             if contains(mode,'ear') % time-align reference if required
                 kr = 2*pi*f*r/c;
                 p = earAlign(kr,az,el,earAz,earEl);

core/toSH_BiMagLS.m

@@ -1,4 +1,4 @@
-function [Hnm,fc,p] = toSH_BiMagLS(H,N,az,el,fs,w,fc,k,r,earAz,earEl)
+function [Hnm,fc,p] = toSH_BiMagLS(H,N,az,el,fs,w,fc,k,r,earAz,earEl,reg_eps)
 % Transform HRTF to SH domain at order N by using ear aligning first [1]
 % and frequency-dependent optimisation above a certain cutoff frequency
 % (MagLS [2,3]). The idea is to align the HRIRs first to reduce the
@@ -83,5 +83,5 @@
 H = H.*exp(-1i*p); % apply correction
 
 %% Then, apply MagLS
-[Hnm,fc] = toSH_MagLS(H,N,az,el,fs,w,fc,k,r);
+[Hnm,fc] = toSH_MagLS(H,N,az,el,fs,w,fc,k,r,reg_eps);
 

core/toSH_FDTA.m

@@ -1,4 +1,4 @@
-function [Hnm,fc] = toSH_FDTA(H,N,az,el,fs,w,fc,r,earAz,earEl)
+function [Hnm,fc] = toSH_FDTA(H,N,az,el,fs,w,fc,r,earAz,earEl,reg_eps)
 % Transform HRTF to SH domain at order N using frequency-dependent 
 % time-alignment [1]. Alignment is performed with the method in [2].
 %
@@ -60,12 +60,6 @@
 H(ind,:,:) = H(ind,:,:).*exp(-1i*p(ind,:,:)); % apply correction
 
 %% Then, get the order-truncated SH-HRTF
-Y = AKsh(N,[],az*180/pi,el*180/pi,'real').';
-if exist('w','var') && ~isempty(w)
-%     Y_inv = 4*pi*w.*Y'; % if integrations weights are provided, use them
-    Y_inv = mult2(4*pi*w,Y'); % use integrations weights if provided
-else
-    Y_inv = pinv(Y); % if not, the pseudoinverse will do just fine
-end
+Y_inv = getYinv(N,az,el,w,reg_eps);
 Hnm = mult3(H,Y_inv);
 

core/toSH_MagLS.m

@@ -1,4 +1,4 @@
-function [Hnm,fc] = toSH_MagLS(H,N,az,el,fs,w,fc,k,r)
+function [Hnm,fc] = toSH_MagLS(H,N,az,el,fs,w,fc,k,r,reg_eps)
 % Transform HRTF to SH domain at order N with the MagLS method [1],
 % following the simplified approach from [2] sec. 4.11.2 and 4.11.3. The
 % main difference of this implementation vs [1] is that the phase is
@@ -56,13 +56,8 @@
 f = linspace(0,fs/2,nfreqs).'; % frequency vector
 
 %% Get the SH matrices
-Y = AKsh(N, [], az*180/pi, el*180/pi, 'real').'; 
-if exist('w','var') && ~isempty(w)
-%     Y_inv = 4*pi*w.*Y'; % if integrations weights are provided, use them
-    Y_inv = mult2(4*pi*w,Y'); % use integrations weights if provided
-else
-    Y_inv = pinv(Y); % if not, the pseudoinverse will do just fine
-end
+Y = getRealSHmatrix(N,az,el);
+Y_inv = getYinv(N,az,el,w,reg_eps);
 
 %% Apply simplified MagLS from [2] sec.4.11.2
 

core/toSH_MagSun.m

@@ -0,0 +1,57 @@
+function [Hnm,fc] = toSH_MagSun(H,N,az,el,fs,w,fc,k,r,reg_eps)
+% Transform HRTF to SH domain at order N with the Magnitude-optimised
+% reconstruction method from Sun's thesis [1].
+%
+% SIMPLE USAGE EXAMPLE:
+%   Hnm = toSH_MagSun(H,15,az,el,48000);
+%
+% INPUT:
+%   H = HRTF up to Nyquist frequency (nfreqs x ndirs x 2 ears)
+%   N = target SH order
+%   az = HRIR azimuth (ndirs x 1) in rad
+%   el = HRIR elevation (ndirs x 1) in rad (0=top, pi/2=front)
+%   fs = sampling frequency in Hz
+%   w = quadrature weights (ndirs x 1); if empty, use pseudoinverse
+%   fc = cutoff frequency in Hz above which phase is "disregarded" in
+%       favour of magnitude; if empty (default), use aliasing frequency
+%   k = length of the transition band in octaves (def=1). E.g.
+%       if fc=623.89 Hz, smooth between 349.65 and 882.31 Hz. If k==0, 
+%       don't smooth.
+%   r = head radius in m (def=0.0875)
+%   
+% OUTPUT:
+%   Hnm = HRTF's SH coefficients (nfreqs x (N+1)^2 x 2 ears)
+%   fc = see above
+%
+% REFERENCES:
+%   [1] Sun, David. "Generation and perception of three-dimensional sound
+%       fields using Higher Order Ambisonics." (2013).
+%
+% AUTHOR: Isaac Engel - isaac.engel(at)imperial.ac.uk
+% September 2021
+
+%% Some parameters
+if ~exist('r','var') || isempty(r)
+    r = 0.0875; % default head radius
+end
+if ~exist('fc','var') || isempty(fc)
+    c = 343; % speed of sound (m/s)
+    fc = N*c/(2*pi*r); % if fc not provided, use aliasing frequency
+end
+if ~exist('k','var') || isempty(k)
+    k = 1; % default smoothing = one octave (half to each side)
+end
+nfreqs = size(H,1);
+f = linspace(0,fs/2,nfreqs).'; % frequency vector
+
+%% Get the SH matrices
+
+Hmag = abs(H);
+Hphase = unwrap(angle(H));
+HphaseAvg = mean(Hphase,2); %HphaseAvg = mult3(Hphase,ones(ndirs))/ndirs;
+alpha = 0.5 + k/log(2) * log(f/fc); % 0 for fc1, 1 for fc2
+alpha = min(max(alpha,0),1); % keep between 0 and 1
+HphaseMod = alpha.*HphaseAvg + (1-alpha).*Hphase;
+Htarg = Hmag.*exp(1i*HphaseMod);
+Y_inv = getYinv(N,az,el,w,reg_eps);
+Hnm = mult3(Htarg,Y_inv); 

core/toSH_SpSub.m

@@ -1,4 +1,4 @@
-function Hnm = toSH_SpSub(H,N,az,el,w,Nmax)
+function Hnm = toSH_SpSub(H,N,az,el,w,Nmax,reg_eps)
 % Transform HRTF to SH domain at order N using spatial subsampling [1].
 % This is equivalent to virtual loudspeaker decoding (first introduced by 
 % [2]) with "mode-matching" (pseudoinverse-based).
@@ -35,13 +35,7 @@
 
 %% Process
 % First, get the SH-HRTF with highest available order
-Ymax = AKsh(Nmax, [], az*180/pi, el*180/pi, 'real').'; % SH coeffs
-if exist('w','var') && ~isempty(w)
-%     Ymax_inv = 4*pi*w.*Ymax'; % if integrations weights are provided, use them
-    Ymax_inv = mult2(4*pi*w,Ymax'); % use integrations weights if provided
-else
-    Ymax_inv = pinv(Ymax); % if not, the pseudoinverse will do just fine
-end
+Ymax_inv = getYinv(Nmax,az,el,w,reg_eps);
 Hnm_max = mult3(H,Ymax_inv);
 
 % Then, resample to an appropriately sized sparse grid (Gauss)
@@ -50,8 +44,7 @@
 Hs = mult3(Hnm_max,Ymax_s); % = cat(3,Hnm_max(:,:,1)*Ymax_s,Hnm_max(:,:,2)*Ymax_s);
 
 % Finally, convert to SH again at the target low order
-Ys = AKsh(N, [], sgrid(:,1)*180/pi, sgrid(:,2)*180/pi, 'real').';
-% Ys_inv = 4*pi*diag(sgrid(:,3))*Ys'; 
-Ys_inv = mult2(4*pi*sgrid(:,3),Ys');
+Ys_inv = getYinv(N,sgrid(:,1),sgrid(:,2));
 Hnm = mult3(Hs,Ys_inv);
+
 

core/toSH_SpSubMod.m

@@ -1,4 +1,4 @@
-function [Hnm,fc] = toSH_SpSubMod(H,N,az,el,fs,w,fc,r,earAz,earEl,Nmax)
+function [Hnm,fc] = toSH_SpSubMod(H,N,az,el,fs,w,fc,r,earAz,earEl,Nmax,reg_eps)
 % Combination of spatial subsampling [1] and frequency-dependent time 
 % alignment [2], as suggested in [3].
 %
@@ -71,4 +71,4 @@
 H(ind,:,:) = H(ind,:,:).*exp(-1i*p(ind,:,:)); % apply correction
 
 %% Then, apply spatial subsampling
-Hnm = toSH_SpSub(H,N,az,el,w,Nmax);
+Hnm = toSH_SpSub(H,N,az,el,w,Nmax,reg_eps);