修改后的MFCC提取

  http://blog.sina.com.cn/u/4840fe2a01000208

//  KMFCC.h: interface for the KMFCC class.
//
//  提取 MFCC 倒谱
//
//  2006年2月8日 柯登峰
//
/**/ //////////////////////////////////////////////////////////////////////
#include < vector >
#include " krfft.h "
#ifndef _KMFCC_KDF_2006_FEB_08_
#define  _KMFCC_KDF_2006_FEB_08_
class  KMFCC  
... {
 KRFFT rfft; // 快速福立叶变换
 float * win;
 float * mel;
 float * real;
 float * imag;
 float * lift;
 int * beg; // 滤波器开始点
 int * end; // 滤波器结束点
 std::vector< std::vector<float> > bank; // 滤波器组
private:
 int m_nCep;
 int m_nFFT;
 int m_nWin;
 int m_nHop;
 int m_nChan;
 int m_nLift;
 int m_nSamRate;
public:
 void MelBank2Mfcc( float * MelBank, float * mfcc );
 void Wav2MelBank( short * wav, float * MelBank );
 void Delta( float * x , float * dx, int nSubDim, int nFrame, int nTotalDim );
 void Wav2Mfcc( short * wav , float * mfcc );
 void CreateMFCC(int nSamRate, int nCep=13, int nWin=1024, int nHop=512, int nChan=35, int nLift=35);
 KMFCC();
 virtual ~KMFCC();
protected:
 void DCT( float * in, float * out, int nChan, int nCep);
 float Frq2Mel( float frq );
 float Mel2Frq( float mel );
 void CreatBank(int nFFT, int nChan,int nSamRate);
} ;
#endif   // _KMFCC_KDF_2006_FEB_08_
 
