diff --git a/SFS_general/direction_vector.m b/SFS_general/direction_vector.m index e2c04e8d..6405d36d 100644 --- a/SFS_general/direction_vector.m +++ b/SFS_general/direction_vector.m @@ -1,7 +1,7 @@ function directions = direction_vector(x1,x2) %DIRECTION_VECTOR return unit vector(s) pointing from x1 to x2 % -% Usage: n = direction_vector(x1,x2) +% Usage: n = direction_vector(x1,[x2]) % % Input parameters: % x1 - starting point(s) [1xn] or [mxn] @@ -11,7 +11,8 @@ % 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. +% 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 @@ -46,15 +47,20 @@ %% ===== Checking of input parameters =================================== -nargmin = 2; +nargmin = 1; nargmax = 2; narginchk(nargmin,nargmax); -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)); +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 diff --git a/SFS_general/retarded_time.m b/SFS_general/retarded_time.m new file mode 100644 index 00000000..494203b5 --- /dev/null +++ b/SFS_general/retarded_time.m @@ -0,0 +1,90 @@ +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, y, z, t, vs, conf) +% +% Input parameters: +% 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) +% +% 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, 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. +% +% 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 = 6; +% nargmax = 6; +% narginchk(nargmin,nargmax); +% isargmatrix(x,y,z,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] 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 = 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 = 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] +% retarded distance +R = (M.*xpara + Delta)./(1-M.^2); +% retarded time +tau = R./c; + +end diff --git a/SFS_helper/isargsecondarysource.m b/SFS_helper/isargsecondarysource.m index 78f7e0f4..5d4706e8 100644 --- a/SFS_helper/isargsecondarysource.m +++ b/SFS_helper/isargsecondarysource.m @@ -54,9 +54,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 917a7783..115049be 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 @@ -99,6 +100,32 @@ % 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*xparallel.^2 + (1-M^2)*|x-xs|.^2 ) + % R = (M*xparallel + R1)./(1-M*2); + % + % see: Ahrens (2012), eq.(5.60) + % + [~,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. %