diff --git a/SFS_general/besselh_derived.m b/SFS_general/besselh_derived.m new file mode 100644 index 00000000..34aca0f8 --- /dev/null +++ b/SFS_general/besselh_derived.m @@ -0,0 +1,65 @@ +function out = besselh_derived(nu,k,z) +% BESSELH_DERIVED derivative of cylindrical hankel function of k kind of order nu, and argument z +% +% Usage: out = besselh_derived(nu,k,z) +% +% Input parameters: +% nu - order of bessel function +% z - argument of bessel function +% +% Output parameters: +% out - value of bessel function at point z +% +% BESSELH_DERIVED(nu,z) derivation of cylindrical hankel function of +% order nu, k kind, and argument z +% +% References: +% (4.1-51) in Ziomek (1995) - "Fundamentals of acoustic field theory +% and space-time signal processing" +% +% see also: besselh + +%***************************************************************************** +% Copyright (c) 2010-2014 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2014 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + + +%% ===== Checking input parameters ======================================= +nargmin = 3; +nargmax = 3; +narginchk(nargmin,nargmax); +isargscalar(nu) +isargscalar(k) +isargnumeric(z) + + +%% ===== Computation ===================================================== +out = 0.5*(besselh(nu-1,k,z) - besselh(nu+1,k,z)); \ No newline at end of file diff --git a/SFS_general/besselj_derived.m b/SFS_general/besselj_derived.m new file mode 100644 index 00000000..5337ed53 --- /dev/null +++ b/SFS_general/besselj_derived.m @@ -0,0 +1,62 @@ +function out = besselj_derived(nu,z) +% BESSELJ_DERIVED derivative of cylindrical bessel function of order nu, and argument z +% +% Usage: out = besselh_derived(nu,z) +% +% Input parameters: +% nu - order of bessel function +% z - argument of bessel function +% +% Output parameters: +% out - value of bessel function at point z +% +% BESSELJ_DERIVED(nu,z) derivation of cylindrical bessel function of +% order nu, and argument z +% +% References: +% (4.1-51) in Ziomek (1995) - "Fundamentals of acoustic field theory +% and space-time signal processing" +% +% see also: besselj + +%***************************************************************************** +% Copyright (c) 2010-2014 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2014 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking input parameters ======================================= +nargmin = 2; +nargmax = 2; +narginchk(nargmin,nargmax); +isargscalar(nu) +isargnumeric(z) + +%% ===== Computation ===================================================== +out = 0.5*(besselj(nu-1,z) - besselj(nu+1,z)); \ No newline at end of file diff --git a/SFS_general/sphbesselh_derived.m b/SFS_general/sphbesselh_derived.m new file mode 100644 index 00000000..d10d98ba --- /dev/null +++ b/SFS_general/sphbesselh_derived.m @@ -0,0 +1,65 @@ +function out = sphbesselh_derived(nu,k,z) +% SPHBESSELH_DERIVED derivative of spherical hankel function of k kind of order nu, and argument z +% +% Usage: out = sphbesselh_derived(nu,k,z) +% +% Input parameters: +% nu - order of bessel function +% z - argument of bessel function +% +% Output parameters: +% out - value of bessel function at point z +% +% SPHBESSELH_DERIVED(nu,z) derivation of spherical hankel function of +% order nu, k kind, and argument z +% +% References: +% (4.1-51) in Ziomek (1995) - "Fundamentals of acoustic field theory +% and space-time signal processing" +% +% see also: sphbesselh + +%***************************************************************************** +% Copyright (c) 2010-2014 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2014 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + + +%% ===== Checking input parameters ======================================= +nargmin = 3; +nargmax = 3; +narginchk(nargmin,nargmax); +isargscalar(nu) +isargscalar(k) +isargnumeric(z) + + +%% ===== Computation ===================================================== +out = 1/(2*nu+1) * (nu*sphbesselh(nu-1,k,z) - (nu+1)*sphbesselh(nu+1,k,z)); \ No newline at end of file diff --git a/SFS_general/sphbesselj_derived.m b/SFS_general/sphbesselj_derived.m new file mode 100644 index 00000000..f71d2ed8 --- /dev/null +++ b/SFS_general/sphbesselj_derived.m @@ -0,0 +1,64 @@ +function out = sphbesselj_derived(nu,z) +% SPHBESSELJ_DERIVED derivative of spherical bessel function of first kind of order nu, and argument z +% +% Usage: out = sphbesselj_derived(nu,z) +% +% Input parameters: +% nu - order of bessel function +% z - argument of bessel function +% +% Output parameters: +% out - value of bessel function at point z +% +% SPHBESSELJ_DERIVED(nu,z) derivation of spherical bessel function of +% order nu, first type, and argument z +% +% References: +% (4.1-51) in Ziomek (1995) - "Fundamentals of acoustic field theory +% and space-time signal processing" +% +% see also: sphbesselj + +%***************************************************************************** +% Copyright (c) 2010-2014 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2014 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + + +%% ===== Checking input parameters ======================================= +nargmin = 2; +nargmax = 2; +narginchk(nargmin,nargmax); +isargscalar(nu) +isargnumeric(z) + + +%% ===== Computation ===================================================== +out = 1/(2*nu+1) * (nu*sphbesselj(nu-1,z) - (nu+1)*sphbesselj(nu+1,z)); \ No newline at end of file diff --git a/SFS_general/sphbessely_derived.m b/SFS_general/sphbessely_derived.m new file mode 100644 index 00000000..685577b3 --- /dev/null +++ b/SFS_general/sphbessely_derived.m @@ -0,0 +1,64 @@ +function out = sphbessely_derived(nu,z) +% SPHBESSELY_DERIVED derivative of spherical bessel function of second kind of order nu, and argument z +% +% Usage: out = sphbessely_derived(nu,z) +% +% Input parameters: +% nu - order of bessel function +% z - argument of bessel function +% +% Output parameters: +% out - value of bessel function at point z +% +% SPHBESSELY_DERIVED(nu,z) derivation of spherical bessel function of +% order nu, second type, and argument z +% +% References: +% (4.1-51) in Ziomek (1995) - "Fundamentals of acoustic field theory +% and space-time signal processing" +% +% see also: sphbesselj + +%***************************************************************************** +% Copyright (c) 2010-2014 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2014 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + + +%% ===== Checking input parameters ======================================= +nargmin = 2; +nargmax = 2; +narginchk(nargmin,nargmax); +isargscalar(nu) +isargnumeric(z) + + +%% ===== Computation ===================================================== +out = 1/(2*nu+1) * (nu*sphbessely(nu-1,z) - (nu+1)*sphbessely(nu+1,z)); \ No newline at end of file diff --git a/SFS_helper/isargsquaredinteger.m b/SFS_helper/isargsquaredinteger.m new file mode 100644 index 00000000..643bf75c --- /dev/null +++ b/SFS_helper/isargsquaredinteger.m @@ -0,0 +1,55 @@ +function isargsquaredinteger(varargin) +%ISARGSQUAREDINTEGER tests if the given arg is a scalar, which is the +%square of an integer +% +% +% Usage: isargsquaredinteger(arg1,arg2,...) +% +% Input options: +% args - list of args +% +% ISARGPOSITIVESCALAR(args) tests if all given args are a positive +% scalar and returns an error otherwise. +% +% See also: isargscalar, isargnegativescalar + +%***************************************************************************** +% Copyright (c) 2010-2015 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2015 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + + +%% ===== Checking for scalar ============================================= +for ii = 1:nargin + if ( ~isnumeric(varargin{ii}) || ~isscalar(varargin{ii}) || ... + mod(sqrt(varargin{ii}), 1) ~= 0 ) + error('%s need to be the square of an integer.',inputname(ii)); + end +end diff --git a/SFS_monochromatic/cht/driving_function_mono_nfchoa_cht_circexp.m b/SFS_monochromatic/cht/driving_function_mono_nfchoa_cht_circexp.m new file mode 100644 index 00000000..3b23045e --- /dev/null +++ b/SFS_monochromatic/cht/driving_function_mono_nfchoa_cht_circexp.m @@ -0,0 +1,112 @@ +function Dm = driving_function_mono_nfchoa_cht_circexp(Pm, f, conf) +%DRIVING_FUNCTION_MONO_NFCHOA_CHT_CIRCEXP computes the circular harmonics +%transform of nfchoa driving functions for a sound field expressed by regular +%circular expansion coefficients. +% +% Usage: D = driving_function_mono_nfchoa_cht_circexp(Pm, f, conf) +% +% Input parameters: +% Pm - regular circular expansion coefficients of virtual +% sound field [n x m] +% f - frequency in Hz [m x 1] or [1 x m] +% conf - configuration struct (see SFS_config) +% +% Output parameters: +% Dm - circular harmonics transform of driving function signal +% [n x m] +% +% DRIVING_FUNCTION_MONO_NFCHOA_CHT_CIRCEXP(Pm, f, conf) returns circular +% harmonics transform of the NFCHOA driving function for a virtual sound +% expressed by regular circular expansion coefficients and the frequency f. + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 3; +nargmax = 3; +narginchk(nargmin,nargmax); +isargmatrix(Pm); +isargvector(f); +isargstruct(conf); +if size(Pm,2) ~= length(f) + error( '%s:number of rows in %s have to match length of %s', ... + upper(mfilename), inputname(1), inputname(2) ); +end + +%% ===== Configuration ================================================== +c = conf.c; +R0 = conf.secondary_sources.size / 2; + +%% ===== Variables ====================================================== +Nce = (size(Pm, 1)-1)/2; + +% frequency depended stuff +omega = 2*pi*row_vector(f); % [1 x Nf] +k = omega./c; % [1 x Nf] +kR0 = k.*R0; % [1 x Nf] + +%% ===== Computation ==================================================== +% Calculate the circular harmonics of driving function +% +% Regular circular expansion of the sound field: +% \~~ oo +% P(x,w) = > P J (kr) . e^(+j m phi) +% /__ m=-oo m m +% +% and 3D free field Green's Function: +% \~~ oo +% G (x,w) = > G J (kr) . e^(+j m phi) +% ls /__ m=-oo m m +% +% with the regular expansion coefficients of Green's Function +% (see Gumerov2004, eq. 3.2.2): +% m -i (2) +% G = ----- . H (kr0) +% n 4 m + +% P +% 1 m +% D = ------ -------- +% m 2pi r0 G +% m + +Dm = zeros(size(Pm)); +l = 0; +for m=-Nce:Nce + l = l+1; + Dm(l,:) = Pm(l,:) ./ besselh(m,2,kR0); +end +% factor from expansion of 2D free field Green's Function +Dm = Dm./(-1j/4); +% normalization due to size of circular array +Dm = Dm./2*pi*R0; diff --git a/SFS_monochromatic/cht/driving_function_mono_nfchoa_cht_sphexp.m b/SFS_monochromatic/cht/driving_function_mono_nfchoa_cht_sphexp.m new file mode 100644 index 00000000..ca0761f2 --- /dev/null +++ b/SFS_monochromatic/cht/driving_function_mono_nfchoa_cht_sphexp.m @@ -0,0 +1,163 @@ +function Dm = driving_function_mono_nfchoa_cht_sphexp(Pnm, f, conf) +%DRIVING_FUNCTION_MONO_NFCHOA_CHT_SPHEXP computes the circular harmonics +%transform of nfchoa driving functions for a sound field expressed by regular +%spherical expansion coefficients. +% +% Usage: D = driving_function_mono_nfchoa_cht_sphexp(Pnm, f, conf) +% +% Input parameters: +% Pnm - regular spherical expansion coefficients of virtual +% sound field [nxm] +% f - frequency in Hz [m x 1] or [1 x m] +% conf - configuration struct (see SFS_config) +% +% Output parameters: +% Dm - regular circular harmonics transform of driving +% function signal +% +% DRIVING_FUNCTION_MONO_NFCHOA_CHT_SPHEXP(Pnm, f, conf) returns circular +% harmonics transform of the NFCHOA driving function for a virtual sound +% expressed by regular spherical expansion coefficients and the frequency f. + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 3; +nargmax = 3; +narginchk(nargmin,nargmax); +isargmatrix(Pnm); +isargsquaredinteger(size(Pnm,1)); +isargvector(f); +isargstruct(conf); +if size(Pnm,2) ~= length(f) + error( '%s:number of rows in %s have to match length of %s', ... + upper(mfilename), inputname(1), inputname(2) ); +end + +%% ===== Configuration ================================================== +c = conf.c; +Xc = conf.secondary_sources.center; +R0 = conf.secondary_sources.size / 2; +xref = conf.xref - Xc; +rref = norm( xref ); % reference radius +thetaref = asin( xref(3)./rref); % reference elevation angle + +%% ===== Variables ====================================================== +Nse = sqrt(size(Pnm,1))-1; + +% frequency depended stuff +omega = 2*pi*row_vector(f); % [1 x Nf] +k = omega./c; % [1 x Nf] +kR0 = k.*R0; % [1 x Nf] +krref = k.*rref; % [1 x Nf] + +%% ===== Computation ==================================================== +% Calculate the circular harmonics of driving function +% +% Regular spherical expansion of the sound field: +% \~~ N \~~ n m m +% P(x,w) = > > P j (kr) . Y (theta, phi) +% /__ n=0 /__ m=-n n n n +% +% and 3D free field Green's Function: +% \~~ oo \~~ n m m +% G (x0,f) = > > G . j (kr) . Y (theta, phi) +% ps /__ n=0 /__ m=-n n n n +% +% with the regular expansion coefficients of Green's Function +% (see Gumerov2004, eq. 3.2.2): +% m (2) -m +% G = -i . k . h (kr0) Y (pi/2, 0) +% n n n + +% initialize empty driving signal +Dm = zeros(2*Nse+1, size(Pnm,2)); +l = 0; +if (rref == 0) + % --- Xref == Xc ------------------------------------------------- + % m + % P + % 1 |m| + % D = ------ -------- + % m 2pi r0 m + % G + % |m| + for m=-Nse:Nse + l = l+1; + v = sphexp_index(m); % n = abs(m) + + Dm(l,:) = Pnm(v,:) ./ ... + ( -1i.*k.* sphbesselh(abs(m),2,kR0).*sphharmonics(abs(m),-m, 0, 0) ); + end +else + % --- Xref ~= Xc -------------------------------------------------- + % + % if the reference position is not in the middle of the array, + % things get a 'little' more complicated + % __ + % \ m m + % /__ P j (k rref) Y (thetaref, 0) + % n 1 l=|m| l l l + % D = ------- ------------------------------------ + % m 2pi r0 __ + % \ m m + % /__ G j (k rref) Y (thetaref, 0) + % l=|m| n l l + % + + % save some intermediate results + hn = zeros(size(Pnm)); + jn = hn; + for n=0:Nse + hn(n+1,:) = sphbesselh(n,2,kR0); % [1 x Nf] + jn(n+1,:) = sphbesselj(n,krref); % [1 x Nf] + end + + for m=-Nse:Nse + Pm = zeros(1,length(f)); + Gm = Pm; + for n=abs(m):Nse + v = sphexp_index(m,n); + % + factor = jn(n+1,:) .* sphharmonics(n, -m, thetaref, 0); + % numerator + Pm = Pm + Pnm(v,:) .* factor; + % denominator + Gm = Gm + (-1i*k) .* hn(n+1,:) .* sphharmonics(n, -m, 0, 0) .* factor; + end + l = l+1; + Dm(v,:) = Pm./Gm; + end +end +% normalization due to size of circular array +Dm = Dm./(2*pi*R0); \ No newline at end of file diff --git a/SFS_monochromatic/cht/sound_field_mono_cht.m b/SFS_monochromatic/cht/sound_field_mono_cht.m new file mode 100644 index 00000000..cfded76c --- /dev/null +++ b/SFS_monochromatic/cht/sound_field_mono_cht.m @@ -0,0 +1,124 @@ +function [P, x, y, z] = sound_field_mono_cht(X,Y,Z,Dm,f,conf) +%SOUND_FIELD_MONO_CHT simulates a sound field reproduced by a circular +%secondary source distribution given with the circular harmonics transform +%of the driving function +% +% Usage: [P, x, y, z] = sound_field_mono_cht(X,Y,Z,Dm,f,conf) +% +% Input parameters: +% X - x-axis / m; single value or [xmin,xmax] or nD-array +% Y - y-axis / m; single value or [ymin,ymax] or nD-array +% Z - z-axis / m; single value or [zmin,zmax] or nD-array +% Dm - circular harmonics transform of nfchoa driving function +% f - frequency in Hz +% conf - configuration struct (see SFS_config) +% +% Output parameters: +% P - resulting soundfield +% +% SOUND_FIELD_MONO_CHT(X,Y,Z,Dm,f,conf) computes the sound field reproduced +% by a continuous, circular secondary source distribution driven by a driving +% function whose circular harmonics transform is given as Dm. +% +% see also: circbasis_mono_grid, sound_field_mono_circexp + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 6; +nargmax = 6; +narginchk(nargmin,nargmax); +isargvector(Dm); +isargnumeric(X,Y,Z); +isargpositivescalar(f); +isargstruct(conf); + +%% ===== Configuration ================================================== +Xc = conf.secondary_sources.center; +R0 = conf.secondary_sources.size / 2; +useplot = conf.plot.useplot; +dimension = conf.dimension; + +%% ===== Variables ====================================================== +% select secondary source type (line source or point source) +switch dimension + case '2D' + exp_func = @(mo) circexp_mono_cht(Dm,mo,f,conf); + sound_field_func = @(x,y,z,Am,mo) sound_field_mono_circexp(x, y, z, Am, ... + mo, f, Xc, conf); + case '2.5D' + exp_func = @(mo) sphexp_mono_cht(Dm,mo,f,conf); + sound_field_func = @(x,y,z,Amn,mo) sound_field_mono_sphexp(x, y, z, Amn, ... + mo, f, Xc, conf); +end + +%% ===== Computation ==================================================== +[x,y,z] = xyz_grid(X,Y,Z,conf); +% find coordinates, which are inside and outside the loudspeaker array +select = sqrt((x-Xc(1)).^2 + (y-Xc(2)).^2 + (z-Xc(3)).^2) <= R0; + +if (numel(x) == 1) x = repmat(x, size(select)); end +if (numel(y) == 1) y = repmat(y, size(select)); end +if (numel(z) == 1) z = repmat(z, size(select)); end + +P = zeros(size(x)); + +% regular (interior) domain inside the loudspeaker array +if any(select(:)) + Pnm = exp_func('R'); + P(select) = sound_field_func(x(select), y(select), z(select), Pnm, 'R'); +end + +% singular (exterior) domain outside the loudspeaker array +if any(~select(:)) + Pnm = exp_func('S'); + if sum( ~select(:) ) == 2 + % this handle cases, where x(~select) would only contain 2 entries, which + % would be interpreted as a range by the sound_field_... function + xtmp = x(~select); + ytmp = y(~select); + ztmp = z(~select); + Ptmp(1) = sound_field_func(xtmp(1), ytmp(1), ztmp(1), Pnm, 'S'); + Ptmp(2) = sound_field_func(xtmp(2), ytmp(2), ztmp(2), Pnm, 'S'); + P(~select) = [Ptmp(1); Ptmp(2)]; + else + P(~select) = sound_field_func(x(~select), y(~select), z(~select), Pnm, 'S'); + end +end + +%% ===== Plotting ======================================================= +if (nargout==0 || useplot) + plot_sound_field(P,X,Y,Z,[],conf); +end + +end \ No newline at end of file diff --git a/SFS_monochromatic/circexp/circbasis_mono.m b/SFS_monochromatic/circexp/circbasis_mono.m new file mode 100644 index 00000000..f1ca4094 --- /dev/null +++ b/SFS_monochromatic/circexp/circbasis_mono.m @@ -0,0 +1,93 @@ +function [Jn,H2n,Yn] = circbasis_mono(r,phi,Nce,k,conf) +%CIRCBASIS_MONO evaluates circular basis functions for given input arguments +% +% Usage: [Jn,H2n,Yn] = circbasis_mono(r,phi,Nce,k,conf) +% +% Input parameters: +% r - distance from z-axis in cylindrical coordinates +% phi - azimuth angle in cylindrical coordinates +% k - wave number +% conf - optional configuration struct (see SFS_config) +% +% Output parameters: +% Jn - cell array of cylindrical bessel functions +% H2n - cell array of cylindrical hankel functions of 2nd kind +% Yn - cell array of cylindrical harmonics +% +% CIRCBASIS_MONO(r, phi, Nce, k, conf) computes cylindrical basis functions +% for the given arguments r and phi. r and phi can be of arbitrary (but same) +% size. Output will be stored in cell arrays (one cell entry for each order) +% of length 2*Nce+1 . Each cell array entry contains a matrix of the same +% size as r and phi. +% +% References: +% Williams (1999) - "Fourier Acoustics", ACADEMIC PRESS +% +% see also: cylbasis_mono_XYZgrid besselj besselh + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 5; +nargmax = 5; +narginchk(nargmin,nargmax); +isargequalsize(r,phi); +isargscalar(k); +isargpositivescalar(Nce); + +%% ===== Configuration ================================================== +showprogress = conf.showprogress; + +%% ===== Computation ==================================================== +kr = k.*r; % argument of bessel functions + +L = 2*Nce + 1; +Jn = cell(L,1); +H2n = cell(L,1); +Yn = cell(L,1); + +l = 0; +for n=-Nce:0 + % negative n + l = l + 1; + Jn{l} = besselj(n,kr); + H2n{l} = Jn{l} - 1j*bessely(n,kr); + Yn{l} = exp(1j*n*phi); + % positive n + Jn{L-l+1} = (-1)^n*Jn{l}; + H2n{L-l+1} = (-1)^n*H2n{l}; + Yn{L-l+1} = conj(Yn{l}); + if showprogress, progress_bar(n+Nce,Nce); end % progress bar +end + +end diff --git a/SFS_monochromatic/circexp/circbasis_mono_grid.m b/SFS_monochromatic/circexp/circbasis_mono_grid.m new file mode 100644 index 00000000..3450be88 --- /dev/null +++ b/SFS_monochromatic/circexp/circbasis_mono_grid.m @@ -0,0 +1,91 @@ +function [jn,h2n,Ynm] = circbasis_mono_grid(X,Y,Z,Nce,f,xq,conf) +%CIRCBASIS_MONO_GRID evaluate spherical basis functions for given grid in +%cartesian coordinates +% +% Usage: [jn,h2n,Ynm] = circbasis_mono_grid(X,Y,Z,Nce,f,xq,conf) +% +% Input parameters: +% X - x-axis / m; single value or [xmin,xmax] or nD-array +% Y - y-axis / m; single value or [ymin,ymax] or nD-array +% Z - z-axis / m; single value or [zmin,zmax] or nD-array +% Nce - maximum order of circular basis functions +% f - frequency in Hz +% xq - center of coordinate system +% conf - configuration struct (see SFS_config) +% +% Output parameters: +% jn - cell array of cylindrical bessel functions +% h2n - cell array of cylindrical hankel functions of 2nd kind +% Ynm - cell array of cylindrical harmonics +% +% CIRCBASIS_MONO_GRID(X,Y,Z,f,xq,conf) computes circular basis functions +% for given grid in cartesian coordinates. This is a wrapper function for +% circbasis_mono. +% For the input of X,Y,Z (DIM as a wildcard) : +% * if DIM is given as single value, the respective dimension is +% squeezed, so that dimensionality of the simulated sound field P is +% decreased by one. +% * if DIM is given as [dimmin, dimmax], a linear grid for the +% respective dimension with a resolution defined in conf.resolution is +% established +% * if DIM is given as n-dimensional array, the other dimensions have +% to be given as n-dimensional arrays of the same size or as a single value. +% Each triple of X,Y,Z is interpreted as an evaluation point in an +% customized grid. For this option, plotting and normalisation is disabled. +% +% see also: circbasis_mono + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 7; +nargmax = 7; +narginchk(nargmin,nargmax); + +isargnumeric(X,Y,Z); +isargpositivescalar(f,Nce); +isargposition(xq); + +%% ===== Computation ==================================================== +[xx,yy] = xyz_grid(X,Y,Z,conf); + +k = 2*pi*f/conf.c; % wavenumber + +% shift coordinates to expansion coordinate +xx = xx-xq(1); +yy = yy-xq(2); +% coordinate transformation +r = sqrt(xx.^2 + yy.^2); +phi = atan2(yy,xx); + +[jn,h2n,Ynm] = circbasis_mono(r,phi,Nce,k,conf); diff --git a/SFS_monochromatic/circexp/circexp_convert_sphexp.m b/SFS_monochromatic/circexp/circexp_convert_sphexp.m new file mode 100644 index 00000000..747874b0 --- /dev/null +++ b/SFS_monochromatic/circexp/circexp_convert_sphexp.m @@ -0,0 +1,65 @@ +function Am = circexp_convert_sphexp(Anm) +%CIRCEXP_CONVERT_SPHEXP converts regular spherical expansion coefficients into +%regular circular expansion coefficients +% +% Usage: Am = circexp_convert_sphexp(Anm) +% +% Input parameters: +% Anm - regular spherical expansion coefficients +% +% Output parameters: +% Am - regular circular expansion coefficients +% +% References: +% Hahn, Winter, Spors (2016) - +% "Local Wave Field Synthesis by Spatial Band-limitation in the +% Circular/Spherical Harmonics Domain", 140th AES Convention + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 1; +nargmax = 1; +narginchk(nargmin,nargmax); +isargvector(Anm); +Nse = sqrt(size(Anm,1))-1; + +%% ===== Computation ==================================================== + +% Implementation of Hahn2016, Eq. (32) +Am = zeros(2*Nse+1,size(Anm,2)); +for m=-Nse:Nse + v = sphexp_index(m); % (n,m) = (abs(m),m); + Am(m+Nse+1,:) = Anm(v,:)./(4*pi.*1j.^(m-abs(m))* ... + sphharmonics(abs(m),-m,0,0)); +end diff --git a/SFS_monochromatic/circexp/circexp_mono_ls.m b/SFS_monochromatic/circexp/circexp_mono_ls.m new file mode 100644 index 00000000..cac46a90 --- /dev/null +++ b/SFS_monochromatic/circexp/circexp_mono_ls.m @@ -0,0 +1,116 @@ +function ABm = circexp_mono_ls(xs,mode,Nce,f,xq,conf) +%CIRCEXP_MONO_LS computes the regular/singular circular expansion of line source +% +% Usage: ABm = circexp_mono_ls(xs,mode,Nce,f,xq,conf) +% +% Input parameters: +% xs - position of line source +% mode - 'R' for regular, 'S' for singular +% Nce - maximum order of circular basis functions +% f - frequency / Hz [1 x Nf] or [Nf x 1] +% xq - optional expansion center / m [1 x 3] +% conf - optional configuration struct (see SFS_config) +% +% Output parameters: +% ABm - regular cylindrical expansion coefficients [2*Nce+1 x Nf] +% +% CIRCEXP_MONO_LS(xs,mode,Nce,f,xq,conf) computes the regular circular +% expansion coefficients for a line source. The expansion will be done +% around the expansion coordinate xq: +% +% Regular Expansion: +% \~~ oo +% p (x,f) = > A R (x-x ) +% ls,R /__ n=-oo n n q +% +% with the expansion coefficients: +% m -i +% A = --- S (x -x ) +% n 4 n s q +% +% Singular Expansion: +% \~~ oo m +% p (x,f) = > B S (x-x ) +% ls,R /__ n=-oo n n q +% +% with the expansion coefficients): +% m -i +% B = --- R (x -x ) +% n 4 n s q +% +% see also: circexp_mono_pw circexp_mono_scatter + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 6; +nargmax = 6; +narginchk(nargmin,nargmax); +isargposition(xs); +isargchar(mode); +isargpositivescalar(Nce); +isargvector(f); +isargposition(xq); + +%% ===== Configuration ================================================== +c = conf.c; + +%% ===== Variables ====================================================== +% convert (xs-xq) into cylindrical coordinates +r = sqrt((xs(1)-xq(1)).^2 + (xs(2)-xq(2)).^2); +phi = atan2(xs(2)-xq(2),xs(1)-xq(1)); + +% frequency +k = 2.*pi.*row_vector(f)./c; +kr = k.*r; + +% select suitable basis function +if strcmp('R', mode) + circbasis = @(nu,z) besselh(nu,2,z); +elseif strcmp('S', mode) + circbasis = @besselj; +else + error('unknown mode:'); +end + +%% ===== Computation ==================================================== +L = 2*Nce+1; +Nf = length(kr); +ABm = zeros(L,Nf); +l = 0; +for n=-Nce:Nce + l = l+1; + ABm(l,:) = -1j/4*circbasis(n,kr).*exp(-1j*n*phi); +end + +end diff --git a/SFS_monochromatic/circexp/circexp_mono_pw.m b/SFS_monochromatic/circexp/circexp_mono_pw.m new file mode 100644 index 00000000..b4c31e0a --- /dev/null +++ b/SFS_monochromatic/circexp/circexp_mono_pw.m @@ -0,0 +1,92 @@ +function Am = circexp_mono_pw(npw,Nce,f,xq,conf) +%CIRCEXP_MONO_PW computes regular circular expansion coefficients of plane wave +% +% Usage: Am = circexp_mono_pw(npw,Nce,f,xq,conf) +% +% Input parameters: +% nk - propagation direction of plane wave +% Nce - maximum order of circular basis functions +% f - frequency / Hz [1 x Nf] or [Nf x 1] +% xq - expansion center +% conf - configuration struct (see SFS_config) +% +% Output parameters: +% Am - regular cylindrical expansion coefficients [2*Nce+1 x Nf] +% +% CIRCEXP_MONO_PW(npw,Nce,f,xq,conf) computes the regular circular +% expansion coefficients for a plane wave. The expansion will be done a +% round the expansion coordinate xq: +% +% \~~ oo +% p (x,f) = > A R (x-x ) +% pw /__ n=-oo n n q +% +% with the cylyndrical expansion coefficients: +% +% n +% A = 4pi i exp (-i*n*phi ) +% n pw +% +% see also: circexp_mono_ls circexp_mono_scatter + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 5; +nargmax = 5; +narginchk(nargmin,nargmax); +isargposition(npw); +isargvector(f); +isargposition(xq); + +%% ===== Configuration ================================================== +c = conf.c; + +%% ===== Variables ====================================================== +% convert nk into cylindrical coordinates +phi = atan2(npw(2),npw(1)); +% delay of plane wave to reference point +npw = npw./vector_norm(npw,2); +delay = 2*pi*row_vector(f)/c*(npw*xq.'); + +%% ===== Computation ==================================================== +L = 2*Nce+1; +Nf = length(delay); +Am = zeros(L,Nf); +l = 0; +for n=-Nce:Nce + l = l+1; + Am(l,:) = (1j)^(-n).*exp(-1j*n*phi).*exp(-1j*delay); +end + +end diff --git a/SFS_monochromatic/circexp/circexp_mono_scatter.m b/SFS_monochromatic/circexp/circexp_mono_scatter.m new file mode 100644 index 00000000..0fe1cb10 --- /dev/null +++ b/SFS_monochromatic/circexp/circexp_mono_scatter.m @@ -0,0 +1,115 @@ +function Bm = circexp_mono_scatter(Am,R,sigma,f,conf) +%CIRCEXP_MONO_SCATTER computes the singular circular expansion coefficients of +%cylinder-scattered field +% +% Usage: Bm = circexp_mono_scatter(Al,R,sigma,f,conf) +% +% Input parameters: +% Am - regular circular expansion of incident field [n x Nf] +% R - radius of cylinder / m +% sigma - complex admittance of scatterer +% f - frequency / Hz [1 x Nf] or [Nf x 1] +% conf - configuration struct (see SFS_config) +% +% Output parameters: +% Bm - singular circular expansion coefficients of +% scattered field [n x Nf] +% +% CIRCEXP_MONO_SCATTER(Am,R,sigma,f,conf) computes the singular +% cylindrical expansion coefficients of a field resulting from a scattering +% of an incident field at a cylindrical. Incident field is descriped by +% regular expansion coefficients (expansion center is expected to be at the +% center of the cylinder xq): +% +% \~~ oo +% p (x,f) = > A R (x-x ) +% ind /__ n=0 n n q +% +% The scattered field is descriped by singular expansion coefficients, +% expanded around the center of the cylinder x0. +% +% \~~ oo +% p (x,f) = > B S (x-x ) +% sca /__ n=0 n n q +% +% Due to the boundary conditions on the surface of the cylinder the +% coefficients are related by: +% +% k . J' (kR) + sigma . H (kR) +% n n +% B = - ------------------------------ . A +% n k . H' (kR) + sigma . H (kR) n +% n n +% +% where k = 2*pi*f/c. +% +% see also: circexp_mono_pw circexp_mono_ls + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 5; +nargmax = 5; +narginchk(nargmin,nargmax); +isargmatrix(Am); +isargscalar(sigma); +isargpositivescalar(R); +isargvector(f); + +%% ===== Configuration ================================================== +c = conf.c; + +%% ===== Variables ====================================================== +% frequency +k = 2.*pi.*row_vector(f)./c; +kR = k.*R; +% select suitable transformation function +if isinf(sigma) + T = @(x) -besselj(x,kR)./besselh(x,2,kR); +elseif sigma == 0 + T = @(x) -besselj_derived(x,kR)./besselh_derived(x,2,kR); +else + T = @(x) -(k.*besselj_derived(x,kR)+sigma.*besselj(x,kR)) ... + ./(k.*besselh_derived(x,2,kR)+sigma.*besselh(x,2,kR)); +end + +%% ===== Computation ==================================================== +Nce = (size(Am,1)-1)/2; +Bm = zeros(size(Am)); +l = 0; +for n=-Nce:Nce + l = l+1; + Bm(l,:) = T(n).*Am(l,:); +end + +end diff --git a/SFS_monochromatic/circexp/circexp_mono_timereverse.m b/SFS_monochromatic/circexp/circexp_mono_timereverse.m new file mode 100644 index 00000000..d3cfddd0 --- /dev/null +++ b/SFS_monochromatic/circexp/circexp_mono_timereverse.m @@ -0,0 +1,60 @@ +function Pm = circexp_mono_timereverse(Pm) +%CIRCEXP_MONO_TIMEREVERSE computes coefficients of a time reversed sound field +% +% Usage: Bnm = circexp_mono_timereverse(Pm) +% +% Input parameters: +% Pm - circular expansion coefficients [n x Nf] +% +% Output parameters: +% Pm - time reversed circexp expansion coefficients [n x Nf] +% +% CIRCEXP_MONO_TIMEREVERSE(Pm) + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTPILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 1; +nargmax = 1; +narginchk(nargmin,nargmax); +isargmatrix(Pm); +if mod(size(Pm,1)-1,2) ~= 0 + error('%s: Number of rows of %s has be to odd',upper(mfilename), ... + inputname(Pm)); +end + +%% ===== Computation ==================================================== +Nce = (size(Pm,1)-1) / 2; +Pm = conj(Pm(end:-1:1,:)).*(-1).^(-Nce:Nce).'; + +end diff --git a/SFS_monochromatic/circexp/circexp_mono_translation.m b/SFS_monochromatic/circexp/circexp_mono_translation.m new file mode 100644 index 00000000..0231e578 --- /dev/null +++ b/SFS_monochromatic/circexp/circexp_mono_translation.m @@ -0,0 +1,112 @@ +function [EF,EFm] = circexp_mono_translation(xt,mode,Nce,f,conf) +%CIRCEXP_MONO_TRANSLATION compute circular translation coefficients +%(multipole re-expansion) +% +% Usage: [EF,EFm] = circexp_mono_translation(xt,mode,Nce,f,conf) +% +% Input parameters: +% xt - translatory shift [1x3] / m +% mode - 'RS' for regular-to-singular reexpansion +% 'RR' for regular-to-regular reexpansion +% 'SR' for singular-to-regular reexpansion +% 'SS' for singular-to-singular reexpansion +% f - frequency / Hz +% conf - optional configuration struct (see SFS_config) +% +% Output parameters: +% EF - circular re-expansion coefficients for t +% EFm - circular re-expansion coefficients for -t +% +% CIRCEXP_MONO_TRANSLATION(t,mode,f,conf) computes the circular re-expansion +% coefficients to perform as translatory shift of circular basis function. +% Multipole Re-expansion computes the circular basis function for a shifted +% coordinate system (x+t) based on the original basis functions for (x). +% +% \~~ inf +% E (x + t) = > (E|F) (xt) F (x) +% n /__ l=-inf l,n l +% +% where {E,F} = {R,S}. R denotes the regular circular basis function, while +% S symbolizes the singular circular basis function. Note that (S|S) and +% (S|R) are equivalent to (R|R) and (R|S), respectively. +% +% see also: circexp_mono_ps, circexp_mono_pw + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 5; +nargmax = 5; +narginchk(nargmin,nargmax); +isargposition(xt); +isargchar(mode); +isargpositivescalar(Nce,f); + +%% ===== Configuration ================================================== +c = conf.c; + +%% ===== Variables ====================================================== +% convert t into spherical coordinates +rt = norm(xt(1:2)); +phit = atan2(xt(2),xt(1)); + +% frequency +k = 2*pi*f/c; +kr = k*rt; + +L = 2*Nce+1; +EF = zeros(L,L); +EFm = EF; + +% select suitable basis function +if strcmp('RR',mode) || strcmp('SS',mode) + circbasis = @(nu) besselj(nu,kr); +elseif strcmp('SR',mode) || strcmp('RS',mode) + circbasis = @(nu) besselh(nu,2,kr); +else + error('unknown mode:'); +end + +%% ===== Computation ==================================================== +s = 0; +for n=-Nce:Nce + s = s+1; + l = 0; + for m=-Nce:Nce + l = l+1; + EF(s,l) = circbasis(m-n) .* exp(+1j.*(m-n).*phit); + EFm(s,l) = EF(s,l) .* (-1)^(n-m); + end +end + +end diff --git a/SFS_monochromatic/circexp/circexp_truncation_order.m b/SFS_monochromatic/circexp/circexp_truncation_order.m new file mode 100644 index 00000000..5e29a94f --- /dev/null +++ b/SFS_monochromatic/circexp/circexp_truncation_order.m @@ -0,0 +1,73 @@ +function Nce = circexp_truncation_order(r,f,nmse,conf) +%CIRCEXP_TRUNCATION_ORDER computes truncation order for circular expansion +%coefficients of an arbitrary sound field +% +% Usage: Nce = circexp_truncation_order(r,f,nmse,conf) +% +% Input parameters: +% r - max 2D distance from expansion center / m +% f - frequency / Hz +% nmse - maximum bound for normalized mean squared error +% conf - optional configuration struct (see SFS_config) +% +% Output parameters: +% Nce - Maximum order for circular expansion +% +% CIRCEXP_TRUNCATION_ORDER(r,f,epsilon,conf) yields the order up to which +% a the circular expansion coefficients of an arbitrary sound field have +% be summed up. For a given frequency and maximum radius the normalized +% truncation mean squared error is below the specified error bound (nmse). +% +% References: +% Kennedy et al. (2007) - "Intrinsic Limits of Dimensionality and +% Richness in Random Multipath Fields", +% IEEE Transactions on Signal Processing + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 4; +nargmax = 4; +narginchk(nargmin,nargmax); +isargpositivescalar(f,r,nmse); + +%% ===== Configuration ================================================== +c = conf.c; + +%% ===== Computation ==================================================== +% See Kennedy et al. (eq. 36) +lambda = c/f; % wave length +delta = max(0, ceil(0.5*log(0.0093/nmse))); +Nce = ceil(pi*r*exp(1)/lambda) + delta; + +end diff --git a/SFS_monochromatic/circexp/sound_field_mono_circbasis.m b/SFS_monochromatic/circexp/sound_field_mono_circbasis.m new file mode 100644 index 00000000..904e620b --- /dev/null +++ b/SFS_monochromatic/circexp/sound_field_mono_circbasis.m @@ -0,0 +1,76 @@ +function P = sound_field_mono_circbasis(Pm,JH2m,Ym) +%SOUND_FIELD_MONO_CIRCBASIS simulates a sound field expressed with +%circular basis functions +% +% Usage: P = sound_field_mono_circbasis(Pm,JH2m,Ym) +% +% Input parameters: +% Pm - regular/singular circular expansion coefficients +% JH2m - cell array of cylindrical bessel/hankel(2nd kind) +% functions +% Ym - cell array of cylindrical harmonics (exponential +% functions) +% +% Output parameters: +% P - resulting soundfield +% +% SOUND_FIELD_MONO_CIRCBASIS(Pm,JH2m,Ym) +% +% see also: circbasis_mono_grid sound_field_mono_circexp + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 3; +nargmax = 3; +narginchk(nargmin,nargmax); +isargvector(Pm); +isargequallength(Pm, Ym); +if length(Ym) ~= length(JH2m) + error('%s, length(%s) has to be length(%s)!',upper(mfilename), ... + inputname(2),inputname(3)); +end + +%% ===== Variables ====================================================== +Nce = length(JH2m) - 1; + +%% ===== Computation ==================================================== +P = zeros(size(JH2m{1})); + +l = 0; +for n=0:Nce + l=l+1; + P = P + Pm(l).*(JH2m{l}.*Ym{l}); +end + +end diff --git a/SFS_monochromatic/circexp/sound_field_mono_circexp.m b/SFS_monochromatic/circexp/sound_field_mono_circexp.m new file mode 100644 index 00000000..405abbc8 --- /dev/null +++ b/SFS_monochromatic/circexp/sound_field_mono_circexp.m @@ -0,0 +1,83 @@ +function [P,x,y,z] = sound_field_mono_circexp(X,Y,Z,Pm,mode,f,xq,conf) +%SOUND_FIELD_MONO_CIRCEXP yields a sound field given with regular/singular +%spherical expansion coefficients +% +% Usage: [P,x,y,z] = sound_field_mono_circexp(X,Y,Z,Pm,mode,f,xq,conf) +% +% Input parameters: +% X - x-axis / m; single value or [xmin,xmax] or nD-array +% Y - y-axis / m; single value or [ymin,ymax] or nD-array +% Z - z-axis / m; single value or [zmin,zmax] or nD-array +% Pm - regular/singular spherical expansion coefficients +% mode - 'R' for regular, 'S' for singular +% f - frequency in Hz +% xq - expansion center coordinates +% conf - configuration struct (see SFS_config) +% +% Output parameters: +% P - resulting sound field +% +% SOUND_FIELD_MONO_CIRCEXP(X,Y,Z,Pm,mode,f,xq,conf) +% +% see also: circbasis_mono_grid, sound_field_mono_circbasis + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 8; +nargmax = 8; +narginchk(nargmin,nargmax); +isargvector(Pm); +if mod(size(Pm,1),2) ~= 0 + error('%s: Number of rows of %s has be to odd',upper(mfilename), ... + inputname(4)); +end +isargnumeric(X,Y,Z); +isargpositivescalar(f); +isargposition(xq); + +%% ===== Variables ====================================================== +Nce = ( size(Pm, 1) - 1 ) / 2; + +%% ===== Computation ==================================================== +switch mode + case 'R' + [fn,~,Ym,x,y,z] = circbasis_mono_grid(X,Y,Z,Nce,f,xq,conf); + case 'S' + [~,fn,Ym,x,y,z] = circbasis_mono_grid(X,Y,Z,Nce,f,xq,conf); + otherwise + error('%s: %s is an unknown mode!',upper(mfilename),mode); +end +P = sound_field_mono_circbasis(Pm,fn,Ym); + +end diff --git a/SFS_monochromatic/driving_functions_mono/driving_function_mono_nfchoa_circexp.m b/SFS_monochromatic/driving_functions_mono/driving_function_mono_nfchoa_circexp.m new file mode 100644 index 00000000..43a38840 --- /dev/null +++ b/SFS_monochromatic/driving_functions_mono/driving_function_mono_nfchoa_circexp.m @@ -0,0 +1,107 @@ +function D = driving_function_mono_nfchoa_circexp(x0, Pm,f,conf) +%computes the nfchoa driving functions for a sound field expressed by regular +%spherical expansion coefficients. +% +% Usage: D = driving_function_mono_nfchoa_circexp(x0,Pm,f,conf) +% +% Input parameters: +% x0 - position of the secondary sources / m [N0x3] +% Pm - regular circular expansion coefficients of sound field +% [N x Nf] +% f - frequency / Hz [Nf x 1] or [1 x Nf] +% conf - configuration struct (see SFS_config) +% +% Output parameters: +% D - driving function signal [nx1] +% +% DRIVING_FUNCTION_MONO_NFCHOA_CIRCEXP(x0,Pm,f,conf) returns the NFCHOA +% driving signals for the given secondary sources x0, the virtual sound ex- +% pressed by regular circular expansion coefficients Pm, and the frequency f. +% +% see also: driving_function_mono_wfs_circexp +% driving_function_mono_nfchoa_sphexp + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 4; +nargmax = 4; +narginchk(nargmin,nargmax); +isargmatrix(Pm,x0); +isargvector(f); +isargstruct(conf); + +%% ===== Configuration ================================================== +c = conf.c; +dimension = conf.dimension; +Xc = conf.secondary_sources.center; +R0 = conf.secondary_sources.size/2; % radius of SSD + +%% ===== Computation ==================================================== +Nce = (size(Pm, 1)-1)/2; +% Calculate the driving function in time-frequency domain + +% secondary source positions +x00 = bsxfun(@minus,x0(:,1:3),Xc); +[phi0, r0] = cart2pol(x00(:,1),x00(:,2)); + +% frequency depended stuff +omega = 2*pi*row_vector(f); % [1 x Nf] +k = omega./c; % [1 x Nf] +kr0 = r0 * k; % [N0 x Nf] + +% initialize empty driving signal +D = zeros(size(kr0)); % [N0 x Nf] + +% Calculate the driving function in time-frequency domain +switch dimension + case '2D' + % === 2-Dimensional ================================================== + % Circular Harmonics of Driving Function + Dm = driving_function_mono_nfchoa_cht_circexp(Pm,f,conf); + + l = 0; + for m=-Nce:Nce + l = l+1; + % second line is for compensating difference between r0 and R0 + D = D + bsxfun(@times, Dm(l,:), exp(1i.*m.*phi0)).* ... + bsxfun(@rdivide, besselh(m,2,kr0), besselh(m,2,k*R0)); + end + case {'2.5D', '3D'} + % convert circular expansion to spherical expansion + Pmn = sphexp_convert_circexp(Pm); + % compute driving function for spherical expansion + D = driving_function_mono_nfchoa_sphexp(x0,Pmn,f,conf); + otherwise + error('%s: the dimension %s is unknown.',upper(mfilename),dimension); +end diff --git a/SFS_monochromatic/driving_functions_mono/driving_function_mono_nfchoa_sphexp.m b/SFS_monochromatic/driving_functions_mono/driving_function_mono_nfchoa_sphexp.m new file mode 100644 index 00000000..6186dba1 --- /dev/null +++ b/SFS_monochromatic/driving_functions_mono/driving_function_mono_nfchoa_sphexp.m @@ -0,0 +1,145 @@ +function D = driving_function_mono_nfchoa_sphexp(x0, Pnm,f,conf) +%computes the nfchoa driving functions for a sound field expressed by regular +%spherical expansion coefficients. +% +% Usage: D = driving_function_mono_nfchoa_sphexp(x0,Pnm,f,conf) +% +% Input parameters: +% x0 - position of the secondary sources / m [N0x3] +% Pnm - regular spherical expansion coefficients of sound field +% [N x Nf] +% f - frequency / Hz [Nf x 1] or [1 x Nf] +% conf - configuration struct (see SFS_config) +% +% Output parameters: +% D - driving function signal [nx1] +% +% DRIVING_FUNCTION_MONO_NFCHOA_SPHEXP(x0,Pnm,f,conf) returns NFCHOA +% driving signals for the given secondary sources, the virtual sound expressed +% by regular spherical expansion coefficients and the frequency f. +% +% see also: driving_function_mono_wfs_sphexp + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 4; +nargmax = 4; +narginchk(nargmin,nargmax); +isargmatrix(Pnm,x0); +isargvector(f); +isargstruct(conf); +if mod(sqrt(size(Pnm, 1)),1) ~= 0 + error(['%s: number of columns of Pnm (%s) is not the square of an', ... + 'integer.'], upper(mfilename), sqrt(size(Pnm, 1))); +end + +%% ===== Configuration ================================================== +c = conf.c; +dimension = conf.dimension; +Xc = conf.secondary_sources.center; +R0 = conf.secondary_sources.size/2; + +%% ===== Variables ====================================================== +Nse = sqrt(size(Pnm, 1))-1; + +% secondary source positions +x00 = bsxfun(@minus,x0(:,1:3),Xc); +[phi0, theta0, r0] = cart2sph(x00(:,1),x00(:,2),x00(:,3)); + +% frequency depended stuff +omega = 2*pi*row_vector(f); % [1 x Nf] +k = omega./c; % [1 x Nf] +kr0 = r0 * k; % [N0 x Nf] + +% initialize empty driving signal +D = zeros(size(kr0)); % [N0 x Nf] + +%% ===== Computation ==================================================== +% Calculate the driving function in time-frequency domain + +% Regular spherical expansion of the sound field: +% \~~ N \~~ n m m +% P(x,w) = > > P j (kr) . Y (theta, phi) +% /__ n=0 /__ m=-n n n n +% +% and 3D free field Green's Function: +% \~~ oo \~~ n m m +% G (x0,f) = > > G . j (kr) . Y (theta, phi) +% ps /__ n=0 /__ m=-n n n n +% +% with the regular expansion coefficients of reen's Function +% (see Gumerov2004, eq. 3.2.2): +% m (2) -m +% G = -i . k . h (kr0) Y (pi/2, 0) +% n n n + +if strcmp('2D',dimension) + % === 2-Dimensional ================================================== + % Convert spherical Expansion to circular expansion (approximation) + Pm = circexp_convert_sphexp(Pnm); + % Compute driving function for circular expansion + D = driving_function_mono_nfchoa_circexp(x0,Pm,f,conf); + % +elseif strcmp('2.5D',dimension) + % === 2.5-Dimensional ================================================ + % Spherical Harmonics of Driving Function + Dm = driving_function_mono_nfchoa_cht_sphexp(Pnm,f,conf); + + l = 0; + for m=-Nse:Nse + l=l+1; + % second line is for compensating possible difference between r0 and R0 + D = D + bsxfun(@times, Dm(l,:), exp(1i.*m.*phi0)).* ... + bsxfun(@rdivide, sphbesselh(abs(m),2,kr0), sphbesselh(abs(m),2,k*R0)); + end + % +elseif strcmp('3D',dimension) + % === 3-Dimensional ================================================== + % Spherical Harmonics of Driving Function + Dnm = driving_function_mono_nfchoa_sht_sphexp(Pnm,f,conf); + + % Inverse Spherical Harmonics Transform + for n=0:Nse + Dn = 0; + for m=-n:n + v = sphexp_index(m, n); + Dn = Dn + bsxfun(@times, Dnm(v,:), sphharmonics(n,m,theta0,phi0)); + end + % compensating possible difference between r0 and R0 + D = D + Dn .* bsxfun(@rdivide, sphbesselh(n,2,kr0), sphbesselh(n,2,k*R0)); + end + % +else + error('%s: the dimension %s is unknown.',upper(mfilename),dimension); +end diff --git a/SFS_monochromatic/driving_functions_mono/driving_function_mono_wfs_circexp.m b/SFS_monochromatic/driving_functions_mono/driving_function_mono_wfs_circexp.m new file mode 100644 index 00000000..3c94697c --- /dev/null +++ b/SFS_monochromatic/driving_functions_mono/driving_function_mono_wfs_circexp.m @@ -0,0 +1,170 @@ +function D = driving_function_mono_wfs_circexp(x0,n0,Pm,mode,f,xq,conf) +%computes the wfs driving functions for a sound field expressed by cylindrical +%expansion coefficients. +% +% Usage: D = driving_function_mono_wfs_circexp(x0,n0,Pm,mode,f,xq,conf) +% +% Input parameters: +% x0 - position of the secondary sources / m [nx3] +% n0 - directions of the secondary sources / m [nx3] +% Pm - circular expansion coefficients of sound field +% mode - 'R' for regular expansion, 'S' for singular expansion +% f - frequency / Hz +% xq - expansion center coordinates / m [1x3] +% conf - configuration struct (see SFS_config) +% +% Output parameters: +% D - driving function signal [nx1] +% +% DRIVING_FUNCTION_MONO_WFS_CIRCEXP(x0,n0,Pm,mode,f,xq,conf) +% +% see also: driving_function_mono_wfs_sphexpS +% +% References: +% Spors2011 - "Local Sound Field Synthesis by Virtual Acoustic Scattering +% and Time-Reversal" (AES131) +% Gumerov,Duraiswami (2004) - "Fast Multipole Methods for the Helmholtz +% Equation in three Dimensions", ELSEVIER + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 7; +nargmax = 7; +narginchk(nargmin,nargmax); +isargmatrix(x0,n0); +isargpositivescalar(f); +isargstruct(conf); +isargposition(xq); + +%% ===== Configuration ================================================== +xref = conf.xref; +c = conf.c; +dimension = conf.dimension; +driving_functions = conf.driving_functions; + +%% ===== Computation ==================================================== +Nce = (size(Pm, 1)-1)/2; + +% apply shift to the center of expansion xq +x = x0(:,1)-xq(1); +y = x0(:,2)-xq(2); + +% conversion to cylindrical coordinates +r0 = sqrt(x.^2 + y.^2); +phi0 = atan2(y,x); + +% frequency depended stuff +omega = 2*pi*f; +k = omega/c; +kr0 = k.*r0; + +% gradient in spherical coordinates +Gradr = zeros(size(x0,1),1); +Gradphi = zeros(size(x0,1),1); +Gradz = zeros(size(x0,1),1); + +% directional weights for conversion spherical gradient into carthesian +% coordinates + point product with normal vector n0 (directional derivative +% in cartesian coordinates) +Sn0r = cos(phi0).*n0(:,1)... + + sin(phi0).*n0(:,2); +Sn0phi = -sin(phi0).*n0(:,1)... + + cos(phi0).*n0(:,2); + +% select suitable basis function +if strcmp('R', mode) + circbasis = @besselj; + circbasis_derived = @besselj_derived; +elseif strcmp('S', mode) + circbasis = @(nu,z) besselh(nu,2,z); + circbasis_derived = @(nu,z) besselh_derived(nu,2,z); +else + error('unknown mode:'); +end + +% Calculate the driving function in time-frequency domain + +% indexing the expansion coefficients +l = 0; + +switch dimension + case {'2D', '2.5D', '3D'} + + % === 2-,2.5-,3-Dimensional ============================================ + switch driving_functions + case 'default' + % --- SFS Toolbox ------------------------------------------------ + % d + % D(x0, w) = -2 --- P(x0,w) + % d n + % with cylindrical expansion of the sound field: + % \~~ oo + % P(x,w) = > B F (x-xq) + % /__ n=-oo n n + % + % where F = {R,S}. + % + % regular cylindrical basis functions: + % + % R (x) = J (kr) . exp(j n phi)) + % n n + % singular cylindrical basis functions + % (2) + % S (x) = H (kr) . exp(j n phi) + % n n + + for n=-Nce:Nce + l = l + 1; + cn_prime = k.*circbasis_derived(n,kr0); + cn = circbasis(n,kr0); + Yn = exp(1j.*n.*phi0); + Gradr = Gradr + ( Pm(l).*cn_prime.*Yn ); + Gradphi = Gradphi + 1./r0.*( Pm(l).*cn.*1j.*n.*Yn ); + end + % directional gradient + D = -2*( Sn0r.*Gradr + Sn0phi.*Gradphi); + + % 2.5D correction factor + if strcmp('2.5D', dimension) + xref = repmat(xref,[size(x0,1) 1]); + g0 = sqrt( 2*pi*vector_norm(xref-x0, 2)./(1j.*k)); + D = D.*g0; + end + otherwise + error(['%s: %s, this type of driving function is not implemented ', ... + 'for 2D/2.5D/3D.'],upper(mfilename),driving_functions); + end + otherwise + error('%s: the dimension %s is unknown.',upper(mfilename),dimension); +end diff --git a/SFS_monochromatic/driving_functions_mono/driving_function_mono_wfs_sphexp.m b/SFS_monochromatic/driving_functions_mono/driving_function_mono_wfs_sphexp.m new file mode 100644 index 00000000..f1692f66 --- /dev/null +++ b/SFS_monochromatic/driving_functions_mono/driving_function_mono_wfs_sphexp.m @@ -0,0 +1,187 @@ +function D = driving_function_mono_wfs_sphexp(x0,n0,Pnm,mode,f,xq,conf) +%computes the wfs driving functions for a sound field expressed by spherical +%expansion coefficients. +% +% Usage: D = driving_function_mono_wfs_sphexp(x0,n0,Pnm,mode,f,xq,conf) +% +% Input parameters: +% x0 - position of the secondary sources / m [nx3] +% n0 - directions of the secondary sources / m [nx3] +% Pnm - singular spherical expansion coefficients of sound field +% mode - 'R' for regular expansion, 'S' for singular expansion +% f - frequency in Hz +% xq - expansion center coordinates / m [1x3] +% conf - configuration struct (see SFS_config) +% +% Output parameters: +% D - driving function signal [nx1] +% +% DRIVING_FUNCTION_MONO_WFS_SPHEXP(x0,n0,Pnm,mode,f,xq,conf) +% +% see also: driving_function_mono_wfs_cylexp +% +% References: +% Spors2011 - "Local Sound Field Synthesis by Virtual Acoustic Scattering +% and Time-Reversal" (AES131) +% Gumerov,Duraiswami (2004) - "Fast Multipole Methods for the Helmholtz +% Equation in three Dimensions", ELSEVIER + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 7; +nargmax = 7; +narginchk(nargmin,nargmax); +isargmatrix(x0,n0); +isargvector(Pnm); +isargpositivescalar(f); +isargchar(mode); +isargposition(xq); +isargstruct(conf); + +%% ===== Configuration ================================================== +c = conf.c; +dimension = conf.dimension; +driving_functions = conf.driving_functions; + +%% ===== Variables ====================================================== +Nse = sqrt(size(Pnm, 1))-1; + +% apply shift to the center of expansion xq +x = x0(:,1)-xq(1); +y = x0(:,2)-xq(2); +z = x0(:,3)-xq(3); + +% conversion to spherical coordinates +r0 = sqrt(x.^2 + y.^2 + z.^2); +phi0 = atan2(y,x); +theta0 = asin(z./r0); + +% frequency depended stuff +omega = 2*pi*f; +k = omega/c; +kr = k.*r0; + +% gradient in spherical coordinates +Gradr = zeros(size(x0,1),1); +Gradphi = zeros(size(x0,1),1); +Gradtheta = zeros(size(x0,1),1); + +% directional weights for conversion spherical gradient into carthesian +% coordinates + point product with normal vector n0 (directional derivative +% in cartesian coordinates) +Sn0r = cos(theta0).*cos(phi0).*n0(:,1)... + + cos(theta0).*sin(phi0).*n0(:,2)... + + sin(theta0) .*n0(:,3); +Sn0phi = -sin(phi0) .*n0(:,1)... + + cos(phi0) .*n0(:,2); +Sn0theta = sin(theta0).*cos(phi0).*n0(:,1)... + + sin(theta0).*sin(phi0).*n0(:,2)... + - cos(theta0) .*n0(:,3); + +% select suitable basis function +if strcmp('R', mode) + sphbasis = @sphbesselj; + sphbasis_derived = @sphbesselj_derived; +elseif strcmp('S', mode) + sphbasis = @(nu,z) sphbesselh(nu,2,z); + sphbasis_derived = @(nu,z) sphbesselh_derived(nu,2,z); +else + error('unknown mode:'); +end + +%% ===== Computation ==================================================== +% Calculate the driving function in time-frequency domain + +% indexing the expansion coefficients +l = 0; + +switch dimension + case {'3D'} + + if (strcmp('default',driving_functions)) + % --- SFS Toolbox ------------------------------------------------ + % + % d + % D(x0, w) = -2 ------ P(x0,w) + % d n0 + % with regular/singular spherical expansion of the sound field: + % \~~ N \~~ n m m + % P(x,w) = > > B F (x-xq) + % /__ n=0 /__ m=-n n n + % + % where F = {R,S}. + % + % regular spherical basis functions: + % m m + % R (x) = j (kr) . Y (theta, phi) + % n n n + % singular spherical basis functions: + % m (2) m + % S (x) = h (kr) . Y (theta, phi) + % n n n + + for n=0:Nse + cn_prime = k.*sphbasis_derived(n,kr); + cn = sphbasis(n, kr); + for m=-n:n + l = l + 1; + Ynm = sphharmonics(n,m, theta0, phi0); + Gradr = Gradr + Pnm(l).*cn_prime.*Ynm; + Gradphi = Gradphi + 1./r0.*Pnm(l).*cn.*1j.*m.*Ynm; + Gradtheta = Gradtheta + 1./(r0.*cos(theta0).^2).*cn.*Pnm(l).*... + ( (n+1).*sin(theta0).*Ynm - ... + sqrt((2*n+1)./(2*n+3).*((n+1).^2-m^2)).* ... + sphharmonics(n+1,m, theta0, phi0) ); + end + end + % directional gradient + D = -2*( Sn0r.*Gradr + Sn0phi.*Gradphi + Sn0theta.*Gradtheta ); + else + error(['%s: %s, this type of driving function is not implemented ', ... + 'for 3D'],upper(mfilename),driving_functions); + end + + case {'2D', '2.5D'} + if (strcmp('default',driving_functions)) + % approximate spherical expansion with circular expansion + Pm = circexp_convert_sphexp(Pnm); + % compute driving function for circular expansion + D = driving_function_mono_wfs_circexp(x0,n0,Pm,mode,f,xq,conf); + else + error(['%s: %s, this type of driving function is not implemented ', ... + 'for 2D/2.5D'],upper(mfilename),driving_functions); + end + otherwise + error('%s: the dimension %s is unknown.',upper(mfilename),dimension); +end diff --git a/SFS_monochromatic/driving_function_mono_sdm_kx.m b/SFS_monochromatic/kx/driving_function_mono_sdm_kx.m similarity index 100% rename from SFS_monochromatic/driving_function_mono_sdm_kx.m rename to SFS_monochromatic/kx/driving_function_mono_sdm_kx.m diff --git a/SFS_monochromatic/driving_functions_mono/driving_function_mono_sdm_kx_fs.m b/SFS_monochromatic/kx/driving_function_mono_sdm_kx_fs.m similarity index 100% rename from SFS_monochromatic/driving_functions_mono/driving_function_mono_sdm_kx_fs.m rename to SFS_monochromatic/kx/driving_function_mono_sdm_kx_fs.m diff --git a/SFS_monochromatic/driving_functions_mono/driving_function_mono_sdm_kx_ls.m b/SFS_monochromatic/kx/driving_function_mono_sdm_kx_ls.m similarity index 100% rename from SFS_monochromatic/driving_functions_mono/driving_function_mono_sdm_kx_ls.m rename to SFS_monochromatic/kx/driving_function_mono_sdm_kx_ls.m diff --git a/SFS_monochromatic/driving_functions_mono/driving_function_mono_sdm_kx_ps.m b/SFS_monochromatic/kx/driving_function_mono_sdm_kx_ps.m similarity index 100% rename from SFS_monochromatic/driving_functions_mono/driving_function_mono_sdm_kx_ps.m rename to SFS_monochromatic/kx/driving_function_mono_sdm_kx_ps.m diff --git a/SFS_monochromatic/driving_functions_mono/driving_function_mono_sdm_kx_pw.m b/SFS_monochromatic/kx/driving_function_mono_sdm_kx_pw.m similarity index 100% rename from SFS_monochromatic/driving_functions_mono/driving_function_mono_sdm_kx_pw.m rename to SFS_monochromatic/kx/driving_function_mono_sdm_kx_pw.m diff --git a/SFS_monochromatic/sound_field_mono_sdm_kx.m b/SFS_monochromatic/kx/sound_field_mono_sdm_kx.m similarity index 100% rename from SFS_monochromatic/sound_field_mono_sdm_kx.m rename to SFS_monochromatic/kx/sound_field_mono_sdm_kx.m diff --git a/SFS_monochromatic/sht/driving_function_mono_nfchoa_sht.m b/SFS_monochromatic/sht/driving_function_mono_nfchoa_sht.m new file mode 100644 index 00000000..af8529e5 --- /dev/null +++ b/SFS_monochromatic/sht/driving_function_mono_nfchoa_sht.m @@ -0,0 +1,87 @@ +function Dnm = driving_function_mono_nfchoa_sht(xs, src, Nse, f, conf) +%DRIVING_FUNCTION_MONO_NFCHOA_SHT computes the spherical harmonics transform of +%nfchoa driving functions for a specified source model +% +% Usage: D = driving_function_mono_nfchoa_sht(xs, src, Nse, f, conf) +% +% Input parameters: +% xs - position of virtual source or direction of plane +% wave [1 x 3] / m +% Nse - maximum order of spherical basis functions +% f - frequency [m x 1] or [1 x m] / Hz +% conf - configuration struct (see SFS_config) +% +% Output parameters: +% Dnm - regular spherical harmonics transform of driving +% function signal [n x m] +% +% DRIVING_FUNCTION_MONO_NFCHOA_SHT(xs, src, Nse, f, conf) returns spherical +% harmonics transform of the NFCHOA driving function with maximum order Nse +% for a specified source model. +% +% see also: driving_function_mono_nfchoa_sht_sphexp + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 5; +nargmax = 5; +narginchk(nargmin,nargmax); +isargposition(xs); +isargchar(src); +isargpositivescalar(Nse); +isargvector(f); +isargstruct(conf); + +%% ===== Computation ==================================================== +% Get SHT of driving signals +if strcmp('pw',src) + % === Plane wave ===================================================== + % Direction of plane wave + nk = xs / norm(xs); + % SHT of Driving signal + Dnm = driving_function_mono_nfchoa_sht_pw(nk, Nse, f,conf); + +elseif strcmp('ps',src) + % === Point source =================================================== + % SHT of Driving signal + Dnm = driving_function_mono_nfchoa_sht_ps(xs, Nse, f,conf); + +elseif strcmp('ls',src) + % === Line source ==================================================== + % SHT of Driving signal + Dnm = driving_function_mono_nfchoa_sht_ls(xs, Nse, f,conf); + +else + error('%s: %s is not a known source type.',upper(mfilename),src); +end diff --git a/SFS_monochromatic/sht/driving_function_mono_nfchoa_sht_ls.m b/SFS_monochromatic/sht/driving_function_mono_nfchoa_sht_ls.m new file mode 100644 index 00000000..74cb232b --- /dev/null +++ b/SFS_monochromatic/sht/driving_function_mono_nfchoa_sht_ls.m @@ -0,0 +1,71 @@ +function Dnm = driving_function_mono_nfchoa_sht_ls(xs, Nse, f, conf) +%DRIVING_FUNCTION_MONO_NFCHOA_SHT_LS computes the spherical harmonics transform +%of nfchoa driving functions for a line source. +% +% Usage: D = driving_function_mono_nfchoa_sht_ls(xs, Nse, f, conf) +% +% Input parameters: +% xs - position of line source [1 x 3] / m +% Nse - maximum order of spherical basis functions +% f - frequency [m x 1] or [1 x m] / Hz +% conf - configuration struct (see SFS_config) +% +% Output parameters: +% Dnm - regular spherical harmonics transform of driving +% function signal [n x m] +% +% DRIVING_FUNCTION_MONO_NFCHOA_SHT_LS(xs, Nse, f, conf) returns spherical +% harmonics transform of the NFCHOA driving function with maximum order Nse +% for a virtual line source at xs. +% +% see also: sphexp_mono_ls driving_function_mono_nfchoa_sht_sphexp + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 4; +nargmax = 4; +narginchk(nargmin,nargmax); +isargposition(xs); +isargpositivescalar(Nse); +isargvector(f); +isargstruct(conf); + +%% ===== Configuration ================================================== +X0 = conf.secondary_sources.center; + +%% ===== Computation ==================================================== +% Calculate spherical expansion coefficients of point source +Pnm = sphexp_mono_ls(xs, 'R', Nse, f, X0, conf); +% Calculate spherical harmonics driving function +Dnm = driving_function_mono_nfchoa_sht_sphexp(Pnm, f, conf); diff --git a/SFS_monochromatic/sht/driving_function_mono_nfchoa_sht_ps.m b/SFS_monochromatic/sht/driving_function_mono_nfchoa_sht_ps.m new file mode 100644 index 00000000..e2b3af8e --- /dev/null +++ b/SFS_monochromatic/sht/driving_function_mono_nfchoa_sht_ps.m @@ -0,0 +1,71 @@ +function Dnm = driving_function_mono_nfchoa_sht_ps(xs, Nse, f, conf) +%DRIVING_FUNCTION_MONO_NFCHOA_SHT_PS computes the spherical harmonics +%transform of nfchoa driving functions for a point source. +% +% Usage: D = driving_function_mono_nfchoa_sht_ps(xs, Nse, f, conf) +% +% Input parameters: +% xs - position of point source [1 x 3] / m +% Nse - maximum order of spherical basis functions +% f - frequency [m x 1] or [1 x m] / Hz +% conf - configuration struct (see SFS_config) +% +% Output parameters: +% Dnm - regular spherical harmonics transform of driving +% function signal [n x m] +% +% DRIVING_FUNCTION_MONO_NFCHOA_SHT_PS(xs, Nse, f, conf) returns spherical +% harmonics transform of the NFCHOA driving function with maximum order Nse +% for a virtual point source at xs. +% +% see also: sphexp_mono_ps driving_function_mono_nfchoa_sht_sphexp + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 4; +nargmax = 4; +narginchk(nargmin,nargmax); +isargposition(xs); +isargpositivescalar(Nse); +isargvector(f); +isargstruct(conf); + +%% ===== Configuration ================================================== +X0 = conf.secondary_sources.center; + +%% ===== Computation ==================================================== +% Calculate spherical expansion coefficients of point source +Pnm = sphexp_mono_ps(xs, 'R', Nse, f, X0, conf); +% Calculate spherical harmonics driving function +Dnm = driving_function_mono_nfchoa_sht_sphexp(Pnm, f, conf); diff --git a/SFS_monochromatic/sht/driving_function_mono_nfchoa_sht_pw.m b/SFS_monochromatic/sht/driving_function_mono_nfchoa_sht_pw.m new file mode 100644 index 00000000..b79139bf --- /dev/null +++ b/SFS_monochromatic/sht/driving_function_mono_nfchoa_sht_pw.m @@ -0,0 +1,71 @@ +function Dnm = driving_function_mono_nfchoa_sht_pw(npw, Nse, f, conf) +%DRIVING_FUNCTION_MONO_NFCHOA_SHT_PW computes the spherical harmonics +%transform of nfchoa driving functions for a point source. +% +% Usage: D = driving_function_mono_nfchoa_sht_pw(npw, Nse, f, conf) +% +% Input parameters: +% npw - unit vector propagation direction of plane wave +% Nse - maximum order of spherical basis functions +% f - frequency [m x 1] or [1 x m] / Hz +% conf - configuration struct (see SFS_config) +% +% Output parameters: +% Dnm - regular spherical harmonics transform of driving +% function signal [n x m] +% +% DRIVING_FUNCTION_MONO_NFCHOA_SHT_PW(npw, Nse, f, conf) returns spherical +% harmonics transform of the NFCHOA driving function with maximum order Nse +% for a virtual plane wave with propagation direction npw. +% +% see also: sphexp_mono_pw driving_function_mono_nfchoa_sht_sphexp + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 4; +nargmax = 4; +narginchk(nargmin,nargmax); +isargposition(npw); +isargpositivescalar(Nse); +isargvector(f); +isargstruct(conf); + +%% ===== Configuration ================================================== +X0 = conf.secondary_sources.center; + +%% ===== Computation ==================================================== +% Calculate spherical expansion coefficients of point source +Pnm = sphexp_mono_pw(npw, Nse, f, X0, conf); +% Calculate spherical harmonics driving function +Dnm = driving_function_mono_nfchoa_sht_sphexp(Pnm, f, conf); diff --git a/SFS_monochromatic/sht/driving_function_mono_nfchoa_sht_sphexp.m b/SFS_monochromatic/sht/driving_function_mono_nfchoa_sht_sphexp.m new file mode 100644 index 00000000..e67b109a --- /dev/null +++ b/SFS_monochromatic/sht/driving_function_mono_nfchoa_sht_sphexp.m @@ -0,0 +1,112 @@ +function Dnm = driving_function_mono_nfchoa_sht_sphexp(Pnm, f, conf) +%DRIVING_FUNCTION_MONO_NFCHOA_SHT_SPHEXP computes the spherical harmonics +%transform of nfchoa driving functions for a sound field expressed by regular +%spherical expansion coefficients. +% +% Usage: D = driving_function_mono_nfchoa_sht_sphexp(Pnm, f, conf) +% +% Input parameters: +% Pnm - regular spherical expansion coefficients of virtual +% sound field [nxm] +% f - frequency in Hz [m x 1] or [1 x m] +% conf - configuration struct (see SFS_config) +% +% Output parameters: +% Dnm - regular spherical harmonics transform of driving +% function signal [n x m] +% +% DRIVING_FUNCTION_MONO_NFCHOA_SHT_SPHEXP(Pnm, f, conf) returns spherical +% harmonics transform of the NFCHOA driving function for a virtual sound +% expressed by regular spherical expansion coefficients and the frequency f. + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 3; +nargmax = 3; +narginchk(nargmin,nargmax); +isargmatrix(Pnm); +isargsquaredinteger(size(Pnm,1)); +isargvector(f); +isargstruct(conf); +if size(Pnm,2) ~= length(f) + error( '%s:number of rows in %s have to match length of %s', ... + upper(mfilename), inputname(1), inputname(2) ); +end + +%% ===== Configuration ================================================== +c = conf.c; +R0 = conf.secondary_sources.size / 2; + +%% ===== Variables ====================================================== +Nse = sqrt(size(Pnm,1))-1; + +% frequency depended stuff +omega = 2*pi*row_vector(f); % [1 x Nf] +k = omega./c; % [1 x Nf] +kR0 = k.*R0; % [1 x Nf] + +%% ===== Computation ==================================================== +% Calculate the spherical harmonics of driving function +% +% Regular spherical expansion of the sound field: +% \~~ N \~~ n m m +% P(x,w) = > > P j (kr) . Y (theta, phi) +% /__ n=0 /__ m=-n n n n +% +% and 3D free field Green's Function: +% \~~ oo \~~ n m m +% G (x0,f) = > > G . j (kr) . Y (theta, phi) +% ps /__ n=0 /__ m=-n n n n +% +% with the regular expansion coefficients of Green's Function +% (see Gumerov2004, eq. 3.2.2): +% m (2) -m +% G = -i . k . h (kr0) Y (pi/2, 0) +% n n n + +% m +% P +% m n +% D = --------------------- +% n (2) +% -j k r0^2 h (k r0) +% n + +Dnm = zeros(size(Pnm)); +for n=0:Nse + v = sphexp_index(-n:n, n); % m=-n:n + Dnm(v,:) = bsxfun(@rdivide, Pnm(v,:), sphbesselh(n,2,kR0)); +end +% order independent factor +Dnm = bsxfun(@rdivide, Dnm, -1j.*(R0.^2)*k); diff --git a/SFS_monochromatic/sht/driving_function_mono_wfs_sht_sphexp.m b/SFS_monochromatic/sht/driving_function_mono_wfs_sht_sphexp.m new file mode 100644 index 00000000..c9820d04 --- /dev/null +++ b/SFS_monochromatic/sht/driving_function_mono_wfs_sht_sphexp.m @@ -0,0 +1,119 @@ +function Dnm = driving_function_mono_wfs_sht_sphexp(Pnm, mode, f, conf) +%DRIVING_FUNCTION_MONO_WFS_SHT_SPHEXP computes spherical harmonics transform +%the wfs driving functions for a sound field expressed by spherical expansion +%coefficients. +% +% Usage: Dnm = driving_function_mono_wfs_sht_sphexp(Pnm, mode, f, conf) +% +% Input parameters: +% Pnm - regular/singular spherical expansion coefficients of +% sound field +% mode - 'R' for regular expansion, 'S' for singular expansion +% f - frequency / Hz +% conf - configuration struct (see SFS_config) +% +% Output parameters: +% Dnm - regular spherical harmonics transform of driving +% function signal [n x m] +% +% DRIVING_FUNCTION_MONO_WFS_SHT_SPHEXP(Pnm, mode, f, conf) + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 4; +nargmax = 4; +narginchk(nargmin,nargmax); +isargmatrix(Pnm); +isargvector(f); +isargchar(mode); +isargstruct(conf); + +%% ===== Configuration ================================================== +c = conf.c; +driving_functions = conf.driving_functions; +r0 = conf.secondary_sources.size / 2; + +%% ===== Variables ====================================================== +Nse = sqrt(size(Pnm, 1))-1; + +% frequency depended stuff +k = 2.*pi.*f(:)./c; +kr0 = k.*r0; + +% select suitable basis function +if strcmp('R', mode) + sphbasis_derived = @(nu) k.*sphbesselj_derived(nu, kr0); +elseif strcmp('S', mode) + sphbasis_derived = @(nu) k.*sphbesselh_derived(nu,2,kr0); +else + error('%s: unknown mode (%s)!', upper(mfilename), mode); +end + +%% ===== Computation ==================================================== +% Calculate the driving function in time-frequency domain + +if (strcmp('default',driving_functions)) + % --- SFS Toolbox ------------------------------------------------ + % + % d + % D(x0, w) = ---- Ps(x0,w) + % d r0 + % with regular/singular spherical expansion of the sound field: + % \~~ N \~~ n m m + % P(x,w) = > > B F (x) + % /__ n=0 /__ m=-n n n + % + % where F = {R,S}. + % + % regular spherical basis functions: + % m m + % R (x) = j (kr) . Y (theta, phi) + % n n n + % singular spherical basis functions: + % m (2) m + % S (x) = h (kr) . Y (theta, phi) + % n n n + + Dnm = zeros( size(Pnm) ); + for n=0:Nse + v = sphexp_index(-n:n, n); % m=-n:n + + Dnm(v,:) = Pnm(v,:) .* repmat( sphbasis_derived(n), 2*n+1, 1); + end + Dnm = -2.*Dnm; + +else + error(['%s: %s, this type of driving function is not implemented ', ... + 'for a 2.5D point source.'], upper(mfilename), driving_functions); +end diff --git a/SFS_monochromatic/sht/sound_field_mono_sht.m b/SFS_monochromatic/sht/sound_field_mono_sht.m new file mode 100644 index 00000000..1aa4e6b8 --- /dev/null +++ b/SFS_monochromatic/sht/sound_field_mono_sht.m @@ -0,0 +1,112 @@ +function [P, x, y, z] = sound_field_mono_sht(X,Y,Z,Dnm,f,conf) +%SOUND_FIELD_MONO_SHT simulates a sound field given with the spherical +%harmonics transform of the driving function +% +% Usage: [P, x, y, z] = sound_field_mono_sht(X,Y,Z,Dnm,f,conf) +% +% Input parameters: +% X - x-axis / m; single value or [xmin,xmax] or nD-array +% Y - y-axis / m; single value or [ymin,ymax] or nD-array +% Z - z-axis / m; single value or [zmin,zmax] or nD-array +% Dnm - spherical harmonics transform of nfchoa driving function +% f - frequency in Hz +% conf - configuration struct (see SFS_config) +% +% Output parameters: +% P - resulting soundfield +% +% SOUND_FIELD_MONO_SHT(X,Y,Z,Dnm,f,conf) +% +% see also: sphbasis_mono_grid, sound_field_mono_sphexp + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 6; +nargmax = 6; +narginchk(nargmin,nargmax); +isargvector(Dnm); +isargsquaredinteger(length(Dnm)); +isargnumeric(X,Y,Z); +isargpositivescalar(f); +isargstruct(conf); + +%% ===== Configuration ================================================== +Xc = conf.secondary_sources.center; +r0 = conf.secondary_sources.size / 2; +useplot = conf.plot.useplot; + +%% ===== Computation ==================================================== +[x,y,z] = xyz_grid(X,Y,Z,conf); +% find coordinates, which are inside and outside the loudspeaker array +select = sqrt((x-Xc(1)).^2 + (y-Xc(2)).^2 + (z-Xc(3)).^2) <= r0; + +if (numel(x) == 1) x = repmat(x, size(select)); end +if (numel(y) == 1) y = repmat(y, size(select)); end +if (numel(z) == 1) z = repmat(z, size(select)); end + +P = zeros(size(x)); + +% regular (interior) domain inside the loudspeaker array +if any(select(:)) + Pnm = sphexp_mono_sht(Dnm,'R',f,conf); + P(select) = sound_field_mono_sphexp(x(select), y(select), z(select), ... + Pnm, 'R', f, Xc,conf); +end + +% singular (exterior) domain outside the loudspeaker array +if any(~select(:)) + Pnm = sphexp_mono_sht(Dnm,'S',f,conf); + if sum( ~select(:) ) == 2 + % this handle cases, where x(~select) would only contain 2 entries, which + % would be interpreted as a range by the sound_field_... function + xtmp = x(~select); + ytmp = y(~select); + ztmp = z(~select); + Ptmp(1) = sound_field_mono_sphexp(xtmp(1), ytmp(1), ztmp(1), ... + Pnm, 'S', f, Xc,conf); + Ptmp(2) = sound_field_mono_sphexp(xtmp(2), ytmp(2), ztmp(2), ... + Pnm, 'S', f, Xc,conf); + P(~select) = [Ptmp(1); Ptmp(2)]; + else + P(~select) = sound_field_mono_sphexp(x(~select), y(~select), z(~select), ... + Pnm, 'S', f, Xc,conf); + end +end + +%% ===== Plotting ======================================================= +if (nargout==0 || useplot) + plot_sound_field(P,X,Y,Z,[],conf); +end + +end \ No newline at end of file diff --git a/SFS_monochromatic/sphexp/sound_field_mono_sphbasis.m b/SFS_monochromatic/sphexp/sound_field_mono_sphbasis.m new file mode 100644 index 00000000..1ae8bc77 --- /dev/null +++ b/SFS_monochromatic/sphexp/sound_field_mono_sphbasis.m @@ -0,0 +1,78 @@ +function P = sound_field_mono_sphbasis(ABnm,jh2n,Ynm) +%SOUND_FIELD_MONO_SPHBASIS simulates a sound field expressed with spherical +%basis functions +% +% Usage: P = sound_field_mono_sphbasis(AB,jh2n,Ynm) +% +% Input parameters: +% ABnm - regular/singular spherical expansion coefficients +% jh2n - cell array of spherical bessel/hankel(2nd kind) functions +% Ynm - cell array of spherical harmonics +% +% Output parameters: +% P - resulting soundfield +% +% SOUND_FIELD_MONO_SPHBASIS(ABnm,jh2n,Ynm) +% +% see also: sphbasis_mono_grid sound_field_mono_sphexp + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 3; +nargmax = 3; +narginchk(nargmin,nargmax); +isargvector(ABnm); +L = size(ABnm,1); +isargsquaredinteger(L); +if length(Ynm) ~= length(jh2n)^2 + error('%s: length(Y) has to be equal length(jh2n)^2!',upper(mfilename)); +end +if length(Ynm) < length(ABnm) + error('%s: length(Y) has to larger equal length(ABnm)!',upper(mfilename)); +end + +%% ===== Variables ====================================================== +Nse = sqrt(L) - 1; + +%% ===== Computation ==================================================== +P = zeros(size(jh2n{1})); + +% spherical basic functions +l = 0; +for n=0:Nse + for m=-n:n + l=l+1; + P = P + ABnm(l)*(jh2n{n+1}.*Ynm{l}); + end +end diff --git a/SFS_monochromatic/sphexp/sound_field_mono_sphexp.m b/SFS_monochromatic/sphexp/sound_field_mono_sphexp.m new file mode 100644 index 00000000..76302c48 --- /dev/null +++ b/SFS_monochromatic/sphexp/sound_field_mono_sphexp.m @@ -0,0 +1,79 @@ +function [P,x,y,z] = sound_field_mono_sphexp(X,Y,Z,ABnm,mode,f,xq,conf) +%SOUND_FIELD_MONO_SPHEXPR simulates a sound field given with regular/singular +%spherical expansion coefficients +% +% Usage: [P,x,y,z] = sound_field_mono_sphexp(X,Y,Z,ABnm,mode,f,xq,conf) +% +% Input parameters: +% X - x-axis / m; single value or [xmin,xmax] or nD-array +% Y - y-axis / m; single value or [ymin,ymax] or nD-array +% Z - z-axis / m; single value or [zmin,zmax] or nD-array +% ABnm - regular/singular spherical expansion coefficients +% mode - 'R' for regular,'S' for singular +% f - frequency in Hz +% xq - expansion center coordinates,default: [0,0,0] +% conf - configuration struct (see SFS_config) +% +% Output parameters: +% P - resulting soundfield +% +% SOUND_FIELD_MONO_SPHEXPR(X,Y,Z,ABnm,mode,f,xq,conf) +% +% see also: sphbasis_mono_grid,sound_field_mono_sphbasis + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 8; +nargmax = 8; +narginchk(nargmin,nargmax); +isargvector(ABnm); +isargsquaredinteger(length(ABnm)); +isargnumeric(X,Y,Z); +isargpositivescalar(f); +isargposition(xq); +isargstruct(conf); + +%% ===== Variables ====================================================== +Nse = sqrt(length(ABnm)) - 1; + +%% ===== Computation ==================================================== +if strcmp('R',mode) + [fn,~,Ynm,x,y,z] = sphbasis_mono_grid(X,Y,Z,Nse,f,xq,conf); +elseif strcmp('S',mode) + [~,fn,Ynm,x,y,z] = sphbasis_mono_grid(X,Y,Z,Nse,f,xq,conf); +else + error('%s: %s is an unknown mode!',upper(mfilename),mode); +end + +P = sound_field_mono_sphbasis(ABnm,fn,Ynm); diff --git a/SFS_monochromatic/sphexp/sphbasis_mono.m b/SFS_monochromatic/sphexp/sphbasis_mono.m new file mode 100644 index 00000000..40c42c68 --- /dev/null +++ b/SFS_monochromatic/sphexp/sphbasis_mono.m @@ -0,0 +1,101 @@ +function [jn,h2n,Ynm] = sphbasis_mono(r,theta,phi,Nse,k,conf) +%SPHBASIS_MONO evaluates spherical basis functions for given input arguments +% +% Usage: [jn,h2n,Ynm] = sphbasis_mono(r,theta,phi,Nse,k,conf) +% +% Input parameters: +% r - distance from origin / m [n1 x n2 x ...] +% theta - elevation angle / rad [n1 x n2 x ...] +% phi - azimuth angle / rad [n1 x n2 x ...] +% Nse - maximum order of spherical basis functions +% k - wave number +% conf - optional configuration struct (see SFS_config) +% +% Output parameters: +% jn - cell array of spherical bessel functions +% h2n - cell array of spherical hankel functions of 2nd kind +% Ynm - cell array of spherical harmonics +% +% SPHBASIS_MONO(r,theta,phi,Nse,k,conf) computes spherical basis functions +% for the given arguments r, theta and phi. r, theta and phi can be of +% arbitrary (but same) size. Output will be stored in cell arrays (one cell +% entry for each order) of length Nse+1 for jn and h2n. For Ynm the lenght +% is (Nse+1).^2. The coefficients of Ynm are stored with the linear index l +% resulting from the order m and the degree n of the spherical harmonics: +% +% m 2 +% Y = Y ; with l = (n+1) - (n - m) +% l n +% +% see also: sphbasis_mono_grid sphharmonics sphbesselj sphbesselh + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 5; +nargmax = 6; +narginchk(nargmin,nargmax); +isargpositivescalar(Nse); +isargequalsize(r,phi,theta); +isargscalar(k); +if nargin. * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 7; +nargmax = 7; +narginchk(nargmin,nargmax); + +isargnumeric(X,Y,Z); +isargpositivescalar(Nse,f); +isargcoord(xq); + +%% ===== Computation ==================================================== +[x,y,z] = xyz_grid(X,Y,Z,conf); + +k = 2*pi*f/conf.c; % wavenumber + +% shift coordinates to expansion coordinate +x = x - xq(1); +y = y - xq(2); +z = z - xq(3); +% coordinate transformation +r = sqrt(x.^2 + y.^2 + z.^2); +phi = atan2(y,x); +theta = asin(z./r); + +[jn,h2n,Ynm] = sphbasis_mono(r,theta,phi,Nse,k,conf); diff --git a/SFS_monochromatic/sphexp/sphexp_access.m b/SFS_monochromatic/sphexp/sphexp_access.m new file mode 100644 index 00000000..6959e531 --- /dev/null +++ b/SFS_monochromatic/sphexp/sphexp_access.m @@ -0,0 +1,92 @@ +function a = sphexp_access(Anm,m1,n1,m2,n2) +%SPHEXP_ACCESS yields array elements of spherical expansion coefficients +% +% Usage: a = sphexp_access(Anm,m1,n1,m2,n2) +% +% Input parameters: +% A - 1D, 2D, 3D array of expansion coefficients (3rd +% dimension for temporal frequency) +% m1 - order of 1st dimension +% n1 - degree of 1st dimension (optional, default = |m1|) +% m2 - order of 2st dimension (optional, default = 0) +% n2 - degree of 2st dimension (optional, default = |m2|) +% +% Output parameters: +% a - coefficients +% +% SPHEXP_ACCESS(Anm,m1,n1,m2,n2) computes one/two indices for accessing +% 1D/2D arrays of spherical expansion coefficients +% A(n1,m1) => A(l1) ; A(n2,m2) => A(l2) +% CAUTION: THIS FUNCTION DOES NOT USE ANY CHECK OF INPUT ARGUMENTS +% +% see also: sphexp_access + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters =================================== +L1 = length(m1); +if nargin < 3 + n1 = abs(m1); +elseif length(n1) < L1 + n1 = repmat(n1,[1 L1]); +elseif length(n1) > L1 + m1 = repmat(m1,[1 length(n2)]); + L1 = length(m2); +end + +if nargin < 4 + m2 = 0; +end + +L2 = length(m2); +if nargin < 5 + n2 = abs(m2); +elseif length(n2) < L2 + n2 = repmat(n2,[1 L2]); +elseif length(n2) > L2 + m2 = repmat(m2,[1 length(n2)]); + L2 = length(m2); +end + +%% ===== Computation ==================================================== +a = zeros(L1,L2,size(Anm,3)); + +s1 = abs(m1) <= n1; +s2 = abs(m2) <= n2; + +if any(s1) && any(s2) + [v,w] = sphexp_index(m1(s1),n1(s1),m2(s2),n2(s2)); + a(s1,s2,:) = Anm(v,w,:); +end + +end diff --git a/SFS_monochromatic/sphexp/sphexp_convert_circexp.m b/SFS_monochromatic/sphexp/sphexp_convert_circexp.m new file mode 100644 index 00000000..9333a12a --- /dev/null +++ b/SFS_monochromatic/sphexp/sphexp_convert_circexp.m @@ -0,0 +1,67 @@ +function Anm = sphexp_convert_circexp(Am) +%SPHEXP_CONVERT_CIRCEXP compute converts regular circular expansion into +%spherical expansion coefficients +% +% Usage: Anm = sphexp_convert_circexp(Am) +% +% Input parameters: +% Am - regular circular expansion coefficients +% +% Output parameters: +% Anm - regular spherical expansion coefficients +% +% References: +% Hahn, Spors (2015) - "Sound Field Synthesis of Virtual Cylindrical +% Waves using Circular and Spherical Loudspeaker +% Arrays", 138th AES Convention + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 1; +nargmax = 1; +narginchk(nargmin,nargmax); +isargvector(Am); +Nce = (length(Am)-1)/2; + +%% ===== Computation ==================================================== + +% Implementation of Hahn2015, Eq. (14) +Anm = zeros((Nce+1).^2,1); +for m=-Nce:Nce + % for theta=0 the legendre polynom is zero if n+m is odd + for n = abs(m):2:Nce + v = sphexp_index(m,n); + Anm(v) = 4*pi.*1j.^(m-n).*sphharmonics(n,-m,0,0).*Am(m+Nce+1); + end +end diff --git a/SFS_monochromatic/sphexp/sphexp_index.m b/SFS_monochromatic/sphexp/sphexp_index.m new file mode 100644 index 00000000..3b7eb9cf --- /dev/null +++ b/SFS_monochromatic/sphexp/sphexp_index.m @@ -0,0 +1,69 @@ +function [l1,l2] = sphexp_index(m1,n1,m2,n2) +%SPHEXP_INDEX calculates index(indices) for the arrays of spherical expansion +%coefficients +% +% Usage: [l1,l2] = sphexp_index(m1,n1,m2,n2) +% +% Input parameters: +% m1 - order of 1st dimension +% n1 - degree of 1st dimension (optional, default = |m1|) +% m2 - order of 2st dimension (optional, default = 0) +% n2 - degree of 2st dimension (optional, default = |m2|) +% +% Output parameters: +% l1 - index for 1st dimension +% l2 - index for 2st dimension +% +% SPHEXP_INDEX(m1,n1,m2,n2) computes one/two indices for accessing +% 1D/2D arrays of spherical expansion coefficients +% A(n1,m1) => A(l1) ; A(n2,m2) => A(l2) +% CAUTION: THIS FUNCTION DOES NOT USE ANY CHECK OF INPUT ARGUMENTS +% +% see also: sphexp_access + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +if nargin < 2 + n1 = abs(m1); +end +if nargin < 3 + m2 = 0; +end +if nargin < 4 + n2 = abs(m2); +end + +%% ===== Computation ==================================================== +l1 = (n1 + 1).^2 - (n1 - m1); +l2 = (n2 + 1).^2 - (n2 - m2); diff --git a/SFS_monochromatic/sphexp/sphexp_mono_cht.m b/SFS_monochromatic/sphexp/sphexp_mono_cht.m new file mode 100644 index 00000000..17fafddf --- /dev/null +++ b/SFS_monochromatic/sphexp/sphexp_mono_cht.m @@ -0,0 +1,82 @@ +function Pnm = sphexp_mono_cht(Dm,mode,f,conf) +%SPHEXP_MONO_CHT yields spherical expansion coefficients of a sound field +%resulting from of a driving function given as a circular harmonics transform +% +% Usage: Pnm = sphexp_mono_cht(Dm,mode,f,conf) +% +% Input parameters: +% Dm - circular harmonics transform of driving function +% mode - 'R' for regular, 'S' for singular +% f - frequency / Hz [Nf x 1] or [1 x Nf] +% conf - configuration struct (see SFS_config) +% +% Output parameters: +% Pnm - spherical expansion coefficients of a synthesized +% sound field reproduced by nfchoa driving function +% +% SPHEXP_MONO_CHT(Dm,mode,f,conf) computes the spherical expansion +% coefficients of sound field reproduced by a circular secondary source +% distribution consisting of SPHERICAL monopoles. The driving function is +% given as its circular harmonics transform. + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 4; +nargmax = 4; +narginchk(nargmin,nargmax); +isargmatrix(Dm); +isargvector(f); +isargchar(mode); +if ~strcmp('R',mode) && ~strcmp('S',mode) + error('%s: unknown mode (%s)!',upper(mfilename),mode); +end +isargstruct(conf); + +%% ===== Configuration ================================================== +R0 = conf.secondary_sources.size / 2; + +%% ===== Variables ====================================================== +Nse = (size(Dm,1)-1)/2; + +%% ===== Computation ==================================================== +% regular/singular spherical expansion of 3D Green's function +Gnm = sphexp_mono_ps([R0,0,0],mode,Nse,f,[0,0,0],conf); +% regular/singular spherical expansion of synthesised sound field +Pnm = zeros(size(Gnm)); +for n=0:Nse + v = sphexp_index(-n:n,n); + Pnm(v,:) = 2*pi*R0*Gnm(v,:).*Dm((-n:n)+Nse+1,:); +end + +end diff --git a/SFS_monochromatic/sphexp/sphexp_mono_ls.m b/SFS_monochromatic/sphexp/sphexp_mono_ls.m new file mode 100644 index 00000000..3e6bfd93 --- /dev/null +++ b/SFS_monochromatic/sphexp/sphexp_mono_ls.m @@ -0,0 +1,75 @@ +function Anm = sphexp_mono_ls(xs,mode,Nse,f,xq,conf) +%SPHEXP_MONO_LS computes the regular/singular spherical expansion of line source +% +% Usage: Anm = sphexp_mono_ls(xs,mode,Nse,f,xq,conf) +% +% Input parameters: +% xs - position of line source +% mode - 'R' for regular, 'S' for singular +% Nse - maximum order of spherical basis functions +% f - frequency +% xq - optional expansion center coordinate +% conf - optional configuration struct (see SFS_config) +% +% Output parameters: +% Al - regular Spherical Expansion Coefficients +% +% SPHEXP_MONO_LS(xs,mode,Nse,f,xq,conf) computes the regular/singular +% spherical expansion coefficients for a point source at xs. The expansion +% will be done around the expansion coordinate xq. +% +% see also: sphexp_mono_ps sphexp_mono_pw sphexp_convert_circexp + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 6; +nargmax = 6; +narginchk(nargmin,nargmax); +isargposition(xs); +isargpositivescalar(f,Nse); +isargposition(xq); +isargstruct(conf); + +%% ===== Computation ==================================================== +% select suitable basis function +if strcmp('R',mode) + Am = circexp_mono_ls(xs,mode,Nse,f,xq,conf); + Anm = sphexp_convert_circexp(Am); +elseif strcmp('S',mode) + to_be_implemented; +else + error('%s: %s, unknown mode',upper(mfilename),mode); +end + +end diff --git a/SFS_monochromatic/sphexp/sphexp_mono_multiscatter.m b/SFS_monochromatic/sphexp/sphexp_mono_multiscatter.m new file mode 100644 index 00000000..da6d11e7 --- /dev/null +++ b/SFS_monochromatic/sphexp/sphexp_mono_multiscatter.m @@ -0,0 +1,101 @@ +function B = sphexp_mono_multiscatter(A,xq,R,sigma,f,conf) +%SPHEXP_MONO_MULTISCATTER compute the singular spherical expansion of a sphere- +%scattered sound field +% +% Usage: B = sphexp_mono_multiscatter(A,xq,R,sigma,f,conf) +% +% Input parameters: +% A - regular spherical expansion coefficients of incident +% sound fields arriving at each sphere [n x Nq] +% xq - positions of the sphere / m [Nq x 3] +% R - radii of the spheres / m [Nq x 1] +% sigma - admittances of the sphere / m [Nq x 1] +% f - frequency / Hz +% +% Output parameters: +% B - regular spherical expansion coefficients of sound fields +% scattered at each sphere [n x Nq] +% +% SPHEXP_MONO_MULTISCATTER(A,xq,R,sigma,f,conf) + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 6; +nargmax = 6; +narginchk(nargmin,nargmax); +isargpositivescalar(f); +isargvector(R,sigma); +isargmatrix(A,xq); +isargstruct(conf); + +L = size(A,1); +isargsquaredinteger(L); + +if length(R) == 1 + R = repmat( R,[1,size(A,2)] ); +elseif length(R) ~= size(A,2) + error('%s: Length of R does not match size of A',upper(mfilename)); +end + +if length(sigma) == 1 + sigma = repmat( sigma,[1,size(A,2)] ); +elseif length(sigma) ~= size(A,2) + error('%s: Length of R does not match size of A',upper(mfilename)); +end + +%% ===== Computation ==================================================== +Nq = size(A,2); +AB = zeros(Nq*L); + +E = ones(L,1); +for qdx=1:Nq + selectq = ((qdx-1)*L+1):(qdx*L); + AB(selectq,selectq) = ... + diag(1./sphexp_mono_scatter(E,R(qdx),sigma(qdx),f,conf)); +end + +for qdx=1:Nq + selectq = ((qdx-1)*L+1):(qdx*L); + for pdx=(qdx+1):Nq + selectp = ((pdx-1)*L+1):(pdx*L); + [SRpq,SRqp] = sphexp_mono_translation(xq(qdx,:)-xq(pdx,:),'SR',f,conf); + AB(selectp,selectq) = -SRpq; + AB(selectq,selectp) = -SRqp; + end +end + +A = reshape(A,[],1); + +B = AB\A; +B = reshape(B,L,Nq); diff --git a/SFS_monochromatic/sphexp/sphexp_mono_ps.m b/SFS_monochromatic/sphexp/sphexp_mono_ps.m new file mode 100644 index 00000000..60698a3d --- /dev/null +++ b/SFS_monochromatic/sphexp/sphexp_mono_ps.m @@ -0,0 +1,134 @@ +function Anm = sphexp_mono_ps(xs,mode,Nse,f,xq,conf) +%SPHEXP_MONO_PS computes regular/singular spherical expansion of point source +% +% Usage: Anm = sphexp_mono_ps(xs,mode,Nse,f,xq,conf) +% +% Input parameters: +% xs - position of point source +% mode - 'R' for regular, 'S' for singular +% Nse - maximum order of spherical basis functions +% f - frequency [Nf x 1] or [1 x Nf] +% xq - optional expansion center coordinate +% conf - optional configuration struct (see SFS_config) +% +% Output parameters: +% Anm - regular Spherical Expansion Coefficients [(Nse+1)^2 x m] +% +% SPHEXP_MONO_PS(xs,mode,f,Nse,xq,conf) computes the regular/singular +% Spherical Expansion Coefficients for a point source at xs. The expansion +% will be done around the expansion coordinate xq: +% +% Regular Expansion: +% \~~ oo \~~ n m m +% p (x,f) = > > A R (x-x ) +% ps,R /__ n=0 /__ m=-n n n q +% +% with the expansion coefficients (Gumerov2004, eq. 3.2.2): +% m -m +% A = -i . k . S (x - x ) +% n n s q +% +% Singular Expansion: +% \~~ oo \~~ n m m +% p (x,f) = > > B S (x-x ) +% ps,S /__ n=0 /__ m=-n n n q +% +% with the expansion coefficients (Gumerov2004, eq. 3.2.2): +% m -m +% B = -i . k . R (x - x ) +% n n s q +% +% The coefficients are stored in linear arrays with index l resulting from +% m and n: +% +% m m 2 +% A = A ; B = B with l = (n+1) - (n - m) +% l n l n +% +% References: +% Gumerov,Duraiswami (2004) - "Fast Multipole Methods for the +% Helmholtz Equation in three +% Dimensions", ELSEVIER +% +% see also: sphexp_mono_ls sphexp_mono_pw + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 6; +nargmax = 6; +narginchk(nargmin,nargmax); +isargposition(xs); +isargpositivescalar(Nse); +isargvector(f); +isargposition(xq); +isargstruct(conf); + +%% ===== Configuration ================================================== +c = conf.c; + +%% ===== Variables ====================================================== +% convert (xs-xq) into spherical coordinates +r = sqrt(sum((xs-xq).^2)); +phi = atan2(xs(2)-xq(2),xs(1)-xq(1)); +theta = asin((xs(3)-xq(3))/r); + +% frequency +k = 2.*pi.*row_vector(f)./c; +kr = k.*r; +Nf = length(kr); + +% select suitable basis function +if strcmp('R',mode) + sphbasis = @(nu) sphbesselh(nu,2,kr); +elseif strcmp('S',mode) + sphbasis = @(nu) sphbesselj(nu,kr); +else + error('unknown mode:'); +end + +%% ===== Computation ==================================================== +Anm = zeros((Nse + 1).^2,Nf); +for n=0:Nse + cn = -1j.*k.*sphbasis(n); + for m=0:n + % spherical harmonics: conj(Y_n^m) = Y_n^-m (Gumerov2004, eq. 2.1.59) + Ynm = sphharmonics(n,m,theta,phi); + % -m + v = sphexp_index(-m,n); + Anm(v,:) = cn.*Ynm; + % +m + v = sphexp_index(m,n); + Anm(v,:) = cn.*conj(Ynm); + end +end diff --git a/SFS_monochromatic/sphexp/sphexp_mono_pw.m b/SFS_monochromatic/sphexp/sphexp_mono_pw.m new file mode 100644 index 00000000..b346463f --- /dev/null +++ b/SFS_monochromatic/sphexp/sphexp_mono_pw.m @@ -0,0 +1,111 @@ +function Anm = sphexp_mono_pw(npw,Nse,f,xq,conf) +%SPHEXP_MONO_PW computes the regular spherical expansion of plane wave +% +% Usage: Anm = sphexp_mono_pw(npw,Nse,f,xq,conf) +% +% Input parameters: +% npw - unit vector propagation direction of plane wave +% Nse - maximum order of spherical basis functions +% f - frequency / Hz [Nf x 1] or [1 x Nf] +% xq - optional expansion coordinate +% conf - optional configuration struct (see SFS_config) +% +% Output parameters: +% Anm - regular Spherical Expansion Coefficients [(Nse+1)^2 x Nf] +% +% SPHEXP_MONO_PW(npw,Nse,f,xq,conf) computes the regular spherical +% expansion coefficients for a plane wave. The expansion will be done around +% the expansion coordinate xq: +% +% \~~ oo \~~ n m m +% p (x,f) = > > A R (x-x ) +% pw /__ n=0 /__ m=-n n n q +% +% with the expansion coefficients (Gumerov, p. 74, eq. 2.3.6): +% +% m -n -m +% A = 4pi i Y (theta , phi ) +% n n pw pw +% +% The coefficients are stored in linear arrays with index l resulting from +% m and n: +% +% m 2 +% A = A ; with l = (n+1) - (n - m) +% l n +% +% References: +% Gumerov,Duraiswami (2004) - "Fast Multipole Methods for the +% Helmholtz Equation in three +% Dimensions", ELSEVIER +% +% see also: sphexp_mono_ls sphexp_mono_ps + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 5; +nargmax = 5; +narginchk(nargmin,nargmax); +isargposition(npw); +isargstruct(conf); +isargposition(xq); +isargpositivescalar(Nse); +isargvector(f); + +%% ===== Configuration ================================================== +c = conf.c; + +%% ===== Computation ==================================================== +% convert npw into spherical coordinates +phi = atan2(npw(2),npw(1)); +theta = asin(npw(3)); + +% phase shift due to expansion coordinate +phase = exp(-1j*2*pi*row_vector(f)/c*(npw*xq.')); +Nf = length(phase); + +Anm = zeros((Nse + 1).^2,Nf); +for n=0:Nse + cn = 4*pi*(1j)^(-n); + for m=0:n + % spherical harmonics: conj(Y_n^m) = Y_n^-m (Gumerov2004, eq. 2.1.59) + Ynm = sphharmonics(n,m,theta,phi); + % -m + v = sphexp_index(-m,n); + Anm(v,:) = cn.*Ynm.*phase; + % +m + v = sphexp_index(m,n); + Anm(v,:) = cn.*conj(Ynm).*phase; + end +end diff --git a/SFS_monochromatic/sphexp/sphexp_mono_scatter.m b/SFS_monochromatic/sphexp/sphexp_mono_scatter.m new file mode 100644 index 00000000..0abf739f --- /dev/null +++ b/SFS_monochromatic/sphexp/sphexp_mono_scatter.m @@ -0,0 +1,124 @@ +function Bnm = sphexp_mono_scatter(Anm,R,sigma,f,conf) +%SPHEXP_MONO_SCATTER compute the singular spherical expansion of a sphere- +%scattered sound field +% +% Usage: Bnm = sphexp_mono_scatter(Anm,R,sigma,f,conf) +% +% Input parameters: +% Anm - regular spherical expansion of incident sound field +% R - radius of sphere / m +% sigma - admittance of sphere (from 0(sound soft) to inf(sound hard)) +% f - frequency / Hz [Nf x 1] or [1 x Nf] +% +% Output parameters: +% Bnm - singular spherical expansion coefficients of scattered +% field +% +% SPHEXP_MONO_SCATTER(Anm,R,sigma,f,conf) computes the singular spherical +% expansion coefficients of a field resulting from a scattering of an incident +% field at a sphere. Incident field is descriped by regular expansion +% coefficients (expansion center is expected to be at the center of the sphere +% xq): +% +% \~~ oo \~~ n m m +% P (x,f) = > > A R (x-x ) +% ind /__ n=0 /__ m=-n n n q +% +% The scattered field is descriped by singular expansion coefficients, +% expanded around the center of the sphere xq. +% +% \~~ oo \~~ n m m +% P (x,f) = > > B S (x-x ) +% sca /__ n=0 /__ m=-n n n q +% +% Due to the boundary conditions on the surface of the sphere the +% coefficients are related by: +% +% k . j' (kR) + sigma . j (kR) +% m n n m +% B = - ------------------------------ . A +% n k . h' (kR) + sigma . h (kR) n +% n n +% +% where k = 2*pi*f/c. The coefficients are stored in linear arrays with +% index l resulting from m and n: +% +% m m 2 +% A = A ; B = B with l = (n+1) - (n - m) +% l n l n +% +% References: +% Gumerov,Duraiswami (2004) - "Fast Multipole Methods for the +% Helmholtz Equation in three +% Dimensions", ELSEVIER +% +% see also: sphexp_mono_ps, sphexp_mono_pw, sphexp_mono_timereverse + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 5; +nargmax = 5; +narginchk(nargmin,nargmax); +isargvector(Anm); +L = length(Anm); +isargsquaredinteger(L); +isargscalar(sigma); +isargpositivescalar(R); +isargvector(f); +isargstruct(conf); + +%% ===== Configuration ================================================== +c = conf.c; + +%% ===== Variables ====================================================== +% frequency +k = 2.*pi.*row_vector(f)./c; +kR = k.*R; + +% select suitable transformation function +if isinf(sigma) + T = @(x) -sphbesselj(x,kR)./sphbesselh(x,2,kR); +elseif sigma == 0 + T = @(x) -sphbesselj_derived(x,kR)./sphbesselh_derived(x,2,kR); +else + T = @(x) -(k.*sphbesselj_derived(x,kR)+sigma.*sphbesselj(x,kR)) ... + ./(k.*sphbesselh_derived(x,2,kR)+sigma.*sphbesselh(x,2,kR)); +end + +%% ===== Computation ==================================================== +Bnm = zeros(size(Anm)); +for n=0:(sqrt(L) - 1) + v = sphexp_index(-n:n,n); + Bnm(v,:) = T(n).*Anm(v,:); +end diff --git a/SFS_monochromatic/sphexp/sphexp_mono_sht.m b/SFS_monochromatic/sphexp/sphexp_mono_sht.m new file mode 100644 index 00000000..89751155 --- /dev/null +++ b/SFS_monochromatic/sphexp/sphexp_mono_sht.m @@ -0,0 +1,99 @@ +function Pnm = sphexp_mono_sht(Dnm,mode,f,conf) +%SPHEXP_MONO_NFCHOA_SHT yields spherical expansion coefficients of a sound field +%resulting from of a driving function given as a spherical harmonics transform +% +% Usage: Pnm = sphexp_mono_sht(Dnm,mode,f,conf) +% +% Input parameters: +% Dnm - spherical harmonics transform of nfchoa driving function +% mode - 'R' for regular, 'S' for singular +% f - frequency / Hz [Nf x 1] or [1 x Nf] +% conf - optional configuration struct (see SFS_config) +% +% Output parameters: +% Pnm - spherical expansion coefficients of a sound field +% reproduced by nfchoa driving function +% +% SPHEXP_MONO_SHT(Dnm,mode,f,conf) + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 4; +nargmax = 4; +narginchk(nargmin,nargmax); +isargmatrix(Dnm); +isargsquaredinteger(size(Dnm,1)); +isargvector(f); +isargchar(mode); +if ~strcmp('R',mode) && ~strcmp('S',mode) + error('%s: unknown mode (%s)!',upper(mfilename),mode); +end +isargstruct(conf); + +%% ===== Configuration ================================================== +r0 = conf.secondary_sources.size / 2; +dimension = conf.dimension; + +%% ===== Variables ====================================================== +Nse = sqrt(size(Dnm,1)) - 1; + +%% ===== Computation ==================================================== +if strcmp('2D',dimension) + % === 2-Dimensional ================================================== + + error('%s: 2D not supported.',upper(mfilename)); + +elseif strcmp('2.5D',dimension) + % === 2.5-Dimensional ================================================ + + % regular/singular expansion of 3D Green's function + Gnm = 2*pi*r0*sphexp_mono_ps([r0,0,0],mode,Nse,f,[0,0,0],conf); +elseif strcmp('3D',dimension) + % === 3-Dimensional ================================================== + + % regular/singular expansion of 3D Green's function + Gnm = sphexp_mono_ps([0,0,r0],mode,Nse,f,[0,0,0],conf); + + for n=0:Nse + v = sphexp_index(-n:n,n); + w = sphexp_index(0,n); + Gnm(v,:) = 2*pi*r0^2*sqrt(4*pi / (2*n+1))*Gnm(w,:); + end +else + error('%s: the dimension %s is unknown.',upper(mfilename),dimension); +end + +Pnm = Gnm .* Dnm; + +end diff --git a/SFS_monochromatic/sphexp/sphexp_mono_timereverse.m b/SFS_monochromatic/sphexp/sphexp_mono_timereverse.m new file mode 100644 index 00000000..2fb1949c --- /dev/null +++ b/SFS_monochromatic/sphexp/sphexp_mono_timereverse.m @@ -0,0 +1,58 @@ +function ABnm = sphexp_mono_timereverse(ABnm) +%SPHEXP_MONO_TIMEREVERSE compute coefficients of a time reversed sound field +% +% Usage: ABnm = sphexp_mono_timereverse(ABnm) +% +% Input parameters: +% ABnm - spherical expansion coefficients of sound field +% +% Output parameters: +% ABnm - spherical expansion coefficients of time reversed +% sound field + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 1; +nargmax = 1; +narginchk(nargmin,nargmax); +isargvector(ABnm); +L = length(ABnm); +isargsquaredinteger(L); + +%% ===== Computation ==================================================== + +for n=0:(sqrt(L)-1) + v = sphexp_index(-n:n,n); + ABnm(v) = conj(sphexp_access(ABnm,n:-1:-n,n)); +end diff --git a/SFS_monochromatic/sphexp/sphexp_mono_translation.m b/SFS_monochromatic/sphexp/sphexp_mono_translation.m new file mode 100644 index 00000000..1a3f01cb --- /dev/null +++ b/SFS_monochromatic/sphexp/sphexp_mono_translation.m @@ -0,0 +1,264 @@ +function [EF,EFm] = sphexp_mono_translation(t,mode,Nse,f,conf) +%SPHEXP_MONO_TRANSLATION computes spherical translation coefficients +%(multipole re-expansion) +% +% Usage: [EF,EFm] = sphexp_mono_translation(t,mode,Nse,f,conf) +% +% Input parameters: +% t - translatory shift [1x3] / m +% mode - 'RS' for regular-to-singular reexpansion +% 'RR' for regular-to-regular reexpansion +% 'SR' for singular-to-regular reexpansion +% 'SS' for singular-to-singular reexpansion +% Nse - maximum order of spherical basis functions +% f - frequency [1 x Nf] / Hz +% conf - configuration struct (see SFS_config) +% +% Output parameters: +% EF - spherical re-expansion coefficients for t +% EFm - spherical re-expansion coefficients for -t +% +% SPHEXP_MONO_TRANSLATION(t,mode,Nse,f,conf) computes the spherical +% re-expansion coefficients to perform as translatory shift of spherical basis +% function. Multipole Re-expansion computes the spherical basis function for a +% shifted coordinate system (x+t) based on the original basis functions for +% (x). +% +% m \~~ inf \~~ l s,m s +% E (x + t) = > > (E|F) (t) F (x) +% n /__ l=0 /__ s=-l l,n l +% +% where {E,F} = {R,S}. R denotes the regular spherical basis function, while +% S symbolizes the singular spherical basis function. Note that (S|S) and +% (S|R) are respectively equivalent to (R|R) and (R|S). The reexpansion +% coefficients can seperated into sectorial (n = abs|m| and/or l = abs|s|) +% and tesseral (else) coefficients. Latter will only be computed, if +% conf.dimensions == '3D'. +% +% References: +% Gumerov,Duraiswami (2004) - "Fast Multipole Methods for the +% Helmholtz Equation in three +% Dimensions", ELSEVIER +% +% see also: sphexp_mono_ps, sphexp_mono_pw + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 5; +nargmax = 5; +narginchk(nargmin,nargmax); +isargposition(t); +isargchar(mode); +isargpositivescalar(Nse); +isargvector(f); +isargstruct(conf); + +%% ===== Configuration ================================================== +showprogress = conf.showprogress; +dimension = conf.dimension; + +%% ===== Variables ====================================================== +% convert t into spherical coordinates +r = sqrt(sum(t.^2)); +phi = atan2(t(2),t(1)); +theta = asin(t(3)./r); + +% frequency +Nf = length(f); +k = 2*pi*f/conf.c; +kr = reshape(k,1,1,Nf)*r; + +L = (2*Nse + 1).^2; +S = zeros(L,L,Nf); + +% Auxilary Coefficients for Computations +[a,b] = sphexp_translation_auxiliary(2*Nse,conf); + +% select suitable basis function +if strcmp('RR',mode) || strcmp('SS',mode) + sphbasis = @(nu) sphbesselj(nu,kr); +elseif strcmp('SR',mode) || strcmp('RS',mode) + sphbasis = @(nu) sphbesselh(nu,2,kr); +else + error('unknown mode:'); +end + +%% ===== Sectorial Coefficients ========================================= +% if translation vector's z-coordinate is zero, many coefficients are zero. +% This is to save some computation time +if t(3) == 0 + inc = 2; +else + inc = 1; +end + +% for n=0, m=0 (Gumerov2004, eq. 3.2.5) +% +% s, 0 s, 0 ___ l -s +% (S|R) = (R|S) = \|4pi (-1) S (t) +% l, 0 l, 0 l +% +% s, 0 s, 0 ___ l -s +% (S|S) = (R|R) = \|4pi (-1) R (t) +% l, 0 l, 0 l +% +% for s=0, l=0 +% +% 0, n 0, m ___ m +% (S|R) = (R|S) = \|4pi S (t) +% 0, m 0, n n +% +% 0, n 0, m ___ m +% (S|S) = (R|R) = \|4pi R (t) +% 0, m 0, n n +% +for l=0:2*Nse % n=l + Hn = sqrt(4*pi)*sphbasis(l); % radial basis function (see Variables) + for s=l:-inc:0 % m=s + % spherical harmonics: conj(Y_n^m) = Y_n^-m (Gumerov2004, eq. 2.1.59) + Ynm = sphharmonics(l,s,theta,phi); + + v = sphexp_index(s,l); + S(v,1,:) = (-1).^l*Hn.*conj(Ynm); % s,l, m=0, n=0 + S(1,v,:) = Hn.*Ynm; % s=0,l=0, m, n + + v = sphexp_index(-s,l); + S(v,1,:) = (-1).^l*Hn.*Ynm; % -s,l, m=0, n=0 + S(1,v,:) = Hn.*conj(Ynm); % s=0,l=0, -m, n + end + if showprogress, progress_bar(v,L); end % progress bar +end + +%% ===== Sectorial Coefficients ========================================== +% for n=|m| (Gumerov2004, eq. 3.2.78) +% +% -m-1 s,m+1 -s s-1,m s-1 s-1,m +% b (E|F)(t) = b (E|F)(t) - b (E|F)(t) +% m+1 l,|m+1| l l-1,|m| l+1 l+1,|m| +% +% -m-1 s,-m-1 s s+1,-m -s-1 s+1,-m +% b (E|F)(t) = b (E|F)(t) - b (E|F)(t) +% m+1 l,|m+1| l l-1,|m| l+1 l+1,|m| +% +% while {E,F} = {R,S} +% +for m=0:Nse-1 + % NOTE: sphexp_access(b, -m-1) == sphexp_access(b, -m-1, m+1) + bm = 1./sphexp_access(b,-m-1); + for l=1:(2*Nse-m-1) + for s=-l:inc:l + % +m + [v,w] = sphexp_index(s,l,m+1); + S(v,w,:) = bm * ( ... + sphexp_access(b,-s ,l) * sphexp_access(S,s-1,l-1,m) - ... + sphexp_access(b,s-1,l+1) * sphexp_access(S,s-1,l+1,m)... + ); + % -m + [v,w] = sphexp_index(s,l,-m-1); + S(v,w,:) = bm * ( ... + sphexp_access(b,s ,l) * sphexp_access(S,s+1,l-1,-m) - ... + sphexp_access(b,-s-1,l+1) * sphexp_access(S,s+1,l+1,-m)... + ); + end + end + if showprogress, progress_bar(m,Nse); end % progress bar +end + +% for l=|s| using symmetry relation (Gumerov2004, eq. 3.2.49 +% +% s,m |s|+n -m,-s +% (E|F)(t) = (-1) (E|F)(t) +% |s|,n n,|s| +% +% while {E,F} = {R,S} +% +for l=1:Nse + for n=inc:(2*Nse-l) + s = [-l,l]; + m = (-n+inc):inc:(n-inc); + [v,w] = sphexp_index( s,l,m,n); + S(v,w,:) = (-1).^(l+n)*sphexp_access(S,-m,n,-s,l).'; + end +end + +%% ===== Tesseral Coefficients ========================================== +% TODO: there is a bug somewhere in here +if strcmp('3D', dimension) + for m=-Nse:Nse + for s=-Nse:Nse + + lowerbound = Nse - max(abs(m),abs(s)); + % left propagation + for n=(0:lowerbound-1)+abs(m) + amn1 = 1./sphexp_access(a,m,n); + amn2 = sphexp_access(a,m,n-1); + for l=(abs(s)-abs(m)+n+1):(2*Nse-n-1) + [v, w] = sphexp_index(m, n+1, s, l); + S(v,w,:) = amn1 * ( ... + amn2 * sphexp_access(S,m,n-1,s,l) - ... + sphexp_access(a,s,l) * sphexp_access(S,m,n ,s,l+1) + ... + sphexp_access(a,s,l-1) * sphexp_access(S,m,n ,s,l-1) ... + ); + end + end + end + % up propagation + for l=(0:lowerbound-1)+abs(s) + asl1 = 1./sphexp_access(a,s,l); + asl2 = sphexp_access(a,s,l-1); + for n=(abs(m)-abs(s)+l+2):(2*Nse-l-1) + [v, w] = sphexp_index(s,l+1, m, n); + S(v,w,:) = asl1 * ( ... + asl2 * sphexp_access(S,s,l-1,m,n) - ... + sphexp_access(a,m,n) * sphexp_access(S,s,l ,m,n+1) + ... + sphexp_access(a,m,n-1) * sphexp_access(S,s,l ,m,n-1)... + ); + end + end + if showprogress, progress_bar(m+Nse,2*Nse); end % progress bar + end +end +%% ====== Final Calculation Steps ======================================= +L = (Nse + 1)^2; +EF = S(1:L,1:L,:); % (E|F)(t) +% SR(-t) +EFm = zeros(L,L,Nf); +for n=0:Nse + for l=0:Nse + m = -n:inc:n; + s = -l:inc:l; + [v,w] = sphexp_index(s,l,m,n); + EFm(v,w,:) = (-1).^(l+n)*sphexp_access(EF,s,l,m,n); + end +end diff --git a/SFS_monochromatic/sphexp/sphexp_translation_auxiliary.m b/SFS_monochromatic/sphexp/sphexp_translation_auxiliary.m new file mode 100644 index 00000000..ca85a7b1 --- /dev/null +++ b/SFS_monochromatic/sphexp/sphexp_translation_auxiliary.m @@ -0,0 +1,114 @@ +function [a,b] = sphexp_translation_auxiliary(Nse,conf) +%SPHEXP_TRANSLATION_AUXILIARY yields the auxiliary coefficients for the +%recursive calculation of spherical translation coefficients +% +% Usage: [a,b] = sphexp_translation_auxiliary(Nse,conf) +% +% Input parameters: +% Nse - maximum degree of the auxiliary functions (optional) +% conf - configuration struct (see SFS_config) +% +% Output parameters: +% a - auxiliary coefficients for the calculation of tesseral +% spherical translation coefficients +% b - auxiliary coefficients for the calculation of sectorial +% spherical translation coefficients +% +% SPHEXP_TRANSLATION_AUXILIARY(Nse,conf) computes the auxiliary coefficients +% for the calculation of tesseral spherical translation coefficients +% (Gumerov2004, eq. 2.2.8 and 3.2.65): +% +% +------------------+ +% m |(n+1+|m|)(n+1-|m|) _ +% a = |------------------ for |m| < n +% n \| (2n+1)(2n+3) +% +% 0 else +% +% Auxiliary coefficients for the calculation of sectorial spherical +% translation coefficients (Gumerov2004, eq. 2.2.10 and 3.2.78/79): +% +% +------------+ +% |(n-m)(n-m-1) _ _ +% |------------ for 0 < m < n +% \|(2n-1)(2n+1) +% +% +------------+ +% m |(n-m)(n-m-1) _ +% b = - |------------ for -n < m < 0 +% n \|(2n-1)(2n+1) +% +% 0 else +% +% The coefficients are stored in linear arrays with index l resulting from +% m and n: +% +% m m 2 +% a = a ; b = a with l = (n+1) - (n - m) +% l n l n +% +% References: +% Gumerov,Duraiswami (2004) - "Fast Multipole Methods for the +% Helmholtz Equation in three +% Dimensions", ELSEVIER +% +% see also: sphexp_mono_translation + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 2; +nargmax = 2; +narginchk(nargmin,nargmax); +isargpositivescalar(Nse); +isargstruct(conf); + +%% ===== Configuration ================================================== +showprogress = conf.showprogress; + +%% ===== Computation ==================================================== +L = (Nse + 1).^2; +a = zeros(L,1); +b = zeros(L,1); +for n=0:Nse + a_denum = (2*n+1)*(2*n+3); + b_denum = (2*n-1)*(2*n+1); + + m = -n:n; + v = sphexp_index(m,n); + + a(v) = sqrt( (n+1+abs(m)).*(n+1-abs(m))./a_denum ); + b(v) = ( 1-(m<0)*2 ) .* sqrt( (n-m-1).*(n-m)./b_denum ); + + if showprogress, progress_bar(v(end),L); end % progress bar +end diff --git a/SFS_monochromatic/sphexp/sphexp_truncate.m b/SFS_monochromatic/sphexp/sphexp_truncate.m new file mode 100644 index 00000000..6721dde3 --- /dev/null +++ b/SFS_monochromatic/sphexp/sphexp_truncate.m @@ -0,0 +1,82 @@ +function Anm = sphexp_truncate(Pnm,N,M,Mshift) +%SPHEXP_TRUNCATE truncates spherical expansion by setting remaining +%coefficients to zero +% +% Usage: Anm = sphexp_truncate(Pnm,N,M,Mshift) +% +% Input parameters: +% Pnm - 1D array of spherical expansion coefficients [n x Nf] +% N - maximum degree of spherical expansion +% M - maximum order of spherical expansion, default: N +% Mshift - shift for asymmetric trunction with respect to order, +% default: 0 +% +% Output parameters: +% Anm - 1D array of bandlimited spherical expansion +% coefficients [n x Nf] +% +% SPHEXP_TRUNCATE(Pnm,N,M,Mshift) sets coefficients belonging to an degree +% n higher than N or an order m exceeding (-M+Mshift:+M+Mshift) to zero. +% +% see also: sphexp_truncation_order + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 2; +nargmax = 4; +narginchk(nargmin,nargmax); +isargmatrix(Pnm); +isargpositivescalar(N); +switch nargin + case 2 + M = N; + Mshift = 0; + case 3 + isargpositivescalar(M); + Mshift = 0; + case 4 + isargscalar(Mshift); +end +%% ===== Variable ======================================================= +Nse = min(sqrt(size(Pnm,1))-1,N); + +%% ===== Computation ==================================================== +Anm = zeros(size(Pnm)); + +for m=max(-M+Mshift,-Nse):min(M+Mshift,Nse) + v = sphexp_index(m,abs(m):Nse); + Anm(v,:) = Pnm(v,:); +end + +end diff --git a/SFS_monochromatic/sphexp/sphexp_truncation_order.m b/SFS_monochromatic/sphexp/sphexp_truncation_order.m new file mode 100644 index 00000000..7e874274 --- /dev/null +++ b/SFS_monochromatic/sphexp/sphexp_truncation_order.m @@ -0,0 +1,77 @@ +function Nse = sphexp_truncation_order(r,f,nmse,conf) +%SPHEXP_TRUNCATION_ORDER yields the bound of summation for a spherical expansion +%of an arbitrary sound field +% +% Usage: Nse = sphexp_truncation_order(r,f,nmse,conf) +% +% Input parameters: +% r - max 3D distance from expansion center / m +% f - frequency / Hz +% nmse - maximum bound for normalized mean squared error +% conf - configuration struct (see SFS_config) +% +% Output parameters: +% Nse - maximum order for spherical expansion +% +% SPHEXP_TRUNCATION_ORDER(r,f,nmse,conf) yields the order up to which +% a the spherical expansion coefficients of an arbitrary sound field have +% be summed up. For a given frequency (f) the normalized truncation mean +% squared error is below the specified error bound (nmse) at any point inside +% a spherical volume around the expansion center with a radius r. Note, that +% this approximation holds for any sound field and might heavily over-estimate +% the needed order in specific cases. +% +% References: +% Kennedy et al. (2007) - "Intrinsic Limits of Dimensionality and +% Richness in Random Multipath Fields", +% IEEE Transactions on Signal Processing +% +% see also: sphexp_truncation + +%***************************************************************************** +% Copyright (c) 2010-2016 Quality & Usability Lab, together with * +% Assessment of IP-based Applications * +% Telekom Innovation Laboratories, TU Berlin * +% Ernst-Reuter-Platz 7, 10587 Berlin, Germany * +% * +% Copyright (c) 2013-2016 Institut fuer Nachrichtentechnik * +% Universitaet Rostock * +% Richard-Wagner-Strasse 31, 18119 Rostock * +% * +% This file is part of the Sound Field Synthesis-Toolbox (SFS). * +% * +% The SFS is free software: you can redistribute it and/or modify it under * +% the terms of the GNU General Public License as published by the Free * +% Software Foundation, either version 3 of the License, or (at your option) * +% any later version. * +% * +% The SFS is distributed in the hope that it will be useful, but WITHOUT ANY * +% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * +% FOR A PARTICULAR PURPOSE. * +% See the GNU General Public License for more details. * +% * +% You should have received a copy of the GNU General Public License along * +% with this program. If not, see . * +% * +% The SFS is a toolbox for Matlab/Octave to simulate and investigate sound * +% field synthesis methods like wave field synthesis or higher order * +% ambisonics. * +% * +% http://github.com/sfstoolbox/sfs sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +nargmin = 4; +nargmax = 4; +narginchk(nargmin,nargmax); +isargpositivescalar(f,r,nmse); +isargstruct(conf); + +%% ===== Configuration ================================================== +c = conf.c; + +%% ===== Computation ==================================================== +% See Kennedy et al. (eq. 42/43) +lambda = c/f; % wave length +delta = max(0,ceil(log(0.67848/nmse))); % nmse dependent term +Nse = ceil(pi*r*exp(1)/lambda) + delta; diff --git a/SFS_plotting/plot_scatterer.m b/SFS_plotting/plot_scatterer.m new file mode 100644 index 00000000..f418807d --- /dev/null +++ b/SFS_plotting/plot_scatterer.m @@ -0,0 +1,9 @@ +function plot_scatterer(xq,R) + +hold on; +for idx=length(R) + rectangle('Position',[xq(idx,1)-R(idx),xq(idx,2)-R(idx), 2*R(idx), 2*R(idx)],... + 'Curvature',[1 1], ... + 'Linewidth',2); +end +hold off; diff --git a/SFS_scattering/cylexpS_mono_multiscatter.m b/SFS_scattering/cylexpS_mono_multiscatter.m new file mode 100644 index 00000000..80513f1d --- /dev/null +++ b/SFS_scattering/cylexpS_mono_multiscatter.m @@ -0,0 +1,57 @@ +function B = cylexpS_mono_multiscatter(A, xq, R, sigma, f, conf) + +%% ===== Checking of input parameters ================================== +nargmin = 5; +nargmax = 6; +narginchk(nargmin,nargmax); +isargmatrix(A,xq); +isargpositivescalar(f); +isargvector(R, sigma); +if nargin