www.pudn.com > asr.rar > MELCEPNormFreq.M, change:2004-11-06,size:6270b


function [c,dStd]=melcepNormFreq(s,fs,w,nc,p,n,inc,fl,fh) 
%MELCEPSTNormFreq Calculate the mel cepstrum of a signal C=(S,FS,W,NC,P,N,INC,FL,FH) 
% mfcc  coefficient are normalizated through all the frame in frequency(CMS), and weight(/), no ceptrum 
% energy coefficient is normalizated by Max. 
% dStd is weight coefficient to watch by plot 
 
% Simple use: c=melcepst(s,fs)	% calculate mel cepstrum with 12 coefs, 256 sample frames 
%				  c=melcepst(s,fs,'e0dD') % include log energy, 0th cepstral coef, delta and delta-delta coefs 
% 
% Inputs: 
%     s	 speech signal 
%     fs  sample rate in Hz (default 11025) 
%     nc  number of cepstral coefficients excluding 0'th coefficient (default 12) 
%     n   length of frame (default power of 2 <30 ms)) 
%     p   number of filters in filterbank (default floor(3*log(fs)) ) 
%     inc frame increment (default n/2) 
%     fl  low end of the lowest filter as a fraction of fs (default = 0) 
%     fh  high end of highest filter as a fraction of fs (default = 0.5) 
% 
%		w   any sensible combination of the following: 
% 
%				'R'  rectangular window in time domain 
%				'N'	Hanning window in time domain 
%				'M'	Hamming window in time domain (default) 
% 
%		      't'  triangular shaped filters in mel domain (default) 
%		      'n'  hanning shaped filters in mel domain 
%		      'm'  hamming shaped filters in mel domain 
% 
%				'p'	filters act in the power domain 
%				'a'	filters act in the absolute magnitude domain (default) 
% 
%			   '0'  include 0'th order cepstral coefficient 
%				'e'  include log energy 
%				'd'	include delta coefficients (dc/dt) 
%				'D'	include delta-delta coefficients (d^2c/dt^2) 
% 
%		      'z'  highest and lowest filters taper down to zero (default) 
%		      'y'  lowest filter remains at 1 down to 0 frequency and 
%			   	  highest filter remains at 1 up to nyquist freqency 
% 
%		       If 'ty' or 'ny' is specified, the total power in the fft is preserved. 
% 
% Outputs:	c     mel cepstrum output: one frame per row 
% 
 
 
%      Copyright (C) Mike Brookes 1997 
% 
%      Last modified Thu Jun 15 09:14:48 2000 
% 
%   VOICEBOX is a MATLAB toolbox for speech processing. Home page is at 
%   http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 
%  
% Modified by LXL  
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   This program is free software; you can redistribute it and/or modify 
%   it under the terms of the GNU General Public License as published by 
%   the Free Software Foundation; either version 2 of the License, or 
%   (at your option) any later version. 
% 
%   This program is distributed in the hope that it will be useful, 
%   but WITHOUT ANY WARRANTY; without even the implied warranty of 
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
%   GNU General Public License for more details. 
% 
%   You can obtain a copy of the GNU General Public License from 
%   ftp://prep.ai.mit.edu/pub/gnu/COPYING-2.0 or by writing to 
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
if nargin<2 fs=11025; end 
if nargin<3 w='M'; end 
if nargin<4 nc=12; end 
if nargin<5 p=floor(3*log(fs)); end 
if nargin<6 n=pow2(floor(log2(0.03*fs))); end 
if nargin<9 
   fh=0.5;    
   if nargin<8 
     fl=0; 
     if nargin<7 
        inc=floor(n/2); 
     end 
  end 
end 
 
if any(w=='R') 
   z=enframe(s,n,inc); 
elseif any (w=='N') 
   z=enframe(s,hanning(n),inc); 
else 
   z=enframe(s,hamming(n),inc); 
