function oy = schroeder (x,tn,tg,td,cg,cd,ag,ad,k)

%     Schroeder reverb function 
%     The structure is:  y = schroeder(x,tn,tg,td,cg,cd,ag,ad,k)
%
%
%     where x = the input signal
%     Input signal array
%
%     tn = the number of taps of the tapped delay line
%     Number of early reflections e.g. 4
%
%     tg = a vector which contains the gain of each tap of the tapped delay line
%     Early reflection coefficients (linear not dB) e.g. [0.7, 0.5, 0.8, 0.6] for -3.1dB,-6dB -1.9dB & -4.4dB  
%
%     td = a vector which contains the delay length of each tap of the tapped delay line
%     Early reflection delays (in samples not time) [59, 100, 195, 225] for 3.7, 6.3, 12.2 14.1 ms reflections @ FS=16k
%
%     cg = a vector of length 4 which contains the gain of each of the comb filters (should be less than 1 for stability)
%     Comb filter (recursive echo) coefficients (linear not dB) e.g. [0.7, 0.6, 0.8, 0.7] for -3.1dB,-4.4dB -1.9dB & -3.1dB    
%
%     cd = a vector of length 4 which contains the delay of each of the comb filters
%     Comb filter (recursive echo) delays (in samples not time) e.g. [400, 531, 640, 734] for 25 33.2 40 45.9 ms reflections @ FS=16k
%
%     ag = the gain of the allpass filters (this should be less than 1 for stability)
%     All Pass filter (more recursive echoes) coefficient (linear not dB) e.g. 0.7 for -3.1 dB reflections
%
%     ad = a vector of length 2 which contains the delay of each of the allpass filters
%     All Pass filter (more recursive echoes) delays (in samples not time) e.g. [520, 495] for 32.5 & 30.9 ms reflections
%
%     k = the gain factor of the direct signal added to the cascaded allpass filters
%     Direct signal level (linear not dB) e.g. 0.5 for -6dB
%
%     y = the output signal
%     Output signal array to be played back
%
%     Example function call. 
%     y = schroeder(x,4,[0.7,0.5,0.8,0.6],[59, 100, 195,225],[0.7,0.6,0.8,0.7],[400, 531, 640, 734],0.7,[520,495],0.5);
%
%
% Gautham J. Mysore - gauthamjm@yahoo.com / Edited by Luis Miranda
%


%send the input signal througth the tapped delay line
ty = tapped_delay_line(x,tn,tg,td);


%increase the delay lengths so that the total delays of the allpass filters are greater than that of the tapped delay line
ad = ad + td(tn);

% send the input to each of the 4 comb filters separately

%ncomb is the number of comb filters. Later, User can define the number of
%the filters by using GUI or data file
ncomb = 4;

for i = 1:ncomb
    outcomb{i} = fbcomb(x,cg(i),cd(i));
    len_outcomb(i) = length(outcomb{i}(:,1));
    
    if i == ncomb
        max_outcomb = max(len_outcomb);
    end
end

for j = 1:ncomb
    outcomb{j} = [(outcomb{j}(:,1));zeros(max_outcomb - length(outcomb{j}(:,1)),1)];

    if j == 1
        %initialize apinput value
        apinput = zeros(length(outcomb{j}(:,1)),1);    
        apinput_temp = zeros(length(outcomb{j}(:,1)),1);
    end
    
    apinput_temp(:,1) = (outcomb{j}(:,1));
    
    apinput = apinput + apinput_temp;
    
    %sum the ouptut of the 4 comb filters
    if j == ncomb
        apinput = (apinput / 4) * ag;
    end
end

% send the output of the comb filters to the allpass filters
y = allpass(apinput,1,ad(1));
y = allpass(y,ag,ad(2));
y = allpass(y,ag,ad(3));

%max value between tappled delay output, output from all pass filter, and
%input signal x.
max_out = max([length(ty), length(y), length(x)]);

%find the total output
ty = [ty; zeros(max_out-length(ty),2)];
y = [y; zeros(max_out-length(y),2)];
x = [x; zeros(max_out-length(x),2)];

%oy = ty + y + k.*x;
oy(:,1) = ty(:,1) + y + k.*x(:,1);
oy(:,2) = ty(:,2) + y + k.*x(:,2);


% normalize the output signal
oy(:,1) = 0.99.*(oy(:,1) ./ max(abs(oy(:,1))));
oy(:,2) = 0.99.*(oy(:,2) ./ max(abs(oy(:,2))));

oy = [oy(:,1) oy(:,2)]
