%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% order俊 蝶弗 LPC spectrum %
% by whitjch 06/02/02       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

news = fopen('ieaou.pcm' ,'rb');
x = fread(news , 'short')/2^15;
fclose(news);
x=x(1:160);                                       % 30 msec
LENGTH=length(x) ;                                % length of x[n]
fs=8000;                                          % sampling rate 
n=0:1/fs:(LENGTH-1)/fs;   

% original speech signal %
subplot(4,2,1),plot(n*1000, x),hold on
xlabel('Time[msec]');    ylabel('Amplitude');
title('Input Signal')

% WIN_LEN=160;                                    % window length 160 (20msec)
N=256;

% short-time spectrum %
f=0:fs/N:fs/2-1;                                  % half frequency
X=abs(fft(x,N));
X_LOG=20*log10(X);                                % LOG (db)
subplot(4,2,2),plot(f/1000,X_LOG(1:N/2)),hold on
xlabel('frequency[kHz]'),ylabel('LOG(db)')
title('Short-Time Spectrum')

% LPC spectrum %
order=[4 8 12 16 20];                                % order
[a,g]=lpc(x,order(1));                               % predictor coefficients
EST_X=freqz(1,a,N/2);
EST_X_LOG=20*log10(abs(EST_X));
subplot(4,2,3),plot(f/1000,EST_X_LOG(1:N/2))
title('LPC spectrum (p=4)')

[a,g]=lpc(x,order(2));                               % predictor coefficients
EST_X=freqz(1,a,N/2);
EST_X_LOG=20*log10(abs(EST_X));
subplot(4,2,4),plot(f/1000,EST_X_LOG(1:N/2))
title('LPC spectrum (p=8)')

[a,g]=lpc(x,order(3));                               % predictor coefficients
EST_X=freqz(1,a,N/2);
EST_X_LOG=20*log10(abs(EST_X));
subplot(4,2,5),plot(f/1000,EST_X_LOG(1:N/2))
title('LPC spectrum (p=12)')

[a,g]=lpc(x,order(4));                               % predictor coefficients
EST_X=freqz(1,a,N/2);
EST_X_LOG=20*log10(abs(EST_X));
subplot(4,2,6),plot(f/1000,EST_X_LOG(1:N/2))
title('LPC spectrum (p=16)')

[a,g]=lpc(x,order(5));                               % predictor coefficients
EST_X=freqz(1,a,N/2);
EST_X_LOG=20*log10(abs(EST_X));
subplot(4,2,7),plot(f/1000,EST_X_LOG(1:N/2))
title('LPC spectrum (p=20)')