New linear interpolation method for HRTFs working in the frequency do…

…main

SFS_config.m

@@ -337,6 +337,22 @@
 % will be done between the two or three nearest HRTFs.
 conf.ir.useinterpolation = true; % boolean
 %
+% You can choose between the following interpolation methods:
+%   'simple'      - Interpolation in the time domain performed samplewise. This
+%                   does not heed the times of arrival of the impulse responses.
+%   'freqdomain'  - Interpolation in the frequency domain performed separately
+%                   for magnitude and phase.
+%                   This method cannot work properly if there is too much noise in
+%                   the phase information at low frequencies which is often the
+%                   case for measured HRTFs. Low frequencies can be corrected
+%                   according to theory, see e.g. the corrected KEMAR HRTFs published
+%                   at http://github.com/spatialaudio/lf-corrected-kemar-hrtfs.
+%                   The implementation of this method suffers from circular shifting,
+%                   see test_interpolation_methods.m in the validation folder. For
+%                   typical HRIRs with leading and trailing zeros, the error is
+%                   negligible.
+conf.ir.interpolationmethod = 'simple';
+%
 % If you have HRIRs in the form of the SimpleFreeFieldHRIR SOFA convention, zeros
 % are padded at the beginning of every impulse response corresponding to their
 % measurement distance. If you know that your measured HRIRs already have a

SFS_general/interpolation.m

@@ -1,5 +1,5 @@
 function B = interpolation(A,X,x)
-%INTPOLATION interpolates the given data A(X) at the point x
+%INTERPOLATION interpolates the given data A(X) at the point x
 %
 %   Usage: B = interpolation(A,X,x)
 %

SFS_general/secondary_source_diameter.m

@@ -84,7 +84,7 @@
     end
     % Find source1 :=  source with largest distance from origin
     [~,idx1] = max(vector_norm(x0(:,1:3),2));
-    % Find source2 := source with maximum distace to source1
+    % Find source2 := source with maximum distance to source1
     [diam,idx2] = max(vector_norm(x0(:,1:3) - ...
         repmat(x0(idx1,1:3),[size(x0,1),1]),2));
     % Center is half-way between source1 and source2

SFS_general/secondary_source_positions.m

@@ -213,7 +213,7 @@
     [~,theta] = cart2sph(x0(:,1),x0(:,2),x0(:,3)); % get elevation
     x0(:,7) = x0(:,7) .* cos(theta);
 elseif strcmp('custom',geometry)
-    % Custom geometry definedy by conf.secondary_sources.x0.
+    % Custom geometry defined by conf.secondary_sources.x0.
     % This could be in the form of a n x 7 matrix, where n is the number of
     % secondary sources or as a SOFA file/struct.
     if ischar(conf.secondary_sources.x0) || isstruct(conf.secondary_sources.x0)

SFS_ir/interpolate_ir.m

@@ -1,7 +1,7 @@
 function [ir_new,x0_new] = interpolate_ir(ir,x0,xs,conf)
 %INTERPOLATE_IR interpolates three given IRs for the given angle
 %
-%   Usage: [ir,x0] = interpolate_ir(ir,x0,xs,conf)
+%   Usage: [ir_new,x0_new] = interpolate_ir(ir,x0,xs,conf)
 %
 %   Input parameters:
 %       ir      - matrix containing impulse responses in the form [M C N], where
@@ -14,16 +14,30 @@
 %       conf    - configuration struct (see SFS_config)
 %
 %   Output parameters:
-%       ir      - impulse response for the given position [1 C N]
-%       x0      - position corresponding to the returned impulse response
+%       ir_new  - impulse response for the given position [1 C N]
+%       x0_new  - position corresponding to the returned impulse response
 %
 %   INTERPOLATE_IR(ir,x0,xs,conf) interpolates the two to three given impulse
 %   responses from ir with their corresponding angles x0 for the given angles
 %   xs and returns an interpolated impulse response.
-%   Note, that the given parameter are not checked if they have all the correct
+%	For the 1D case, the interpolation method differs depending on the setting
+%	of conf.ir.interpolationmethod:
+%     'simple'      - Interpolation in the time domain performed samplewise.
+%	                  This does not heed the times of arrival of the impulse
+%                     responses.
+%     'freqdomain'  - Interpolation in the frequency domain performed separately
+%                     for magnitude and phase.
+%   Note, that the given parameters are not checked if they have all the correct
 %   dimensions in order to save computational time, because this function could
 %   be called quite often.
 %
