《MFCC特征提取.docx》由会员分享，可在线阅读，更多相关《MFCC特征提取.docx（10页珍藏版）》请在三一办公上搜索。
1、MFCC特征提取/* *MFCC特征提取程序 *读取一个音频文件(.wav)，将根据帧长分割后的每帧2阶MFCC *系数写在输出文件中，以,为间隔 */ #include #include #include #include #include #include #include #include /MFCC的特征选取，最后提取了13维左右最后结果保存在文件中 using namespace std; #define SP_EMPHASIS_FACTOR 0.97f /* 预加重的系数 */ typedef struct _TWavHeader int rId; /标志符 int rLen; /
2、数据大小,包括数据头的大小和音频文件的大小 int wId; /格式类型 int fId; /fmt int fLen; /Sizeof(WAVEFORMATEX) short wFormatTag; /编码格式，包括WAVE_FORMAT_PCM，WAVEFORMAT_ADPCM等 short nChannels; /声道数，单声道为1，双声道为2 int nSamplesPerSec; /采样频率 int nAvgBytesPerSec; /每秒的数据量 short nBlockAlign; /块对齐 short wBitsPerSample; /WAVE文件的采样大小 int dId;
3、/data int wSampleLength; /音频数据的大小 TWavHeader; const int FS=16; /*修改帧长*/ const long FrmLen=1024; /可修改帧长 const unsigned long FFTLen=512; /参与FFT运算的512个数据 const double PI=3.1415926536; const int FiltNum=40; /滤波器组数，一共40组 const int PCEP=13; /最后得到的关于的13个MFCC的系数 double HammingFrmLen; int temp_1; /计算次数的 vect
4、orMFCCcoefficient; static double last=0; /一窗数据最后一个点的值，此点用于预加重 void preemphasis(double* buf, double* result, short FrmLen ); /预加重 void InitHamming; void HammingWindow(double* result,double* data); void compute_fft(double *buffer,vectorcomplex & vecList); void FFT(const unsigned long & ulN, vectorcomp
5、lex & vecList); /FFT的实际程序 void InitFilt(double *FiltCoe1, double *FiltCoe2, int *Num); /初始化滤波器 void CFilt(double *spdata, double *FiltCoe1, double *FiltCoe2, int *Num, double *En,vectorcomplex & vecList); /计算每个滤波器组内的总能量 void MFCC(double *En, double *Cep); /计算MFCC的13个系数 int main / TWavHeader wavehead
6、er; FILE *sourcefile, *MFCCFile; short bufferFrmLen; double dBuffFrmLen; double resultFrmLen; /预加重结果 double dataFrmLen; /加窗后得到的数据 double FiltCoe1FFTLen/2+1; /左系数 double FiltCoe2FFTLen/2+1; /右系数 int NumFFTLen/2+1; /一般而言，每个点会包含在相邻的两个滤波器中，这里是与该点相关的第二个滤波器 double EnFiltNum+1; /频带能量 double CepPCEP;/MFCC结果
7、 /对定义的变量赋初值 temp_1=0; int i = 0; for (i = 0; i FrmLen; i+) bufferi = 0; dBuffi = resulti = datai = 0.0f; for (i = 0; i FFTLen/2 + 1; i+) FiltCoe1i = FiltCoe2i = 0.0f; Numi = 0; for(i = 0; i FiltNum + 1; i+) Eni = 0.0f; vectorcomplex vecList;/FFT计算之后的数据 /if( (sourcefile = fopen( E:testtesta.wav, rb )
8、 = NULL ) / C4996 errno_t err; /* 输入文件位置*/ char infilename30; coutinfilename; if( (err = fopen_s( &sourcefile, infilename, rb ) !=0 ) / Note: fopen is deprecated; consider using fopen_s instead printf( The file test.wav was not openedn ); else printf( The file test.wav was openedn ); /*输出文件位置*/ char
9、 outfilename30; coutoutfilename; if( (err = fopen_s( &MFCCFile, outfilename, wb ) !=0 ) / Note: fopen is deprecated; consider using fopen_s instead printf( The file MFCCResult.txt was not openedn ); else printf( The file MFCCResult.txt was openedn ); /fread(&waveheader,sizeof(struct _TWavHeader),1,s
10、ourcefile); /Reads data from a stream InitHamming;/初始化汉明窗 InitFilt(FiltCoe1,FiltCoe2,Num); /初始化MEL滤波系数 while ( fread(buffer,sizeof(short),FrmLen,sourcefile) = FrmLen ) for (int j = 0; j FrmLen; j+) dBuffj = (double)bufferj;/拿到一帧数据 preemphasis(dBuff,result,FrmLen);/预加重结果存在result里面 HammingWindow(resul
11、t,data); /给一帧数据加窗,存在data里面 compute_fft(data,vecList); CFilt(data, FiltCoe1, FiltCoe2, Num, En,vecList); MFCC(En, Cep); /char endl = n; for(int j = 0; j PCEP; j+) if(j=1) fprintf(MFCCFile, %.1f, Cepj); temp_1+; /if ( j = PCEP -1) / fprintf(MFCCFile, %f, Cepj); /temp_1+; vecList.clear; /Moves the file
12、 pointer to a specified location. fseek(sourcefile, -FrmLen/2, SEEK_CUR);/考虑到帧移，每次移动半帧 /fprintf(MFCCFile,%d,temp_1); /写总个数 fclose(sourcefile); /*关闭文件*/ fclose(MFCCFile); /*关闭写MFCC结果文件*/ int length = MFCCcoefficient.size; for(int i=0; i length; +i) cout MFCCcoefficienti ; if( (i+1) % 13 = 0 ) cout en
13、dl; /getchar; return 0; /预加重 void preemphasis(double* buf, double* result, short FrmLen) int i; result0 = buf0 - SP_EMPHASIS_FACTOR * last; for(i=1;iFrmLen;i+) resulti = bufi - SP_EMPHASIS_FACTOR * bufi-1; last = buf(FrmLen-1)/2; /假设每次移半帧 /汉明窗初始化 void InitHamming double twopi; int i; twopi=2*PI; for
14、( i=0;iFrmLen;i+) Hammingi=(double)(0.54-0.46*cos(i*twopi/(double)(FrmLen-1); /给一帧数据加窗 void HammingWindow(double* result,double* data) int i; for(i=0;iFrmLen;i+) datai=resulti*Hammingi; void compute_fft(double *data,vectorcomplex & vecList) for(int i=0;iFFTLen;+i) if(iFrmLen) complex temp(datai); ve
15、cList.push_back(temp); else complex temp(0); /如果计算的FFT长度大于窗长，则不足的部分用零填充。得到的效果差不多 vecList.push_back(temp); FFT(FFTLen,vecList); void FFT(const unsigned long & ulN, vectorcomplex & vecList) /得到指数，这个指数实际上指出了计算FFT时内部的循环次数 unsigned long ulPower = 0; /指数 unsigned long ulN1 = ulN - 1; /ulN1=511 while(ulN1
16、0) ulPower+; ulN1 /= 2; /反序，因为FFT计算后的结果次序不是顺序的，需要反序来调整。可以在FFT实质部分计算之前先调整，也可以在结果 /计算出来后再调整。本程序中是先调整，再计算FFT实质部分 bitset bsIndex; /二进制容器 unsigned long ulIndex; /反转后的序号 unsigned long ulK; for(unsigned long long p = 0; p ulN; p+) ulIndex = 0; ulK = 1; bsIndex = bitset(p); for(unsigned long j = 0; j p) /只有
17、大于时，才调整，否则又调整回去了 complex c = vecListp; vecListp = vecListulIndex; vecListulIndex = c; /计算旋转因子 vectorcomplex vecW; for(unsigned long i = 0; i ulN / 2; i+) vecW.push_back(complex(cos(2 * i * PI / ulN) , -1 * sin(2 * i * PI / ulN); /计算FFT unsigned long ulGroupLength = 1; /段的长度 unsigned long ulHalfLengt
18、h = 0; /段长度的一半 unsigned long ulGroupCount = 0; /段的数量 complex cw; /WH(x) complex c1; /G(x) + WH(x) complex c2; /G(x) - WH(x) for(unsigned long b = 0; b ulPower; b+) ulHalfLength = ulGroupLength; ulGroupLength *= 2; for(unsigned long j = 0; j ulN; j += ulGroupLength) for(unsigned long k = 0; k ulHalfL
19、ength; k+) cw = vecWk * ulN / ulGroupLength * vecListj + k + ulHalfLength; c1 = vecListj + k + cw; c2 = vecListj + k - cw; vecListj + k = c1; vecListj + k + ulHalfLength = c2; /* 设置滤波器参数 输入参数：无 输出参数：*FiltCoe1-三角形滤波器左边的系数 *FiltCoe2-三角形滤波器右边的系数 *Num -决定每个点属于哪一个滤波器 */ void InitFilt(double *FiltCoe1, do
20、uble *FiltCoe2, int *Num) int i,k; double Freq; double FiltFreqFiltNum+2; /40个滤波器，故有42各滤波器端点。每一个滤波器的左右端点分别是前一个及后一个滤波器的中心频率所在的点 double BWFiltNum+1; /带宽，即每个相邻端点之间的频率跨度 double low = (double)( 400.0 / 3.0 ); /* 滤波器组的最低频率，即第一个端点值 */ short lin = 13; /* 1000Hz以前的13个滤波器是线性的分布的 */ double lin_spacing = (doubl
21、e)( 200.0 / 3.0 ); /* 相邻滤波器中心的距离为66.6Hz */ short log = 27; /* 1000Hz以后是27个对数线性分布的滤波器 */ double log_spacing = 1.0711703f; /* 相邻滤波器左半边宽度的比值 */ for ( i=0; ilin; i+) FiltFreqi = low + i * lin_spacing; for ( i=lin; ilin+log+2; i+) FiltFreqi = FiltFreqlin-1 * (double)pow( log_spacing, i - lin + 1 ); for (
22、 i=0; iFiltNum+1; i+) BWi = FiltFreqi+1 - FiltFreqi; for(i = 0 ; i= FFTLen/2 ; i+ ) Numi = 0; bool bFindFilt = false; for(i = 0 ; i = FFTLen/2 ; i+) Freq = FS * 1000.0F * i / (double)(FFTLen); bFindFilt = false; for(k = 0; k = FiltFreqk & Freq = FiltFreqk+1) bFindFilt = true; if(k = FiltNum) FiltCoe
23、1i=0.0F; else FiltCoe1i = (Freq - FiltFreqk) / (double)(BWk) * 2.0f / (BWk + BWk+1); if(k = 0) FiltCoe2i = 0.0F; else FiltCoe2i = (FiltFreqk+1 - Freq) / (double)(BWk) * 2.0f / (BWk + BWk-1); Numi = k; /当k=FiltNum时，它为第FiltNum个滤波器，实际上它并不存在。这里只是为了计算方便，假设有第FiltNum个滤波器存在。 /但其实这并不影响结果 break; if (!bFindFil
24、t) Numi = 0; /这时，该点不属于任何滤波器，因为其左右系数皆为0，所以可以假定它属于某个滤波器，而不会影响结果。这里我 /将其设为第一个滤波器。 FiltCoe1i=0.0F; FiltCoe2i=0.0F; /* 根据滤波器参数计算频带能量 输入参数：*spdata -预处理之后的一帧语音信号 *FiltCoe1-三角形滤波器左边的系数 *FiltCoe2-三角形滤波器右边的系数 *Num -决定每个点属于哪一个滤波器 输出参数：*En -输出对数频带能量 */ /把属于某一频带的能量全部加起来了 /CFilt(data, FiltCoe1, FiltCoe2, Num, En,
25、vecList); veclist ： FFT计算出的结果 Num:决定每个点属于哪一个滤波器 void CFilt(double *spdata, double *FiltCoe1, double *FiltCoe2, int *Num, double *En,vectorcomplex & vecList) double temp=0; int id, id1, id2; for(id = 0; id FiltNum ; id+) Enid=0.0F; for(id = 0 ; id 0) & (id1 FiltNum) id2 = id1-1; Enid1 = Enid1 + FiltC
26、oe1id * temp; Enid2 = Enid2 + FiltCoe2id * temp; for(id = 0 ; id FiltNum ; id+) if (Enid != 0) Enid=(double)log10(Enid); /* 计算MFCC系数 输入参数：*En -对数频带能量 */ void MFCC(double* En, double* Cep) int idcep, iden; / double Cep13; for(idcep = 0 ; idcep PCEP ; idcep+) Cepidcep = 0.0f; for(iden = 0 ; iden FiltN
27、um ; iden+) /离散余弦变换 if(iden = 0) Cepidcep = Cepidcep + Eniden * (double)cos(idcep * (iden+0.5f) * PI/(FiltNum) * 10.0f * sqrt(1/(double)FiltNum); else Cepidcep = Cepidcep + Eniden * (double)cos(idcep * (iden+0.5f) * PI/(FiltNum) * 10.0f * sqrt(2/(double)FiltNum); MFCCcoefficient.push_back(Cepidcep);
