From 74c7ad29db98781b00d381364ef845e194449f3e Mon Sep 17 00:00:00 2001 From: Fiete Winter Date: Tue, 6 Nov 2018 11:28:52 +0100 Subject: [PATCH 1/3] initial commit for aliasing model --- SFS_aliasing/aliasing_extended_control.m | 81 ++++++++++++++ SFS_aliasing/aliasing_extended_modal.m | 81 ++++++++++++++ SFS_aliasing/aliasing_frequency.m | 103 ++++++++++++++++++ .../local_wavenumber_vector.m | 74 ++++--------- SFS_aliasing/minmax_kt_circle.m | 54 +++++++++ SFS_start.m | 2 + 6 files changed, 343 insertions(+), 52 deletions(-) create mode 100644 SFS_aliasing/aliasing_extended_control.m create mode 100644 SFS_aliasing/aliasing_extended_modal.m create mode 100644 SFS_aliasing/aliasing_frequency.m rename SFS_general/aliasing_frequency.m => SFS_aliasing/local_wavenumber_vector.m (52%) create mode 100644 SFS_aliasing/minmax_kt_circle.m diff --git a/SFS_aliasing/aliasing_extended_control.m b/SFS_aliasing/aliasing_extended_control.m new file mode 100644 index 00000000..fc83f30f --- /dev/null +++ b/SFS_aliasing/aliasing_extended_control.m @@ -0,0 +1,81 @@ +function f = aliasing_extended_control(x0, kSx0, x, minmax_kGt_fun, minmax_kSt_fun, conf) +%ALIASING_EXTENDED_CONTROL aliasing frequency for an extended listening area +%with an defined control area where the sound field synthesis is prioritized +% +% Usage: f = aliasing_extended_control(x0, kSx0, x, minmax_kGt_fun, minmax_kSt_fun, conf) +% +% Input options: +% x0 - position, direction, and sampling distance of +% secondary sources [N0x7] / m +% kSx0 - normalised local wavenumber vector kS(x0) +% of virtual sound field at x0 [N0x3] +% x - position for which aliasing frequency is calculated +% [Nx3] +% minmax_kGt_fun - function handle to determine the extremal value of +% the tangential component of k_G(x-x0) +% [kGtmin, kGtmax] = minmax_kGt_fun(x0,x) +% minmax_kSt_fun - function handle to determine the extremal value of +% the tangential component of k_S(x0) +% [kStmin, kStmax] = minmax_kSt_fun(x0) +% conf - configuration struct (see SFS_config) +% +% Output parameters: +% f - aliasing frequency [Nx1] +% + +%***************************************************************************** +% The MIT License (MIT) * +% * +% Copyright (c) 2010-2018 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 * +%***************************************************************************** + + +phik = cart2pol(kSx0(:,1),kSx0(:,2)); % azimuth angle of kS(x0) +phin0 = cart2pol(x0(:,4),x0(:,5)); % azimuth angle of normal vector n0(x0) + +% secondary source selection +select = cos(phin0 - phik) >= 0; +x0 = x0(select,:); +% kSt(x0) (tangential component of kS(x0) ) +kSt = sin(phin0(select) - phik(select)); + +% mininum and maximum values of kSt(x_0) +[kStmin, kStmax] = minmax_kSt_fun(x0); +select = kSt >= kStmin & kSt <= kStmax; +x0 = x0(select,:); +kSt = kSt(select); + +% sampling distance +deltax0 = abs(x0(:,7)); + +f = inf(size(x,1), 1); +for xdx = 1:size(x,1); + % mininum and maximum values of k_Gt(x - x_0) + % (tangential component of k_G(x-x0)) + [kGtmin, kGtmax] = minmax_kGt_fun(x0,x(xdx,:)); + % aliasing frequency for x + f(xdx) = conf.c./max(deltax0.*max(abs(kSt-kGtmin),abs(kSt-kGtmax))); +end diff --git a/SFS_aliasing/aliasing_extended_modal.m b/SFS_aliasing/aliasing_extended_modal.m new file mode 100644 index 00000000..64f17d96 --- /dev/null +++ b/SFS_aliasing/aliasing_extended_modal.m @@ -0,0 +1,81 @@ +function f = aliasing_extended_modal(x0, kSx0, x, minmax_kGt_fun, xc, M, conf) +%ALIASING_EXTENDED_MODAL aliasing frequency for an extended listening area for +%an circular control area at xc with R=M/k where synthesis is focused on. +% +% Usage: f = aliasing_extended_control(x0, kSx0, x, minmax_kGt_fun, minmax_kSt_fun, conf) +% +% Input options: +% x0 - position, direction, and sampling distance of +% secondary sources [N0x7] / m +% kSx0 - normalised local wavenumber vector of virtual sound +% field [N0x3] +% x - position for which aliasing frequency is calculated +% [Nx3] +% minmax_kGt_fun - function handle to determine the extremal value of +% the tangential component of k_G(x-x0) +% [kGtmin, kGtmax] = minmax_kGt_fun(x0) +% xc - center of circular control area +% M - modal order which defines the radius R=M/k +% conf - configuration struct (see SFS_config) +% +% Output parameters: +% f - aliasing frequency [Nx1] +% + +%***************************************************************************** +% The MIT License (MIT) * +% * +% Copyright (c) 2010-2018 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 * +%***************************************************************************** + +phik = cart2pol(kSx0(:,1),kSx0(:,2)); % azimuth angle of kS(x0) +phin0 = cart2pol(x0(:,4),x0(:,5)); % azimuth angle of normal vector n0(x0) + +% secondary source selection +select = cos(phin0 - phik) >= 0; +x0 = x0(select,:); +% k_St(x0) (tangential component of k_S(x0) ) +kSt = sin(phin0(select) - phik(select)); + +% sampling distance +deltax0 = abs(x0(:,7)); + +f = inf(size(x,1), 1); +for xdx = 1:size(x,1); + % mininum and maximum values of k_Gt(x - x_0) + % (tangential component of k_G(x-x0)) + [kGtmin, kGtmax] = minmax_kGt_fun(x0,x(xdx,:)); + % aliasing frequency for x0 + f0 = conf.c./(deltax0.*max(abs(kSt-kGtmin),abs(kSt-kGtmax))); + % radius of local region at f0 + rM = M.*conf.c./(2.*pi.*f0); + % mininum and maximum values of kSt(x_0) at f0 + [kStmin, kStmax] = minmax_kt_circle(x0, xc, rM); + select = kSt > kStmin & kSt < kStmax; + if sum(select) ~= 0 + f(xdx) = min(f0(select)); + end +end diff --git a/SFS_aliasing/aliasing_frequency.m b/SFS_aliasing/aliasing_frequency.m new file mode 100644 index 00000000..ef183e29 --- /dev/null +++ b/SFS_aliasing/aliasing_frequency.m @@ -0,0 +1,103 @@ +function f = aliasing_frequency(method, xs, src, area, x, varargin) +%ALIASING_FREQUENCY Summary of this function goes here +% Detailed explanation goes here + +%***************************************************************************** +% The MIT License (MIT) * +% * +% Copyright (c) 2010-2018 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 * +%***************************************************************************** + +switch area + case 'position' + conf = varargin{1}; + % the position x is a circular area with radius 0 + minmax_kGt_fun = @(x0p,xp) minmax_kt_circle(x0p,xp,0); % handle + case {'circle', 'circular'} + Rl = varargin{1}; % radius of circular area + conf = varargin{2}; + minmax_kGt_fun = @(x0p,xp) minmax_kt_circle(x0p,xp,Rl); % handle + case 'ellipse' + to_be_implemented(mfilename); + otherwise + error('%s: unknown type (%s) for listening area', upper(mfilename), area); +end + +N0 = conf.secondary_sources.number; +switch conf.secondary_sources.geometry + case {'circle', 'circular'} + R = conf.secondary_sources.size / 2; + deltax = 2*pi*R/N0; + case 'linear' + L = conf.secondary_sources.size; + deltax = L/(N0-1); + otherwise + error('%s: unsupported geometry (%s) for secondary sources', ... + upper(mfilename), conf.secondary_sources.geometry); +end + +% densely sampled SSD for calculating the aliasing frequency +conf.secondary_sources.number = 2048; +x0 = secondary_source_positions(conf); +x0 = secondary_source_selection(x0,xs,src); % a_S(x0) >= 0 +x0(:,7) = deltax; % correct sampling + +% local wavenumber of virtual sound field at x0 +kSx0 = local_wavenumber_vector(x0(:,1:3), xs, src); + +switch method + case {'wfs', 'WFS'} + % infinite control area + minmax_kSt_fun = @(x0p) minmax_kt_circle(x0p,[0,0,0],Inf); % handle + f = aliasing_extended_control(x0, kSx0, x, minmax_kGt_fun, ... + minmax_kSt_fun, conf); + case {'nfchoa', 'HOA'} + % circular control area at array center with Rl = M/k + M = nfchoa_order(N0,conf); + xc = conf.secondary_sources.center; + f = aliasing_extended_modal(x0, kSx0, x, minmax_kGt_fun, xc, M, conf); + case {'localwfs-sbl', 'LWFS-SBL'} + % circular control area at xref with Rl = M/k + xc = conf.xref; + M = conf.localwfs_sbl.order; + f = aliasing_extended_modal(x0, kSx0, x, minmax_kGt_fun, xc, M, conf); + case {'localwfs-vss', 'LWFS-VSS'} + switch conf.localwfs_vss.geometry + case {'circle', 'circular'} + xc = conf.localwfs_vss.center; + Rc = conf.localwfs_vss.size / 2; + minmax_kSt_fun = @(x0p) minmax_kt_circle(x0p,xc,Rc); % handle + otherwise + error('%s: unsupported geometry (%s) for LWFS-VSS', ... + upper(mfilename), conf.localwfs_vss.geometry); + end + f = aliasing_extended_control(x0, kSx0, x, minmax_kGt_fun, ... + minmax_kSt_fun, conf); + otherwise + error('%s: unsupported SFS method (%s)', ... + upper(mfilename), method) +end + diff --git a/SFS_general/aliasing_frequency.m b/SFS_aliasing/local_wavenumber_vector.m similarity index 52% rename from SFS_general/aliasing_frequency.m rename to SFS_aliasing/local_wavenumber_vector.m index c80648ef..abfefcca 100644 --- a/SFS_general/aliasing_frequency.m +++ b/SFS_aliasing/local_wavenumber_vector.m @@ -1,36 +1,5 @@ -function [fal,dx0] = aliasing_frequency(x0,conf) -%ALIASING_FREQUENCY aliasing frequency for the given secondary sources +function kvec = local_wavenumber_vector(x, xs, src) % -% Usage: [fal,dx0] = aliasing_frequency([x0],conf) -% -% Input parameters: -% x0 - secondary sources / m -% conf - configuration struct (see SFS_config) -% -% Output parameters: -% fal - aliasing frequency / Hz -% dx0 - mean distance between secondary sources / m -% -% ALIASING_FREQUENCY(x0,conf) returns the aliasing frequency for the given -% secondary sources. First the mean distance dx0 between the secondary sources -% is calculated, afterwards the aliasing frequency is calculated after Spors -% (2009) as fal = c/(2*dx0). If no secondary sources x0 are provided, they are -% first calculated by calling secondary_source_positions(). -% For a calculation that includes the dependency on the listener position have -% a look at Start (1997). -% -% See also: sound_field_mono_wfs, secondary_source_positions, -% secondary_source_distance -% -% References: -% Spors and Ahrens (2009) - "Spatial sampling artifacts of wave field -% synthesis for the reproduction of virtual point sources", 126th -% Convention of the Audio Engineering Society, Paper 7744, -% http://www.aes.org/e-lib/browse.cfm?elib=14940 -% -% Start (1997) - "Direct Sound Enhancement by Wave Field Synthesis", -% PhD thesis, TU Delft, -% http://resolver.tudelft.nl/uuid:c80d5b58-67d3-4d84-9e73-390cd30bde0d %***************************************************************************** % The MIT License (MIT) * @@ -61,28 +30,29 @@ % http://sfstoolbox.org sfstoolbox@gmail.com * %***************************************************************************** - -%% ===== Checking of input parameters ==================================== -nargmin = 1; -nargmax = 2; +%% ===== Checking of input parameters =================================== +nargmin = 3; +nargmax = 3; narginchk(nargmin,nargmax); -if nargin 1 | -sqrt(1-rho.^2) > klt; +kGtmin(select) = -1; +kGtmin(~select) = klt(~select).*sqrt(1-rho(~select).^2) ... + - kln(~select).*rho(~select); + +kGtmax = zeros(size(xc,1),1); +select = rho > 1 | sqrt(1-rho.^2) < klt; +kGtmax(select) = +1; +kGtmax(~select) = klt(~select).*sqrt(1-rho(~select).^2) ... + + kln(~select).*rho(~select); diff --git a/SFS_start.m b/SFS_start.m index 98587c13..c24a399c 100644 --- a/SFS_start.m +++ b/SFS_start.m @@ -63,6 +63,7 @@ function SFS_start(printbanner) % Add the base path and the needed sub-directories if exist('addpath') addpath(basepath); + addpath([basepath,'/SFS_aliasing']); addpath([basepath,'/SFS_analysis']); addpath([basepath,'/SFS_binaural_synthesis']); addpath([basepath,'/SFS_general']); @@ -86,6 +87,7 @@ function SFS_start(printbanner) end else path(path,basepath); + path(path,[basepath,'/SFS_aliasing']); path(path,[basepath,'/SFS_analysis']); path(path,[basepath,'/SFS_binaural_synthesis']); path(path,[basepath,'/SFS_general']); From e6a096598afea98ebf3491372e72bc8b06aa6387 Mon Sep 17 00:00:00 2001 From: Fiete Winter Date: Thu, 15 Nov 2018 17:44:42 +0100 Subject: [PATCH 2/3] add support for discretised virtual secondary source distribution and discretised plane wave decomposition --- SFS_aliasing/aliasing_extended_sbl.m | 161 ++++++++++++++++++ ...nded_control.m => aliasing_extended_vss.m} | 79 ++++++--- SFS_aliasing/aliasing_frequency.m | 80 ++++++--- ...sing_extended_modal.m => minmax_kt_line.m} | 101 ++++++----- 4 files changed, 333 insertions(+), 88 deletions(-) create mode 100644 SFS_aliasing/aliasing_extended_sbl.m rename SFS_aliasing/{aliasing_extended_control.m => aliasing_extended_vss.m} (59%) rename SFS_aliasing/{aliasing_extended_modal.m => minmax_kt_line.m} (51%) diff --git a/SFS_aliasing/aliasing_extended_sbl.m b/SFS_aliasing/aliasing_extended_sbl.m new file mode 100644 index 00000000..a7b0407b --- /dev/null +++ b/SFS_aliasing/aliasing_extended_sbl.m @@ -0,0 +1,161 @@ +function f = aliasing_extended_sbl(x0, kS, x, minmax_kGt_fun, xc, M, Npw, conf) +%ALIASING_EXTENDED_SBL aliasing frequency for an extended listening area for +%an circular control area at xc with R=M/k where synthesis is focused on. +% +% Usage: f = aliasing_extended_control(x0, kSx0, x, minmax_kGt_fun, minmax_kSt_fun, conf) +% +% Input options: +% x0 - position, direction, and sampling distance of +% secondary sources [N0x7] / m +% kSx0 - normalised local wavenumber vector of virtual sound +% field [N0x3] +% x - position for which aliasing frequency is calculated +% [Nx3] +% minmax_kGt_fun - function handle to determine the extremal value of +% the tangential component of k_G(x-x0) +% [kGtmin, kGtmax] = minmax_kGt_fun(x0) +% xc - center of circular control area +% M - modal order which defines the radius R=M/k +% Npw - +% conf - configuration struct (see SFS_config) +% +% Output parameters: +% f - aliasing frequency [Nx1] +% + +%***************************************************************************** +% The MIT License (MIT) * +% * +% Copyright (c) 2010-2018 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 * +%***************************************************************************** + +% shift coordinates to expansion centre xc +x = bsxfun(@minus, x, xc); +x0(:,1:3) = bsxfun(@minus, x0(:,1:3), xc); + +phin0 = cart2pol(x0(:,4),x0(:,5)); % azimuth angle of normal vector n0(x0) +deltax0 = abs(x0(:,7)); % sampling distance delta_x0(x0) +phik = cart2pol(kS(:,1),kS(:,2)); % azimuth angle of kS(x0) + +select = cos(phin0 - phik) >= 0; % secondary source selection +gdx = find(select); + +% quantities for selected secondary sources +x0select = x0(select,:); +[phix0select, rx0select] = cart2pol(x0select(:,1),x0select(:,2)); +deltax0select = deltax0(select); % sampling distance +kSselect = kS(select,:); +phikselect = phik(select); +kStselect = sin(phin0(select) - phikselect); % tangential component of kS(x0) + +% wavelength at which secondary source x0 turn inactive due to spatial +% bandwidth limitation +lambdaMselect = 2.*pi/M.*rx0select.*abs(sin(phix0select - phikselect)); + +f = inf(size(x,1), 1); +if Npw < M + f(:) = 0; + return +end + +eps = 0.01; + +for xdx = 1:size(x,1); + + vectorxx0 = bsxfun(@minus, x(xdx,:), x0(:,1:3)); % vector x - x0 + [phixx0,rxx0] = cart2pol(vectorxx0(:,1),vectorxx0(:,2)); % polar coordinates + % mininum and maximum values of k_Gt(x - x_0) + % (tangential component of k_G(x-x0)) + [kGtmin, kGtmax] = minmax_kGt_fun(x0,x(xdx,:)); + kGtminselect = kGtmin(gdx); + kGtmaxselect = kGtmax(gdx); + + lambda = 0; + + % aliasing wavelength for eta=1 and zeta=0 + lambda10 = deltax0select.* ... + max(abs(kStselect-kGtminselect),abs(kStselect-kGtmaxselect)); + % select all wavelengths, that above the limit of the spatial bandwidth + % limitation + select = lambda10 >= lambdaMselect; + if any(select) + lambda = max(lambda, max(lambda10(select))); + end + + if isinf(Npw) + f(xdx) = conf.c./lambda; + continue; + end + + % aliasing wavelength for eta=0 and zeta=1 + lambda01 = 2.*pi./Npw.*rxx0(gdx).*abs(sin(phixx0(gdx)-phikselect)); + % select all wavelengths, that above the limit of the spatial bandwidth + % limitation + select = lambda01 >= lambdaMselect; + if any(select) + lambda = max(lambda, max(lambda01(select))); + end + + % aliasing wavelength for eta=1 and zeta=1 + for sdx=1:size(x0select,1); % interate over all active x0' + + % secondary source selection for ALL x0 with a single kS(x0') + active = x0(:,4:6)*kSselect(sdx,:).' >= 0; + + % vector x0 - x0' + vectorx0x0p = bsxfun(@minus, x0(active,1:3), x0select(sdx,1:3)); + % polar angle and radius of vector (x0-x0') + [phix0x0p,rx0x0p] = cart2pol(vectorx0x0p(:,1),vectorx0x0p(:,2)); + + % range for aliasing wavelength at x0' caused by Discrete SSD (eta = 1) + lambda10x0p1 = ... + abs(sin(phin0(active) - phikselect(sdx))-kGtmin(active)) ... + .*deltax0(active); + lambda10x0p2 = ... + abs(sin(phin0(active) - phikselect(sdx))-kGtmax(active)) ... + .*deltax0(active); + lambda10x0pmin = min(lambda10x0p1,lambda10x0p2); + lambda10x0pmax = max(lambda10x0p1,lambda10x0p2); + + % aliasing wavelength at x0' caused by Plane Wave Decomposition (zeta = 1) + lambda01x0p = 2.*pi/Npw.*rx0x0p.*abs(sin(phix0x0p - phikselect(sdx))); + + % select x0' where PWD wavelength is in the range of SSD wavelength + select = lambda01x0p < (1+eps).*lambda10x0pmax & ... + lambda10x0pmin < (1+eps).*lambda01x0p; + if any(select) + lambda11 = lambda01x0p(select); + % select all wavelengths, that above the limit of the spatial bandwidth + % limitation + select = lambda11 >= lambdaMselect(sdx); + if any(select) + lambda = max(lambda, max(lambda11(select))); + end + end + end + + f(xdx) = conf.c./lambda; +end diff --git a/SFS_aliasing/aliasing_extended_control.m b/SFS_aliasing/aliasing_extended_vss.m similarity index 59% rename from SFS_aliasing/aliasing_extended_control.m rename to SFS_aliasing/aliasing_extended_vss.m index fc83f30f..6b708c2a 100644 --- a/SFS_aliasing/aliasing_extended_control.m +++ b/SFS_aliasing/aliasing_extended_vss.m @@ -1,4 +1,4 @@ -function f = aliasing_extended_control(x0, kSx0, x, minmax_kGt_fun, minmax_kSt_fun, conf) +function f = aliasing_extended_vss(x0, xv, kSxv, x, minmax_kGt_fun, conf) %ALIASING_EXTENDED_CONTROL aliasing frequency for an extended listening area %with an defined control area where the sound field synthesis is prioritized % @@ -52,30 +52,71 @@ % http://sfstoolbox.org sfstoolbox@gmail.com * %***************************************************************************** - -phik = cart2pol(kSx0(:,1),kSx0(:,2)); % azimuth angle of kS(x0) +phik = cart2pol(kSxv(:,1),kSxv(:,2)); % azimuth angle of kS(xv) +phinv = cart2pol(xv(:,4),xv(:,5)); % azimuth angle of normal vector nv(xv) phin0 = cart2pol(x0(:,4),x0(:,5)); % azimuth angle of normal vector n0(x0) -% secondary source selection -select = cos(phin0 - phik) >= 0; -x0 = x0(select,:); -% kSt(x0) (tangential component of kS(x0) ) -kSt = sin(phin0(select) - phik(select)); - -% mininum and maximum values of kSt(x_0) -[kStmin, kStmax] = minmax_kSt_fun(x0); -select = kSt >= kStmin & kSt <= kStmax; -x0 = x0(select,:); -kSt = kSt(select); +% selection of virtual secondary source +select = cos(phinv - phik) >= 0; +xv = xv(select,:); +% k_S,tv(xv) (tangential component of kS(xv) at xv ) +kStv = sin(phinv(select) - phik(select)); % sampling distance deltax0 = abs(x0(:,7)); +deltaxv = abs(xv(:,7)); f = inf(size(x,1), 1); + +eps = 0.01; for xdx = 1:size(x,1); - % mininum and maximum values of k_Gt(x - x_0) - % (tangential component of k_G(x-x0)) - [kGtmin, kGtmax] = minmax_kGt_fun(x0,x(xdx,:)); - % aliasing frequency for x - f(xdx) = conf.c./max(deltax0.*max(abs(kSt-kGtmin),abs(kSt-kGtmax))); + + % mininum and maximum values of k_G,tv(x-xv) + % (tangential component of k_G(x-xv) at xv) + [kGtvmin, kGtvmax] = minmax_kGt_fun(xv,x(xdx,:)); + lambda = max(deltaxv.*max(abs(kStv-kGtvmin),abs(kStv-kGtvmax))); + + % mininum and maximum values of k_G,t0(x-x0) + % (tangential component of k_G(x-x0) at x0) + [kGt0min, kGt0max] = minmax_kGt_fun(x0,x(xdx,:)); + + for vdx=1:size(xv,1) % interate over all active xv + + vectorxvx0 = bsxfun(@minus, xv(vdx,1:3), x0(:,1:3)); % vector xv - x0 + + % secondary source selection for ALL x0 and a single xv + active = vectorxvx0*xv(vdx,4:6).' >= 0; + + % polar angle and radius of vector (xv-x0) + phixvx0 = cart2pol(vectorxvx0(active,1),vectorxvx0(active,2)); + + % k_FS,tv(xv-x0) (tangential component of k_FS(xv-x0) at xv) + kFStv = sin(phinv(vdx) - phixvx0); + % k_FS,t0(xv-x0) (tangential component of k_FS(xv-x0) at x0) + kFSt0 = sin(phin0(active) - phixvx0); + + % range for aliasing wavelength at x0 caused by discrete SSD (eta = 1) + lambda10_1 = abs(kFSt0 - kGt0min(active)).*deltax0(active); + lambda10_2 = abs(kFSt0 - kGt0max(active)).*deltax0(active); + + lambda10min = min(lambda10_1,lambda10_2); + lambda10max = max(lambda10_1,lambda10_2); + + % aliasing wavelength at xv caused by discrete virtual SSD (zeta = 1) + lambda01 = abs((kStv(vdx) - kFStv).*deltaxv(vdx)); + + % select the x0 where the virtual SSD wavelength is in the range of SSD + % wavelength + select = lambda01 < (1+eps).*lambda10max & lambda10min < (1+eps).*lambda01; + if any(select) + lambda = max(lambda, max(lambda01(select))); + end + + select = lambda01 < eps; + if any(select) + lambda = max(lambda, max(lambda10max(select))); + end + end + + f(xdx) = conf.c./lambda; end diff --git a/SFS_aliasing/aliasing_frequency.m b/SFS_aliasing/aliasing_frequency.m index ef183e29..e5b7c31d 100644 --- a/SFS_aliasing/aliasing_frequency.m +++ b/SFS_aliasing/aliasing_frequency.m @@ -31,71 +31,105 @@ % http://sfstoolbox.org sfstoolbox@gmail.com * %***************************************************************************** +%% ===== Checking input parameters ======================================= +nargmin = 6; +nargmax = 7; +narginchk(nargmin,nargmax); +isargchar(method, area, src); +isargxs(xs); +ismatrix(x); + switch area case 'position' - conf = varargin{1}; % the position x is a circular area with radius 0 minmax_kGt_fun = @(x0p,xp) minmax_kt_circle(x0p,xp,0); % handle + conf = varargin{1}; + case 'line' + to_be_implemented(mfilename); case {'circle', 'circular'} Rl = varargin{1}; % radius of circular area - conf = varargin{2}; + isargpositivescalar(Rl); minmax_kGt_fun = @(x0p,xp) minmax_kt_circle(x0p,xp,Rl); % handle + conf = varargin{2}; case 'ellipse' to_be_implemented(mfilename); otherwise error('%s: unknown type (%s) for listening area', upper(mfilename), area); end +isargstruct(conf); + +%% ===== Computation ==================================================== +% prepare virtual secondary sources N0 = conf.secondary_sources.number; switch conf.secondary_sources.geometry case {'circle', 'circular'} R = conf.secondary_sources.size / 2; - deltax = 2*pi*R/N0; + deltax0 = 2*pi*R/N0; case 'linear' L = conf.secondary_sources.size; - deltax = L/(N0-1); + deltax0 = L/(N0-1); otherwise error('%s: unsupported geometry (%s) for secondary sources', ... upper(mfilename), conf.secondary_sources.geometry); end % densely sampled SSD for calculating the aliasing frequency -conf.secondary_sources.number = 2048; +conf.secondary_sources.number = 256; x0 = secondary_source_positions(conf); -x0 = secondary_source_selection(x0,xs,src); % a_S(x0) >= 0 -x0(:,7) = deltax; % correct sampling +x0(:,7) = deltax0; % correct sampling -% local wavenumber of virtual sound field at x0 +% normalised local wavenumber vector of virtual sound field kS(x0) kSx0 = local_wavenumber_vector(x0(:,1:3), xs, src); switch method case {'wfs', 'WFS'} - % infinite control area - minmax_kSt_fun = @(x0p) minmax_kt_circle(x0p,[0,0,0],Inf); % handle - f = aliasing_extended_control(x0, kSx0, x, minmax_kGt_fun, ... - minmax_kSt_fun, conf); + M = Inf; % WFS does not have any bandwidth limitation + xc = [0,0,0]; % no effect since M = inf + Npw = Inf; % WFS does not involve a discretised Plane Wave Decomposition + f = aliasing_extended_sbl(x0,kSx0,x,minmax_kGt_fun,xc,M,Npw,conf); case {'nfchoa', 'HOA'} % circular control area at array center with Rl = M/k - M = nfchoa_order(N0,conf); - xc = conf.secondary_sources.center; - f = aliasing_extended_modal(x0, kSx0, x, minmax_kGt_fun, xc, M, conf); + M = nfchoa_order(N0,conf); % + xc = conf.secondary_sources.center; % expansion around the array center + Npw = Inf; % HOA does not involve a discretised Plane Wave Decomposition + f = aliasing_extended_sbl(x0,kSx0,x,minmax_kGt_fun,xc,M,Npw,conf); case {'localwfs-sbl', 'LWFS-SBL'} % circular control area at xref with Rl = M/k - xc = conf.xref; M = conf.localwfs_sbl.order; - f = aliasing_extended_modal(x0, kSx0, x, minmax_kGt_fun, xc, M, conf); + xc = conf.xref; + Npw = conf.localwfs_sbl.Npw; + f = aliasing_extended_sbl(x0,kSx0,x,minmax_kGt_fun,xc,M,Npw,conf); case {'localwfs-vss', 'LWFS-VSS'} + + % prepare virtual secondary sources + Nv = conf.localwfs_vss.number; switch conf.localwfs_vss.geometry case {'circle', 'circular'} - xc = conf.localwfs_vss.center; - Rc = conf.localwfs_vss.size / 2; - minmax_kSt_fun = @(x0p) minmax_kt_circle(x0p,xc,Rc); % handle + Rv = conf.localwfs_vss.size / 2; + deltaxv = 2*pi*Rv/Nv; + case 'linear' + Lv = conf.localwfs_vss.size; + deltaxv = Lv/(Nv-1); otherwise - error('%s: unsupported geometry (%s) for LWFS-VSS', ... + error('%s: unsupported geometry (%s) for virtual secondary sources', ... upper(mfilename), conf.localwfs_vss.geometry); end - f = aliasing_extended_control(x0, kSx0, x, minmax_kGt_fun, ... - minmax_kSt_fun, conf); + + % densely sampled virtual SSD for calculating the aliasing frequency + virtualconf = conf; + virtualconf.secondary_sources.size = conf.localwfs_vss.size; + virtualconf.secondary_sources.center = conf.localwfs_vss.center; + virtualconf.secondary_sources.geometry = conf.localwfs_vss.geometry; + virtualconf.secondary_sources.number = 256; + + xv = secondary_source_positions(virtualconf); + xv(:,7) = deltaxv; + + % local wavenumber of virtual sound field at xv + kSxv = local_wavenumber_vector(xv(:,1:3), xs, src); + + f = aliasing_extended_vss(x0,xv,kSxv,x,minmax_kGt_fun,conf); otherwise error('%s: unsupported SFS method (%s)', ... upper(mfilename), method) diff --git a/SFS_aliasing/aliasing_extended_modal.m b/SFS_aliasing/minmax_kt_line.m similarity index 51% rename from SFS_aliasing/aliasing_extended_modal.m rename to SFS_aliasing/minmax_kt_line.m index 64f17d96..25967bf5 100644 --- a/SFS_aliasing/aliasing_extended_modal.m +++ b/SFS_aliasing/minmax_kt_line.m @@ -1,25 +1,4 @@ -function f = aliasing_extended_modal(x0, kSx0, x, minmax_kGt_fun, xc, M, conf) -%ALIASING_EXTENDED_MODAL aliasing frequency for an extended listening area for -%an circular control area at xc with R=M/k where synthesis is focused on. -% -% Usage: f = aliasing_extended_control(x0, kSx0, x, minmax_kGt_fun, minmax_kSt_fun, conf) -% -% Input options: -% x0 - position, direction, and sampling distance of -% secondary sources [N0x7] / m -% kSx0 - normalised local wavenumber vector of virtual sound -% field [N0x3] -% x - position for which aliasing frequency is calculated -% [Nx3] -% minmax_kGt_fun - function handle to determine the extremal value of -% the tangential component of k_G(x-x0) -% [kGtmin, kGtmax] = minmax_kGt_fun(x0) -% xc - center of circular control area -% M - modal order which defines the radius R=M/k -% conf - configuration struct (see SFS_config) -% -% Output parameters: -% f - aliasing frequency [Nx1] +function [kGtmin, kGtmax] = minmax_kt_line(x0, xl, Ll, alphal) % %***************************************************************************** @@ -51,31 +30,61 @@ % http://sfstoolbox.org sfstoolbox@gmail.com * %***************************************************************************** -phik = cart2pol(kSx0(:,1),kSx0(:,2)); % azimuth angle of kS(x0) -phin0 = cart2pol(x0(:,4),x0(:,5)); % azimuth angle of normal vector n0(x0) +%% ===== Main ============================================================ + +% intersection of: +% g1: x = xl + gamma * nl -Ll/2 <= gamma <= Ll/2 +% g2 (n0-axis): 0 = (x - x0)^T * n0 +% +% -> gamma = ((x0 - xl)^T * n0) / (nl^T * n0) + +xl = bsxfun(@minus, xl, x0(:,1:3)); % shift xl about x0 +nl = [cos(alphal), sin(alphal), 0]; % direction vector of line +n0 = x0(:,4:6); +phin0 = cart2pol(n0(:,1),n0(:,2)); -% secondary source selection -select = cos(phin0 - phik) >= 0; -x0 = x0(select,:); -% k_St(x0) (tangential component of k_S(x0) ) -kSt = sin(phin0(select) - phik(select)); +gamma = -(xl(:,1).*n0(:,1) + xl(:,2).*n0(:,2))./(n0*nl.'); -% sampling distance -deltax0 = abs(x0(:,7)); +select = gamma <= 0 | gamma > Ll/2; +gamma1(select) = +Ll/2; +gamma1(~select) = gamma; -f = inf(size(x,1), 1); -for xdx = 1:size(x,1); - % mininum and maximum values of k_Gt(x - x_0) - % (tangential component of k_G(x-x0)) - [kGtmin, kGtmax] = minmax_kGt_fun(x0,x(xdx,:)); - % aliasing frequency for x0 - f0 = conf.c./(deltax0.*max(abs(kSt-kGtmin),abs(kSt-kGtmax))); - % radius of local region at f0 - rM = M.*conf.c./(2.*pi.*f0); - % mininum and maximum values of kSt(x_0) at f0 - [kStmin, kStmax] = minmax_kt_circle(x0, xc, rM); - select = kSt > kStmin & kSt < kStmax; - if sum(select) ~= 0 - f(xdx) = min(f0(select)); - end +select = gamma < -Ll/2 | gamma > 0; +gamma2(select) = -Ll/2; +gamma2(~select) = gamma; + +select = ~isinf(gamma1); +if any(select) + x1(select,:) = bsxfun(@plus, xl(select,:), gamma1(select).*nl); +elseif any(~select) + x1(~select,:) = sign(gamma1(~select))*nl; end + +select = ~isinf(gamma2); +if any(select) + x2(select,:) = bsxfun(@plus, xl(select,:), gamma2(select).*nl); +elseif any(~select) + x2(~select,:) = sign(gamma2(~select))*nl; +end + +if abs(x1) < 1e-10 + x1 = x2; + degenerated = true; +else + degenerated = false; +end +if abs(x2) < 1e-10 + x2 = x1; + if degenerated + error('degenerated case'); + end +end + +phi1 = cart2pol(x1(:,1),x1(:,2)); +phi2 = cart2pol(x2(:,1),x2(:,2)); + +kGt1 = sin(phin0 - phi1); +kGt2 = sin(phin0 - phi2); + +kGtmin = min(kGt1,kGt2); +kGtmax = max(kGt1,kGt2); From ea8cd4ec7b7fa4b5cd41c319d91b522efaa0c630 Mon Sep 17 00:00:00 2001 From: Fiete Winter Date: Thu, 15 Nov 2018 17:48:04 +0100 Subject: [PATCH 3/3] add test function for aliasing frequency --- validation/test_aliasing.m | 225 +++++++++++++++++++++++++++++++++++++ 1 file changed, 225 insertions(+) create mode 100644 validation/test_aliasing.m diff --git a/validation/test_aliasing.m b/validation/test_aliasing.m new file mode 100644 index 00000000..5de09851 --- /dev/null +++ b/validation/test_aliasing.m @@ -0,0 +1,225 @@ +function status = test_aliasing(modus) +%TEST_ALIASING tests the prediction of aliasing frequency +% +% Usage: status = test_driving_functions(modus) +% +% Input parameters: +% modus - 0: numerical +% 1: visual +% +% Output parameters: +% status - true or false + +%***************************************************************************** +% The MIT License (MIT) * +% * +% Copyright (c) 2010-2018 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 * +%***************************************************************************** + + +% TODO: add mode to save data as reference data +status = false; + + +%% ===== Checking of input parameters =================================== +nargmin = 1; +nargmax = 1; +narginchk(nargmin,nargmax); + + +%% ===== Configuration =================================================== +conf = SFS_config; +conf.plot.usenormalisation = false; +conf.usetapwin = false; + +conf.nfchoa.order = 25; + +conf.localwfs_sbl.order = 25; +conf.localwfs_sbl.Npw = 1024; +conf.localwfs_sbl.fc = 500; + +conf.localwfs_vss.geometry = 'circular'; +conf.localwfs_vss.size = 1.0; +conf.localwfs_vss.number = 16; +conf.localwfs_vss.usetapwin = false; +conf.localwfs_vss.consider_secondary_sources = false; +conf.localwfs_vss.consider_target_field = false; + +N0gt = 256; +Nvgt = 256; +f = 1500; + +% test scenarios +scenarios = { ... + 'WFS' , 'linear', 'pw', [ 0.5 -1.0 0.0]; ... + 'WFS' , 'linear', 'ps', [ 0.0 3.0 0.0]; ... + 'WFS' , 'linear', 'fs', [ 0.0 -1.0 0.0 0.0 -1.0 0.0]; ... + 'LWFS-SBL' , 'linear', 'pw', [ 0.5 -1.0 0.0]; ... + 'LWFS-SBL' , 'linear', 'ps', [ 0.0 1.0 0.0]; ... + 'LWFS-VSS' , 'linear', 'pw', [ 0.5 -1.0 0.0]; ... + 'LWFS-VSS' , 'linear', 'ps', [ 0.0 1.0 0.0]; ... + 'WFS' , 'circular', 'pw', [ 0.5 0.5 0.0]; ... + 'WFS' , 'circular', 'ps', [ 0.0 2.5 0.0]; ... + 'WFS' , 'circular', 'fs', [ 0.0 0.5 0.0 0.0 -1.0 0.0]; ... + 'LWFS-SBL' , 'circular', 'pw', [ 0.5 0.5 0.0]; ... + 'LWFS-SBL' , 'circular', 'ps', [ 0.0 2.5 0.0]; ... + 'LWFS-VSS' , 'circular', 'pw', [ 0.5 0.5 0.0]; ... + 'LWFS-VSS' , 'circular', 'ps', [ 0.0 2.5 0.0]; ... + 'HOA' , 'circular', 'pw', [ 0.5 0.5 0.0]; ... + 'HOA' , 'circular', 'ps', [ 0.0 2.5 0.0]; ... +}; + +% Start testing +for ii=1:size(scenarios) + + method = scenarios{ii,1}; + conf.secondary_sources.geometry = scenarios{ii,2}; + + % set scenario + switch scenarios{ii,2} + case 'linear' + X = [-1.55 1.55]; + Y = [-2.85 0.15]; + Z = 0; + conf.xref = [0.5 -1.5 0]; + conf.secondary_sources.size = 3; + conf.secondary_sources.number = 16; + case {'circular', 'circle'} + X = [-1.55 1.55]; + Y = [-1.55 1.55]; + Z = 0; + conf.xref = [0.5, -0.75, 0]; + conf.secondary_sources.size = 3; + conf.secondary_sources.number = 36; + end + xref = conf.xref; + conf.localwfs_vss.center = conf.xref; + + + % evaluation coordinates for aliasing frequency + conf.resolution = 51; + [xS,yS,~] = xyz_grid(X,Y,Z,conf); + xvec = [xS(:), yS(:)]; + xvec(:,3) = 0; + + src = scenarios{ii,3}; + if strcmp('fs', src) + srcgt = 'ps'; + else + srcgt = src; + end + + xs = scenarios{ii,4}; + + fS = aliasing_frequency(method, xs, src, 'position', xvec, conf); + fS = reshape(fS, size(xS)); + + if modus + suffix = sprintf(... + 'virtual %s, %s array (L = %d)\n%s', src, ... + conf.secondary_sources.geometry, conf.secondary_sources.number, method); + + % reference amplitude + g = sound_field_mono(xref(1),xref(2),xref(3),[xs(1:3),0,1,0,1],srcgt,1,... + f,conf); + + % compute synthesised sound field and aliasing error + conf.resolution = 150; + gtconf = conf; + gtconf.secondary_sources.number = N0gt; + gtconf.localwfs_vss.number = Nvgt; + switch method + case 'WFS' + P = sound_field_mono_wfs(X,Y,Z,xs,src,f,conf); + Pgt = sound_field_mono_wfs(X,Y,Z,xs,src,f,gtconf); + case 'HOA' + P = sound_field_mono_nfchoa(X,Y,Z,xs,src,f,conf); + Pgt = sound_field_mono_nfchoa(X,Y,Z,xs,src,f,gtconf); + suffix = sprintf('%s (M = %d)', suffix, conf.nfchoa.order); + case 'LWFS-SBL' + P = sound_field_mono_localwfs_sbl(X,Y,Z,xs,src,f,conf); + Pgt = sound_field_mono_localwfs_sbl(X,Y,Z,xs,src,f,gtconf); + suffix = sprintf('%s (M = %d, Npw = %d)', suffix, ... + conf.localwfs_sbl.order, conf.localwfs_sbl.Npw); + case 'LWFS-VSS' + P = sound_field_mono_localwfs_vss(X,Y,Z,xs,src,f,conf); + Pgt = sound_field_mono_localwfs_vss(X,Y,Z,xs,src,f,gtconf); + + end + + x0 = secondary_source_positions(conf); + + % plot aliasing frequency as function of x + figure; + subplot(2,2,1); + imagesc(X,Y,fS) + hold on + contour(xS,yS,fS,[f,f],'k'); + hold off + set(gca, 'YDir', 'normal', 'CLim', [f-1000 f+1000]); + title(sprintf('aliasing frequency\n %s', suffix)); + + % plot synthesised sound field + subplot(2,2,2); + imagesc(X,Y,real(P)./abs(g)); + hold on + contour(xS,yS,fS,[f,f],'k'); + hold off + set(gca, 'YDir', 'normal', 'CLim', [-1 1]); + title(sprintf('synthesised sound field (f=%2.2f Hz)\n %s', f, suffix)); + + % plot ground truth sound field + subplot(2,2,3); + imagesc(X,Y,real(Pgt)./abs(g)); + set(gca, 'YDir', 'normal', 'CLim', [-1 1]); + title(sprintf('ground truth sound field (f=%2.2f Hz)\n %s', f, suffix)); + + % plot aliasing error + subplot(2,2,4); + imagesc(X,Y,db((P-Pgt)./Pgt)); + hold on + contour(xS,yS,fS,[f,f],'k'); + hold off + set(gca, 'YDir', 'normal', 'CLim', [-40 10]); + title(sprintf('aliasing error (f=%2.2f Hz)\n %s', f, suffix)); + + for idx = 1:4 + subplot(2,2,idx); + xlabel('x / m'); + ylabel('y / m'); + xlim(X); + ylim(Y); + hold on + draw_loudspeakers(x0,[1 1 0], conf); + hold off + set_colorbar(conf); + colorbar; + end + end +end + + +status = true;