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_vss.m b/SFS_aliasing/aliasing_extended_vss.m new file mode 100644 index 00000000..6b708c2a --- /dev/null +++ b/SFS_aliasing/aliasing_extended_vss.m @@ -0,0 +1,122 @@ +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 +% +% 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(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) + +% 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_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 new file mode 100644 index 00000000..e5b7c31d --- /dev/null +++ b/SFS_aliasing/aliasing_frequency.m @@ -0,0 +1,137 @@ +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 * +%***************************************************************************** + +%% ===== Checking input parameters ======================================= +nargmin = 6; +nargmax = 7; +narginchk(nargmin,nargmax); +isargchar(method, area, src); +isargxs(xs); +ismatrix(x); + +switch area + case 'position' + % 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 + 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; + deltax0 = 2*pi*R/N0; + case 'linear' + L = conf.secondary_sources.size; + 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 = 256; +x0 = secondary_source_positions(conf); +x0(:,7) = deltax0; % correct sampling + +% normalised local wavenumber vector of virtual sound field kS(x0) +kSx0 = local_wavenumber_vector(x0(:,1:3), xs, src); + +switch method + case {'wfs', 'WFS'} + 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; % 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 + M = conf.localwfs_sbl.order; + 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'} + 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 virtual secondary sources', ... + upper(mfilename), conf.localwfs_vss.geometry); + end + + % 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) +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_aliasing/minmax_kt_line.m b/SFS_aliasing/minmax_kt_line.m new file mode 100644 index 00000000..25967bf5 --- /dev/null +++ b/SFS_aliasing/minmax_kt_line.m @@ -0,0 +1,90 @@ +function [kGtmin, kGtmax] = minmax_kt_line(x0, xl, Ll, alphal) +% + +%***************************************************************************** +% 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 * +%***************************************************************************** + +%% ===== 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)); + +gamma = -(xl(:,1).*n0(:,1) + xl(:,2).*n0(:,2))./(n0*nl.'); + +select = gamma <= 0 | gamma > Ll/2; +gamma1(select) = +Ll/2; +gamma1(~select) = gamma; + +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); 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']); 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;