From 1f36d96686ae187ea096aefdc55d080f072b3cf2 Mon Sep 17 00:00:00 2001 From: VeraE Date: Thu, 23 Jun 2016 17:42:24 +0200 Subject: [PATCH 01/17] add conf field for interpolation method --- SFS_config.m | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/SFS_config.m b/SFS_config.m index 583eeaff..5f15e88f 100644 --- a/SFS_config.m +++ b/SFS_config.m @@ -331,7 +331,14 @@ % disabled the HRTF returned by a nearest neighbour search is used instead. % Depending on the geometry of the measured HRTF data set, the interpolation % will be done between the two or three nearest HRTFs. -conf.ir.useinterpolation = true; % boolean +conf.ir.useinterpolation = false; % boolean +% +% You can choose between the following interpolation methods: +% 'simple' - Interpolation in the time domain performed samplewise. This +% does not heed the time of arrival of the HRTFs. +% 'freqdomain' - Interpolation in the frequency domain performed separately +% for magnitude and phase of the HRTF. +conf.ir.interpolationmethod = 'freqdomain'; % % 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 From 845d1a6a0ffa103e635fb969831dd26407f016c1 Mon Sep 17 00:00:00 2001 From: VeraE Date: Fri, 24 Jun 2016 14:05:57 +0200 Subject: [PATCH 02/17] add interpolation of magnituge and phase in frequency domain --- SFS_ir/interpolate_ir.m | 37 +++++++++++++++++++++++++++++++++++-- 1 file changed, 35 insertions(+), 2 deletions(-) diff --git a/SFS_ir/interpolate_ir.m b/SFS_ir/interpolate_ir.m index 8769f2f2..a4829490 100644 --- a/SFS_ir/interpolate_ir.m +++ b/SFS_ir/interpolate_ir.m @@ -64,6 +64,14 @@ %% ===== Configuration ================================================== useinterpolation = conf.ir.useinterpolation; +% Check for old configuration +if ~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 +96,33 @@ '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' + TF = fft(ir,[],3); + % Calculate interpolation for magnitude and phase separately + mag = abs(TF); + pha = unwrap(angle(TF),[],3); + % Calculate interpolation only for the first half of the spectrum + idx_half = floor(size(mag,3)/2)+1; + mag_new(1,1,:) = interpolation(squeeze(mag(1:2,1,1:idx_half))',... + x0(:,1:2),xs); + mag_new(1,2,:) = interpolation(squeeze(mag(1:2,2,1:idx_half))',... + x0(:,1:2),xs); + pha_new(1,1,:) = interpolation(squeeze(pha(1:2,1,1:idx_half))',... + x0(:,1:2),xs); + pha_new(1,2,:) = interpolation(squeeze(pha(1:2,2,1:idx_half))',... + x0(:,1:2),xs); + ir_new(1,1,:) = ifft(squeeze(mag_new(1,1,:)... + .*exp(1i*pha_new(1,1,:))),size(mag,3),'symmetric'); + ir_new(1,2,:) = ifft(squeeze(mag_new(1,2,:)... + .*exp(1i*pha_new(1,1,:))),size(mag,3),'symmetric'); + otherwise + error('%s: %s is an unknown interpolation method.',... + upper(mfilename),interpolationmethod); + end else % --- 2D interpolation --- warning('SFS:irs_intpol3D',... From b368184c97a39e7dd29d7b687bb51a162787f13b Mon Sep 17 00:00:00 2001 From: VeraE Date: Fri, 24 Jun 2016 14:23:05 +0200 Subject: [PATCH 03/17] refine backwards compatibility check --- SFS_ir/interpolate_ir.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SFS_ir/interpolate_ir.m b/SFS_ir/interpolate_ir.m index a4829490..29132ac0 100644 --- a/SFS_ir/interpolate_ir.m +++ b/SFS_ir/interpolate_ir.m @@ -65,7 +65,7 @@ %% ===== Configuration ================================================== useinterpolation = conf.ir.useinterpolation; % Check for old configuration -if ~isfield(conf.ir,'interpolationmethod') +if useinterpolation && ~isfield(conf.ir,'interpolationmethod') warning('SFS:irs_intpolmethod',... 'no interpolation method provided, will use method ''simple''.'); interpolationmethod = 'simple'; From a1205a681a02386417b81c2af3d4f99dfa925a3d Mon Sep 17 00:00:00 2001 From: VeraE Date: Mon, 27 Jun 2016 15:26:21 +0200 Subject: [PATCH 04/17] change conf.ir.useinterpolation in SFS_config back to true --- SFS_config.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SFS_config.m b/SFS_config.m index 5f15e88f..7a3e88ff 100644 --- a/SFS_config.m +++ b/SFS_config.m @@ -331,7 +331,7 @@ % disabled the HRTF returned by a nearest neighbour search is used instead. % Depending on the geometry of the measured HRTF data set, the interpolation % will be done between the two or three nearest HRTFs. -conf.ir.useinterpolation = false; % boolean +conf.ir.useinterpolation = true; % boolean % % You can choose between the following interpolation methods: % 'simple' - Interpolation in the time domain performed samplewise. This From 9abc582124675858909be3cea73408eb411721ce Mon Sep 17 00:00:00 2001 From: VeraE Date: Tue, 5 Jul 2016 17:46:41 +0200 Subject: [PATCH 05/17] correct interpolation to avoid phase aliasing --- SFS_ir/interpolate_ir.m | 33 ++++++++++++++++++--------------- 1 file changed, 18 insertions(+), 15 deletions(-) diff --git a/SFS_ir/interpolate_ir.m b/SFS_ir/interpolate_ir.m index 29132ac0..28dfe836 100644 --- a/SFS_ir/interpolate_ir.m +++ b/SFS_ir/interpolate_ir.m @@ -101,24 +101,27 @@ 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' - TF = fft(ir,[],3); - % Calculate interpolation for magnitude and phase separately + % Upsample to avoid phase aliasing in unwrapping of phase + TF = fft(ir,4*size(ir,3),3); + % Magnitude and phase will be interpolated separately mag = abs(TF); pha = unwrap(angle(TF),[],3); % Calculate interpolation only for the first half of the spectrum - idx_half = floor(size(mag,3)/2)+1; - mag_new(1,1,:) = interpolation(squeeze(mag(1:2,1,1:idx_half))',... - x0(:,1:2),xs); - mag_new(1,2,:) = interpolation(squeeze(mag(1:2,2,1:idx_half))',... - x0(:,1:2),xs); - pha_new(1,1,:) = interpolation(squeeze(pha(1:2,1,1:idx_half))',... - x0(:,1:2),xs); - pha_new(1,2,:) = interpolation(squeeze(pha(1:2,2,1:idx_half))',... - x0(:,1:2),xs); - ir_new(1,1,:) = ifft(squeeze(mag_new(1,1,:)... - .*exp(1i*pha_new(1,1,:))),size(mag,3),'symmetric'); - ir_new(1,2,:) = ifft(squeeze(mag_new(1,2,:)... - .*exp(1i*pha_new(1,1,:))),size(mag,3),'symmetric'); + % and only for original bins + idx_half = floor(size(TF,3)/2)+1; + mag_new(1,:) = interpolation(... + squeeze(mag(1:2,1,1:4:idx_half))',x0(:,1:2),xs); + mag_new(2,:) = interpolation(... + squeeze(mag(1:2,2,1:4:idx_half))',x0(:,1:2),xs); + pha_new(1,:) = interpolation(... + squeeze(pha(1:2,1,1:4:idx_half))',x0(:,1:2),xs); + pha_new(2,:) = interpolation(... + squeeze(pha(1:2,2,1:4:idx_half))',x0(:,1:2),xs); + % Calculate interpolated impulse response from magnitude and phase + ir_new(1,1,:) = ifft(mag_new(1,:)... + .*exp(1i*pha_new(1,:)),size(ir,3),'symmetric'); + ir_new(1,2,:) = ifft(mag_new(2,:)... + .*exp(1i*pha_new(2,:)),size(ir,3),'symmetric'); otherwise error('%s: %s is an unknown interpolation method.',... upper(mfilename),interpolationmethod); From 0a372d9b999e1d0d68143d8c990947f54e4c1c62 Mon Sep 17 00:00:00 2001 From: VeraE Date: Tue, 5 Jul 2016 18:03:24 +0200 Subject: [PATCH 06/17] add references --- SFS_ir/interpolate_ir.m | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/SFS_ir/interpolate_ir.m b/SFS_ir/interpolate_ir.m index 28dfe836..bb0f801b 100644 --- a/SFS_ir/interpolate_ir.m +++ b/SFS_ir/interpolate_ir.m @@ -24,6 +24,13 @@ % 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 %***************************************************************************** @@ -101,6 +108,8 @@ 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 From 7d8022ade0c42dd484f311d7b684edbe22ddf104 Mon Sep 17 00:00:00 2001 From: Vera Erbes Date: Tue, 5 Jul 2016 18:08:59 +0200 Subject: [PATCH 07/17] remove additional tabs --- SFS_ir/interpolate_ir.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SFS_ir/interpolate_ir.m b/SFS_ir/interpolate_ir.m index bb0f801b..e63615bf 100644 --- a/SFS_ir/interpolate_ir.m +++ b/SFS_ir/interpolate_ir.m @@ -108,8 +108,8 @@ 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) - % + % 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 From e748d62920758724d13b97d52c969ce349844a7a Mon Sep 17 00:00:00 2001 From: VeraE Date: Tue, 5 Jul 2016 18:25:59 +0200 Subject: [PATCH 08/17] align comment --- SFS_ir/interpolate_ir.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SFS_ir/interpolate_ir.m b/SFS_ir/interpolate_ir.m index e63615bf..161ccc54 100644 --- a/SFS_ir/interpolate_ir.m +++ b/SFS_ir/interpolate_ir.m @@ -108,8 +108,8 @@ 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) - % + % 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 From 289d15420b31b8fa33831ef90c7f1ca155b71a18 Mon Sep 17 00:00:00 2001 From: Hagen Wierstorf Date: Wed, 6 Jul 2016 10:21:31 +0200 Subject: [PATCH 09/17] Clean up --- SFS_ir/interpolate_ir.m | 60 ++++++++++++++++++++--------------------- 1 file changed, 30 insertions(+), 30 deletions(-) diff --git a/SFS_ir/interpolate_ir.m b/SFS_ir/interpolate_ir.m index 161ccc54..b87a10d1 100644 --- a/SFS_ir/interpolate_ir.m +++ b/SFS_ir/interpolate_ir.m @@ -104,36 +104,36 @@ deg(x0(1,1)), deg(x0(2,1)), ... deg(x0(1,2)), deg(x0(2,2))); 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 - mag = abs(TF); - pha = 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; - mag_new(1,:) = interpolation(... - squeeze(mag(1:2,1,1:4:idx_half))',x0(:,1:2),xs); - mag_new(2,:) = interpolation(... - squeeze(mag(1:2,2,1:4:idx_half))',x0(:,1:2),xs); - pha_new(1,:) = interpolation(... - squeeze(pha(1:2,1,1:4:idx_half))',x0(:,1:2),xs); - pha_new(2,:) = interpolation(... - squeeze(pha(1:2,2,1:4:idx_half))',x0(:,1:2),xs); - % Calculate interpolated impulse response from magnitude and phase - ir_new(1,1,:) = ifft(mag_new(1,:)... - .*exp(1i*pha_new(1,:)),size(ir,3),'symmetric'); - ir_new(1,2,:) = ifft(mag_new(2,:)... - .*exp(1i*pha_new(2,:)),size(ir,3),'symmetric'); - otherwise - error('%s: %s is an unknown interpolation method.',... - upper(mfilename),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 --- From ac5d9c64173bf52d5e7b48b359d7c185b312fdf0 Mon Sep 17 00:00:00 2001 From: VeraE Date: Wed, 6 Jul 2016 11:02:58 +0200 Subject: [PATCH 10/17] update header comment --- SFS_ir/interpolate_ir.m | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/SFS_ir/interpolate_ir.m b/SFS_ir/interpolate_ir.m index b87a10d1..40a215dc 100644 --- a/SFS_ir/interpolate_ir.m +++ b/SFS_ir/interpolate_ir.m @@ -20,6 +20,12 @@ % 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. +% 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 time of arrival of the HRTFs. +% 'freqdomain' - Interpolation in the frequency domain performed separately +% for magnitude and phase of the HRTF. % Note, that the given parameter are not checked if they have all the correct % dimensions in order to save computational time, because this function could % be called quite often. From 14ecc3cbf71fd888f0c538c0c1e84318d94dcd4b Mon Sep 17 00:00:00 2001 From: VeraE Date: Wed, 20 Jul 2016 16:16:39 +0200 Subject: [PATCH 11/17] update comments to explain new interpolation method --- SFS_config.m | 15 ++++++++++++--- SFS_ir/interpolate_ir.m | 11 ++++++----- 2 files changed, 18 insertions(+), 8 deletions(-) diff --git a/SFS_config.m b/SFS_config.m index ab73bce9..9c7d74ea 100644 --- a/SFS_config.m +++ b/SFS_config.m @@ -338,10 +338,19 @@ % % You can choose between the following interpolation methods: % 'simple' - Interpolation in the time domain performed samplewise. This -% does not heed the time of arrival of the HRTFs. +% does not heed the times of arrival of the impulse responses. % 'freqdomain' - Interpolation in the frequency domain performed separately -% for magnitude and phase of the HRTF. -conf.ir.interpolationmethod = 'freqdomain'; +% 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 diff --git a/SFS_ir/interpolate_ir.m b/SFS_ir/interpolate_ir.m index 40a215dc..52820aa1 100644 --- a/SFS_ir/interpolate_ir.m +++ b/SFS_ir/interpolate_ir.m @@ -14,8 +14,8 @@ % 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 @@ -23,10 +23,11 @@ % 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 time of arrival of the HRTFs. +% This does not heed the times of arrival of the impulse +% responses. % 'freqdomain' - Interpolation in the frequency domain performed separately -% for magnitude and phase of the HRTF. -% Note, that the given parameter are not checked if they have all the correct +% 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. % From dcc0d0e2e7caa2d04bdb21bdfa4513bef4feca99 Mon Sep 17 00:00:00 2001 From: Hagen Wierstorf Date: Mon, 1 Aug 2016 10:51:25 +0200 Subject: [PATCH 12/17] Update comment --- SFS_ir/interpolate_ir.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SFS_ir/interpolate_ir.m b/SFS_ir/interpolate_ir.m index 52820aa1..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 From 3f36877b370f58a10ee20acda75e6017feb370bc Mon Sep 17 00:00:00 2001 From: VeraE Date: Mon, 8 Aug 2016 18:38:45 +0200 Subject: [PATCH 13/17] add test function for interpolation methods --- validation/test_interpolation_methods.m | 252 ++++++++++++++++++++++++ 1 file changed, 252 insertions(+) create mode 100644 validation/test_interpolation_methods.m diff --git a/validation/test_interpolation_methods.m b/validation/test_interpolation_methods.m new file mode 100644 index 00000000..2fc9476f --- /dev/null +++ b/validation/test_interpolation_methods.m @@ -0,0 +1,252 @@ +function status = test_interpolation_methods(modus) +%TEST_INTERPOLATION_METHOD 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 +hrtf_file = 'KEMAR_HRTFs_lfcorr.sofa'; +if ~exist(hrtf_file,'file') + disp('Download'); + url = ['https://github.com/spatialaudio/lf-corrected-kemar-hrtfs/', ... + 'blob/master/' hrtf_file '?raw=true']; + basepath = get_sfs_path(); + file_path = [basepath '/data/HRTFs/']; + download_file(url,[file_path hrtf_file]); + addpath(genpath(file_path)) +end +% Load HRTF data set +hrtf = SOFAload(hrtf_file); + + +%% ===== Interpolate shifted Dirac impulses ============================== +% interpolation points +x0 = [1 3; 0 0; 0 0]; +% target point +xs = [2; 0; 0]; +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,:)']; +xs_close = hrtf.SourcePosition(idx1,:)'; +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,:)']; +xs_dist = hrtf.SourcePosition(idx30,:)'; +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; From e6f7bd51aeb7dcde910b42917ee051c1c25bc67b Mon Sep 17 00:00:00 2001 From: VeraE Date: Mon, 8 Aug 2016 18:46:17 +0200 Subject: [PATCH 14/17] typos --- SFS_general/interpolation.m | 2 +- SFS_general/secondary_source_diameter.m | 2 +- SFS_general/secondary_source_positions.m | 2 +- SFS_ssr/ssr_brs_wfs.m | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) 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_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: From 7456634ab634e7b4d522a47ebcefa88d75565ff1 Mon Sep 17 00:00:00 2001 From: VeraE Date: Mon, 8 Aug 2016 19:00:02 +0200 Subject: [PATCH 15/17] switch degree and radians for nicer warning in test function --- validation/test_interpolation_methods.m | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/validation/test_interpolation_methods.m b/validation/test_interpolation_methods.m index 2fc9476f..e31103df 100644 --- a/validation/test_interpolation_methods.m +++ b/validation/test_interpolation_methods.m @@ -83,9 +83,9 @@ %% ===== Interpolate shifted Dirac impulses ============================== % interpolation points -x0 = [1 3; 0 0; 0 0]; +x0 = [1 3; 0 0; 0 0]/180*pi; % target point -xs = [2; 0; 0]; +xs = [2; 0; 0]/180*pi; N = 50; %length impulse responses % 1. Interpolate between Dirac impulses with one sample in between @@ -171,8 +171,8 @@ 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,:)']; -xs_close = hrtf.SourcePosition(idx1,:)'; +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,:,:); @@ -185,8 +185,8 @@ 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,:)']; -xs_dist = hrtf.SourcePosition(idx30,:)'; +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,:,:); From 1c922755aa1111ebc9f0f7d8996f8c67f2eaa5de Mon Sep 17 00:00:00 2001 From: Hagen Wierstorf Date: Tue, 9 Aug 2016 12:44:49 +0200 Subject: [PATCH 16/17] Update automatic file download --- validation/test_interpolation_methods.m | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/validation/test_interpolation_methods.m b/validation/test_interpolation_methods.m index e31103df..e3fcf8ef 100644 --- a/validation/test_interpolation_methods.m +++ b/validation/test_interpolation_methods.m @@ -67,15 +67,14 @@ %% ===== Main ============================================================ % Check if HRTF data set with low frequency correction is available, % download otherwise -hrtf_file = 'KEMAR_HRTFs_lfcorr.sofa'; +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_file '?raw=true']; - basepath = get_sfs_path(); - file_path = [basepath '/data/HRTFs/']; - download_file(url,[file_path hrtf_file]); - addpath(genpath(file_path)) + 'blob/master/' hrtf '?raw=true']; + download_file(url,hrtf_file); end % Load HRTF data set hrtf = SOFAload(hrtf_file); From f60d1d403a3519ef152ca63ed1aca94e20e32e34 Mon Sep 17 00:00:00 2001 From: Hagen Wierstorf Date: Tue, 9 Aug 2016 12:46:20 +0200 Subject: [PATCH 17/17] Fix typo in comment --- validation/test_interpolation_methods.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/validation/test_interpolation_methods.m b/validation/test_interpolation_methods.m index e3fcf8ef..cc29ea4e 100644 --- a/validation/test_interpolation_methods.m +++ b/validation/test_interpolation_methods.m @@ -1,5 +1,5 @@ function status = test_interpolation_methods(modus) -%TEST_INTERPOLATION_METHOD demonstrates the difference between interpolation +%TEST_INTERPOLATION_METHODS demonstrates the difference between interpolation %methods % % Usage: status = test_interpolation_methods(modus) @@ -10,7 +10,7 @@ % % 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