% Very simple speech vocoding routine

% Input is x, a column vector
%
% Parameters to be specified:
% 
%	fs - sampling frequency (Hz)
% 	noverlap - number of samples of overlap of successive windows
%	nfft - fft size for window
%
%
% Last modified 7/28/98 - rms


% Parameter values for Angelo Ripepi's experiment

fs = 11025;	% in Hz
f0 = 0; % monotone excitation frequency, in Hz
fnyq = fs/2;	% Nyquist frequency
nfft = 256;
noverlap = nfft/2;
window = hamming(256);
ncoeffs = 14;		% Number of LPC coefficients used


% reading in input
load welcome;		% loads in a variable called x from MAT file welcome.mat

% Windowing signal ... columns are frames

[xspec,f,t] = SPECPREP(x,nfft,fs,window,noverlap);


% Computing excitation function p, which is either pulses with frequency
% f0, or random noise (if f0 ==0)

if f0 == 0
	p = randn(1,length(x));
else
	N0 = round(fs/f0);	% discretized sample period
	p = [1:length(x)];
	p = mod(p,N0)==0;	% Pulse train with period N0	
end

[pspec,f,t] = SPECPREP(p,nfft,fs,window,noverlap);


% Computing LPC coefficients
% Note that input is split into even and odd frames, to compensate for
% 50 percent overlap of frames

odd = [1:2:size(xspec,2)];	% index vector for even columns
even = [2:2:size(xspec,2)];	% index vector for even columns

xodd = xspec(:,odd);
xeven = xspec(:,even);

yodd = zeros(size(xodd));
yeven = zeros(size(xeven));

[aodd,godd] = lpc(xodd,ncoeffs);
aodd = aodd';

[aeven,geven] = lpc(xeven,ncoeffs);
aeven = aeven';

b = 1;					% All pole filter => no numerator coeffs

ziodd = zeros(size(aodd,1)-1,size(aodd,2));	% Initial condition vector
zieven = zeros(size(aeven,1)-1,size(aeven,2));

% Computation of output on a frame-by-frame basis

% Different indexing needed for first frame to keep subscripts in bounds
	
	[yodd(:,1) ziodd(:,1)]  = filter(b,aodd(:,1),pspec(:,1),zeros(size(aodd,1)-1,1));
	yodd(:,1) = godd(1)*yodd(:,1);

	[yeven(:,1) zieven(:,1)]  = filter(b,aeven(:,1),pspec(:,1),zeros(size(aeven,1)-1,1));
	yeven(:,1) = geven(1)*yeven(:,1);

for l = 2:size(xodd,2),		% General filtering loop for odd frames
	[yodd(:,l) ziodd(:,l)]  = filter(b,aodd(:,l),pspec(:,l),ziodd(:,l-1));
	yodd(:,l) = godd(l)*yodd(:,l);
end
	
for l = 2:size(xeven,2),	% General filtering loop for even frames
	[yeven(:,l) zieven(:,l)]  = filter(b,aeven(:,l),pspec(:,l),zieven(:,l-1));
	yeven(:,l) = geven(l)*yeven(:,l);
end
	

% Now recombine even and odd frames
	
youtodd = reshape(yodd,size(yodd,1)*size(yodd,2),1);
youteven = reshape(yeven,size(yeven,1)*size(yeven,2),1);

if length(youtodd)==length(youteven)	% Is the number of frames even?
	youtodd = [youtodd' zeros(noverlap,1)']';	% Stuff zeros in front and in back
	youteven = [zeros(noverlap,1)' youteven']';

else								% If odd number of frames, only stuff evens
	youteven = [zeros(noverlap,1)' youteven' zeros(noverlap,1)']';

end
	

	
yout = youtodd+youteven;	
yout = yout/max(yout);
sound(yout,11025)