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