Chinaunix首页 | 论坛 | 博客
  • 博客访问: 453769
  • 博文数量: 73
  • 博客积分: 1170
  • 博客等级: 少尉
  • 技术积分: 720
  • 用 户 组: 普通用户
  • 注册时间: 2012-02-20 15:48
文章分类

全部博文(73)

文章存档

2013年(9)

2012年(64)

我的朋友

分类: 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].image = X1[i + p].image + X1[i + p + bfsize / 2].image * W[i * (1<    
    X2[i + p + bfsize / 2].real = X1[i + p].real - X1[i + p + bfsize / 2].real * W[i * (1<    X2[i + p + bfsize / 2].image = X1[i + p].image - X1[i + p + bfsize / 2].image * 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特征量的总长度
}

阅读(8986) | 评论(2) | 转发(1) |
给主人留下些什么吧!~~

foreverhst2014-10-31 14:36:43

Access Violation

foreverhst2014-10-31 14:34:43

您好,程序为什么会出现栈溢出问题