MFCC特征提取(可用程序)

/******************************************************MFCC特征提取程序*读取一个音频文件(.wav)，将根据帧长分割后的每帧2阶MFCC*系数写在输出文件中，以,为间隔*******************************************************/#includeiostream#includefstream#includestring#includecstdio#includecmath#includevector#includecomplex#includebitset//MFCC的特征选取，最后提取了13维左右最后结果保存在文件中usingnamespacestd;#defineSP_EMPHASIS_FACTOR0.97f/*预加重的系数*/typedefstruct_TWavHeader{intrId;//标志符（RIFF）intrLen;//数据大小,包括数据头的大小和音频文件的大小intwId;//格式类型（WAVE）intfId;//fmtintfLen;//Sizeof(WAVEFORMATEX)shortwFormatTag;//编码格式，包括WAVE_FORMAT_PCM，WAVEFORMAT_ADPCM等shortnChannels;//声道数，单声道为1，双声道为2intnSamplesPerSec;//采样频率intnAvgBytesPerSec;//每秒的数据量shortnBlockAlign;//块对齐shortwBitsPerSample;//WAVE文件的采样大小intdId;//dataintwSampleLength;//音频数据的大小}TWavHeader;constintFS=16;/*修改帧长*/constlongFrmLen=1024;//可修改帧长constunsignedlongFFTLen=512;//参与FFT运算的512个数据constdoublePI=3.1415926536;constintFiltNum=40;//滤波器组数，一共40组constintPCEP=13;//最后得到的关于的13个MFCC的系数doubleHamming[FrmLen];inttemp_1;//计算次数的vectordoubleMFCCcoefficient;staticdoublelast=0;//一窗数据最后一个点的值，此点用于预加重voidpreemphasis(double*buf,double*result,shortFrmLen);//预加重voidInitHamming();voidHammingWindow(double*result,double*data);voidcompute_fft(double*buffer,vectorcomplexdouble&vecList);voidFFT(constunsignedlong&ulN,vectorcomplexdouble&vecList);//FFT的实际程序voidInitFilt(double*FiltCoe1,double*FiltCoe2,int*Num);//初始化滤波器voidCFilt(double*spdata,double*FiltCoe1,double*FiltCoe2,int*Num,double*En,vectorcomplexdouble&vecList);//计算每个滤波器组内的总能量voidMFCC(double*En,double*Cep);//计算MFCC的13个系数intmain(){//TWavHeaderwaveheader;FILE*sourcefile,*MFCCFile;shortbuffer[FrmLen];doubledBuff[FrmLen];doubleresult[FrmLen];//预加重结果doubledata[FrmLen];//加窗后得到的数据doubleFiltCoe1[FFTLen/2+1];//左系数doubleFiltCoe2[FFTLen/2+1];//右系数intNum[FFTLen/2+1];//一般而言，每个点会包含在相邻的两个滤波器中，这里是与该点相关的第二个滤波器doubleEn[FiltNum+1];//频带能量doubleCep[PCEP];//MFCC结果//对定义的变量赋初值temp_1=0;inti=0;for(i=0;iFrmLen;i++){buffer[i]=0;dBuff[i]=result[i]=data[i]=0.0f;}for(i=0;iFFTLen/2+1;i++){FiltCoe1[i]=FiltCoe2[i]=0.0f;Num[i]=0;}for(i=0;iFiltNum+1;i++){En[i]=0.0f;}vectorcomplexdoublevecList;//FFT计算之后的数据//if((sourcefile=fopen(E:\test\test\a.wav,rb))==NULL)//C4996errno_terr;/*输入文件位置*/charinfilename[30];cout请输入音频文件路径(需要加入后缀名*.wav)\n;cininfilename;if((err=fopen_s(&sourcefile,infilename,rb))!=0)//Note:fopenisdeprecated;considerusingfopen_sinsteadprintf(Thefile'test.wav'wasnotopened\n);elseprintf(Thefile'test.wav'wasopened\n);/*输出文件位置*/charoutfilename[30];cout请输入MFCC输出文件路径(需要加入后缀名*.txt)\n;cinoutfilename;if((err=fopen_s(&MFCCFile,outfilename,wb))!=0)//Note:fopenisdeprecated;considerusingfopen_sinsteadprintf(Thefile'MFCCResult.txt'wasnotopened\n);elseprintf(Thefile'MFCCResult.txt'wasopened\n);//fread(&waveheader,sizeof(struct_TWavHeader),1,sourcefile);//ReadsdatafromastreamInitHamming();//初始化汉明窗InitFilt(FiltCoe1,FiltCoe2,Num);//初始化MEL滤波系数while(fread(buffer,sizeof(short),FrmLen,sourcefile)==FrmLen){for(intj=0;jFrmLen;j++)dBuff[j]=(double)buffer[j];//拿到一帧数据preemphasis(dBuff,result,FrmLen);//预加重结果存在result里面HammingWindow(result,data);//给一帧数据加窗,存在data里面compute_fft(data,vecList);CFilt(data,FiltCoe1,FiltCoe2,Num,En,vecList);MFCC(En,Cep);//charendl='\n';for(intj=0;jPCEP;j++){if(j==1){fprintf(MFCCFile,%.1f,,Cep[j]);temp_1++;}//if(j==PCEP-1)//{fprintf(MFCCFile,%f,,Cep[j]);//temp_1++;}}vecList.clear();//Movesthefilepointertoaspecifiedlocation.fseek(sourcefile,-FrmLen/2,SEEK_CUR);//考虑到帧移，每次移动半帧}//fprintf(MFCCFile,%d,,temp_1);//写总个数fclose(sourcefile);/*关闭文件*/fclose(MFCCFile);/*关闭写MFCC结果文件*/intlength=MFCCcoefficient.size();for(inti=0;ilength;++i){coutMFCCcoefficient[i]'';if((i+1)%13==0)coutendl;}//getchar();return0;}//预加重voidpreemphasis(double*buf,double*result,shortFrmLen){inti;result[0]=buf[0]-SP_EMPHASIS_FACTOR*last;for(i=1;iFrmLen;i++){result[i]=buf[i]-SP_EMPHASIS_FACTOR*buf[i-1];}last=buf[(FrmLen-1)/2];//假设每次移半帧}//汉明窗初始化voidInitHamming(){doubletwopi;inti;twopi=2*PI;for(i=0;iFrmLen;i++){Hamming[i]=(double)(0.54-0.46*cos(i*twopi/(double)(FrmLen-1)));}}//给一帧数据加窗voidHammingWindow(double*result,double*data){inti;for(i=0;iFrmLen;i++){data[i]=result[i]*Hamming[i];}}voidcompute_fft(double*data,vectorcomplexdouble&vecList){for(inti=0;iFFTLen;++i){if(iFrmLen){complexdoubletemp(data[i]);vecList.push_back(temp);}else{complexdoubletemp(0);//如果计算的FFT长度大于窗长，则不足的部分用零填充。得到的效果差不多vecList.push_back(temp);}}FFT(FFTLen,vecList);}voidFFT(constunsignedlong&ulN,vectorcomplexdouble&vecList){//得到指数，这个指数实际上指出了计算FFT时内部的循环次数unsignedlongulPower=0;//指数unsignedlongulN1=ulN-1;//ulN1=511while(ulN10){ulPower++;ulN1/=2;}//反序，因为FFT计算后的结果次序不是顺序的，需要反序来调整。可以在FFT实质部分计算之前先调整，也可以在结果//计算出来后再调整。本程序中是先调整，再计算FFT实质部分bitsetsizeof(unsignedlong)*8bsIndex;//二进制容器unsignedlongulIndex;//反转后的序号unsignedlongulK;for(unsignedlonglongp=0;pulN;p++){ulIndex=0;ulK=1;bsIndex=bitsetsizeof(unsignedlong)*8(p);for(unsignedlongj=0;julPower;j++){ulIndex+=bsIndex.test(ulPower-j-1)?ulK:0;ulK*=2;}if(ulIndexp)//只有大于时，才调整，否则又调整回去了{complexdoublec=vecList[p];vecList[p]=vec