From a4f0824897515ba2759a70a36eb554d1aa3765e2 Mon Sep 17 00:00:00 2001 From: Fiete Winter Date: Fri, 11 Mar 2016 10:42:55 +0100 Subject: [PATCH 1/5] uniformly moving point source --- SFS_helper/isargsecondarysource.m | 10 ++++---- SFS_monochromatic/greens_function_mono.m | 31 ++++++++++++++++++++++++ 2 files changed, 36 insertions(+), 5 deletions(-) diff --git a/SFS_helper/isargsecondarysource.m b/SFS_helper/isargsecondarysource.m index d7555e54..9ac51303 100644 --- a/SFS_helper/isargsecondarysource.m +++ b/SFS_helper/isargsecondarysource.m @@ -57,9 +57,9 @@ function isargsecondarysource(varargin) 'and weights.'], ... inputname(ii)); end - if ~all(abs(vector_norm(x0(:,4:6),2)-1)<1e-10) - error(['The norm of the direction of %s is not 1. ', ... - 'You can use the direction_vector function to get ', ... - 'correct direction values.'],inputname(ii)); - end +% if ~all(abs(vector_norm(x0(:,4:6),2)-1)<1e-10) +% error(['The norm of the direction of %s is not 1. ', ... +% 'You can use the direction_vector function to get ', ... +% 'correct direction values.'],inputname(ii)); +% end end diff --git a/SFS_monochromatic/greens_function_mono.m b/SFS_monochromatic/greens_function_mono.m index a140b1e1..955a555c 100644 --- a/SFS_monochromatic/greens_function_mono.m +++ b/SFS_monochromatic/greens_function_mono.m @@ -25,6 +25,7 @@ % References: % H. Wierstorf, J. Ahrens, F. Winter, F. Schultz, S. Spors (2015) - % "Theory of Sound Field Synthesis" +% J. Ahrens (2012) - "Analytic Methods of Sound Field Synthesis" % % See also: sound_field_mono @@ -102,6 +103,36 @@ % G = 1/(4*pi) .* (1i*omega/c + 1./r) .* scalar./r.^2 .* exp(-1i*omega/c.*r); +elseif strcmp('mps',src) + % Source model for a uniformly moving point source: retarded 3D Green's + % function. + % + % 1 e^(+-i w/c R) + % G(x-xs,vs,w) = --- -------------- + % 4pi R_1 + % + % vs = v * ns % velocity vector + % M = v/c % Mach number + % xparallel = (x-xs)*ns % component in direction of movement + % R_1 = sqrt( M^2*xorth.^2 + (1-M^2)*|x-xs|.^2 ) + % R = (M*xorth + R1)./(1-M*2); + % + % see: Ahrens (2012), eq.(5.60) + % + v = norm(xs(4:6)); % velocity of sound source + nxs = xs(4:6)./v; % direction of movement + xs = xs(1:3); % + M = v/c; % + % shift coordinates + x = x-xs(1); + y = y-xs(2); + z = z-xs(3); + % component of x in direciton of movement: scalar = x*nxs + xparallel = nxs(1).*x + nxs(2).*y + nxs(3).*z; + + R1 = sqrt( M^2.*xparallel.^2 + (1-M^2)*(x.^2 + y.^2 + z.^2) ); + Rplus = (M*xparallel + R1)./(1-M*2); + G = 1/(4*pi) * exp(-1i*omega/c.*Rplus)./R1; elseif strcmp('ls',src) % Source model for a line source: 2D Green's function. % From 8977aac1fb8dca90b3e2d04983aece19710309a6 Mon Sep 17 00:00:00 2001 From: Fiete Winter Date: Tue, 10 May 2016 20:28:58 +0200 Subject: [PATCH 2/5] fix bug --- SFS_monochromatic/greens_function_mono.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SFS_monochromatic/greens_function_mono.m b/SFS_monochromatic/greens_function_mono.m index 955a555c..5da7a4dc 100644 --- a/SFS_monochromatic/greens_function_mono.m +++ b/SFS_monochromatic/greens_function_mono.m @@ -131,7 +131,7 @@ xparallel = nxs(1).*x + nxs(2).*y + nxs(3).*z; R1 = sqrt( M^2.*xparallel.^2 + (1-M^2)*(x.^2 + y.^2 + z.^2) ); - Rplus = (M*xparallel + R1)./(1-M*2); + Rplus = (M*xparallel + R1)./(1-M.^2); G = 1/(4*pi) * exp(-1i*omega/c.*Rplus)./R1; elseif strcmp('ls',src) % Source model for a line source: 2D Green's function. From 81102b70e0eb10d9017d98c763246752950a4295 Mon Sep 17 00:00:00 2001 From: Hagen Wierstorf Date: Thu, 26 May 2016 12:21:34 +0200 Subject: [PATCH 3/5] Add possibility of one input arg to direction_vector() --- SFS_general/direction_vector.m | 148 ++++++++++++++++++--------------- 1 file changed, 82 insertions(+), 66 deletions(-) diff --git a/SFS_general/direction_vector.m b/SFS_general/direction_vector.m index b8ebd4ed..6405d36d 100644 --- a/SFS_general/direction_vector.m +++ b/SFS_general/direction_vector.m @@ -1,66 +1,82 @@ -function n = direction_vector(x1,x2) -%DIRECTION_VECTOR returns a unit vector pointing from x1 to x2 -% -% Usage: n = direction_vector(x1,x2) -% -% Input parameters: -% x1 - starting point -% x2 - ending point -% -% Output parameters: -% n - unit vector pointing in the direction from x1 to x2 -% -% DIRECTION_VECTOR(x1,x2) calculates the unit vector pointing from the -% n-dimensional point x1 to the n-dimensional point x2. The vectors x1 -% and x2 can each be stored in a matrix containing m different direction -% vectors. -% -% See also: secondary_source_positions - -%***************************************************************************** -% 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); -if size(x1)~=size(x2) - error('%s: x1 and x2 had to have the same size.',upper(mfilename)); -end - - -%% ==== Main ============================================================= -n = zeros(size(x1)); -for ii=1:size(x1,1) - n(ii,:) = (x2(ii,:)-x1(ii,:)) / norm(x2(ii,:)-x1(ii,:)); -end +function directions = direction_vector(x1,x2) +%DIRECTION_VECTOR return unit vector(s) pointing from x1 to x2 +% +% Usage: n = direction_vector(x1,[x2]) +% +% Input parameters: +% x1 - starting point(s) [1xn] or [mxn] +% x2 - ending point(s) [1xn] or [mxn] +% +% Output parameters: +% n - unit vector(s) pointing in the direction(s) from x1 to x2 +% +% DIRECTION_VECTOR(x1,x2) calculates the unit vectors pointing from +% n-dimensional points x1 to the n-dimensional points x2. If only x1 is +% provided the direction vectors from zero to x1 are calculated. +% +% See also: secondary_source_positions + +%***************************************************************************** +% The MIT License (MIT) * +% * +% Copyright (c) 2010-2016 SFS Toolbox Developers * +% * +% Permission is hereby granted, free of charge, to any person obtaining a * +% copy of this software and associated documentation files (the "Software"), * +% to deal in the Software without restriction, including without limitation * +% the rights to use, copy, modify, merge, publish, distribute, sublicense, * +% and/or sell copies of the Software, and to permit persons to whom the * +% Software is furnished to do so, subject to the following conditions: * +% * +% The above copyright notice and this permission notice shall be included in * +% all copies or substantial portions of the Software. * +% * +% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR * +% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * +% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL * +% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER * +% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING * +% FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER * +% DEALINGS IN THE SOFTWARE. * +% * +% The SFS Toolbox allows to simulate and investigate sound field synthesis * +% methods like wave field synthesis or higher order ambisonics. * +% * +% http://sfstoolbox.org sfstoolbox@gmail.com * +%***************************************************************************** + + +%% ===== Checking of input parameters =================================== +nargmin = 1; +nargmax = 2; +narginchk(nargmin,nargmax); +if nargin==1 + x2 = x1; + x1 = zeros(1,size(x2,2)); +else + if size(x1,2)~=size(x2,2) + error('%s: x1 and x2 need to have the same dimension.',upper(mfilename)); + end + if size(x1,1)~=size(x2,1) && ~(size(x1,1)==1 | size(x2,1)==1) + error(['%s: x1 and x2 need to have the same size, or one needs to ', ... + 'be a vector.'],upper(mfilename)); + end +end + + +%% ==== Main ============================================================= +% Made both matrices the same size +m1 = size(x1,1); +m2 = size(x2,1); +m = max(m1,m2); +n = size(x1,2); +if m1>m2 + x2 = repmat(x2,[m 1]); +elseif m2>m1 + x1 = repmat(x1,[m 1]); +end +% Calculate direction vectors +directions = zeros([m n]); +for ii=1:size(x1,1) + directions(ii,:) = (x2(ii,:)-x1(ii,:)) / norm(x2(ii,:)-x1(ii,:)); +end From 801bdae0e3820add02aeb35f8169b63fc739c895 Mon Sep 17 00:00:00 2001 From: Fiete Winter Date: Fri, 27 May 2016 17:25:27 +0200 Subject: [PATCH 4/5] optimized function for retarded time --- SFS_general/retarded_time.m | 86 +++++++++++++++++++++++++++++++++++++ 1 file changed, 86 insertions(+) create mode 100644 SFS_general/retarded_time.m diff --git a/SFS_general/retarded_time.m b/SFS_general/retarded_time.m new file mode 100644 index 00000000..936bf5f1 --- /dev/null +++ b/SFS_general/retarded_time.m @@ -0,0 +1,86 @@ +function [tau, R, Delta] = retarded_time(x, t, vs, conf) +% RETARDED_TIME computes the retarded time for a uniformly moving point source +% +% Usage: [tau, R, Delta] = retarded_time(x, t, vs, conf) +% +% Input parameters: +% x - position of observer / m [nx3] +% t - absolute time / s [nxm] or [1xm] +% vs - velocity vector of sound source / (m/s) [1x3] or [nx3] +% conf - configuration struct (see SFS_config) +% +% Output parameters: +% tau - retarded time between sound source and observer / s [nxm] +% R - retarded distance (R = tau*c) / m [nxm] +% Delta - auxilary distance / m [nxm] +% +% RETARDED_TIMES(x, t, vs, conf) the retarded time for a point source +% uniformly moving along a line. Hereby, vs defines its direction +% and its velocity. The position of the sound source at t=0 is the coordinate +% origin. +% +% See also: + +%***************************************************************************** +% The MIT License (MIT) * +% * +% Copyright (c) 2010-2016 SFS Toolbox Team * +% * +% Permission is hereby granted, free of charge, to any person obtaining a * +% copy of this software and associated documentation files (the "Software"), * +% to deal in the Software without restriction, including without limitation * +% the rights to use, copy, modify, merge, publish, distribute, sublicense, * +% and/or sell copies of the Software, and to permit persons to whom the * +% Software is furnished to do so, subject to the following conditions: * +% * +% The above copyright notice and this permission notice shall be included in * +% all copies or substantial portions of the Software. * +% * +% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR * +% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * +% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL * +% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER * +% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING * +% FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER * +% DEALINGS IN THE SOFTWARE. * +% * +% The SFS Toolbox allows to simulate and investigate sound field synthesis * +% methods like wave field synthesis or higher order ambisonics. * +% * +% http://sfstoolbox.org sfstoolbox@gmail.com * +%***************************************************************************** + +%% ===== Checking of input parameters ================================== +% nargmin = 4; +% nargmax = 4; +% narginchk(nargmin,nargmax); +% isargmatrix(x,t,vs); +% isargstruct(conf); + + +%% ===== Configuration ================================================== +c = conf.c; + + +%% ===== Computation ==================================================== + +v = vector_norm(vs,2); % velocity of sound source [n x 1] or [1 x 1] +nvs = bsxfun(@rdivide,vs,v); % direction of movement [n x 3] or [1 x 3] +M = v/c; % Mach number [n x 1] or [1 x 1] + +% vector between observer and point source d(t) = x - vs*t +% [n x m x 3] +d = bsxfun(@minus, reshape(x,[],1,3), bsxfun(@times,t, reshape(vs,[],1,3))); +% instant distance between sound source and observer position +r = vector_norm(d, 3); % [n x m] +% component of x - v*t parallel to the direction of movement: (x-v*t)^T * nvs +% [n x m] +xpara = sum( bsxfun(@times, d, reshape(nvs,[],1,3)), 3); +% auxiliary distance +Delta = sqrt( M.^2.*xpara.^2 + (1-M.^2).*(r.^2) ); % [n x 1] +% retarded distance +R = (M.*xpara + Delta)./(1-M.^2); +% retarded time +tau = R./c; + +end From b51c43da8087ebffd7276d803680ea34a2436b05 Mon Sep 17 00:00:00 2001 From: Fiete Winter Date: Fri, 27 May 2016 18:12:24 +0200 Subject: [PATCH 5/5] improve retarded time function --- SFS_general/retarded_time.m | 30 ++++++++++++++---------- SFS_monochromatic/greens_function_mono.m | 30 ++++++++++-------------- 2 files changed, 30 insertions(+), 30 deletions(-) diff --git a/SFS_general/retarded_time.m b/SFS_general/retarded_time.m index 936bf5f1..494203b5 100644 --- a/SFS_general/retarded_time.m +++ b/SFS_general/retarded_time.m @@ -1,10 +1,12 @@ -function [tau, R, Delta] = retarded_time(x, t, vs, conf) +function [tau, R, Delta] = retarded_time(x, y, z, t, vs, conf) % RETARDED_TIME computes the retarded time for a uniformly moving point source % -% Usage: [tau, R, Delta] = retarded_time(x, t, vs, conf) +% Usage: [tau, R, Delta] = retarded_time(x, y, z, t, vs, conf) % % Input parameters: -% x - position of observer / m [nx3] +% x - x coordinate position of observer / m [nx1] or [1x1] +% y - y coordinate position of observer / m [nx1] or [1x1] +% z - z coordinate position of observer / m [nx1] or [1x1] % t - absolute time / s [nxm] or [1xm] % vs - velocity vector of sound source / (m/s) [1x3] or [nx3] % conf - configuration struct (see SFS_config) @@ -14,7 +16,7 @@ % R - retarded distance (R = tau*c) / m [nxm] % Delta - auxilary distance / m [nxm] % -% RETARDED_TIMES(x, t, vs, conf) the retarded time for a point source +% RETARDED_TIMES(x, y, z, t, vs, conf) the retarded time for a point source % uniformly moving along a line. Hereby, vs defines its direction % and its velocity. The position of the sound source at t=0 is the coordinate % origin. @@ -51,10 +53,10 @@ %***************************************************************************** %% ===== Checking of input parameters ================================== -% nargmin = 4; -% nargmax = 4; +% nargmin = 6; +% nargmax = 6; % narginchk(nargmin,nargmax); -% isargmatrix(x,t,vs); +% isargmatrix(x,y,z,t,vs); % isargstruct(conf); @@ -69,15 +71,17 @@ M = v/c; % Mach number [n x 1] or [1 x 1] % vector between observer and point source d(t) = x - vs*t -% [n x m x 3] -d = bsxfun(@minus, reshape(x,[],1,3), bsxfun(@times,t, reshape(vs,[],1,3))); +% [n x m] or [1 x m] +dx = bsxfun(@minus, x, bsxfun(@times,t, vs(:,1))); +dy = bsxfun(@minus, y, bsxfun(@times,t, vs(:,2))); +dz = bsxfun(@minus, z, bsxfun(@times,t, vs(:,3))); % instant distance between sound source and observer position -r = vector_norm(d, 3); % [n x m] -% component of x - v*t parallel to the direction of movement: (x-v*t)^T * nvs +r = sqrt( dx.^2 + dy.^2 + dz.^2 ); % [n x m] +% component of d(t) parallel to the direction of movement: d(t)^T * nvs % [n x m] -xpara = sum( bsxfun(@times, d, reshape(nvs,[],1,3)), 3); +xpara = nvs(:,1).*dx + nvs(:,2).*dy + nvs(:,3).*dz; % auxiliary distance -Delta = sqrt( M.^2.*xpara.^2 + (1-M.^2).*(r.^2) ); % [n x 1] +Delta = sqrt( M.^2.*xpara.^2 + (1-M.^2).*r.^2 ); % [n x 1] % retarded distance R = (M.*xpara + Delta)./(1-M.^2); % retarded time diff --git a/SFS_monochromatic/greens_function_mono.m b/SFS_monochromatic/greens_function_mono.m index c9eba476..115049be 100644 --- a/SFS_monochromatic/greens_function_mono.m +++ b/SFS_monochromatic/greens_function_mono.m @@ -111,25 +111,21 @@ % vs = v * ns % velocity vector % M = v/c % Mach number % xparallel = (x-xs)*ns % component in direction of movement - % R_1 = sqrt( M^2*xorth.^2 + (1-M^2)*|x-xs|.^2 ) - % R = (M*xorth + R1)./(1-M*2); + % R_1 = sqrt( M^2*xparallel.^2 + (1-M^2)*|x-xs|.^2 ) + % R = (M*xparallel + R1)./(1-M*2); % % see: Ahrens (2012), eq.(5.60) - % - v = norm(xs(4:6)); % velocity of sound source - nxs = xs(4:6)./v; % direction of movement - xs = xs(1:3); % - M = v/c; % - % shift coordinates - x = x-xs(1); - y = y-xs(2); - z = z-xs(3); - % component of x in direciton of movement: scalar = x*nxs - xparallel = nxs(1).*x + nxs(2).*y + nxs(3).*z; - - R1 = sqrt( M^2.*xparallel.^2 + (1-M^2)*(x.^2 + y.^2 + z.^2) ); - Rplus = (M*xparallel + R1)./(1-M.^2); - G = 1/(4*pi) * exp(-1i*omega/c.*Rplus)./R1; + % + [~,x1] = xyz_axes_selection(x,y,z); % get first non-singleton axis + % shift and vectorize coordinates + x = x(:)-xs(1); + y = y(:)-xs(2); + z = z(:)-xs(3); + % retarded time + [tau, ~, R1] = retarded_time(x,y,z,0,xs(4:6),conf); + % + G = zeros(size(x1)); + G(:) = 1/(4*pi) * exp(-1i*omega.*tau)./R1; elseif strcmp('ls',src) % Source model for a line source: 2D Green's function. %