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