www.pudn.com > LPC.rar > proclpc.m
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 &amt; 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)';