%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This work is supplementary material for the book %
% %
% Jens Ahrens, Analytic Methods of Sound Field Synthesis, Springer-Verlag %
% Berlin Heidelberg, 2012, http://dx.doi.org/10.1007/978-3-642-25743-8 %
% %
% It has been downloaded from http://soundfieldsynthesis.org and is %
% licensed under a Creative Commons Attribution-NonCommercial-ShareAlike %
% 3.0 Unported License. Please cite the book appropriately if you use %
% these materials in your own work. %
% %
% (c) 2012 by Jens Ahrens %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
% This script shows how to simulate WFS of a moving source with a
% trajectory that is not parallel to the x-axis. The trick is to rotate the
% underlying spatial grid as well as the secondary source distribution, do
% the simulation, and then rotate the coordinate system. The result is a
% secondary source distribution that is parallel to the x-axis but the
% source trajectory is not. This example is not contained in the book.
% See http://www.soundfieldsynthesis.org/moving-sound-sources.
v = 240;
theta = 30 / 180 * pi; % rotation angle of source trajectory
f = 500;
t = 0;
c = 343;
omega = 2*pi*f;
M = v / c;
dx_0 = 0.1;
d_ref = 1;
% create spatial grid
X = linspace( -3, 3, 500 );
Y = linspace( -1, 5, 500 );
[ x, y ] = meshgrid( X, Y );
z = 0;
% rotate spatial grid in opposite direction compared to the source
% trajectory (Trick!)
x_rot = cos( -theta ) .* x - sin( -theta ) .* y;
y_rot = sin( -theta ) .* x + cos( -theta ) .* y;
% initialize s
s = zeros( size( x_rot ) );
x_0_prime = -5 : dx_0 : 5; % this moves along the secondary source distribution
y_0_prime = ones( size( x_0_prime ) );
% rotate secondary source distribution (the source still moves along x-axis)
x_0s = cos( -theta ) * x_0_prime - sin( -theta ) * y_0_prime;
y_0s = sin( -theta ) * x_0_prime + cos( -theta ) * y_0_prime;
z_0 = 0;
% loop over secondary sources (Eq. (5.69))
for index = 1 : length( x_0_prime )
x_0 = x_0s( index );
y_0 = y_0s( index );
disp( [ 'Processing secondary source at x_0 = ' num2str( x_0, '%0.1f' ) ' m.' ] );
r_0 = sqrt( ( x_rot - x_0 ).^2 + ( y_rot - y_0 ).^2 + ( z - z_0 ).^2 );
% Eq. (5.58) and (5.69)
Delta = sqrt( ( x_0 - v * ( t - r_0/c ) ).^2 + ( y_0.^2 + z_0.^2 ) .* ( 1 - M^2 ) );
% Eq. (5.57) and (5.69)
tau = ( M .* ( x_0 - v * ( t - r_0/c ) ) + Delta ) / ( c * ( 1 - M^2 ) );
% Eq. (5.65)
dx = ( ( x_0 - v*t ) ./ Delta.^2 + 1 / ( c * ( 1 - M^2 ) ) * ( M + ( x_0 - v*t ) ./ Delta ) * 1i*omega ) .* exp( 1i .* omega .* ( t - r_0/c - tau ) ) ./ Delta;
% Eq. (5.66)
dy = y_0 ./ Delta .* ( ( 1 - M^2 ) ./ Delta + 1/c * 1i*omega ) .* exp( 1i .* omega .* ( t - r_0/c - tau ) ) ./ Delta;
% Eq. (3.93)
d = sqrt( 2*pi*d_ref / ( 1i * omega/c ) ) .* 2 .* ( cos( -theta + pi/2 ) * dx + sin( -theta + pi/2 ) * dy );
% Eq. (5.69)
s = s + 1 / (4*pi) .* d ./ r_0;
end
% normalize
s = s ./ abs( s( end/2, end/2 ) );
% plot inherently rotates sound field by -theta (now, the source moves
% along a trajectory that is not parallel to the secondary source
% distribution
figure;
imagesc( X, Y, real( s ), [ -1.5 1.5 ] );
set( gca, 'YDir','normal' );
axis square;
colormap jet;
hold on;
% plot secondary sources
plot( ( -3 : dx_0 : 3 ), 1, 'kx' )
% plot source trajectory
plot( [ -cos( theta ) cos( theta ) ] * 5, [ -sin( theta ) sin( theta ) ] * 5, 'Color', [ 1 1 1 ], 'Linestyle', '--', 'Linewidth', 2 )
% plot position of virtual source
plot( 0, 0, 'k.' )
hold off;
xlabel( '$x$ (m)', 'Interpreter', 'latex' );
ylabel( '$y$ (m)', 'Interpreter', 'latex' );
set( gcf, 'Color', [ 1, 1, 1 ] );