开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 多媒体 > 语音合成与识别 > 查看源码
FE_spectrum.cpp
资源名称：speech_analysis.rar　[点击查看]
上传用户：italyroyal
上传日期：2013-05-06
资源大小：473k
文件大小：10k
源码类别：

语音合成与识别

开发平台：

Visual C++

FE_spectrum.cpp：源码内容
  1. ///////////////////////////////////////////////////////////////////////////////
  2. // This is a part of the Feature program.
  3. // Version: 1.0
  4. // Date: February 22, 2003
  5. // Programmer: Oh-Wook Kwon
  6. // Copyright(c) 2003 Oh-Wook Kwon. All rights reserved. owkwon@ucsd.edu
  7. ///////////////////////////////////////////////////////////////////////////////
  9. #include "StdAfx.h" #include "FE_feature.h"
  11. int Fe::fft_spectrum_basic(short *sample, int frameSize, float *spectrum, int fftSize)
  12. {
  13. _fft_spectrum_basic(sample, frameSize, spectrum, fftSize, 0, 0);
  14. return 1;
  15. }
  18. int Fe::_fft_spectrum_basic(short *sample, int frameSize, float *spectrum, int fftSize, int cep_smooth, int cepFilterLen)
  19. {
  20. int i;
  21. vector<float> frameA(frameSize);
  22. vector<float> power(fftSize);
  24. for(i=0;i<frameSize;i++) frameA[i]=sample[i];
  25. preemphasize(&frameA[0], frameSize, m_emphFac);
  26. m_window.Windowing(&frameA[0], frameSize, WIN_HAMMING);
  27. int uiLogFFTSize = (int)(log((double)fftSize)/log((double)2)+0.5);
  29. compute_spectrum(&frameA[0], &power[0], frameSize, uiLogFFTSize);
  31. float log10db=(float)10/log(10);
  32. for(i=0;i<(int)fftSize;i++)
  33. spectrum[i] = log10db*LogE(power[i]);
  35. if(cep_smooth) {
  36. vector<float> a(fftSize);
  37. vector<float> b(fftSize);
  39. for(i=0;i<fftSize;i++){
  40. a[i] = spectrum[i];
  41. b[i] = 0;
  42. }
  43. PRFFT_NEW( &a[0], &b[0], uiLogFFTSize, fftSize, -1);
  44. if(cepFilterLen <= 0) cepFilterLen = 1;
  45. for (i=cepFilterLen ; i < fftSize-(cepFilterLen-1); i++){
  46. a[i] = b[i] = 0;
  47. }
  48. for(i=fftSize-1;i>fftSize/2;i--){
  49. a[i] = a[fftSize-i];
  50. b[i] = b[fftSize-i];
  51. }
  52.      PRFFT_NEW( &a[0], &b[0], uiLogFFTSize, fftSize, 1) ;
  53. for (i=0 ; i < m_fftSize; i++) 
  54. spectrum[i] = a[i];
  55. }
  56. return 1;
  57. }
  60. int Fe::fft_cepstrum_basic(short *sample, int frameSize, float *fft_cep, int ceporder, int fftSize)
  61. {
  62. vector<float> frameA(frameSize);
  63. preprocessing(sample, frameSize, &frameA[0]);
  64. m_window.Windowing(&frameA[0], frameSize, WIN_HAMMING);
  65. _fft_cepstrum_basic(&frameA[0],frameSize,fft_cep,ceporder,fftSize);
  66. cepstral_window(fft_cep, ceporder, m_lifter);
  68. /* owkwon: To make dynamic range consistent with other kinds of cepstral coefficients */
  69. /* Needs further study */
  70. for (int i=0;i<=ceporder;i++) fft_cep[i] /= 8;
  71. return 1;
  72. }
  75. int Fe::_fft_cepstrum_basic(float *sample, int frameSize, float *fft_cep, int ceporder, int fftSize)
  76. {
  77. int i;
  78. int uiLogFFTSize = (int)(log((double)fftSize)/log((double)2)+0.5);
  79. vector<float> a(fftSize);
  80. vector<float> b(fftSize);
  81. for(i=0;i<frameSize;i++){ a[i] = sample[i]; b[i]=0; }
  82. for(i=frameSize;i<fftSize;i++) a[i] = b[i] = 0;
  83. PRFFT_NEW( &a[0], &b[0], uiLogFFTSize, frameSize, 1);
  84. for(i=0;i<fftSize;i++){
  85. a[i] = (float)(0.5*LogE(a[i]*a[i]+b[i]*b[i]));
  86. b[i] = 0;
  87. }
  88. PRFFT_NEW( &a[0], &b[0], uiLogFFTSize, fftSize, -1);
  89. for(i=0;i<=ceporder;i++) fft_cep[i] = a[i];
  91. return 1;
  92. }
  95. int Fe::filterbank_basic(short *sample, int frameSize, float *filter_bank, int fborder, int fftSize)
  96. {
  97. vector<float> frameA(frameSize);
  98. preprocessing(sample, frameSize, &frameA[0]);
  99. m_window.Windowing(&frameA[0], frameSize, WIN_HAMMING);
  100. _filterbank_basic(&frameA[0], frameSize, filter_bank, fborder, fftSize, m_cepSmooth, (int)(m_sampleRate/MAX_PITCH_FREQ));
  101. /* Convert to dB scale */
  102. float log10db=(float)(10/log(10));
  103. for(int i=0;i<fborder;i++) filter_bank[i]*=log10db;
  104. return 1;
  105. }
  108. int Fe::_filterbank_basic(float *sample, int frameSize, float *filter_bank, int fborder, int fftSize, int cep_smooth, int cepFilterLen)
  109. {
  110. if(m_MelFB.size() != fborder || fftSize != m_MelFBfftSize){
  111. m_MelFBfftSize=fftSize;
  112. MfccInitMelFilterBanks ((float)64.0, (float)m_sampleRate, fftSize, fborder);
  113. }
  115. int i;
  116. int uiLogFFTSize = (int)(log((double)fftSize)/log((double)2)+0.5);
  117. vector<float> a(fftSize);
  118. vector<float> b(fftSize);
  119. for(i=0;i<frameSize;i++){ a[i] = sample[i]; b[i]=0; }
  120. for(i=frameSize;i<fftSize;i++) a[i] = b[i] = 0;
  121. PRFFT_NEW( &a[0], &b[0], uiLogFFTSize, frameSize, 1);
  122. for(i=0;i<fftSize;i++){
  123. a[i] = a[i]*a[i] + b[i]*b[i];
  124. b[i] = 0.0;
  125. }
  126. MfccMelFilterBank (&a[0], fborder, &filter_bank[0], 1);
  127. for (i = 0; i < fborder; i++) filter_bank[i] = 0.5*LogE(filter_bank[i]);
  129. return(1);
  130. }
  133. int Fe::compute_spectrum(float *input, float *spectrum, int winlength, int log2length) 
  134. {
  135.   int i, log2length_i;
  136.   int npoints, npoints2;
  138.   npoints  = (int) (pow(2.0,(float) log2length) + 0.5);
  139.   npoints2 = npoints / 2;
  141.   vector<FeComplex> complex_in(npoints);
  142.   for(i=0;i<winlength;i++){
  143.     complex_in[i].m_re = input[i];
  144.     complex_in[i].m_im = (float)0.0;
  145.   }
  146.   for(i = winlength; i < npoints; i++){
  147.     complex_in[i].m_re = (float)0.0;
  148.     complex_in[i].m_im = (float)0.0;
  149.   }
  151.   log2length_i = (int)log2length;
  152.   FAST_new(&complex_in[0],log2length_i);
  153.   
  154.   spectrum[0] = complex_in[0].m_re * complex_in[0].m_re;
  155.   spectrum[npoints2] = complex_in[npoints2].m_re * complex_in[npoints2].m_re;
  157.   /*
  158.     Only the first half of the spectrum[] array is filled with data.
  159.     The second half is just a reflection of the first half.
  160.   */
  161.   for(i=1;i<npoints2;i++){
  162.     spectrum[i] = complex_in[i].m_re * complex_in[i].m_re + complex_in[i].m_im * complex_in[i].m_im;
  163.   }
  164.   for(i=npoints-1;i>npoints2;i--){
  165.   spectrum[i] = spectrum[npoints-i];
  166.   }
  168.   return(0);
  169. }
  172. int Fe::smooth_spectrum(float *spectrum, int pointsN)
  173. {
  174. vector<float> tmp(pointsN);
  175. tmp[0]=spectrum[0];
  176. tmp[1]=(2*spectrum[1]+spectrum[0]+spectrum[2])/(float)4;
  177. for(int i=2;i<pointsN-2;i++){
  178. tmp[i]=(float)(4*spectrum[i]+2*spectrum[i-1]+2*spectrum[i+1]+spectrum[i-2]+spectrum[i+2])/(float)10;
  179. }
  180. tmp[pointsN-2]=(2*spectrum[pointsN-2]+spectrum[pointsN-3]+spectrum[pointsN-1])/(float)4;
  181. tmp[pointsN-1]=spectrum[pointsN-1];
  183. memmove(spectrum, &tmp[0], sizeof(float)*pointsN);
  184. return(1);
  185. }
  188. /*========================================================================
  189. ** Discrete Fourier analysis routine
  190. ** from IEEE Programs for Digital Signal Processing
  191. ** G. D. Bergland and M. T. Dolan, original authors
  192. ** Translated from the FORTRAN with some changes by Paul Kube
  193. ** Modified to return the power spectrum by Chuck Wooters
  194. ========================================================================*/
  196. /**************************************************************************
  197. FAST_new - In-place radix 2 decimation in time FFT
  198. Requires pointer to complex array, x and power of 2 size of FFT, m
  199. (size of FFT = 2**m).  Places FFT output on top of input FeComplex array.
  200. void FAST_new(FeComplex *x, int m)
  201. *************************************************************************/
  203. void Fe::FAST_new(FeComplex *x, int m)
  204. {
  205.     FeComplex u,temp,tm;
  206.     FeComplex *xi,*xip,*xj,*wptr;
  207.     int i,j,k,l,le,windex;
  208.     double arg,w_real,w_imag,wrecur_real,wrecur_imag,wtemp_real;
  210. if(m == 0) return;       /* if m=0 then done */
  211. int n = 1 << m; /* length of fft */
  212.     if(n != m_fftW.size()) {
  213. m_fftW.resize(n); 
  214.         le = n/2;
  216. for(int nn=0;nn < le-1; nn++)
  217. m_fftW[nn].m_re = m_fftW[nn].m_im = (float)0.0;
  219.         arg = 4.0*atan(1.0)/le;
  220.         wrecur_real = w_real = cos(arg);
  221.         wrecur_imag = w_imag = -sin(arg);
  222.         xj = &m_fftW[0];
  223.         for (j = 1 ; j < le ; j++) {
  224.             xj->m_re = (float)wrecur_real;
  225.             xj->m_im = (float)wrecur_imag;
  226.             xj++;
  227.             wtemp_real = wrecur_real*w_real - wrecur_imag*w_imag;
  228.             wrecur_imag = wrecur_real*w_imag + wrecur_imag*w_real;
  229.             wrecur_real = wtemp_real;
  230.         }
  231.     }
  233.     le = n;
  234.     windex = 1;
  235.     for (l = 0 ; l < m ; l++) {
  236.         le = le/2;
  238.         for(i = 0 ; i < n ; i = i + 2*le) {
  239.             xi = x + i;
  240.             xip = xi + le;
  241.             temp.m_re = xi->m_re + xip->m_re;
  242.             temp.m_im = xi->m_im + xip->m_im;
  243.             xip->m_re = xi->m_re - xip->m_re;
  244.             xip->m_im = xi->m_im - xip->m_im;
  245.             *xi = temp;
  246.         }
  248.         wptr = &m_fftW[0] + windex - 1;
  249.         for (j = 1 ; j < le ; j++) {
  250.             u = *wptr;
  251.             for (i = j ; i < n ; i = i + 2*le) {
  252.                 xi = x + i;
  253.                 xip = xi + le;
  254.                 temp.m_re = xi->m_re + xip->m_re;
  255.                 temp.m_im = xi->m_im + xip->m_im;
  256.                 tm.m_re = xi->m_re - xip->m_re;
  257.                 tm.m_im = xi->m_im - xip->m_im;
  258.                 xip->m_re = tm.m_re*u.m_re - tm.m_im*u.m_im;
  259.                 xip->m_im = tm.m_re*u.m_im + tm.m_im*u.m_re;
  260.                 *xi = temp;
  261.             }
  262.             wptr = wptr + windex;
  263.         }
  264.         windex = 2*windex;
  265.     }
  267.     j = 0;
  268.     for (i = 1 ; i < (n-1) ; i++) {
  269.         k = n/2;
  270.         while(k <= j) {
  271.             j = j - k;
  272.             k = k/2;
  273.         }
  274.         j = j + k;
  275.         if (i < j) {
  276.             xi = x + i;
  277.             xj = x + j;
  278.             temp = *xj;
  280.             *xj = *xi;
  281.             *xi = temp;
  282.         }
  283.     }
  284. }
  287. /* pruned FFT */
  288. void PRFFT_NEW(float *a, float *b, int m, int n_pts, int iff)
  289. {
  290. int l,lmx,lix,lm,li,j1,j2,ii,jj,nv2,nm1,k,lmxy,n;
  291. float scl,s,c,arg,T_a,T_b;
  293. n = 1 << m;
  294. for( l=1 ; l<=m ; l++ ) {
  295. lmx = 1 << (m - l) ;
  296. lix = 2 * lmx ;
  297. scl = 2 * (float)M_PI/(float)lix ;
  299. if(lmx < n_pts) lmxy = lmx ;
  300. else lmxy = n_pts ;
  301. for( lm = 1 ; lm <= lmxy ; lm++ ) {
  302. arg=((float)(lm-1)*scl) ;
  303. c = (float)cos(arg) ;
  304. s = iff * (float)sin(arg) ;
  306. for( li = lix ; li <= n ; li += lix ) {
  307. j1 = li - lix + lm ;
  308. j2 = j1 + lmx ;
  309. if(lmxy != n_pts ) {
  310. T_a=a[j1-1] - a[j2-1] ;
  311. /* T_a : real part */
  312. T_b=b[j1-1] - b[j2-1] ;
  313. /* T_b : imaginary part */
  314. a[j1-1] = a[j1-1] + a[j2-1] ;
  315. b[j1-1] = b[j1-1] + b[j2-1] ;
  316. a[j2-1] = T_a*c + T_b*s ;
  317. b[j2-1] = T_b*c - T_a*s ;
  318. } else{
  319. a[j2-1] = a[j1-1]*c + b[j1-1]*s ;
  320. b[j2-1] = b[j1-1]*c - a[j1-1]*s ;
  321. }
  322. }
  323. }
  324. }
  325. nv2 = n/2 ;
  326. nm1 = n - 1 ;
  327. jj = 1 ;
  329. for( ii = 1 ; ii <= nm1 ;) {
  330. if( ii < jj ) {
  331. T_a = a[jj-1] ;
  332. T_b = b[jj-1] ;
  333. a[jj-1] = a[ii-1] ;
  334. b[jj-1] = b[ii-1] ;
  335. a[ii-1] = T_a ;
  336. b[ii-1] = T_b ;
  337. }
  338. k = nv2 ;
  339. while( k < jj ) {
  340. jj = jj - k ;
  341. k = k/2 ;
  342. }
  343. jj = jj + k ;
  344. ii = ii + 1 ;
  345. }
  346. if(iff == (-1)){
  347. for( l=0 ; l<n ; l++ ) {
  348. a[l]/=(float)n;
  349. b[l]/=(float)n;
  350. }
  351. }
  352. }