//  KMFCC.cpp: implementation of the KMFCC class.
//
//  提取 MFCC 倒谱
//
//  2005年2月8日 柯登峰
//
/**/ //////////////////////////////////////////////////////////////////////
#include  " stdafx.h "
#include  " KMFCC.h "
#include  < math.h >
/**/ //////////////////////////////////////////////////////////////////////
//  Construction/Destruction
/**/ //////////////////////////////////////////////////////////////////////
KMFCC::KMFCC():
m_nCep( 0 ),
m_nWin( 0 ),
m_nHop( 0 ),
m_nChan( 0 ),
m_nLift( 0 ),
m_nFFT( 0 ),
real(NULL),
imag(NULL),
mel(NULL),
win(NULL),
beg(NULL),
end(NULL),
lift(NULL)
... {
}
KMFCC:: ~ KMFCC()
... {
 delete [] lift; lift = NULL;
 delete [] real; real = NULL;
 delete [] imag; imag = NULL;
 delete [] mel; mel = NULL;
 delete [] win; win = NULL;
 delete [] beg; beg = NULL;
 delete [] end; end = NULL;
}
void  KMFCC::CreateMFCC( int  nSamRate,  int  nCep,  int  nWin,  int  nHop,  int  nChan,  int  nLift)
... {
 const float pi = (float)3.14159265358979;
 
 m_nFFT = (int)(log(nWin-1.0)/log(2.0)+1);
 m_nFFT = (int)pow(2,m_nFFT);
 
 m_nCep = nCep;
 m_nWin = nWin;
 m_nHop = nHop;
 m_nChan = nChan;
 m_nLift = nLift;
 m_nSamRate = nSamRate;
 
 CreatBank(m_nFFT,nChan,nSamRate);
 
 delete [] lift; lift = new float[nCep+1];
 delete [] real; real = new float[m_nFFT];
 delete [] imag; imag = new float[m_nFFT];
 delete [] mel; mel = new float[nChan];
 delete [] win; win = new float[nWin];
 
 for(int i=0;i<nWin;i++)
  win[i] = 0.54 - 0.46*cos(2*pi*(i+1)/(nWin+1));
 for(i=0;i<nCep;i++)
  lift[i] = (1.0 + 0.5*nLift*sin(pi*(i+1)/nLift))/(nLift+1.0);
}
void  KMFCC::CreatBank( int  nFFT,  int  nChan,  int  nSamRate)
... {
 int i,j;
 float dm,s,c,e,m,sum;
 float pi = (float)3.14159265359;
 
 float * cen = new float[nChan+2];
 
 bank.resize(nChan);
 for(i=0;i<nChan;i++)...{
  bank[i].resize(nFFT);
 }
 
 delete []beg; beg = new int[nChan];
 delete []end; end = new int[nChan];
 
 dm = Frq2Mel(nSamRate/2.0)/(nChan+1);
 for(i=0;i<=nChan+1;i++)...{
  cen[i] = Mel2Frq(dm*i)/nSamRate*nFFT;
 }
 
 for(i=0;i<nChan;i++)...{
  s = cen[i]; // 开始
  c = cen[i+1]; // 中央
  e = cen[i+2]; // 结束
  
  beg[i] = int(s+1);
  end[i] = int(e-1);
  
  sum = 0.0;
  for(j=beg[i];j<=c;j++)...{
   m = j-s;
   sum += m;
   bank[i][j] = m;
  }
  for(;j<=end[i];j++)...{
   m = e-j;
   sum += m;
   bank[i][j] = m;
  }
  for(j=beg[i];j<=end[i];j++)...{
   bank[i][j] /= sum;
  }
 }
 
 delete [] cen;
}
float  KMFCC::Frq2Mel( float  frq) ... { return log(1.0+frq/700.0)*1127.01048; }
float  KMFCC::Mel2Frq( float  mel) ... { return 700.0*(exp(mel/1127.01048)-1.0); }
void  KMFCC::Wav2Mfcc( short   * wav,  float   *  mfcc)
... {
 Wav2MelBank(wav,mel);
 MelBank2Mfcc(mel,mfcc);
}
void  KMFCC::DCT( float   * in ,  float   * out , int  nChan,  int  nCep)
... {
 int i,k;
 float m,pi = (float)3.14159265359;
 float w = 0.0;
 
 for(k=1;k<nCep;k++)...{ // 第 0 维放能量
  m = 0.0;
  for(i=0;i<nChan;i++)...{
   w = cos(pi*k*(i+0.5)/(nChan));
   m += w * in[i];
  }
  out[k] = m;
 }
}
void  KMFCC::Delta( float   * x,  float   * dx,  int  nDim,  int  nFrame,  int  nHop)
... {
 int i,j,k,t1,t2;
 float s;
 
 for(j=0;j<nDim;j++)...{
  for(i=0;i<nFrame;i++)...{
   s = 0.0; 
   t1=i; 
   t2=i;
   for(k=1;k<=2;k++)...{
    if(t1>0) t1--;
    if(t2<nFrame-1) t2++;
    s += k*(x[t2*nHop+j]-x[t1*nHop+j]);
   }
   dx[i*nHop+j] = s/3; // k = 1,2  so sum(1,2)=3
  }
 }
}
void  KMFCC::Wav2MelBank( short   * wav,  float   * MelBank)
... {
 int i,j;
 float m;
 
 // 零均值化
 m=0; 
 for(i=0;i<m_nWin;i++) m+=wav[i];
 m/=m_nWin; 
 
 // 高频预加重
 real[0]=wav[0]*win[0];
 for(i=1;i<m_nWin;i++)
  real[i] = (wav[i]-0.97*wav[i-1]-0.03*m)*win[i];
 for(i=m_nWin;i<m_nFFT;i++)
  real[i] = 0.0;
 
 // 傅立叶谱
 rfft.doRFFT(real,imag,m_nFFT,0);
 
 // 功率谱
 for(i=0;i<m_nFFT;i++)...{
  real[i] = (real[i]*real[i]+imag[i]*imag[i]);
 }
 
 // 滤波器组
 for(i=0;i<m_nChan;i++)...{
  m = 0.0;
  for(j=beg[i];j<=end[i];j++)...{
   m += bank[i][j] * real[j];
  }
  MelBank[i] = m;
 }
}
void  KMFCC::MelBank2Mfcc( float   * MelBank,  float   * mfcc)
... {
 int i;
 float m;
 float eps=(float)(1e-0);
 
 // 波动能量
 m = 0.0; for(i=1;i<m_nChan;i++) m += MelBank[i];
 if(m<eps) m = eps;
 mfcc[0] = log(m); // 第 0 维放能量
 
 // 对数能量
 for(i=0;i<m_nChan;i++)...{
  if(MelBank[i]<eps) MelBank[i] = eps;
  MelBank[i] = log(MelBank[i]);
 }
 
 // 倒谱
 DCT(MelBank,mfcc,m_nChan,m_nCep);
 
 // 提升
 for(i=1;i<m_nCep;i++)
  mfcc[i] *= lift[i];
}

 
//  KFE.h: interface for the KFE class.
//
//  特征提取 0.013x实时
//
//  2005年12月28日 柯登峰
//
/**/ //////////////////////////////////////////////////////////////////////
#ifndef _FE_KDF_2005_DEC_28_
#define  _FE_KDF_2005_DEC_28_
#include < windows.h >
#include  " krfft.h "
#include  " kmfcc.h "
struct  FEHEADER ... {
 int  nFsKHz;   // 采样频率（千赫兹）
 int  nFrmMs;   // 帧长（毫秒）
 int  nSftMs;   // 帧移（毫秒）
 int  nMaxFrm;  // 最大帧数
 int  nFtrDim;  // 特征维数
 char szFtrInfo[256]; // 特征信息  
 char szVersion[256]; // 版本信息
} ;
class  KFE  
... {
public:
 BOOL KCloseFE( );
 BOOL KRunFE( short * pWav, int nLen, float *& feature, int & nFrame  );
 BOOL KOpenFE( char * dir , FEHEADER & hdr );
 KFE();
 virtual ~KFE();
private:
 KMFCC mfcc;
 float * m_mean;
 float * m_bank;
 float * m_feature;
 int m_nFrame;
 int m_nWin;
 int m_nHop;
 int m_nDim;
 int m_nChan;
 int m_nMaxFrm; 
} ;
#endif   // _FE_KDF_2005_DEC_28_

