分类： C/C++
2012-06-28 12:39:01
//ASR.h
#include
#include
#define MAXDATA (256*400) //一般采样数据大小,语音文件的数据不能大于该数据
#define SFREMQ (8000) //采样数据的采样频率8khz
typedef struct WaveStruck//wav数据结构
{
//data head
struct HEAD{
char cRiffFlag[4];
int nFileLen;
char cWaveFlag[4];//WAV文件标志
char cFmtFlag[4];
int cTransition;
short nFormatTag;
short nChannels;
int nSamplesPerSec;//采样频率,mfcc为8khz
int nAvgBytesperSec;
short nBlockAlign;
short nBitNumPerSample;//样本数据位数，mfcc为12bit
} head;
//data block
struct BLOCK{
char cDataFlag[4];//数据标志符(data)
int nAudioLength;//采样数据总数
} block;
} WAVE;
typedef struct //定义实数结构
{
double real;
double image;
}complex;
//获取wav数据样本数据,sample数组的长度至少为MAXDATA,函数返回读入的数据的长度
int getWaveData(char* file,double *sample);
int copyData(double *src,double *dest,int num);//实数数据复制
//所有mfcc操作都是基于每帧操作的,len为总长度,unitsize为采样点组成的观察单位,256||512
//emp_v为欲加重因子,hanming_v为汉明窗因子,fft_v为傅里叶变化级数
//filter_v为mel滤波器个数,dct_v为dct变换因子
//函数返回mfcc识别后特征参数数组，数组长度存储在第0个元素中
double* MFCC(double *sample,int len ,int unitsize=256,double emp_v=0.97,double hanming_v=0.46,int filter_v=20,int dct=12);
//MFCC欲加重,len为数组sample的长度，factor为加重因子
bool _mfcc_preEmphasize(double *sample,int len,double factor=0.9);
//框化数据即分帧，同事在原始采样数据中插入重叠区数据
//函数返回总的帧数,和分帧后的数据sample
int _mfcc_Frame(double *sample,int unitsize,int len);
//加汉明窗,factor为汉明窗因子
bool _mfcc_HanmingWindow(double *sample,int unitsize,int frames,double hanming_factor=0.46);
//快速傅里叶变换，输入参数sample为实数，无虚数部分,len为sample的总大小，也为傅里叶变换点数
//函数返回时complex为fft变换后的频率能量复数数值,总长度为len,频率从0-SFREMQ(8khz)
complex* _mfcc_FFT(double *sample,int len);
//功率计算,返回计算后的功率实数powersample，总长度为len,频率从0-SFREMQ(8khz)
bool _mfcc_Power(complex* fftdata,double *powersample,int len);
//三角滤波,len为powersample数字长度,num为mei滤波个数
//函数返回mel频率下的滤波器内的对数能量double*(长度为num)
double* _mfcc_filter(double *powersample,int len,int num);
//离散余弦变换,lnpower(长度为len)为滤波器内积对数能量,num与mel滤波器一致,dctnum为dct变换级数
//函数返回时double*的长度为dctnum,代表dct变换后的倒频参数
double* _mfcc_DCT(double *lnpower,int len,int dctnum);
//求音框1的能量+音框2的能量+....+音框frames的能量+dct参数
//函数返回double *,其长度为(dctlen+1)*frames
double* _mfc_Logenergy(double *powersample,int unitsize,int frames,double * dctdata,int dctlen);
//差量倒频谱参数,函数返回double*，数组长度为len*2
double* _mfcc_differential(double *logenergy,int len);
//asp.cpp
#include "stdafx.h"
#include "ASR.h"
//获取wav数据样本数据
int getWaveData(char* file,double *sample)
{
WAVE wave[1];
FILE * f;
f = fopen(file,"rb");
if(!f)
{
printf("Cannot open %s for reading\n",file);
return -1;
}
//读取wav文件头并且分析
fread(wave,1,sizeof(wave),f);
if(wave[0].head.cWaveFlag[0]=='W'&&wave[0].head.cWaveFlag[1]=='A'
&&wave[0].head.cWaveFlag[2]=='V'&&wave[0].head.cWaveFlag[3]=='E')//判断是否是wav文件
{
printf("It's not .wav file\n" );
return -1;
}
if(wave[0].head.nSamplesPerSec!=SFREMQ&&wave[1].head.nBitNumPerSample!=12)//判断是否采样频率是8khz,12bit量化
{
printf("It's not 8khz and 12 bit\n" );
return -1;
}
if(wave[0].block.nAudioLength>MAXDATA/2)//wav文件不能太大,为sample长度的一半
{
printf("wav file is to long\n" );
return -1;
}
//读取采样数据
fread(sample,sizeof(char),wave[0].block.nAudioLength,f);
fclose(f);
return wave[0].block.nAudioLength;
}
int copyData(double *src,double *dest,int num)
{
if(src==NULL||dest==NULL)
{
printf("src and dest are NULL\n");
return -1;
}
for(int i=0;i
*(dest++)=*(src++);
}
return i;
}
//开始mfc操作
//所有mfcc操作都是基于每帧操作的,unit为采样点组成的观察单位,256||512
//函数返回double*长度存储在第0个元素中
double* MFCC(double *sample,int len,int unitsize,double emp_v,double hanming_v,int filter_v,int dct_v)
{
double *mfccA=NULL;
//欲加重
_mfcc_preEmphasize(sample,len,emp_v);
//框化，在原始采样数据中添加重叠区数据
int frames=_mfcc_Frame(sample,unitsize,len);
//加汉明窗
_mfcc_HanmingWindow(sample,unitsize,frames,hanming_v);
//fft变化
complex *fftdata=_mfcc_FFT(sample,unitsize*frames);
//频谱的功率sample
_mfcc_Power(fftdata,sample,len);
if(fftdata!=NULL)
delete fftdata;
//保存频谱能量
double *power=new double(len);
copyData(sample,power,len);
//滤波，得到filter_v个滤波器内的对数能量
double *filterpower=_mfcc_filter(sample,len,filter_v);
//dct变换,输入filterpower为滤波器内积对数能量，输出lnpower为dct倒频谱参数
double *dctpara=_mfcc_DCT(filterpower,filter_v,dct_v);
if(filterpower!=NULL)
delete filterpower;
//log对数能量logenergy,长度为frames+dct_v
double *logenergy=_mfc_Logenergy(power,unitsize,frames,dctpara,dct_v);
//保存对数能量logenergy
mfccA=new double((frames+dct_v)*3+1);
mfccA[0]=(frames+dct_v)*3+1;//第0位存储数组的总长度
copyData(&logenergy[0],&mfccA[1],frames+dct_v);
if(power!=NULL)
delete power;
if(dctpara!=NULL)
delete dctpara;
//差量倒频谱运算,返回double*长度为(frames+dct_v)*2
double *diff=_mfcc_differential(logenergy,frames+dct_v);
copyData(&diff[0],&mfccA[frames+dct_v+1],(frames+dct_v)*2);
if(logenergy!=NULL)
delete logenergy;
if(diff!=NULL)
delete diff;
return mfccA;
}
bool _mfcc_preEmphasize(double *sample,int len,double factor)//MFCC欲加重,factor为加重因子
{
//欲加重公式s2(n) = s(n) - a*s(n-1)
for(int i=1;i
sample[i]=sample[i]-factor*sample[i-1];
}
return true;
}
//数据框化即分帧，插入重叠区的数据到sample中，返回帧的总数
int _mfcc_Frame(double *sample,int unitsize,int len)
{
double *mirrordata=new double[len/2];
copyData(&sample[0],&mirrordata[0],len/2);
double temp[512];
if(unitsize>512)
{
printf("unitsize of mfcc is biger than 512");
return false;
}
//因为要操作附加的重叠区域，故分成三部分,第一部分数据不变
//开始和结尾部分不包含附加重叠区的操作
//中间部分，包括附件重叠区域变换
int m=unitsize/2;//重叠区域的单位大小为正常区的一半
int frames=2;//加上重叠区后的帧数
int index=unitsize;//当前添加后的数据的排列序号
for(int i=unitsize;i
//重叠区域数据
copyData(&mirrordata[i-m/2],&temp[0],m);
//添加重叠区的数据
copyData(&temp[0],&sample[index],m);//写回数据
index+=m;
frames+=1;
//原采样数据
copyData(&mirrordata[i],&temp[0],unitsize);
copyData(&temp[0],&sample[index],unitsize);//写回数据
index+=unitsize;
frames+=1;
}
//框化的最后一部分为尾数部分,少于或者等于unitsize大小
if(len%unitsize==0)//刚好等于unitsize大小
{
copyData(&mirrordata[i],&temp[0],unitsize);
copyData(&temp[0],&sample[index],unitsize);//写回数据
}else{
copyData(&mirrordata[i],&sample[index],len%unitsize);
copyData(&temp[0],&sample[index],len%unitsize);//写回数据
}
frames+=1;
delete[] mirrordata;
return frames;
}
//加汉明窗,factor为汉明窗因子
bool _mfcc_HanmingWindow(double *sample,int unitsize,int frames,double hanming_factor)
{
//S'(n) = S(n)*W(n)，W(n, a) = (1 - a) - a cos(2pn/(N-1))，0≦n≦N-1，N为窗口大小
if(hanming_factor>=1||hanming_factor<=0)
{
printf("hanming hanming_factor 0-1\n");
return false;
}
if(unitsize>512)
{
printf("unitsize of mfcc is biger than 512");
return false;
}
int windowsize;//框化后的窗口大小
double HanmingWindow[512];//汉明窗口乘法因子
double fDataTemp;
for(int i=0;i
if((i+1)%2==0)
{
windowsize=unitsize/2;//重叠区为正常区的一半
}else
windowsize=unitsize;
for(int j=0; j
fDataTemp = (double) 2.0 * 3.141592653589793*j/(windowsize-1);
HanmingWindow[j] = (double)(1-hanming_factor)-hanming_factor*cos(fDataTemp);
sample[i*frames+j]=sample[i*frames+j]*HanmingWindow[j];
}
}
return true;
}
//快速傅里叶变换，输入参数sample为实数，无虚数部分,len为sample的总大小，也为傅里叶变换点数
//函数返回时complex为fft变换后的频率能量复数数值,总长度为len,频率从0-SFREMQ(8khz)
complex* _mfcc_FFT(double *sample,int len)
{
//实数转为复数
complex *TD=new complex[len];
complex *FD=new complex[len];
int i;
for(i=0;i
TD[i].real=sample[i];
TD[i].image=0;
}
//fft计算
LONG count; // Fourier变换点数
int j,k; // 循环变量
int bfsize,p; // 中间变量
double angle; // 角度
complex *W,*X1,*X2,*X;
count =len; // 计算Fourier变换点数为数组的长度
int r=count>>1;//变换级数
W = new complex[count / 2];
X1 = new complex[count];
X2 = new complex[count]; // 分配运算所需存储器
// 计算加权系数（旋转因子w的i次幂表)
for(i = 0; i < count / 2; i++)
{
angle = -i * 3.141592653589793 * 2 / count;
W[ i ].real=cos(angle);
W[ i ].image=sin(angle);
}
// 将时域点写入X1
memcpy(X1, TD, sizeof(complex) * count);
// 采用蝶形算法进行快速Fourier变换
for(k = 0; k < r; k++ )
{
for(j = 0; j < (1 << k); j++ )
{
bfsize = 1 << (r-k);
for(i = 0; i < bfsize / 2; i++ )
{
p = j * bfsize;
X2[i + p].real = X1[i + p].real + X1[i + p + bfsize / 2].real * W[i * (1<
X2[i + p + bfsize / 2].real = X1[i + p].real - X1[i + p + bfsize / 2].real * W[i * (1<
}
X = X1;
X1 = X2;
X2 = X;
}
// 重新排序，转为频率f(0-count)的值
for(j = 0; j < count; j++)
{
p = 0;
for(i = 0; i < r; i++ )
{
if (j&(1<
p =1<<(r-i-1);
}
}
FD[j]=X1[p];
}
// 释放内存
delete[] TD;
delete W;
delete X1;
delete X2;
return FD;
}
//对应频谱的功率
bool _mfcc_Power(complex* fftdata,double *powersample,int len)
{
for(int i=0;i
powersample[i]=sqrt(fftdata[i].real*fftdata[i].real+fftdata[i].image*fftdata[i].image);
}
return true;
}
//三角滤波,len为powersample数字长度,num为mei滤波个数
//函数返回mel频率下的对数能量double*(长度为num)
double* _mfcc_filter(double *powersample,int len,int filterNum)
{
double* melf=new double(len);//频率--->mel频率
double* lnpower=new double(filterNum);//对数能量
memset(lnpower,0,filterNum);
double f=0;
for(int i=0;i
f=i*SFREMQ/len;//求出每个数组的相对应的频率
melf[i]=2595*log10(1+f/700.0);//由频率求出mel频率
}
double filterFW=melf[len-1]/(2*filterNum);//mel频率滤波器的半宽带
double filterpara;//三角滤波参数
int j=1;//滤波器计数
for(i=0;i
j=(int)ceil(melf[i]/filterFW);//ceil大于x的最小整数,求出滤波器序号为1-n
if(j>filterNum)//跳过最右边半个滤波器数据
break;
filterpara=1+(melf[i]-j*filterFW)/filterFW; //斜率为1，上升
powersample[i]=filterpara*powersample[i];
lnpower[j-1]+=powersample[i];//累计一个滤波器内的能量
}
for(i=0;i
if(melf[i]
continue;
}
j=(int)floor(melf[i]/filterFW);//floor小于x的最大整数
filterpara=1-(melf[i]-j*filterFW)/filterFW; //斜率为1，下降
powersample[i]=filterpara*powersample[i];
lnpower[j-1]+=powersample[i];//累计一个滤波器内的能量
}
//滤波器能量求对数,得到类似正态分布
//对数能量（定义为讯号的平方和，再取以 10 为底的对数值，再乘以 10）
for(j=0;j
lnpower[j]=log10(lnpower[j])*10;
}
delete[] melf;
return lnpower;//返回对数能量
}
//离散余弦变换,lnpower(长度为len)为滤波器内积对数能量,num与mel滤波器一致,dctnum为dct变换级数
//函数返回时lnpower的长度为dctnum,代表dct变换后的倒频参数
double* _mfcc_DCT(double *lnpower,int len,int dctnum)
{
double *temp=new double(dctnum);
memset(&lnpower[0],0,len);
for(int i=1;i<=dctnum;i++)
{
for(int j=1;j<=len;j++)
{
temp[i-1]+=lnpower[j-1]*cos(i*(j-0.5)*3.141592653589793/len);
}
}
return temp;
}
//求音框1的能量+音框2的能量+....+音框frames的能量+dct参数
//函数返回double *,其长度为(dctlen+frames)
double* _mfc_Logenergy(double *powersample,int unitsize,int frames,double * dctdata,int dctlen)
{
double *logenergy=new double((frames+dctlen));
for(int i=0;i
logenergy[i]=0;
for(int j=0;j
logenergy[i]+=powersample[i*unitsize+j];
}
logenergy[i]=log10(logenergy[0])*10;//取10的对数
}
copyData(&dctdata[0],&logenergy[i],dctlen);
return logenergy;
}
//差量倒频谱参数,函数返回double*，数组长度为len*2
double* _mfcc_differential(double *logenergy,int len)
{
int M=2;//差量倒频谱参数
double suma;
double sumb;
double *diff=new double(len*2);
//差量计算,△Cm(t) = [Sum(g=-M, M, Cm(t+g)g] / [Sum(g=-M, M, g^2], 这里 M 的值一般是取 2 或 3
copyData(&logenergy[0],&diff[0],len);
int i,j;
for(i=1;i<=len;i++)
{
suma=0;
sumb=0;
for(j=-M;j<=M;j++)
{
suma+=diff[i-1+j]*j;
sumb+=j*j;
}
diff[i-1]=suma/sumb;
}
//差差量计算
copyData(&diff[0],&diff[len],len);
for(i=len+1;i<=2*len;i++)
{
suma=0;
sumb=0;
for(j=-M;j<=M;j++)
{
suma+=diff[i-1+j]*j;
sumb+=j*j;
}
diff[i-1]=suma/sumb;
}
return diff;
}
void main()
{
double sample[MAXDATA]={0};//wave文件采样数据存放点
int len=getWaveData("1.wav",sample);
double *mfcca=MFCC(sample,len);
len=mfcca[0];//mfcc特征量的总长度
}
