www.pudn.com > LPCC.rar > LPCC.m

function [parameter, yPreEmp]= wave2lpc(y, fs, FP) 
 
 a=nargin; 
if nargin<1,  selfdemo;  return; end 
if nargin<2; fs=8000; end 
if nargin<3, FP=mfccParamSet(fs); end 
 
y=double(y);	% Convert to double 
y=y-mean(y);	% Shift to zero mean, "the value is very samll" 
 
%% ====== Step 1: pre-emphasis.   
yPreEmp = filter([1, -0.96875], 1, y); 
 
%% ====== Step 2: frame blocking. 
framedY = buffer2(y, FP.frameSize, FP.overlap); 
 
 
%% ====== Step 3: Add hamming window. 
for k = 1:size(framedY, 2), 
	% ====== Step 3: hamming window. 
	Wframe(:,k)  = hamming(FP.frameSize).*framedY(:,k); 
 
     
%% ====== Step4: LPC coeffinient. 
    order=10; 
    % autocorrelation, Rx     
    Rx=LPC(Wframe(:,k), order); 
    Rx_max=Rx(1); 
     
    % 利用第回求解 LPC 係數 
    % 初始值設定 
    alpha(order+1, order+1)= 0; 
    Emin = Rx(1); 
    if Rx(1) ==0            Emin=1.0e-30;    end 
    
    for L= 1:order 
            sum= 0; 
            for i=1 : L -1 
                    sum= sum + alpha(i, L-1) * Rx(L-i); 
            end 
            k=(Rx(L)-sum)/Emin; 
            alpha(L, L)= k;           
  
            for i =1: L-1 
                    alpha(i,L)=alpha(i, L-1) - k*alpha(L-i, L-1); 
            end 
            Emin= Emin*(1-k*k); 
    end 
    a=alpha(order, :);            
end 
 
     
%% ====== Step 4: fast fourier transform. 
for i = 1:size(framedY, 2), 
	fftMag = abs(fft(Wframe)); 
	halfIndex = floor((FP.frameSize+1)/2); 
	fftMag = fftMag(1:halfIndex);   
	fftMag = interp1(1:halfIndex,fftMag,1:1/FP.alpha:halfIndex)';	% VTLN 
	fftMag = [fftMag;zeros(halfIndex-length(fftMag),1)]; 
     
%% ====== Step 5: triangular bandpass filter. 
%	tbfCoef = triBandFilter(fftMag, FP.tbfNum, fstart, fcenter, fstop); 
	tbfCoef = triBandFilter(fftMag, FP.tbfNum, filterBankParam); 
	% ====== Step 6: cosine transform. (Using DCT to get L order mel-scale-cepstrum parameters.) 
	mfcc = melCepstrum(FP.cepsNum, FP.tbfNum, tbfCoef); 
	parameter = [parameter mfcc']; 
end 
 
%% ====== Add energy 
if (FP.useEnergy==1) 
	energy = sum(framedY.^2)/FP.frameSize; 
	logEnergy = 10*log10(eps+energy); 
	parameter = [parameter; logEnergy]; 
end 
 
%% ====== Compute delta energy and delta cepstrum 
% with delta is better for telephone digit recognition HMM 
if (FP.useDelta>=1) 
	deltaWindow = 2; 
	paraDelta = deltaFunction(deltaWindow, parameter); 
	parameter = [parameter; paraDelta]; 
end 
if (FP.useDelta==2) 
	paraDeltaDelta = deltaFunction(deltaWindow, paraDelta); 
	parameter = [parameter; paraDeltaDelta]; 
end 
 
%% ====== Subfunction ====== 
% === Self demo 
function selfdemo 
waveFile='test.wav'; 
[y, fs]=wavread(waveFile); 
FP=mfccParamSet(fs); 
FP.useDelta=0; 
mfcc0= feval(mfilename, y, fs, FP); 
fprintf('No. of extracted frames = %d\n', size(mfcc0, 2)); 
subplot(3,1,1); surf(mfcc0); box on; axis tight; title(sprintf('MFCC of "%s"', waveFile)); 
FP.useDelta=1; 
mfcc1=feval(mfilename, y, fs, FP); 
subplot(3,1,2); surf(mfcc1); box on; axis tight; title(sprintf('MFCC of "%s"', waveFile)); 
FP.useDelta=2; 
mfcc2=feval(mfilename, y, fs, FP); 
subplot(3,1,3); surf(mfcc2); box on; axis tight; title(sprintf('MFCC of "%s"', waveFile)); 
 
%% === Triangular band-pass filters 
function tbfCoef = triBandFilter(fftMag, P, filterBankParam) 
fstart=filterBankParam(1,:); 
fcenter=filterBankParam(2,:); 
fstop=filterBankParam(3,:); 
% Triangular bandpass filter. 
for i=1:P 
   for j = fstart(i):fcenter(i), 
      filtmag(j) = (j-fstart(i))/(fcenter(i)-fstart(i)); 
   end 
   for j = fcenter(i)+1:fstop(i), 
      filtmag(j) = 1-(j-fcenter(i))/(fstop(i)-fcenter(i)); 
   end 
   tbfCoef(i) = sum(fftMag(fstart(i):fstop(i)).*filtmag(fstart(i):fstop(i))'); 
end 
tbfCoef=log(eps+tbfCoef.^2); 
 
%% === TBF coefficients to MFCC 
function mfcc = melCepstrum(L, P, tbfCoef) 
% DCT to find MFCC 
for i = 1:L 
   coef = cos((pi/P)*i*(linspace(1,P,P)-0.5))'; 
   mfcc(i) = sum(coef.*tbfCoef'); 
end 
 
%% === Delta function 
function parameter = deltaFunction(deltaWindow,parameter) 
% compute delta cepstrum and delta log energy. 
rows  = size(parameter,1); 
cols  = size(parameter,2); 
%temp  = [zeros(rows,deltaWindow) parameter zeros(rows,deltaWindow)]; 
temp  = [parameter(:,1)*ones(1,deltaWindow) parameter parameter(:,end)*ones(1,deltaWindow)]; 
temp2 = zeros(rows,cols); 
denominator = sum([1:deltaWindow].^2)*2; 
for i = 1+deltaWindow : cols+deltaWindow, 
   subtrahend = 0; 
   minuend    = 0; 
   for j = 1 : deltaWindow, 
      subtrahend = subtrahend + temp(:,i+j)*j; 
      minuend = minuend + temp(:,i-j)*(-j); 
   end; 
   temp2(:,i-deltaWindow) = (subtrahend + minuend)/denominator; 
end; 
parameter = temp2; 
 
 
function auto_correlaction=LPC(Sample, order) 
len = length(Sample); 
for k=1:order 
                %矩陣內積 
                R(k,1)=dot(Sample(k:end), Sample (1:len-k+1) );               
end 
auto_correlaction=R;