%  This scripts reads in a wave file of a neutral voiced phoneme, high pass filters
%  to detrend the data and remove room noise, computes linear predictive
%  coefficients (LPCs) for the segment and plots the prediction error and
%  reconstructed signal with FIR and IIR implementations of the filter.
%  The pitch is estimated from the error signal and a simple synthesised
%  signal is created to imitate real speech at a different pitch.
%  A pole-zero diagram is create for the IIR reconstruction of the signal
%  the formants are identified and the PSD and Spectrum from the LPCs is
%  plotted.
%
%    Written by Kevin D. Donohue (donohue=@engr.uky.edu) July 2002

%  
winlens = 50;  %PSD window length in milliseconds
[y,fs] = wavread('aaa3.wav');  %  Read in wavefile
winlen = winlens*fs/1000;
[cb,ca] = butter(5,2*100/fs,'high');  %  Filter to remove LF recording noise
yf = filtfilt(cb,ca,y);
[a,er] = lpc(yf,10);  %  Compute LPC coefficent with model order 10

predy = filter(a,1,yf);  %  Compute prediction error with all zero filter

kd=1;   %  Starting figure number
figure(kd) ; plot(predy); title('Prediction error'); xlabel('Samples'); ylabel('Amplitude')

recon = filter(1,a,predy);  %  Compute reconstructed signal from error and all-pole filter
figure(kd+1)

%  Plot reconstructed signal
plot(recon,'b')

hold on
%  Plot with orginal delayed by a unit so it does not entirely overlap
%  the perfectly reconstructed signal 
plot(yf(2:end),'r')
hold off
xlabel('Samples'); ylabel('Amplitude')
title('Reconstructed Signal (blue) and Original (red)')

%  By examining a the error sequence,
%  generate a simple impulse sequence to simulate its period
%  (about 103 sample period)
g = [];
for k=1:150
    g = [g, 1, zeros(1,103)];
end
%  Run simulated error sequence through all pole filter
sim = filter(1,a,g);
%soundsc([(sim')/std(sim); zeros(fix(fs)*1,1); yf/std(yf)],fs)

%  Plot pole zero diagram
figure(kd+2)
r  = (roots(a))
w = [0:.001:2*pi];
plot(real(r),imag(r),'xr',real(exp(j*w)),imag(exp(j*w)),'b')
title('Pole diagram of vocal tract filter')
xlabel('Real'); ylabel('Imaginary')

%  Find resonant frequencies corresponding to poles 
froots = (fs/2)*angle(r)/pi;
nf = find(froots > 0 & froots < fs/2);  %  Find those corresponding to complex conjugate poles
figure(kd+3)
%  Examine average specturm with formant frequencies
[pd,f] = pwelch(yf,hamming(winlen),fix(winlen/2),2*winlen,fs);  
dbspec = 20*log10(pd);
mxp = max(dbspec);  %  Find max and min points for graphing verticle lines
mnp = min(dbspec);
plot(f,dbspec,'b')  %  Plot PSD
hold
% Over lines on plot where formant frequencies were estimated from LPCs
for k=1:length(nf)
plot([froots(nf(k)), froots(nf(k))], [mnp(1), mxp(1)], 'k--')
end
hold off
title('PSD plot with formant frequencies (Black broken lines)')
xlabel('Hertz')
ylabel('dB')
%  Get spectrum from the AR (LPC) parameters
[hz,fz] = freqz(1, a, 1024, fs);
figure(kd+4)
plot(fz,abs(hz))
title('Spectrum Generated by LPCs')
xlabel('Hertz')
ylabel('Amplitude')