% MATLAB/Octave script to simulate the sound field generated by a (virtual)
% endfire loudspeaker array
%
%
% Background and results have been published in the paper:
%
% Spors, Sascha and Wierstorf, Hagen, A Virtual Enfire Louspeaker Array for
% the Generation of Sound Beams, AIA-DAGA 2013, March 2013, Meran, Italy. 
%
% PLEASE refer to this paper when using the script for your research.
%
%
% The script is based on the Sound Field Synthesis toolbox which you can
% obtain here: http://dev.qu.tu-berlin.de/projects/sfs-toolbox
%
% Sascha Spors, Sascha.Spors@uni-rostock.de
% 14.8.2013


%*****************************************************************************
% Copyright (c) 2013      Institut fuer Nachrichtentechnik                   *
%                         Universitaet Rostock                               *
%                         Richard-Wagner-Strasse 31, 18119 Rostock           *
%                                                                            *
%                                                                            *
% This 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 script 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.                       *
%*****************************************************************************

clear all
addpath('/Users/spors/Documents/src/SFS');      % path to SFS toolbox

%% --- set parameters for simulation ---
virtual_sources = 1;        % virtual (focused) sources for endfire array
delayandsum = 0;            % delay-and-sum beamfilter design

c = 343;
f = 1000;                     % frequency of simulated sound field

fh=1715*4;                  % upper frequency limit of endfire array design
N=32;                       % number of endfire sources

dy0=0.3;                    % distance of first endfire source from x-axis
dx0=0;                      % distance of first endfire source from y-axis
theta=90/180*pi;            % angle of endfire array

beta=1e-6;                  % regularization factor
thetad=0;                   % steering angle used for design of endfire beamforming

L = 10;                     % length of loudspeaker array
sdx0 = 0.10;                % distance of secondary sources



%% configuration for SFS toolbox
SFS_start;
conf = SFS_config;  


conf.secondary_sources.geometry = 'linear';
conf.secondary_sources.size = L;   % array length for WFS
conf.secondary_sources.number = conf.secondary_sources.size / sdx0;

conf.driving_functions = 'spors2009eq6';

conf.xref = [0 3 0];
conf.useplot = 0;
conf.usetapwin = true;
conf.tapwinlen = 0.6;

conf.plot.resolution = 150;
Xlim = 2.5;                 % spatial limits of simulated sound field
Ylim = 5;


%% --- variables for script  ---
vdx0 = c/(2*fh);                                    % distance of virtual endfire sources
vx0=[zeros(1,N)' dy0+(1:N)'.*vdx0 zeros(1,N)']';    % positions of virtual sources
vxd=[0 1 0];                                        % direction of focused sources



%% --- design of beamfilter ---

if(delayandsum)
    % delay-and-sum design
    dt = vx0(2,:)/c;
    F = 1/N * exp(-1i*2*pi*f*dt);
else
    % optimal beamfilter [Boone et al., JAES 57(5), 2009]
    W = exp(1i*2*pi*f/c*vx0(2,:)' .* cos(thetad));
    Szz=zeros(N,N);
    for n=1:N
        for m=1:N
            Szz(n,m) = sinc(2*f/c*(vx0(2,n)-vx0(2,m)));
        end
    end
    F = W' * inv(Szz + beta.*eye(N)) / (W'*inv(Szz + beta.*eye(N))*W);
end

%% --- (rotate and shift) positions of virtual sources ---
for n=1:N
   vx0(:,n) = vdx0*n * [cos(theta) sin(theta) 0] + [dx0 dy0 0];
end


%% --- compute synthesized sound field ---
for n=1:N
    
    if(virtual_sources)
        
        [P,x,y,z]=sound_field_mono_wfs([-Xlim,Xlim],[-0.15,Ylim],0,[vx0(:,n)' vxd],'fs',f,conf);
        x0 = secondary_source_positions(conf);
        
        % normalize level of focused source
        idy = find(y>=vx0(2,n),1,'first');
        idx = find(x>=vx0(1,n),1,'first');
        level(n)=abs(P(idy,idx));
        P = P ./ level(n);
        
    else
        
        [P,x,y,z] = sound_field_mono_point_source([-Xlim,Xlim],[-0.15,Ylim],0,vx0(:,n)',f,conf);
        x0=[];
    
    end
        
    % apply beamfilters
    P = P .* ( ones(size(P)) * F(n));
    
    % sum up contributions from individual endfire sources
    if(n==1)
        P2 = P;
    else
        P2 = P2 + P;
    end
end
P = P2;



%% --- show results ---
conf.xref = 1.5*[cos(theta) sin(theta)] + [dx0 0];
P = norm_sound_field_at_xref(P,x,y,0,conf);

% synthesized sound field
conf.plot.usedb = 0;
conf.usenormalisation = false;
plot_sound_field(P,x,y,z,x0,conf);
axis([-Xlim Xlim -0.15 Ylim]);
title('Syntehsized sound field');
colormap jet;

% level of synthesized sound field
conf.plot.usedb = 1;
conf.usenormalisation = false;
plot_sound_field(P,x,y,z,x0,conf);
hold on
[c,h]=contour(x,y,20*log10(abs(P)),[-20,-20],'LineColor',[0 0 0],'LineWidth',1);
clabel(c,h);
hold off
axis([-Xlim Xlim -0.15 Ylim]);
title('Level of synthesized sound field (dB)');

colormap jet;