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

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

语音合成与识别

开发平台：

Visual C++

FE_enhance.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"
  10. #include <stdio.h>
  11. #include <stdarg.h>
  12. #include <stdlib.h>
  13. #include <string.h>
  14. #include <math.h>
  15. #include <assert.h>
  16. #include "FE_feature.h"
  19. /*------------------------------------------
  20. * Definition of Noise Reduction Parameters 
  21. *-----------------------------------------*/
  23. /* general */
  24. #ifndef DEFAULT_SAMPLING_RATE
  25. #define DEFAULT_SAMPLING_RATE        16000
  26. #endif
  28. /* Wiener filter design */
  29. #define NR_NB_FRAME_THRESHOLD_NSE      100
  30. #define NR_LAMBDA_NSE                 0.99
  31. #define NR_EPS                       1e-16 /* owkwon: originally 1e-16 */
  32. #define NR_BETA                       0.98
  33. #define NR_ETA_TH                     0.01 /* owkwon: reduced to -40 dB for kWaves, originally 0.079432823 corresponding to -22 dB */
  34. #define NR_STARTING_FREQ                 0   /* 0Hz  */
  36. /* Spectral subtraction design */
  37. #define NR_SS_TH                      0.01
  39. /* VAD for noise estimation (VADNest) */
  40. #define NR_NB_FRAME_THRESHOLD_LTE       10
  41. #define NR_LAMBDA_LTE                 0.97
  42. #define NR_SNR_THRESHOLD_UPD_LTE        20
  43. #define NR_ENERGY_FLOOR                 80
  44. #define NR_MIN_FRAME                    10
  45. #define NR_SNR_THRESHOLD_VAD            24 /* owkwon: increased, originally 15 */
  46. #define NR_MIN_SPEECH_FRAME_HANGOVER     4
  47. #define NR_HANGOVER                      5
  49. #define NR_SNR_THRESHOLD_NO_UPD_LTE     80 /* owkwon: there sometimes exist only a few silence frames less than NR_MIN_FRAME */
  52. /*----------------------
  53.  * Constant Definitions
  54.  *----------------------*/
  55. #ifndef M_PI
  56. #define M_PI        ((float)3.14159265358979323846)
  57. #endif
  60. /*----------------------
  61.  * Macro Definitions
  62.  *----------------------*/
  63. #define NR_T(z)          ((z)>0 ? (z) : 0)
  64. #define NR_SQ(z)                 ((z)*(z))
  65. #define Round(x)    ((int)((x)+0.5))
  68. /*---------------------
  69.  * Function Prototypes
  70.  *---------------------*/
  71. extern void PRFFT_NEW(float *a, float *b, int m, int n_pts, int iff);
  74. int Fe::EnhanceMain(const char *infile, const char *outfile, int sampleRate, int isWiener)
  75. {
  76. FILE *fpin=fopen(infile,"rb");
  77. if(fpin==0){
  78. fprintf(stderr, "File open error (%s)n",infile);
  79. assert(0);
  80. return 0;
  81. }
  82. fseek(fpin,0L,SEEK_END);
  83. int fsize = ftell(fpin);
  84. rewind(fpin);
  85. int sampleN=fsize/sizeof(short);
  86. vector<short> data(sampleN);
  87. fread(&data[0], sizeof(short), sampleN, fpin);
  88. fclose(fpin);
  90. enhance_basic(&data[0], sampleN, sampleRate, isWiener);
  92. FILE *fpout=fopen(outfile,"wb");
  93. if(fpout==0){
  94. fprintf(stderr, "File open error (%s)n",outfile);
  95. assert(0);
  96. return 0;
  97. }
  98. fwrite(&data[0], sizeof(short), sampleN, fpout);
  99. fclose(fpout);
  100. return 1;
  101. }
  104. void Fe::enhance_basic(short *sample, int sampleN, int samplingRate, int isWiener)
  105. {
  106. int t;
  107. Wiener wiener;
  108. wiener.Init(samplingRate, isWiener);
  109. float out[NR_MAX_WIN_SIZE];
  110. wiener.InitNewUtterance(0);
  111. int frameN=sampleN/wiener.m_shiftSize;
  112. int outN=0;
  113. for(t=0; ;t++){
  114. int i;
  115. int n=my_min(wiener.m_shiftSize,my_max(0,sampleN-t*wiener.m_shiftSize));
  116. FeReturnCode status;
  117. if(n>=wiener.m_shiftSize){
  118. status=wiener.OneFrame(sample+t*wiener.m_shiftSize, n, &out[0], t);
  119. }
  120. else{ /* feed in dummy frames */
  121. short tmpA[NR_MAX_WIN_SIZE];
  122. for(i=0;i<n;i++) tmpA[i]=sample[t*wiener.m_shiftSize+i];
  123. for(i=n;i<wiener.m_shiftSize;i++) tmpA[i]=0;
  124. status=wiener.OneFrame(tmpA, wiener.m_shiftSize, &out[0], t);
  125. }
  126. if(status==FE_NULL) continue;
  127. int r=my_min(wiener.m_shiftSize,my_max(0,sampleN-outN));
  128. assert(outN+r <= sampleN);
  129. for(i=0;i<r;i++){
  130. sample[outN+i]=(short)out[i];
  131. }
  132. outN+=r;
  133. if(outN>=sampleN) break;
  134. }
  135. #if 0
  136. wiener.SaveInput("test.input.raw",0);
  137. wiener.SaveDenoised("test.output1.raw",0);
  138. #endif
  139. wiener.Close();
  140. }
  143. Wiener::Wiener()
  144. {
  145. m_sampleRate=DEFAULT_SAMPLING_RATE;
  146. m_NumChannels=NR_NUM_CHANNELS;
  147. m_localFrameX=0;
  148. }
  151. Wiener::~Wiener()
  152. {
  153. }
  156. void Wiener::Close()
  157. {
  158. }
  161. int Wiener::GetSample(short *sample, int sampleN)
  162. {
  163. long startX=(m_inputEndX)%NR_BUF_SIZE;
  164. if(startX+sampleN <= NR_BUF_SIZE){
  165. memmove((char*)(m_inputSpeech+startX), sample, 2*(sampleN));
  166. }
  167. else{
  168. memmove((char*)(m_inputSpeech+startX), sample, 2*(NR_BUF_SIZE-startX));
  169. memmove((char*)(m_inputSpeech), sample+(NR_BUF_SIZE-startX), 2*(sampleN-NR_BUF_SIZE+startX));
  170. }
  171. m_inputEndX += sampleN;
  172. return sampleN;
  173. }
  176. int Wiener::Init(int samplingRate, int isWiener)
  177. {
  178. int i, kHz;
  179. m_isWiener=isWiener;
  180. m_sampleRate=samplingRate;
  182. /*-------------------------------------------------------------------*/
  183. /* Set parameters FrameLength, FrameShift and FFTLength              */
  184. /*-------------------------------------------------------------------*/
  185. kHz=(int)(m_sampleRate/1000.0+0.5);
  186. if (kHz == SAMPLING_FREQ_1) { /* 8 kHz */
  187. m_winSize = FRAME_LENGTH_1;
  188. m_shiftSize = FRAME_SHIFT_1;
  189. m_fftSize = FFT_LENGTH_1;
  190. }
  191. else if (kHz == SAMPLING_FREQ_2) { /* 11 kHz */
  192. m_winSize = FRAME_LENGTH_2;
  193. m_shiftSize = FRAME_SHIFT_2;
  194. m_fftSize = FFT_LENGTH_2;
  195. }
  196. else if (kHz == SAMPLING_FREQ_3) { /* 16 kHz */
  197. m_winSize = FRAME_LENGTH_3;
  198. m_shiftSize = FRAME_SHIFT_3;
  199. m_fftSize = FFT_LENGTH_3;
  200. }
  201. else {
  202. float shiftMs=10, winMs=25;
  203. m_winSize = (int)(winMs/1000*m_sampleRate);
  204. m_shiftSize = (int)(shiftMs/1000*m_sampleRate);
  205. m_fftSize = 1; while (m_fftSize<m_winSize) m_fftSize *= 2;
  206. }
  208. if(m_shiftSize>NR_MAX_FRAME_SHIFT){
  209. fprintf(stderr, "ERROR:   Too large frame shift size '%d'!rn", m_shiftSize);
  210. assert(0);
  211. return 0;
  212. }
  214. if(m_isWiener==1){
  215. m_specLength=m_fftSize/4+1;
  216. }
  217. else{
  218. m_specLength=m_fftSize/2+1;
  219. }
  221. m_inputEndX=0;
  222. #ifdef _DEBUG
  223. m_denEndX = 0;
  224. #endif
  225. m_localFrameX=0;
  227. /*-------------------------------------------------------*/
  228. /* Initialize Hanning window and scaling factor          */
  229. /*-------------------------------------------------------*/
  230. InitHanning(m_HanningWin, m_winSize);
  232. if(m_isWiener==1){
  233. m_bufStartX = m_shiftSize-20;
  234. InitHanning(m_HanningWin2, NR_FL);
  235. InitMelFilterBanks ((float)NR_STARTING_FREQ, (float)m_sampleRate, 2*(m_specLength-1), m_NumChannels);
  236. InitMelIDCTMatrix (m_melIdctMatrix, m_NumChannels);
  237. }
  239. /* amplitude scaling factor */
  240. m_scaleFactor=0;
  241. for(i=0;i<m_shiftSize;i++){
  242. float sum=0;
  243. for(int n=0;n+i<m_winSize;n+=m_shiftSize) sum += m_HanningWin[i+n];
  244. if(sum>m_scaleFactor) m_scaleFactor=sum;
  245. }
  247. /* over-subtraction factor */
  248. if(m_isWiener==0){
  249. m_oversubGain = 8;
  250. m_oversubCutoffFreq = 800;
  251. float fc=m_fftSize*m_oversubCutoffFreq/(float)m_sampleRate;
  252. for(i=0;i<m_specLength;i++){
  253. m_oversubFactor[i] = m_oversubGain/(1+i/fc);
  254. }
  255. }
  257.     return 1;
  258. }
  261. FeReturnCode Wiener::InitNewUtterance(const char *fname)
  262. {
  263. int i;
  264. m_inputEndX = 0;
  265. #ifdef _DEBUG
  266. m_denEndX = 0;
  267. #endif
  268. m_localFrameX=0;
  270. /* Wiener filter design */
  271. m_nbFrameX=1;
  273. /* VADNest */
  274. m_nbSpeechFrame=0;
  275. m_meanEn=0;
  276. m_flagVADNest=0;
  277. m_hangOver=0;
  278. m_nbFrameVADNest=0;
  280. for(i=0;i<m_specLength;i++){
  281. m_lastSpectrum2[i]=0;
  282. m_lastSpectrum[i]=0;
  283. m_sqrtNoisePSD[i]=(float)NR_EPS;
  284. m_sqrtDen3PSD[i]=0;
  285. }
  287. for(i=0;i<4*m_shiftSize;i++){
  288. m_buf_in[i]=0;
  289. m_buf_out[i]=0;
  290. }
  292. return FE_OK;
  293. }
  296. FeReturnCode Wiener::OneFrame(short *sample, int sampleN, float *out, int frameX)
  297. {
  298.     int i, k, extraFrameN;
  299. long startX=m_shiftSize*m_localFrameX;
  300. long validX;
  301. k=(m_winSize-m_shiftSize)/m_shiftSize;
  302. extraFrameN=((m_winSize-m_shiftSize)%m_shiftSize == 0 ? k : k+1);
  304. sampleN=GetSample(sample,sampleN);
  305. assert(m_localFrameX-frameX <= 4+extraFrameN);
  307. float si[NR_MAX_WIN_SIZE], so[NR_MAX_WIN_SIZE], den[NR_MAX_WIN_SIZE];
  308. FeReturnCode status;
  310. /* shift down one frame in the first buffer and insert a new frame */
  311. for(i=0;i<3*m_shiftSize;i++) m_buf_in[i]=m_buf_in[i+m_shiftSize];
  312. for(i=0;i<m_shiftSize;i++) m_buf_in[3*m_shiftSize+i]=(float)m_inputSpeech[(startX+i)%NR_BUF_SIZE];
  313. for(i=0;i<m_winSize;i++) si[i]=m_buf_in[i];
  315. if(m_isWiener==1){
  316. status=OneFrameWiener(si, so);
  317. }
  318. else{
  319. status=OneFrameSS(si, so);
  320. }
  321. if(status==FE_NULL){
  322. m_localFrameX++;
  323. return FE_NULL;
  324. }
  326. if(m_isWiener==0){
  327. /* overlap and save */
  328. for(i=0;i<3*m_shiftSize;i++) m_buf_out[i]=m_buf_out[i+m_shiftSize];
  329. for(i=0;i<m_shiftSize;i++) m_buf_out[3*m_shiftSize+i]=0;
  330. for(i=0;i<m_winSize;i++) m_buf_out[i]+=so[i];
  332. /* get the result */
  333. for(i=0;i<m_shiftSize;i++) den[i]=m_buf_out[i]/m_scaleFactor;
  334. }
  335. else{
  336. for(i=0;i<m_shiftSize;i++) den[i]=so[i];
  337. }
  339. /* save the denoised speech at the ring buffer */
  340. for(i=0;i<m_shiftSize;i++){
  341. m_outSpeech[(startX+i)%NR_OUT_BUF_SIZE]=den[i];
  342. }
  343. #ifdef _DEBUG
  344. /* save the denoised speech at the ring buffer */
  345. validX=startX-2*m_shiftSize;
  346. for(i=0;i<m_shiftSize;i++){
  347. m_denSpeech[(validX+i)%NR_BUF_SIZE]=(short)(den[i]);
  348. }
  349. #endif
  351. /* NR has 5 frame delay due to one buffer and extraFrameN delay to make one frame for feature extraction */
  352. if(m_localFrameX<5+extraFrameN){
  353. m_localFrameX++;
  354. return FE_NULL;
  355. }
  357. validX=startX-extraFrameN*m_shiftSize-((m_isWiener==1) ? 2*m_shiftSize : 0);
  358. for(i=0;i<m_winSize;i++) {
  359. out[i]=m_outSpeech[(validX+i)%NR_OUT_BUF_SIZE];
  360. }
  362. #ifdef _DEBUG
  363. m_denEndX += m_shiftSize;
  364. #endif
  365. m_localFrameX++;
  366. return FE_SPEECH;
  368. }
  371. FeReturnCode Wiener::OneFrameWiener(float *si, float *out)
  372. {
  373. VADNest(m_localFrameX, si);
  374. EstimateSpectrum(si, m_spec, m_spec_re, m_spec_im, 1);
  375. ComputeMeanPSD(m_spec, m_lastSpectrum, m_lastSpectrum2, m_flagVADNest, m_sqrtInPSD);
  376. if(m_localFrameX<2) {
  377. int i;
  378. for(i=0;i<m_winSize;i++) out[i]=si[i];
  379. return FE_NULL;
  380. }
  381. DesignWiener(m_localFrameX, m_flagVADNest, m_spec, m_sqrtInPSD, m_sqrtNoisePSD, m_sqrtDen3PSD, m_wienerFilter);
  382. if(m_isWiener==0){
  383. ApplyFilter(m_spec_re, m_spec_im, m_wienerFilter, out);
  384. }
  385. else{
  386. MelFilterBank(m_wienerFilter, m_H2mel);
  387. MelIDCT(m_H2mel, m_hWFmirr);
  388. assert(m_shiftSize-m_bufStartX > (NR_FL-1)/2);
  389. ApplyWiener(si+m_shiftSize-m_bufStartX, m_hWFmirr, m_hWFw, out);
  390. }
  391. return FE_SPEECH;
  392. }
  395. FeReturnCode Wiener::OneFrameSS(float *si, float *out)
  396. {
  397. EstimateSpectrum(si, m_spec, m_spec_re, m_spec_im, 0);
  398. ComputeMeanPSD(m_spec, m_lastSpectrum, m_lastSpectrum2, m_flagVADNest, m_sqrtInPSD);
  399. if(m_localFrameX<2) {
  400. int i;
  401. for(i=0;i<m_winSize;i++) out[i]=si[i];
  402. return FE_NULL;
  403. }
  404. DesignSpecsub(m_localFrameX, m_flagVADNest, m_spec, m_sqrtInPSD, m_sqrtNoisePSD, m_sqrtDen3PSD, m_ssFilter);
  405. ApplyFilter(m_spec_re, m_spec_im, m_ssFilter, out);
  406. return FE_SPEECH;
  407. }
  410. void Wiener::InitHanning (float *win, int len)
  411. {
  412.     int i;
  413.     for (i = 0; i < len; i++){
  414. win[i] = (float)(0.5 - 0.5 * cos (2 * M_PI * (i + 0.5)/len));
  415. }
  416. }
  419. void Wiener::EstimateSpectrum(float *s, float *spectrum, float *re, float *im, int subSample)
  420. {
  421. int i, n, m;
  422. float x[NR_MAX_WIN_SIZE];
  424. /* apply Hanning window */
  425. for(n=0;n<m_winSize;n++){
  426. re[n] = s[n]*m_HanningWin[n];
  427. im[n] = 0;
  428. }
  430. /* pad zeros */
  431. for(n=m_winSize;n<m_fftSize;n++){
  432. re[n]=0;
  433. im[n]=0;
  434. }
  436. /* do FFT */
  437. m = (int) (log10 (m_fftSize) / log10 (2));
  438. PRFFT_NEW(re, im, m, m_fftSize, 1);
  440. /* take power spectrum, P(bin)=|X(bin)|^2 */
  441. x[0] = re[0] * re[0];  /* DC */
  442. for (n = 1; n < m_fftSize / 2; n++) {
  443. x[n] = re[n] * re[n] + im[n] * im[n];
  444. }
  445. x[m_fftSize / 2] = re[m_fftSize / 2] * re[m_fftSize / 2];  /* pi/2 */
  447. if(subSample==1){
  448. /* smooth and subsample spectrum */
  449. assert(m_specLength==m_fftSize/4+1);
  450. for(i=0;i<m_fftSize/4;i++){
  451. spectrum[i]=(x[2*i]+x[2*i+1])/2;
  452. }
  453. spectrum[m_fftSize/4]=x[m_fftSize/2];
  454. }
  455. else{
  456. for(i=0;i<m_fftSize/2+1;i++){
  457. spectrum[i]=x[i];
  458. }
  459. }
  460. }
  463. void Wiener::ComputeMeanPSD(float *spectrum, float *lastSpectrum, float *lastSpectrum2, int flagVADNest, float *sqrtInPSD)
  464. {
  465. int i;
  466. for(i=0;i<m_specLength;i++){
  467. sqrtInPSD[i]=(float)sqrt((double)(lastSpectrum2[i]+lastSpectrum[i]+spectrum[i])/(double)3.0);
  468. lastSpectrum2[i]=lastSpectrum[i];
  469. lastSpectrum[i]=spectrum[i];
  470. }
  471. }
  474. void Wiener::DesignWiener(int t, int flagVADNest, const float *in, const float *sqrtInPSD, float *sqrtNoisePSD, float *sqrtDen3PSD, float *filter)
  475. {
  476. int i;
  477. float lambdaNSE;
  478. float sqrtDenPSD[NR_MAX_SPEC_LENGTH], sqrtDen2PSD[NR_MAX_SPEC_LENGTH];
  479. float sqrtEta[NR_MAX_SPEC_LENGTH], sqrtEta2[NR_MAX_SPEC_LENGTH], h[NR_MAX_SPEC_LENGTH];
  481. /* set the forgetting factor */
  482. if(t<NR_NB_FRAME_THRESHOLD_NSE)
  483. lambdaNSE=1-1/(float)(m_nbFrameX);
  484. else
  485. lambdaNSE=(float)NR_LAMBDA_NSE;
  487. if(t>=2) m_nbFrameX++;
  489. if(flagVADNest==0){ /* this is the result of the previous frame */
  490. for(i=0;i<m_specLength;i++){
  491. /* Update. Pnoise(tn)^0.5 = max(a * Pnoise(tn-1)^0.5 + (1-a) * Pin_PSD(tn)^0.5, EPS) */
  492. sqrtNoisePSD[i]=(float)my_max(lambdaNSE*sqrtNoisePSD[i]+(1-lambdaNSE)*sqrtInPSD[i], NR_EPS);
  493. }
  494. }
  495. else{
  496. /* No change. Pnoise(t)^0.5 = Pnoise(tn)^0.5 */
  497. }
  499. if(t>=2){
  500. /* Be careful about sqrt and square! */
  501. for(i=0;i<m_specLength;i++){
  502. sqrtDenPSD[i]=(float)(NR_BETA*sqrtDen3PSD[i]+(1-NR_BETA)*NR_T(sqrtInPSD[i]-sqrtNoisePSD[i]));
  503. sqrtEta[i]=sqrtDenPSD[i]/sqrtNoisePSD[i];
  504. h[i]=sqrtEta[i]/(1+sqrtEta[i]);
  505. sqrtDen2PSD[i]=h[i]*(sqrtInPSD[i]);
  506. sqrtEta2[i]=my_max(sqrtDen2PSD[i]/sqrtNoisePSD[i], (float)NR_ETA_TH);
  507. filter[i]=sqrtEta2[i]/(1+sqrtEta2[i]);
  508. sqrtDen3PSD[i]=(float)(filter[i]*sqrt(in[i]));
  509. }
  510. }
  511. }
  514. void Wiener::DesignSpecsub(int t, int flagVADNest, const float *in, const float *sqrtInPSD, float *sqrtNoisePSD, float *sqrtDen3PSD, float *filter)
  515. {
  516. int i;
  517. float lambdaNSE;
  518. float s, n;
  519. float globalFac;
  521. /* set the forgetting factor */
  522. if(t<NR_NB_FRAME_THRESHOLD_NSE)
  523. lambdaNSE=1-1/(float)(m_nbFrameX);
  524. else
  525. lambdaNSE=(float)NR_LAMBDA_NSE;
  527. if(t>=2) m_nbFrameX++;
  529. if(flagVADNest==0 && t<NR_NB_FRAME_THRESHOLD_NSE){ /* this is the result of the previous frame */
  530. for(i=0;i<m_specLength;i++){
  531. sqrtNoisePSD[i]=(float)sqrt(my_max(lambdaNSE*NR_SQ(sqrtNoisePSD[i])
  532. +(1-lambdaNSE)*(in[i]-NR_SQ(sqrtDen3PSD[i])), NR_EPS*NR_EPS));
  533. }
  534. }
  535. else{
  536. /* No change. Pnoise(t)^0.5 = Pnoise(tn)^0.5 */
  537. }
  539. if(t<2) return;
  541. /* compute filter with over-subtraction */
  542. s=0, n=0;
  543. for(i=0;i<m_specLength;i++){
  544. s+=NR_SQ(sqrtInPSD[i]); n+=NR_SQ(sqrtNoisePSD[i]);
  545. }
  546. s/=m_specLength; n/=m_specLength;
  547. globalFac=my_max(n/s,(float)NR_SS_TH);
  549. for(i=0;i<m_specLength;i++){
  550. float instSNR = NR_SQ(sqrtNoisePSD[i])/(NR_SQ(sqrtNoisePSD[i])+NR_SQ(sqrtInPSD[i]));
  551. float instSN = NR_SQ(sqrtNoisePSD[i])/NR_SQ(sqrtInPSD[i]+(float)1e-20);
  552. filter[i] = (float)sqrt(my_max(1-(1+m_oversubFactor[i]*instSNR)*instSN, (float)NR_SS_TH*NR_SS_TH*instSN));
  553. sqrtDen3PSD[i]=(float)(filter[i]*sqrt(in[i]));
  554. }
  556. return;
  557. }
  560. void Wiener::VADNest(int t, const float *s)
  561. {
  562. int i;
  563. float lambdaLTE, lambdaLTEhigherE=(float)0.99;
  564. float frameEn=0, mean=0;
  565. int M=2*m_shiftSize; /* changed to get smooth estimate */
  567. if(t<NR_NB_FRAME_THRESHOLD_LTE)
  568. lambdaLTE=1-1/(float)(m_nbFrameVADNest+1);
  569. else
  570. lambdaLTE=(float)NR_LAMBDA_LTE;
  572. for(i=0;i<M;i++) {
  573. frameEn += s[i]*s[i];
  574. mean += s[i];
  575. }
  576. frameEn = frameEn - mean*mean/M; /* owkwon: Use the AC energy only */
  577. frameEn = frameEn*80/(float)M; /* owkwon: normalize to a value for 80 samples */
  579. if(t>=2 && frameEn>0) { /* Start updating after only valid non-zero data */
  580. m_nbFrameVADNest++;
  581. }
  583. frameEn = (float)(0.5+16/log(2)*log((64+frameEn)/64)); /* 16/log(2) ~ 23.08 */
  584. if(frameEn<NR_ENERGY_FLOOR) frameEn=NR_ENERGY_FLOOR; /* to avoid the worst case of all-zero data */
  585. if(m_nbFrameVADNest==1){
  586. m_meanEn=frameEn; /* first non-zero frame */
  587. }
  589. if(frameEn-m_meanEn<NR_SNR_THRESHOLD_NO_UPD_LTE && (frameEn-m_meanEn < NR_SNR_THRESHOLD_UPD_LTE || t < NR_MIN_FRAME)){
  590. if(frameEn<m_meanEn || t<NR_MIN_FRAME)
  591. m_meanEn = m_meanEn + (1-lambdaLTE)*(frameEn-m_meanEn);
  592. else
  593. m_meanEn = m_meanEn + (1-lambdaLTEhigherE)*(frameEn-m_meanEn);
  595. if(m_meanEn<NR_ENERGY_FLOOR) m_meanEn=NR_ENERGY_FLOOR;
  596. }
  598. if(m_nbFrameVADNest>0 && frameEn-m_meanEn>=NR_SNR_THRESHOLD_NO_UPD_LTE){ /* owkwon:  */
  599. m_flagVADNest=1;
  600. m_nbSpeechFrame++;
  601. }
  602. else if(frameEn-m_meanEn > NR_SNR_THRESHOLD_VAD && t>=NR_MIN_FRAME){ /* owkwon: changed the inquality sign to '>' */
  603. m_flagVADNest=1;
  604. m_nbSpeechFrame++;
  605. }
  606. else{
  607. if((m_nbSpeechFrame < NR_MIN_SPEECH_FRAME_HANGOVER) && (m_nbSpeechFrame >= 2)){ /* owkwon: added the condition '&& (m_nbSpeechFrame >= 2)' */
  608. m_hangOver=NR_HANGOVER;
  609. }
  610. m_nbSpeechFrame=0;
  611. if(m_hangOver != 0){
  612. m_hangOver=m_hangOver-1;
  613. m_flagVADNest=1;
  614. }
  615. else{
  616. m_flagVADNest=0;
  617. }
  618. }
  619. }
  622. void Wiener::ApplyFilter(float *re, float *im, float *h, float *out)
  623. {
  624. int i;
  625. re[0] *= h[0]; im[0] = 0;
  626. for(i=1;i<m_fftSize/2;i++){
  627. re[i] *= h[i]; im[i] *= h[i];
  628. re[m_fftSize-i] = re[i]; im[m_fftSize-i] = -im[i];
  629. }
  630. re[m_fftSize/2] = 0;
  632. int m = (int) (log10 (m_fftSize) / log10 (2));
  633. PRFFT_NEW(re,im,m,m_fftSize,-1);
  635. for(i=0;i<m_winSize;i++){
  636. out[i]=re[i];
  637. }
  638. }
  641. #ifdef _DEBUG
  642. int Wiener::SaveInput(const char *fname, int offsetX)
  643. {
  644. int i, sampleN;
  645. int endX=my_max(0,m_inputEndX-offsetX*m_shiftSize);
  646. FILE *fp=fopen(fname,"wb");
  647. if(fp==0){
  648. assert(0);
  649. }
  650. sampleN=my_min(NR_BUF_SIZE,endX);
  651. for(i=endX-sampleN;i<endX;i++){
  652. fwrite((char*)(&m_inputSpeech[i%NR_BUF_SIZE]), sizeof(short), 1, fp);
  653. }
  654. fclose(fp);
  655. return 1;
  656. }
  659. int Wiener::SaveDenoised(const char *fname, int offsetX)
  660. {
  661. int i, sampleN;
  662. FILE *fp=fopen(fname,"wb");
  663. if(fp==0){
  664. assert(0);
  665. }
  666. int endX=m_inputEndX-offsetX*m_shiftSize;
  667. sampleN=my_min(NR_BUF_SIZE,endX);
  668. for(i=endX-sampleN;i<endX;i++){
  669. short s=m_denSpeech[i%NR_BUF_SIZE];
  670. fwrite((char*)&s, sizeof(short), 1, fp);
  671. }
  672. fclose(fp);
  673. return 1;
  674. }
  675. #endif
  678. void Wiener::MelFilterBank(float *h2, float *h2mel)
  679. {
  680. int i,k;
  681. for(k=0;k<=m_NumChannels+1;k++){
  682. WfMelFB *melFB=&m_MelFB[k];
  683. float sum = h2[melFB->m_centerX]; /* The center weight is always 1 */
  684. for(i=melFB->m_lowX;i<melFB->m_centerX;i++){
  685. sum += m_MelWeight[i]*h2[i];
  686. }
  687. for(i=melFB->m_centerX+1;i<=melFB->m_highX;i++){
  688. sum += (1-m_MelWeight[i])*h2[i];
  689. }
  690. h2mel[k]=sum/melFB->m_sumWeight;
  691. }
  692. }
  695. void Wiener::MelIDCT(float *h2mel, float *hWFmirr)
  696. {
  697. int n,k;
  698. float hWF[NR_NUM_CHANNELS+2];
  700. for(n=0;n<=NR_NUM_CHANNELS+1;n++){
  701. float sum=0;
  702. for(k=0;k<=NR_NUM_CHANNELS+1;k++)
  703. sum += h2mel[k]*m_melIdctMatrix[k*(NR_NUM_CHANNELS+2)+n];
  704. hWF[n]=sum;
  705. }
  707. for(n=0;n<=NR_NUM_CHANNELS+1;n++)
  708. hWFmirr[n]=hWF[n];
  709. for(n=NR_NUM_CHANNELS+2;n<=2*(NR_NUM_CHANNELS+1);n++)
  710. hWFmirr[n]=hWF[2*(NR_NUM_CHANNELS+1)+1-n];
  711. }
  714. void Wiener::ApplyWiener(float *s, float *hWFmirr, float *hWFw, float *out)
  715. {
  716. int i,n;
  717. float hWFcaus[2*(NR_NUM_CHANNELS+1)+1];
  718. float hWFtrunc[NR_FL];
  719. int FL2=(NR_FL-1)/2;
  720. float fac=0;
  721. for(n=0;n<=NR_NUM_CHANNELS;n++)
  722. hWFcaus[n]=hWFmirr[n+NR_NUM_CHANNELS+1];
  723. for(n=NR_NUM_CHANNELS+1;n<=2*(NR_NUM_CHANNELS+1);n++)
  724. hWFcaus[n]=hWFmirr[n-NR_NUM_CHANNELS-1];
  726. for(n=0;n<=NR_FL-1;n++)
  727. hWFtrunc[n]=hWFcaus[n+NR_NUM_CHANNELS+1-FL2];
  729. for(n=0;n<=NR_FL-1;n++){
  730. hWFw[n]=m_HanningWin2[n]*hWFtrunc[n];
  731. fac += m_HanningWin2[n];
  732. }
  734. for(n=0;n<m_shiftSize;n++){
  735. float sum=0;
  736. for(i=-FL2; i<=FL2; i++) sum += hWFw[i+FL2]*s[n-i];
  737. out[n]=sum;
  738. }
  739. }
  742. /*---------------------------------------------------------------------------
  743. Initialize data structure for FFT windows (mel filter bank).
  744. Computes start (bin[k-1]), center (bin[k]) and end (bin[k+1]) points
  745. of each mel filter to FFT index.
  746. Input:
  747.    m_sampleRate     Sampling frequency
  748.    m_FFTLength    FFT length
  749.    m_NumChannels  Number of channels
  750. Output:
  751.   Initialized weights and mel filters
  752.  *---------------------------------------------------------------------------*/
  753. void Wiener::InitMelFilterBanks (float startingFrequency, float samplingRate, int fftLength, int numChannels)
  754. {
  755.     int i, k;
  756.     float freq[NR_NUM_CHANNELS+2];
  757. int bin[NR_NUM_CHANNELS+2];
  758. float start_mel, fs_per_2_mel;
  760.     /* Constants for calculation */
  761. freq[0] = startingFrequency;
  762.     start_mel = (float)(2595.0 * log10 (1.0 + startingFrequency / 700.0));
  763. bin[0] = Round(fftLength * freq[0] / samplingRate);
  764. freq[numChannels+1] = (float)(samplingRate / 2.0);
  765.     fs_per_2_mel = (float)(2595.0 * log10 (1.0 + (samplingRate / 2.0) / 700.0));
  766. bin[numChannels+1] = Round(fftLength * freq[numChannels+1] / samplingRate);
  768. /* Calculating mel-scaled frequency and the corresponding FFT-bin */
  769.     /* number for the lower edge of the band                          */
  770. for (k = 1; k <= numChannels; k++) {
  771.         freq[k] = (float)(700 * (pow (10, (start_mel + (float) k / (numChannels + 1) * (fs_per_2_mel - start_mel)) / 2595.0) - 1.0));
  772. bin[k] = Round(fftLength * freq[k] / samplingRate);
  773. }
  775. /* This part is never used to compute MFCC coefficients */
  776. /* but initialized for completeness                     */
  777. for(i = 0; i<bin[0]; i++){
  778. m_MelWeight[i]=0;
  779. }
  780. m_MelWeight[fftLength/2]=1;
  782. /* Initialize low, center, high indices to FFT-bin */
  783. m_MelFB[0].m_lowX=0;
  784. m_MelFB[0].m_centerX=bin[0];
  785. m_MelFB[0].m_highX=bin[1]-1;
  786. for (k = 1; k <= numChannels; k++) {
  787. m_MelFB[k].m_lowX=bin[k-1]+1;
  788. m_MelFB[k].m_centerX=bin[k];
  789. m_MelFB[k].m_highX=bin[k+1]-1;
  790. }
  791. m_MelFB[numChannels+1].m_lowX=bin[numChannels]+1;
  792. m_MelFB[numChannels+1].m_centerX=bin[numChannels+1];
  793. m_MelFB[numChannels+1].m_highX=bin[numChannels+1];
  795. for (k = 1; k <= numChannels+1; k++) {
  796. for(i = bin[k-1]; i<bin[k]; i++){
  797. m_MelWeight[i]=(i-bin[k-1])/(float)(bin[k]-bin[k-1]);
  798. }
  799.     }
  800. m_MelWeight[bin[numChannels+1]]=0;
  802. m_MelFB[0].m_sumWeight=(bin[1]-bin[0]+1)/(float)2;
  803. for (k=1;k<=numChannels;k++){
  804. m_MelFB[k].m_sumWeight=(bin[k+1]-bin[k-1])/(float)2;
  805. }
  806. m_MelFB[numChannels+1].m_sumWeight=(bin[numChannels+1]-bin[numChannels]+1)/(float)2;
  807.     return;
  808. }
  811. /*---------------------------------------------------------------------------
  812. Initializes matrix for IDCT computation (IDCT is implemented
  813. as matrix-vector multiplication). The DCT matrix is of size
  814. (numChannels+2)-by-(numChannels+2).
  815.  *---------------------------------------------------------------------------*/
  816. int Wiener::InitMelIDCTMatrix (float *idctMatrix, int numChannels)
  817. {
  818.     int k, n;
  819. float center[NR_NUM_CHANNELS+2], df[NR_NUM_CHANNELS+2];
  820. float fac=m_sampleRate/(float)(2*(m_specLength-1));
  821. center[0]=0;
  822. center[NR_NUM_CHANNELS+1]=m_sampleRate/(float)2;
  823. for(k=1;k<=NR_NUM_CHANNELS;k++){
  824. int i;
  825. WfMelFB *melFB=&m_MelFB[k];
  826. float sum = (float)melFB->m_centerX;  /* The center weight is always 1 */
  827. for(i=melFB->m_lowX;i<melFB->m_centerX;i++) sum += m_MelWeight[i]*i;
  828. for(i=melFB->m_centerX+1;i<=melFB->m_highX;i++) sum += (1-m_MelWeight[i])*i;
  829. center[k]=sum*fac/melFB->m_sumWeight; /* owkwon: Use normalized W(k,i) */
  830. }
  832. df[0]=(center[1]-center[0])/(float)m_sampleRate;
  833. for(k=1;k<=NR_NUM_CHANNELS;k++){
  834. df[k]=(center[k+1]-center[k-1])/(float)m_sampleRate;
  835. }
  836. df[NR_NUM_CHANNELS+1]=(center[NR_NUM_CHANNELS+1]-center[NR_NUM_CHANNELS])/(float)m_sampleRate;
  838.     for (k = 0; k < numChannels+2; k++){
  839.         for (n = 0; n < numChannels+2; n++)
  840.             idctMatrix[k*(numChannels+2)+n] = (float) cos ((2*M_PI*n*center[k])/(float)m_sampleRate) * df[k];
  841. }
  842. return 1;
  843. }