end 
f=rfft(z.'); 
[m,a,b]=melbankm(p,n,fs,fl,fh,w); 
pw=f(a:b,:).*conj(f(a:b,:)); 
pth=max(pw(:))*1E-6; 
if any(w=='p') 
   y=log(max(m*pw,pth)); 
else 
   ath=sqrt(pth); 
   y=log(max(m*abs(f(a:b,:)),ath)); 
end 
% H(z)=z-1/z ;  fy=y*h(n) 
[nfy,ncy]=size(y); 
for j=1:ncy 
    for i=1:nfy-2 
        fy(i,j)=y(i+2,j)-y(i,j); 
    end 
end     
y=fy; 
% y is normalizated 
y=y'; 
[nfy,ncy]=size(y); 
yMean=mean(y); 
ymMean=zeros(nfy,ncy); 
for j=1:ncy 
    for i=1:nfy 
        ymMean(i,j)=yMean(j); 
    end 
end 
c=y; 
[nf,nc]=size(c); 
cStd=std(c,1); 
dStd=1./cStd; 
sStd=zeros(nf,nc); 
for j=1:nc 
    for i=1:nf 
        sStd(i,j)=cStd(j); 
    end 
end 
c=(c-ymMean)./sStd; 
 
if any(w=='e') 
   % energy normalization    
   egy=sum(pw);  
   egymax=max(egy);  
   egy=egy/egymax; 
   c=[log(egy).' c];    
%   c=[log(sum(pw)).' c]; 
end 
 
% calculate derivative 
if any(w=='D') 
  vf=(4:-1:-4)/60; 
  af=(1:-1:-1)/2; 
  ww=ones(5,1); 
  cx=[c(ww,:); c; c(nf*ww,:)]; 
  ncc=nc; 
  if any(w=='e') 
      ncc=nc+1; 
  end 
  vx=reshape(filter(vf,1,cx(:)),nf+10,ncc); 
  vx(1:8,:)=[]; 
  ax=reshape(filter(af,1,vx(:)),nf+2,ncc); 
  ax(1:2,:)=[]; 
  vx([1 nf+2],:)=[]; 
  % vx is normalizated 
  vxMean=mean(vx); 
  vxStd=std(vx,1); 
  [sx,sy]=size(vx); 
  mMean=zeros(sx,sy); 
  for j=1:sy 
    for i=1:sx 
        mMean(i,j)=vxMean(j); 
    end 
  end 
  sStd=zeros(sx,sy); 
  for j=1:sy 
    for i=1:sx 
        sStd(i,j)=vxStd(j); 
    end 
  end 
  vx=(vx-mMean)./sStd; 
  % ax is normalizated 
  axMean=mean(ax); 
  axStd=std(ax,1); 
  [sx,sy]=size(ax); 
  mMean=zeros(sx,sy); 
  for j=1:sy 
    for i=1:sx 
        mMean(i,j)=axMean(j); 
    end 
  end 
  sStd=zeros(sx,sy); 
  for j=1:sy 
    for i=1:sx 
        sStd(i,j)=axStd(j); 
    end 
  end 
  ax=(ax-mMean)./sStd; 
 
  if any(w=='d') 
     c=[c vx ax]; 
  else 
     c=[c ax]; 
  end 
   
elseif any(w=='d') 
  vf=(4:-1:-4)/60; 
  ww=ones(4,1); 
  cx=[c(ww,:); c; c(nf*ww,:)]; 
  ncc=nc; 
  if any(w=='e') 
      ncc=nc+1; 
  end 
  vx=reshape(filter(vf,1,cx(:)),nf+8,ncc); 
  vx(1:8,:)=[]; 
  % vx is normalizated 
  vxMean=mean(vx); 
  vxStd=std(vx,1); 
  [sx,sy]=size(vx); 
  mMean=zeros(sx,sy); 
  for j=1:sy 
    for i=1:sx 
        mMean(i,j)=vxMean(j); 
    end 
  end 
  sStd=zeros(sx,sy); 
  for j=1:sy 
    for i=1:sx 
        sStd(i,j)=vxStd(j); 
    end 
  end 
  vx=(vx-mMean)./sStd; 
  c=[c vx]; 
end 
  
if nargout<1 
   [nf,nc]=size(c); 
   t=((0:nf-1)*inc+(n-1)/2)/fs; 
   ci=(1:nc)-any(w=='0')-any(w=='e'); 
   imh = imagesc(t,ci,c.'); 
   axis('xy'); 
   xlabel('Time (s)'); 
   ylabel('Mel-cepstrum coefficient'); 
	map = (0:63)'/63; 
	colormap([map map map]); 
	colorbar; 
end