% LPCP: function to compute & plot amplitude and LPC spectrum
% Usage: [a,s,f]=lpcp(wave,rate,window,ncoef);
% wave: input waveform
% rate: sample rate in Hz (default 10000 Hz)
% window options: 'hann', 'hamm', 'kais', or 'rect' (default=hamming)
% ncoef: number of LP coefficients (default 12)
% [a,s,f]: log magnitude (peak-normalized dB), lpc spectrum, frequency in Hz

function [amp,s,freq]=lpcp(x,rate,window,ncoef);

if ~exist('rate','var'), 
    rate=10000; 
end;
if ~exist('window','var'),
    window='hamm'; 
end;
if ~exist('ncoef','var'), 
    ncoef=12; 
end;
x=x(:);  % make x a column vector
n=length(x);

pre = [1 -0.9]; % pre-emphasis filter
x = filter(pre,1,x);

% apply window function
window=lower(window);
if window=='rect',
    x=x.*ones(n,1);
elseif window=='hamm',
    x=x.*hamming(n);
elseif window=='hann',
    x=x.*hanning(n);
elseif window=='kais',
    x=x.*kwind(ones(n,1),n/2-1);
end;
% compute log magnitude spectrum
m=fft(x,n);
no2=round(n/2);
a=20*log10( abs ( m ) / n);  
amp=a(1:no2);
amp=amp-max(amp);
freq=linspace(0,rate/2000,length(amp))';
plot(freq,amp);
hold on;
% [a,g]=lpc(x,ncoef); 
% Compute autocorrelation function 
% nfft=2.^nextpow2(2*(length(x))); 
nfft=2*length(x); 
X=fft(x,nfft);
Pxx=X.*conj(X);
pxx=ifft(Pxx);
rxx=fftshift(real(pxx));
rxx(1)=[];
% rxx=xcorr(x);
% remove rxx at negative lags
rxx(1:length(x)-1)=[]; 
% Setup the normal equations
R = toeplitz(rxx(1:ncoef));
a = inv(R)*rxx(2:(ncoef+1));
a = [1; -a];
s=fft(a,n);
s=-10*log10(real(s.*conj(s)));
s(round(length(s)/2)+1:end)=[];
s=s-max(s)+5; % for convenience- 5 dB offset
plot(freq,s,'r','LineWidth',2);
hold off;
xlabel('Frequency (kHz)');
ylabel('Amplitude (dB)');
axis([0 rate/2000 -90 10]);