function [aCoeff,resid,pitch,G,parcor,stream] = proclpc(data,sr,L,fr,fs,preemp)
% USAGE: [aCoeff,resid,pitch,G,parcor,stream] = proclpc(data,sr,L,fr,fs,preemp)
%
% This function computes the LPC (linear-predictive coding) coefficients that
% describe a speech signal. The LPC coefficients are a short-time measure of
% the speech signal which describe the signal as the output of an all-pole
% filter. This all-pole filter provides a good description of the speech
% articulators; thus LPC analysis is often used in speech recognition and 
% speech coding systems. The LPC parameters are recalculated, by default in
% this implementation, every 20ms.
%
% The results of LPC analysis are a new representation of the signal
% s(n) = G e(n) - sum from 1 to L a(i)s(n-i)
% where s(n) is the original data. a(i) and e(n) are the outputs of the LPC 
% analysis with a(i) representing the LPC model. The e(n) term represents 
% either the speech source's excitation, or the residual: the details of the 
% signal that are not captured by the LPC coefficients. The G factor is a
% gain term.
%
% LPC analysis is performed on a monaural sound vector (data) which has been
% sampled at a sampling rate of "sr". The following optional parameters modify
% the behaviour of this algorithm.
% L - The order of the analysis. There are L+1 LPC coefficients in the output
% array aCoeff for each frame of data. L defaults to 13.
% fr - Frame time increment, in ms. The LPC analysis is done starting every
% fr ms in time. Defaults to 20ms (50 LPC vectors a second)
% fs - Frame size in ms. The LPC analysis is done by windowing the speech
% data with a rectangular window that is fs ms long. Defaults to 30ms
% preemp - This variable is the epsilon in a digital one-zero filter which 
% serves to preemphasize the speech signal and compensate for the 6dB
% per octave rolloff in the radiation function. Defaults to .9378.
%
% The output variables from this function are
% aCoeff - The LPC analysis results, a(i). One column of L numbers for each
% frame of data
% resid - The LPC residual, e(n). One column of sr*fs samples representing
% the excitation or residual of the LPC filter.
% pitch - A frame-by-frame estimate of the pitch of the signal, calculated
% by finding the peak in the residual's autocorrelation for each frame.
% G - The LPC gain for each frame.
% parcor - The parcor coefficients. The parcor coefficients give the ratio
% between adjacent sections in a tubular model of the speech 
% articulators. There are L parcor coefficients for each frame of 
% speech.
% stream - The LPC analysis' residual or excitation signal as one long vector.
% Overlapping frames of the resid output combined into a new one-
% dimensional signal and post-filtered.
%
% The synlpc routine inverts this transform and returns the original speech
% signal.
%
% This code was graciously provided by:
% Delores Etter (University of Colorado, Boulder) and 
% Professor Geoffrey Orsak (Southern Methodist University) 
% It was first published in
% Orsak, G.C. et al. "Collaborative SP education using the Internet and
% MATLAB" IEEE SIGNAL PROCESSING MAGAZINE Nov. 1995. vol.12, no.6, pp.
% 23-32.
% Modified and debugging plots added by Kate Nguyen and Malcolm Slaney
% A more complete set of routines for LPC analysis can be found at
% http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
% (c) 1998 Interval Research Corporation 
if (nargin<3), L = 13; end
if (nargin<4), fr = 20; end
if (nargin<5), fs = 30; end
if (nargin<6), preemp = .9378; end
[row col] = size(data);
if col==1 data=data'; end
nframe = 0; 
msfr = round(sr/1000*fr); % Convert ms to samples
msfs = round(sr/1000*fs); % Convert ms to samples
duration = length(data);
speech = filter([1 -preemp], 1, data)'; % Preemphasize speech
msoverlap = msfs - msfr;
ramp = [0:1/(msoverlap-1):1]'; % Compute part of window

for frameIndex=1:msfr:duration-msfs+1 % frame rate=20ms
frameData = speech(frameIndex:(frameIndex+msfs-1)); % frame size=30ms
nframe = nframe+1;
autoCor = xcorr(frameData); % Compute the cross correlation
autoCorVec = autoCor(msfs+[0:L]);
% Levinson's method
err(1) = autoCorVec(1);
k(1) = 0;
A = [];
for index=1:L
numerator = [1 A.']*autoCorVec(index+1:-1:2);
denominator = -1*err(index);
k(index) = numerator/denominator; % PARCOR coeffs
A = [A+k(index)*flipud(A); k(index)]; 
err(index+1) = (1-k(index)^2)*err(index);
end
aCoeff(:,nframe) = [1; A];
parcor(:,nframe) = k';
% Calculate the filter 
% response
% by evaluating the 
% z-transform
if 0
gain=0;
cft=0:(1/255):1;
for index=1:L
gain = gain + aCoeff(index,nframe)*exp(-i*2*pi*cft).^index;
end
gain = abs(1./gain);
spec(:,nframe) = 20*log10(gain(1:128))';
plot(20*log10(gain));
title(nframe);
drawnow;
end
% Calculate the filter response
% from the filter's impulse
% response (to check above).
if 0
impulseResponse = filter(1, aCoeff(:,nframe), [1 zeros(1,255)]);
freqResp = 20*log10(abs(fft(impulseResponse)));
plot(freqResp);
end
errSig = filter([1 A'],1,frameData); % find excitation noise

G(nframe) = sqrt(err(L+1)); % gain
autoCorErr = xcorr(errSig); % calculate pitch & voicing information
[B,I] = sort(autoCorErr);
num = length(I);
if B(num-1) > .01*B(num)
pitch(nframe) = abs(I(num) - I(num-1));
else
pitch(nframe) = 0;
end
% calculate additional info to improve the compressed sound quality
resid(:,nframe) = errSig/G(nframe);
if(frameIndex==1) % add residual frames using a trapezoidal window
stream = resid(1:msfr,nframe);
else
stream = [stream; 
overlap+resid(1:msoverlap,nframe).*ramp; 
resid(msoverlap+1:msfr,nframe)];
end
if(frameIndex+msfr+msfs-1 > duration)
stream = [stream; resid(msfr+1:msfs,nframe)];
else
overlap = resid(msfr+1:msfs,nframe).*flipud(ramp); 
end 
end
stream = filter(1, [1 -preemp], stream)';