+%	References:
+%		K. Hartung, J. Braasch, S. J. Sterbing (1999) - "Comparison of different
+%		methods for the interpolation of head-related transfer functions".
+%		Proc. of the 16th AES Conf.
+%		K. Itoh (1982) - "Analysis of the phase unwrapping algorithm". Applied
+%		Optics 21(14), 2470
+%
 %   See also: get_ir, interpolation
 
 %*****************************************************************************
@@ -64,6 +78,14 @@
 
 %% ===== Configuration ==================================================
 useinterpolation = conf.ir.useinterpolation;
+% Check for old configuration
+if useinterpolation && ~isfield(conf.ir,'interpolationmethod')
+    warning('SFS:irs_intpolmethod',...
+        'no interpolation method provided, will use method ''simple''.');
+    interpolationmethod = 'simple';
+else
+    interpolationmethod = conf.ir.interpolationmethod;
+end
 % Precision of the wanted angle. If an impulse response within the given
 % precision could be found no interpolation is applied.
 prec = 0.001; % ~ 0.05 deg
@@ -88,8 +110,38 @@
              'and (%.1f,%.1f) deg.'], ...
             deg(x0(1,1)), deg(x0(2,1)), ...
             deg(x0(1,2)), deg(x0(2,2)));
-        ir_new(1,1,:) = interpolation(squeeze(ir(1:2,1,:))',x0(:,1:2),xs);
-        ir_new(1,2,:) = interpolation(squeeze(ir(1:2,2,:))',x0(:,1:2),xs);
+        switch interpolationmethod
+        case 'simple'
+            ir_new(1,1,:) = interpolation(squeeze(ir(1:2,1,:))',x0(:,1:2),xs);
+            ir_new(1,2,:) = interpolation(squeeze(ir(1:2,2,:))',x0(:,1:2),xs);
+        case 'freqdomain'
+            % see Itoh (1982), Hartung et al. (1999)
+            %
+            % Upsample to avoid phase aliasing in unwrapping of phase
+            TF = fft(ir,4*size(ir,3),3);
+            % Magnitude and phase will be interpolated separately
+            magnitude = abs(TF);
+            phase = unwrap(angle(TF),[],3);
+            % Calculate interpolation only for the first half of the spectrum
+            % and only for original bins
+            idx_half = floor(size(TF,3)/2)+1;
+            magnitude_new(1,:) = interpolation(...
+                squeeze(magnitude(1:2,1,1:4:idx_half))',x0(:,1:2),xs);
+            magnitude_new(2,:) = interpolation(...
+                squeeze(magnitude(1:2,2,1:4:idx_half))',x0(:,1:2),xs);
+            phase_new(1,:) = interpolation(...
+                squeeze(phase(1:2,1,1:4:idx_half))',x0(:,1:2),xs);
+            phase_new(2,:) = interpolation(...
+                squeeze(phase(1:2,2,1:4:idx_half))',x0(:,1:2),xs);
+            % Calculate interpolated impulse response from magnitude and phase
+            ir_new(1,1,:) = ifft(magnitude_new(1,:) ...
+                .* exp(1i*phase_new(1,:)),size(ir,3),'symmetric');
+            ir_new(1,2,:) = ifft(magnitude_new(2,:)...
+                .* exp(1i*phase_new(2,:)),size(ir,3),'symmetric');
+        otherwise
+            error('%s: %s is an unknown interpolation method.', ...
+                upper(mfilename),interpolationmethod);
+        end
     else
         % --- 2D interpolation ---
         warning('SFS:irs_intpol3D',...

SFS_ssr/ssr_brs_wfs.m

@@ -11,7 +11,7 @@
 %       src     - source type: 'pw' - plane wave
 %                              'ps' - point source
 %                              'fs' - focused source
-%       irs     - IR data set for the second sources
+%       irs     - IR data set for the secondary sources
 %       conf    - configuration struct (see SFS_config)
 %
 %   Output parameters: