From 7916efa7a4990899960bd46e958be4b374dbc59e Mon Sep 17 00:00:00 2001 From: Vera Erbes Date: Tue, 9 Aug 2016 12:49:58 +0200 Subject: [PATCH] New linear interpolation method for HRTFs working in the frequency domain --- SFS_config.m | 16 ++ SFS_general/interpolation.m | 2 +- SFS_general/secondary_source_diameter.m | 2 +- SFS_general/secondary_source_positions.m | 2 +- SFS_ir/interpolate_ir.m | 64 +++++- SFS_ssr/ssr_brs_wfs.m | 2 +- validation/test_interpolation_methods.m | 251 +++++++++++++++++++++++ 7 files changed, 329 insertions(+), 10 deletions(-) create mode 100644 validation/test_interpolation_methods.m diff --git a/SFS_config.m b/SFS_config.m index 6eee15c8..2de2abc5 100644 --- a/SFS_config.m +++ b/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 diff --git a/SFS_general/interpolation.m b/SFS_general/interpolation.m index c5d78fbd..8d55a6d9 100644 --- a/SFS_general/interpolation.m +++ b/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) % diff --git a/SFS_general/secondary_source_diameter.m b/SFS_general/secondary_source_diameter.m index 69e2560a..96d3416a 100644 --- a/SFS_general/secondary_source_diameter.m +++ b/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 diff --git a/SFS_general/secondary_source_positions.m b/SFS_general/secondary_source_positions.m index b2f246f1..572b97c2 100644 --- a/SFS_general/secondary_source_positions.m +++ b/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) diff --git a/SFS_ir/interpolate_ir.m b/SFS_ir/interpolate_ir.m index 8769f2f2..c6132dd4 100644 --- a/SFS_ir/interpolate_ir.m +++ b/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',... diff --git a/SFS_ssr/ssr_brs_wfs.m b/SFS_ssr/ssr_brs_wfs.m index a586cae7..b2688b4b 100644 --- a/SFS_ssr/ssr_brs_wfs.m +++ b/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: diff --git a/validation/test_interpolation_methods.m b/validation/test_interpolation_methods.m new file mode 100644 index 00000000..cc29ea4e --- /dev/null +++ b/validation/test_interpolation_methods.m @@ -0,0 +1,251 @@ +function status = test_interpolation_methods(modus) +%TEST_INTERPOLATION_METHODS demonstrates the difference between interpolation +%methods +% +% Usage: status = test_interpolation_methods(modus) +% +% Input parameters: +% modus - 0: numerical (not available) +% 1: visual +% +% Output parameters: +% status - true or false +% +% TEST_INTERPOLATION_METHODS(modus) demonstrates the different interpolation +% methods in the interpolate_ir function for shifted Dirac impulses and for +% HRIRs of different directions + +%***************************************************************************** +% The MIT License (MIT) * +% * +% Copyright (c) 2010-2016 SFS Toolbox Developers * +% * +% Permission is hereby granted, free of charge, to any person obtaining a * +% copy of this software and associated documentation files (the "Software"), * +% to deal in the Software without restriction, including without limitation * +% the rights to use, copy, modify, merge, publish, distribute, sublicense, * +% and/or sell copies of the Software, and to permit persons to whom the * +% Software is furnished to do so, subject to the following conditions: * +% * +% The above copyright notice and this permission notice shall be included in * +% all copies or substantial portions of the Software. * +% * +% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR * +% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * +% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL * +% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER * +% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING * +% FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER * +% DEALINGS IN THE SOFTWARE. * +% * +% The SFS Toolbox allows to simulate and investigate sound field synthesis * +% methods like wave field synthesis or higher order ambisonics. * +% * +% http://sfstoolbox.org sfstoolbox@gmail.com * +%***************************************************************************** + + +status = false; + + +%% ===== Checking of input parameters ==================================== +nargmin = 1; +nargmax = 1; +narginchk(nargmin,nargmax); +if ~modus + warning('%s: numerical modus not available.',upper(mfilename)); + status = true; + return; +end + + +%% ===== Settings ======================================================== +conf = struct; +conf.ir.useinterpolation = true; + + +%% ===== Main ============================================================ +% Check if HRTF data set with low frequency correction is available, +% download otherwise +basepath = get_sfs_path(); +hrtf = 'KEMAR_HRTFs_lfcorr.sofa'; +hrtf_file = [basepath '/data/HRTFs/' hrtf]; +if ~exist(hrtf_file,'file') + disp('Download'); + url = ['https://github.com/spatialaudio/lf-corrected-kemar-hrtfs/', ... + 'blob/master/' hrtf '?raw=true']; + download_file(url,hrtf_file); +end +% Load HRTF data set +hrtf = SOFAload(hrtf_file); + + +%% ===== Interpolate shifted Dirac impulses ============================== +% interpolation points +x0 = [1 3; 0 0; 0 0]/180*pi; +% target point +xs = [2; 0; 0]/180*pi; +N = 50; %length impulse responses + +% 1. Interpolate between Dirac impulses with one sample in between +h1 = zeros(1,N); +h1(3) = 1; +h2 = zeros(1,N); +h2(5) = 1; + +irs1(1,1,:) = h1; +irs1(2,1,:) = h2; +irs1(:,2,:) = irs1(:,1,:); %redundant second channel + +conf.ir.interpolationmethod = 'simple'; +h_int1_simple = interpolate_ir(irs1,x0,xs,conf); +conf.ir.interpolationmethod = 'freqdomain'; +h_int1_fd = interpolate_ir(irs1,x0,xs,conf); + +% 2. Interpolate between neighbouring Dirac impulses at impulse reponse start +h1 = zeros(1,N); +h1(3) = 1; +h3 = zeros(1,N); +h3(4) = 1; + +irs2(1,1,:) = h1; +irs2(2,1,:) = h3; +irs2(:,2,:) = irs2(:,1,:); %redundant second channel + +conf.ir.interpolationmethod = 'simple'; +h_int2_simple = interpolate_ir(irs2,x0,xs,conf); +conf.ir.interpolationmethod = 'freqdomain'; +h_int2_fd = interpolate_ir(irs2,x0,xs,conf); + +% 3. Interpolate between neighbouring Dirac impulses in middle of impulse response +h4 = zeros(1,N); +h4(25) = 1; +h5 = zeros(1,N); +h5(26) = 1; + +irs3(1,1,:) = h4; +irs3(2,1,:) = h5; +irs3(:,2,:) = irs3(:,1,:); + +conf.ir.interpolationmethod = 'simple'; +h_int3_simple = interpolate_ir(irs3,x0,xs,conf); +conf.ir.interpolationmethod = 'freqdomain'; +h_int3_fd = interpolate_ir(irs3,x0,xs,conf); + +% Plots +% impulse responses +figure + plot(0:N-1,squeeze(irs1(1,1,:)),'k'), hold on + plot(0:N-1,squeeze(irs1(2,1,:)),'b') + plot(0:N-1,squeeze(h_int1_simple(1,1,:)),'r') + plot(0:N-1,squeeze(h_int1_fd(1,1,:)),'m') + grid + xlabel('samples'), ylabel('amplitude') + legend('h_1','h_2','simple interp','freqdomain interp') + title('Interpolation between Dirac impulses with one sample in between') + +figure + plot(0:N-1,squeeze(irs2(1,1,:)),'k'), hold on + plot(0:N-1,squeeze(irs2(2,1,:)),'b') + plot(0:N-1,squeeze(h_int2_simple(1,1,:)),'r') + plot(0:N-1,squeeze(h_int2_fd(1,1,:)),'m') + grid + xlabel('samples'), ylabel('amplitude') + legend('h_1','h_2','simple interp','freqdomain interp') + title('Interpolation of neighbouring Dirac impulses at impulse response start') + +figure + plot(0:N-1,squeeze(irs3(1,1,:)),'k'), hold on + plot(0:N-1,squeeze(irs3(2,1,:)),'b') + plot(0:N-1,squeeze(h_int3_simple(1,1,:)),'r') + plot(0:N-1,squeeze(h_int3_fd(1,1,:)),'m') + grid + xlabel('samples'), ylabel('amplitude') + legend('h_1','h_2','simple interp','freqdomain interp') + title('Interpolation of neighbouring Dirac impulses in middle of impulse response') + + +%% ===== Interpolate HRIRs =============================================== +% 1. Interpolate between close HRIRs +idx0 = 181; %index for 0° azimuth +idx1 = 182; %index for 1° azimuth +idx2 = 183; %index for 2° azimuth +x0_close = [hrtf.SourcePosition(idx0,:)' hrtf.SourcePosition(idx2,:)']/180*pi; +xs_close = hrtf.SourcePosition(idx1,:)'/180*pi; +hrir_close = [hrtf.Data.IR(idx0,:,:); hrtf.Data.IR(idx2,:,:)]; +hrir_close_ref = hrtf.Data.IR(idx1,:,:); + +conf.ir.interpolationmethod = 'simple'; +hrir_close_simple = interpolate_ir(hrir_close,x0_close,xs_close,conf); +conf.ir.interpolationmethod = 'freqdomain'; +hrir_close_fd = interpolate_ir(hrir_close,x0_close,xs_close,conf); + +% 2. Interpolate between distant HRIRs +idx0 = 181; %index for 0° azimuth +idx30 = 211; %index for 30° azimuth +idx60 = 241; %index for 60° azimuth +x0_dist = [hrtf.SourcePosition(idx0,:)' hrtf.SourcePosition(idx60,:)']/180*pi; +xs_dist = hrtf.SourcePosition(idx30,:)'/180*pi; +hrir_dist = [hrtf.Data.IR(idx0,:,:); hrtf.Data.IR(idx60,:,:)]; +hrir_dist_ref = hrtf.Data.IR(idx30,:,:); + +conf.ir.interpolationmethod = 'simple'; +hrir_dist_simple = interpolate_ir(hrir_dist,x0_dist,xs_dist,conf); +conf.ir.interpolationmethod = 'freqdomain'; +hrir_dist_fd = interpolate_ir(hrir_dist,x0_dist,xs_dist,conf); + +% Plots +% impulse responses +figure + plot(0:hrtf.API.N-1,squeeze(hrir_close(1,1,:)),'k'), hold on + plot(0:hrtf.API.N-1,squeeze(hrir_close(2,1,:)),'b') + plot(0:hrtf.API.N-1,squeeze(hrir_close_ref(1,1,:)),'g') + plot(0:hrtf.API.N-1,squeeze(hrir_close_simple(1,1,:)),'r') + plot(0:hrtf.API.N-1,squeeze(hrir_close_fd(1,1,:)),'m') + grid + xlabel('samples'), ylabel('amplitude') + axis([0 160 -0.6 0.6]) + legend('hrir_1','hrir_2','hrir_{ref}','simple interp','freqdomain interp') + title('Interpolation of close HRIRs') + +figure + plot(0:hrtf.API.N-1,squeeze(hrir_dist(1,1,:)),'k'), hold on + plot(0:hrtf.API.N-1,squeeze(hrir_dist(2,1,:)),'b') + plot(0:hrtf.API.N-1,squeeze(hrir_dist_ref(1,1,:)),'g') + plot(0:hrtf.API.N-1,squeeze(hrir_dist_simple(1,1,:)),'r') + plot(0:hrtf.API.N-1,squeeze(hrir_dist_fd(1,1,:)),'m') + grid + xlabel('samples'), ylabel('amplitude') + axis([0 160 -0.6 0.6]) + legend('hrir_1','hrir_2','hrir_{ref}','simple interp','freqdomain interp') + title('Interpolation of distant HRIRs') + +% magnitude responses +f = (0:hrtf.API.N-1)/hrtf.API.N*hrtf.Data.SamplingRate; +figure + semilogx(f,db(abs(fft(squeeze(hrir_close(1,1,:))))),'k'), hold on + semilogx(f,db(abs(fft(squeeze(hrir_close(2,1,:))))),'b') + semilogx(f,db(abs(fft(squeeze(hrir_close_ref(1,1,:))))),'g') + semilogx(f,db(abs(fft(squeeze(hrir_close_simple(1,1,:))))),'r') + semilogx(f,db(abs(fft(squeeze(hrir_close_fd(1,1,:))))),'m') + grid + xlabel('frequency in Hz'), ylabel('amplitude in dB') + axis([0 hrtf.Data.SamplingRate/2 -60 20]) + legend('hrir_1','hrir_2','hrir_{ref}','simple interp','freqdomain interp',... + 'Location','SW') + title('Interpolation of close HRIRs') + +figure + semilogx(f,db(abs(fft(squeeze(hrir_dist(1,1,:))))),'k'), hold on + semilogx(f,db(abs(fft(squeeze(hrir_dist(2,1,:))))),'b') + semilogx(f,db(abs(fft(squeeze(hrir_dist_ref(1,1,:))))),'g') + semilogx(f,db(abs(fft(squeeze(hrir_dist_simple(1,1,:))))),'r') + semilogx(f,db(abs(fft(squeeze(hrir_dist_fd(1,1,:))))),'m') + grid + xlabel('frequency in Hz'), ylabel('amplitude in dB') + axis([0 hrtf.Data.SamplingRate/2 -60 20]) + legend('hrir_1','hrir_2','hrir_{ref}','simple interp','freqdomain interp',... + 'Location','SW') + title('Interpolation of distant HRIRs') + +status = true;