processReverb(data, stAlgo)

function [stAlgo, status, data] = processReverb(data, stAlgo)

%Version 3.2
%Copyright  2003 Stephen G. McGovern

status = 0;
n = stAlgo.n;
src = stAlgo.src;
mic = stAlgo.mic;
rm = stAlgo.rm;
r = stAlgo.r;
fs = stAlgo.fs

nn=[-n:1:n];                          % Index for the sequence
rms=nn+0.5-0.5*(-1).^nn;              % Part of equations 2,3,& 4
srcs=(-1).^(nn);                      % part of equations 2,3,& 4
xi=[srcs*src(1)+rms*rm(1)-mic(1)];    % Equation 2
yj=[srcs*src(2)+rms*rm(2)-mic(2)];    % Equation 3
zk=[srcs*src(3)+rms*rm(3)-mic(3)];    % Equation 4

[i,j,k]=meshgrid(xi,yj,zk);           % convert vectors to 3D matrices
d=sqrt(i.^2+j.^2+k.^2);               % Equation 5
time=round(fs*d/343)+1;               % Similar to Equation 6

[e,f,g]=meshgrid(nn, nn, nn);         % convert vectors to 3D matrices
c=r.^(abs(e)+abs(f)+abs(g));          % Equation 9
e=c./d;                               % Equivalent to Equation 10

h=full(sparse(time(:),1,e(:)));       % Equivalent to equation 11

Ly=length(data)+length(h)-1;  % Begin of convolution data*h
Ly2=pow2(nextpow2(Ly));    % Find smallest power of 2 that is > Ly
X=fft(data, Ly2);		   % Fast Fourier transform
H=fft(h, Ly2);	           % Fast Fourier transform
Y=X.*H;        	           %
y=real(ifft(Y, Ly2));      % Inverse fast Fourier transform
y=y(1:1:Ly);               % Take just the first N elements
y=y/max(abs(y));
ystart = find(y>0.0001);
y=y([ystart(1):end],:);% Normalize the output
data = y;