function [ y ] = vocoder_reverb(carrier, modul, chan, numband, overlap, fs, delays, rt60, wavout)
% The channel vocoder requires two signals: a carrier and a modulator.
% The vocoder takes the two signals, and using their spectral information, 
% creates a third signal. The Channel Vocoder modulates the carrier signal 
% with the modulation signal. This signal is then routed through a reverb unit. 
% The addition of the reverb not only adds to the creative possibilities 
% of the function, but allows the user to reconstruct the synthesised sound
% in a pre determined "space".
% The vocoder was initially designed by William A. Sethares; 
% the reverb was initially designed by Riitta Väänänen
%
% vocoder_reverb(carrier, modul, chan, numband, overlap, fs, delays, rt60,
% wavout);
% This function has many parameters which include: 
% (1) carrier - specifying the carrier signal (import .wav file - usually a
% square/sawtooth wave)
% (2) modul - specifying the modulator signal (import .wav file - usually a
% voice/instrument)
% (3) chan - number of channels (determines the length of the FFT) (e.g. 512)
% (4) number of bandwidths the signal will be broken into (e.g. 32)
% (5) the overlap of the FFT (a ratio between 0 and 1). (e.g 0.25)
% (6) the sampling frequency (e.g. 44100)
% (7) the number of samples between each delay for the reverberation (e.g. 2191 2971 3253  3307)
% (8) rt60 - time it takes reverb to decay 60dB in seconds (e.g. 1.5)
% (9) the final output wav files name. (must be a '.wav')
 
if numband>chan, error('# bands must be < # channels'), end
[rc, cc]   = size(carrier); if cc>rc, carrier = carrier'; end % Determines the size of the carrier signal
[rm, cm]   = size(modul);   if cm>rm, modul = modul'; end % Determines the size of the modulator signal
st         = min(rc,cc);                         % Are the signals stereo or mono?
if st~= min(rm,cm), error('carrier and modulator must have same number of tracks'); end
len        = min(length(carrier),length(modul)); % find shortest length of tracks
carrier    = carrier(1:len,1:st);                % shorten carrier if needed
modul      = modul(1:len,1:st);                  % shorten modulator if needed
L          = 2*chan;                             % window length/FFT length
w          = hanning(L); if st==2, w=[w w]; end  % window/ stereo window
bands      = 1:round(chan/numband):chan;         % indices for frequency bands     
bands(end) = chan;
voc_out          = zeros(len,st);                      % output vector
ii         = 0;
while ii*L*overlap+L <= len
    ind    = round([1+ii*L*overlap:ii*L*overlap+L]);
    FFTmod = fft( modul(ind,:) .* w );    % window & take FFT of modulator
    FFTcar = fft( carrier(ind,:) .* w );  % window & take FFT of carrier  
    syn    = zeros(chan,st);              % place for synthesized output
    for jj = 1:numband-1                  % for each frequency band
        b        = [bands(jj):bands(jj+1)-1]; % current band
        syn(b,:) = FFTcar(b,:)*diag(mean(abs(FFTmod(b,:))));
    end                                   % take product of spectra
    midval   = FFTmod(1+L/2,:).*FFTcar(1+L/2,:); % midpoint is special
    synfull  = [syn; midval; flipud( conj( syn(2:end,:) ) );]; % + and - frequencies
    timsig   = real( ifft(synfull) );     % invert back to time
    voc_out(ind,:) = voc_out(ind,:) + timsig;         % add back into time waveform   
    ii       = ii+1;
end
voc_out = 0.8*voc_out/max(max(abs(voc_out)));               % normalize output
wavwrite(voc_out,fs,'vocoded_output.wav');

%%%%%%%%%% REVERB CODE START %%%%%%%%%%%%
% allpass filter delays:
 apdelays = [441 713];
 apgains = [0.5 0.5];

% Feedback gain:
Kp = zeros(1,length(delays));
Kp = 10.^((-3.0.*delays)./(rt60*fs))

% Lowpass filter coefficients:
Bp = zeros(1,length(delays));
alpha = 0.25; % the ratio between the reverb times at nyquist and zero frequency
Bp = ones(1,length(delays)) - 2./(1 + Kp.^(1 - 1/alpha))

beta = Kp.*(1-Bp)

% lowpass filter (first order allpole) delays
lp_d = zeros(1,length(delays));

%y = zeros(length(delays),length(x(1,:)));
y = zeros(1,length(voc_out));

% initialization of delay lines
delaylines = zeros(length(delays),max(delays));

% allpass delay lines:
apdelaylines = zeros(2,max(apdelays));

% delay line pointers
dl_p = ones(1,length(delays));

ap_dlp1 = 1;
ap_dlp2 = 1
% temporary variables
temp_out = zeros(length(delays),1)


%length(x(1,:))

ap_out1 = 0;
ap_out2 = 0;

for i=1:length(voc_out),
   %i
   
   for j = 1:length(delays),
     y(i) = y(i) + temp_out(j); % sum the delay line outputs
   
   	% input of the delay line (input x + output fed back to the delay line)
   	%delaylines(j,dl_p(j)) = x(i) + delaylines(rem(dl_p(j),delays(j)) + 1) * beta(j);
 		delaylines(j,dl_p(j)) = voc_out(i) + temp_out(j);
   
		% delay line pointer update 
 		dl_p(j) = mod((dl_p(j)+1),delays(j)) + 1;
		
       % temporary output of each delay line;
       % This is a lowpass filtered and attenuated delay line output
      temp_out(j) = (delaylines(j,dl_p(j)) + Bp(j)*lp_d(j))*beta(j);
      
      % store the delay line output to the lowpass filter delay 
      % for use in the next round
      lp_d(j)=delaylines(j,dl_p(j));

   end
   
   	
      % Allpass filter 1:
   	apdelaylines(1,ap_dlp1) = y(i) + ap_out1 * apgains(1);	   
      
      % delay pointer update
      ap_dlp1 = mod(ap_dlp1+1,apdelays(1)) + 1;
      
      % output computation
      ap_out1 =  apdelaylines(1,ap_dlp1) - apgains(1) * y(i);
      
      %y(i) = ap_out1 + x(i);
      
      
      % Allpass filter 2:
      
      apdelaylines(2,ap_dlp2) = y(i) + ap_out2 * apgains(2);	   
      
      % delay pointer update
      ap_dlp2 = mod(ap_dlp2+1,apdelays(2)) + 1;
      
      % output computation
      ap_out2 =  apdelaylines(2,ap_dlp2) - apgains(2) * y(i);
      
      y(i) = ap_out2 + voc_out(i);
   
wavwrite(y,fs,wavout);
end  