//  KFE.cpp: implementation of the KFE class.
//
//  特征提取
//
//  2005年12月28日 柯登峰
//
/**/ //////////////////////////////////////////////////////////////////////
#include  " stdafx.h "
#include  " KFE.h "
/**/ //////////////////////////////////////////////////////////////////////
//  Construction/Destruction
/**/ //////////////////////////////////////////////////////////////////////
KFE::KFE(): m_feature(NULL),m_mean(NULL),m_bank(NULL)
... {
}
KFE:: ~ KFE()
... {
 delete [] m_feature; m_feature = 0;
 delete [] m_mean; m_mean = 0;
 delete [] m_bank; m_bank = 0;
}
BOOL KFE::KOpenFE( char   *  dir, FEHEADER  &  hdr)
... {
 m_nWin = hdr.nFrmMs*hdr.nFsKHz;
 m_nHop = hdr.nSftMs*hdr.nFsKHz;
 m_nMaxFrm = hdr.nMaxFrm;
 m_nDim = hdr.nFtrDim*3; // MFCC,D,A
 m_nChan = 26;
 mfcc.CreateMFCC(hdr.nFsKHz*1000,hdr.nFtrDim,m_nWin,m_nHop,m_nChan,m_nChan);
 delete [] m_feature; m_feature = new float[m_nMaxFrm*m_nDim];
 delete [] m_bank; m_bank = new float[m_nMaxFrm*m_nChan];
 delete [] m_mean; m_mean = new float[m_nChan];
 
 hdr.nFtrDim = m_nDim;
 strcpy(hdr.szFtrInfo,"MFCC_E_D_A_Z");
 strcpy(hdr.szVersion,"柯登峰 2006年2月8日");
 return TRUE;
}
BOOL KFE::KRunFE( short   * pWav,  int  nLen,  float   *&  feature,  int   & nFrame)
... {
 int i,j,dim = m_nDim/3;
 nFrame = 0;
 // MFCC
 for(i=0,j=0; i+m_nWin<nLen; i+=m_nHop,j+=m_nDim,nFrame++)...{
  mfcc.Wav2Mfcc(pWav+i,m_feature+j);
 }
 // Z
 for(i=0;i<dim;i++)...{
  float mean = 0.0;
  for(j=0;j<nFrame;j++)
   mean += m_feature[j*m_nDim+i];
  mean /= nFrame;
  for(j=0;j<nFrame;j++)
   m_feature[j*m_nDim+i] -= mean;
 }
 // D
 mfcc.Delta(m_feature,m_feature+dim,dim,nFrame,dim*3);
 // A
 mfcc.Delta(m_feature+dim,m_feature+2*dim,dim,nFrame,dim*3);
 feature = m_feature;
 return TRUE;
}
BOOL KFE::KCloseFE()
... {
 return TRUE;
}