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; //标志符（RIFF）

int rLen; //数据大小,包括数据头的大小和音频文件的大小

int wId; //格式类型（"WAVE"）

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; //"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 Hamming[FrmLen];

int temp_1; //计算次数的

vectorMFCCcoefficient;

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,vector >& vecList);

void FFT(const unsigned long & ulN, vector >& vecList); //FFT的实际程序

void InitFilt(double *FiltCoe1, double *FiltCoe2, int *Num); //初始化滤波器

void CFilt(double *spdata, double *FiltCoe1, double *FiltCoe2, int *Num, double

*En,vector >& vecList); //计算每个滤波器组内的总能量

void MFCC(double *En, double *Cep); //计算MFCC的13个系数

int main()

{

// TWavHeader waveheader;

FILE *sourcefile, *MFCCFile;

short buffer[FrmLen];

double dBuff[FrmLen];

double result[FrmLen]; //预加重结果

double data[FrmLen]; //加窗后得到的数据

double FiltCoe1[FFTLen/2+1]; //左系数

double FiltCoe2[FFTLen/2+1]; //右系数

int Num[FFTLen/2+1]; //一般而言，每个点会包含在相邻的两个滤波器中，这里是与该点相关的第二个滤波器

double En[FiltNum+1]; //频带能量

double Cep[PCEP];//MFCC结果

//对定义的变量赋初值

temp_1=0;

int i = 0;

for (i = 0; i < FrmLen; i++)

{

buffer[i] = 0;

dBuff[i] = result[i] = data[i] = 0.0f;

}

for (i = 0; i < FFTLen/2 + 1; i++)

{

FiltCoe1[i] = FiltCoe2[i] = 0.0f;

Num[i] = 0;

}

for(i = 0; i < FiltNum + 1; i++)

{

En[i] = 0.0f;

}

vector> vecList;//FFT计算之后的数据

//if( (sourcefile = fopen( "E:\test\test\a.wav", "rb" )) == NULL ) // C4996

errno_t err;

/* 输入文件位置*/

char infilename[30];

cout<<"请输入音频文件路径(需要加入后缀名*.wav)\n";

cin>>infilename;

if( (err = fopen_s( &sourcefile, infilename, "rb" )) !=0 )

// Note: fopen is deprecated; consider using fopen_s instead

printf( "The file 'test.wav' was not opened\n" );

else

printf( "The file 'test.wav' was opened\n" );

/*输出文件位置*/

char outfilename[30];

cout<<"请输入MFCC输出文件路径(需要加入后缀名*.txt)\n";

cin>>outfilename;

if( (err = fopen_s( &MFCCFile, outfilename, "wb" )) !=0 )

// Note: fopen is deprecated; consider using fopen_s instead

printf( "The file 'MFCCResult.txt' was not opened\n" );

else

printf( "The file 'MFCCResult.txt' was opened\n" );

//fread(&waveheader,sizeof(struct _TWavHeader),1,sourcefile); //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++)

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);

//char endl = '\n';

for(int j = 0; j < PCEP; j++)

{

if(j==1)

{

fprintf(MFCCFile, "%.1f,", Cep[j]);

temp_1++;

}

//if ( j == PCEP -1)

// {fprintf(MFCCFile, "%f,", Cep[j]);

//temp_1++;}

}

vecList.clear();

//Moves the file 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 << MFCCcoefficient[i] << ' ';

if( (i+1) % 13 == 0 )

cout << endl;

}

//getchar();

return 0;

}

//预加重

void preemphasis(double* buf, double* result, short FrmLen)

