LPC vocoder-ucla-working-2002.rar_LPC_LPC vocoder_WORKING_predic

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)';