processRIRverb(data, stAlgo)

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


%Covolution code by Stephen G. McGovern 2003

status = 0;

irdata = [];
h=stAlgo.ir(:,1);
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(:,1), 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
irdata = stAlgo.mix.*y;
irdata(:,1) = irdata+((1-stAlgo.mix)*[data(:,1);zeros(length(y)-length(data),1)]);

if size(stAlgo.ir,2)>1
    chn=1;
    if size(data,2)>1
    chn=2;
    end
    irdata = [irdata(:,1) zeros(length(irdata),1)];
    h=stAlgo.ir(:,2);
    Ly=length(data(:,chn))+length(h)-1;  % Begin of convolution data*h
    Ly2=pow2(nextpow2(Ly));    % Find smallest power of 2 that is > Ly
    X=fft(data(:,chn), 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));
    y=y([ystart(1):end],:);% Normalize the output
    irdata(:,2) = stAlgo.mix.*y+(1-stAlgo.mix).*[data(:,chn);zeros(length(y)-length(data),1)];
end

data = irdata;