{

int i;

result[0] = buf[0] - SP_EMPHASIS_FACTOR * last;

for(i=1;i

{

result[i] = buf[i] - SP_EMPHASIS_FACTOR * buf[i-1];

}

last = buf[(FrmLen-1)/2]; //假设每次移半帧

}

//汉明窗初始化

void InitHamming()

{

double twopi;

int i;

twopi=2*PI;

for( i=0;i

{

Hamming[i]=(double)(0.54-0.46*cos(i*twopi/(double)(FrmLen-1)));

}

}

//给一帧数据加窗

void HammingWindow(double* result,double* data)

{

int i;

for(i=0;i

{

data[i]=result[i]*Hamming[i];

}

}

void compute_fft(double *data,vector >& vecList)

{

for(int i=0;i

{

if(i

{

complex temp(data[i]);

vecList.push_back(temp);

}

else

{

complex temp(0); //如果计算的FFT长度大于窗长，则不足的部分用零填充。得到的效果差不多

vecList.push_back(temp);

}

}

FFT(FFTLen,vecList);

}

void FFT(const unsigned long & ulN, vector >& vecList)

{

//得到指数，这个指数实际上指出了计算FFT时内部的循环次数

unsigned long ulPower = 0; //指数

unsigned long ulN1 = ulN - 1; //ulN1=511

while(ulN1 > 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 < ulPower; j++)

{

ulIndex += bsIndex.test(ulPower - j - 1) ? ulK : 0;

ulK *= 2;

}

if(ulIndex > p) //只有大于时，才调整，否则又调整回去了

{

complex c = vecList[p];

vecList[p] = vecList[ulIndex];

vecList[ulIndex] = c;

}

}

//计算旋转因子

vector > 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 ulHalfLength = 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 < ulHalfLength; k++)

{

cw = vecW[k * ulN / ulGroupLength] * vecList[j + k + ulHalfLength];

c1 = vecList[j + k] + cw;

c2 = vecList[j + k] - cw;

vecList[j + k] = c1;

vecList[j + k + ulHalfLength] = c2;

}

}

}

}

/*

设置滤波器参数

输入参数：无

输出参数：*FiltCoe1---三角形滤波器左边的系数

*FiltCoe2---三角形滤波器右边的系数

*Num ---决定每个点属于哪一个滤波器

*/

void InitFilt(double *FiltCoe1, double *FiltCoe2, int *Num)

{

int i,k;

double Freq;

double FiltFreq[FiltNum+2]; //40个滤波器，故有42各滤波器端点。每一个滤波器的左右端点分别是前一个及后一个滤波器的中心频率所在的点

double BW[FiltNum+1]; //带宽，即每个相邻端点之间的频率跨度

double low = (double)( 400.0 / 3.0 ); /* 滤波器组的最低频率，即第一个端点值*/

short lin = 13; /* 1000Hz以前的13个滤波器是线性的分布的*/

double lin_spacing = (double)( 200.0 / 3.0 ); /* 相邻滤波器中心的距离为66.6Hz */

short log = 27; /* 1000Hz以后是27个对数线性分布的滤波器*/

double log_spacing = 1.0711703f; /* 相邻滤波器左半边宽度的比值*/

for ( i=0; i

{

FiltFreq[i] = low + i * lin_spacing;

}

for ( i=lin; i

{

FiltFreq[i] = FiltFreq[lin-1] * (double)pow( log_spacing, i - lin + 1 );

}

for ( i=0; i

{

BW[i] = FiltFreq[i+1] - FiltFreq[i];

}

for(i = 0 ; i<= FFTLen/2 ; i++ )

{

Num[i] = 0;

}

bool bFindFilt = false;

for(i = 0 ; i <= FFTLen/2 ; i++)

{

Freq = FS * 1000.0F * i / (double)(FFTLen);

bFindFilt = false;

for(k = 0; k <= FiltNum; k++)

{

if(Freq >= FiltFreq[k] && Freq <= FiltFreq[k+1])

{

bFindFilt = true;

if(k == FiltNum)

{

FiltCoe1[i]=0.0F;

}

else

{

FiltCoe1[i] = (Freq - FiltFreq[k]) / (double)(BW[k]) * 2.0f / (BW[k] + BW[k+1]);

}

if(k == 0)

{

FiltCoe2[i] = 0.0F;

}

else

{

FiltCoe2[i] = (FiltFreq[k+1] - Freq) / (double)(BW[k]) * 2.0f / (BW[k] + BW[k-1]);

}

Num[i] = k; //当k==FiltNum时，它为第FiltNum个滤波器，实际上它并不存在。这里只是为了计算方便，假设有第FiltNum个滤波器存在。

//但其实这并不影响结果

break;

}

}

if (!bFindFilt)

{

Num[i] = 0; //这时，该点不属于任何滤波器，因为其左右系数皆为0，所以可以假定它属于某个滤波器，而不会影响结果。这里我

//将其设为第一个滤波器。

FiltCoe1[i]=0.0F;

FiltCoe2[i]=0.0F;

}

}

}

/*

根据滤波器参数计算频带能量

输入参数：*spdata ---预处理之后的一帧语音信号

*FiltCoe1---三角形滤波器左边的系数

*FiltCoe2---三角形滤波器右边的系数

*Num ---决定每个点属于哪一个滤波器

输出参数：*En ---输出对数频带能量

*/

//把属于某一频带的能量全部加起来了

//CFilt(data, FiltCoe1, FiltCoe2, Num, En,vecList); veclist ：FFT计算出的结果Num:决定每个点属于哪一个滤波器

void CFilt(double *spdata, double *FiltCoe1, double *FiltCoe2, int *Num, double

*En,vector >& vecList)

{

double temp=0;

int id, id1, id2;

for(id = 0; id < FiltNum ; id++)

{

En[id]=0.0F;

}

for(id = 0 ; id <= FFTLen/2 ; id++)

{

temp = vecList[id].real() * vecList[id].real() + vecList[id].imag() * vecList[id].imag();

temp=temp / ( (FrmLen/2) * (FrmLen/2) );

id1 = Num[id];

if (id1 == 0)

En[id1] = En[id1] + FiltCoe1[id] * temp;

if (id1 == FiltNum)

En[id1-1] = En[id1-1] + FiltCoe2[id] * temp;

if ((id1 > 0) && (id1 < FiltNum))

{

id2 = id1-1;

En[id1] = En[id1] + FiltCoe1[id] * temp;

En[id2] = En[id2] + FiltCoe2[id] * temp;

}

}

for(id = 0 ; id < FiltNum ; id++)

{

if (En[id] != 0)

En[id]=(double)log10(En[id]);

}

}

/*

计算MFCC系数

输入参数：*En ---对数频带能量

*/

void MFCC(double* En, double* Cep)

{

int idcep, iden;

// double Cep[13];

for(idcep = 0 ; idcep < PCEP ; idcep++)

{

Cep[idcep] = 0.0f;

for(iden = 0 ; iden < FiltNum ; iden++) //离散余弦变换

{

if(iden == 0)

Cep[idcep] = Cep[idcep] + En[iden] * (double)cos(idcep * (iden+0.5f) * PI/(FiltNum)) * 10.0f * sqrt(1/(double)FiltNum);

else

Cep[idcep] = Cep[idcep] + En[iden] * (double)cos(idcep * (iden+0.5f) * PI/(FiltNum)) * 10.0f * sqrt(2/(double)FiltNum);

}

MFCCcoefficient.push_back(Cep[idcep]);

}

}