00001
00023
00024
00025
00026
00027
00028
00029
00030 #include <sent/stddefs.h>
00031 #include <sent/mfcc.h>
00032
00033 #ifdef MFCC_SINCOS_TABLE
00034
00041 static void
00042 make_costbl_hamming(MFCCWork *w, int framesize)
00043 {
00044 int i;
00045 float a;
00046
00047 w->costbl_hamming = (double *)mymalloc(sizeof(double) * framesize);
00048 a = 2.0 * PI / (framesize - 1);
00049 for(i=1;i<=framesize;i++) {
00050
00051 w->costbl_hamming[i-1] = 0.54 - 0.46 * cos(a * (i - 1));
00052 }
00053 w->costbl_hamming_len = framesize;
00054 #ifdef MFCC_TABLE_DEBUG
00055 jlog("Stat: mfcc-core: generated Hamming cos table (%d bytes)\n",
00056 w->costbl_hamming_len * sizeof(double));
00057 #endif
00058 }
00059
00066 static void
00067 make_fft_table(MFCCWork *w, int n)
00068 {
00069 int m;
00070 int me, me1;
00071
00072 w->costbl_fft = (double *)mymalloc(sizeof(double) * n);
00073 w->sintbl_fft = (double *)mymalloc(sizeof(double) * n);
00074 for (m = 1; m <= n; m++) {
00075 me = 1 << m;
00076 me1 = me / 2;
00077 w->costbl_fft[m-1] = cos(PI / me1);
00078 w->sintbl_fft[m-1] = -sin(PI / me1);
00079 }
00080 w->tbllen = n;
00081 #ifdef MFCC_TABLE_DEBUG
00082 jlog("Stat: mfcc-core: generated FFT sin/cos table (%d bytes)\n", w->tbllen * sizeof(double));
00083 #endif
00084 }
00085
00093 static void
00094 make_costbl_makemfcc(MFCCWork *w, int fbank_num, int mfcc_dim)
00095 {
00096 int size;
00097 int i, j, k;
00098 float B, C;
00099
00100 size = fbank_num * mfcc_dim;
00101 w->costbl_makemfcc = (double *)mymalloc(sizeof(double) * size);
00102
00103 B = PI / fbank_num;
00104 k = 0;
00105 for(i=1;i<=mfcc_dim;i++) {
00106 C = i * B;
00107 for(j=1;j<=fbank_num;j++) {
00108 w->costbl_makemfcc[k] = cos(C * (j - 0.5));
00109 k++;
00110 }
00111 }
00112 w->costbl_makemfcc_len = size;
00113 #ifdef MFCC_TABLE_DEBUG
00114 jlog("Stat: mfcc-core: generated MakeMFCC cos table (%d bytes)\n",
00115 w->costbl_makemfcc_len * sizeof(double));
00116 #endif
00117 }
00118
00126 static void
00127 make_sintbl_wcep(MFCCWork *w, int lifter, int mfcc_dim)
00128 {
00129 int i;
00130 float a, b;
00131
00132 w->sintbl_wcep = (double *)mymalloc(sizeof(double) * mfcc_dim);
00133 a = PI / lifter;
00134 b = lifter / 2.0;
00135 for(i=0;i<mfcc_dim;i++) {
00136 w->sintbl_wcep[i] = 1.0 + b * sin((i+1) * a);
00137 }
00138 w->sintbl_wcep_len = mfcc_dim;
00139 #ifdef MFCC_TABLE_DEBUG
00140 jlog("Stat: mfcc-core: generated WeightCepstrum sin table (%d bytes)\n",
00141 w->sintbl_wcep_len * sizeof(double));
00142 #endif
00143 }
00144
00145 #endif
00146
00155 float Mel(int k, float fres)
00156 {
00157 return(1127 * log(1 + (k-1) * fres));
00158 }
00159
00168 void
00169 InitFBank(MFCCWork *w, Value *para)
00170 {
00171 float mlo, mhi, ms, melk;
00172 int k, chan, maxChan, nv2;
00173
00174
00175 w->fb.fftN = 2; w->fb.n = 1;
00176 while(para->framesize > w->fb.fftN){
00177 w->fb.fftN *= 2; w->fb.n++;
00178 }
00179
00180 nv2 = w->fb.fftN / 2;
00181 w->fb.fres = 1.0E7 / (para->smp_period * w->fb.fftN * 700.0);
00182 maxChan = para->fbank_num + 1;
00183 w->fb.klo = 2; w->fb.khi = nv2;
00184 mlo = 0; mhi = Mel(nv2 + 1, w->fb.fres);
00185
00186
00187 if (para->lopass >= 0) {
00188 mlo = 1127*log(1+(float)para->lopass/700.0);
00189 w->fb.klo = ((float)para->lopass * para->smp_period * 1.0e-7 * w->fb.fftN) + 2.5;
00190 if (w->fb.klo<2) w->fb.klo = 2;
00191 }
00192
00193 if (para->hipass >= 0) {
00194 mhi = 1127*log(1+(float)para->hipass/700.0);
00195 w->fb.khi = ((float)para->hipass * para->smp_period * 1.0e-7 * w->fb.fftN) + 0.5;
00196 if (w->fb.khi>nv2) w->fb.khi = nv2;
00197 }
00198
00199
00200 w->fb.cf = (float *)mymalloc((maxChan + 1) * sizeof(float));
00201 ms = mhi - mlo;
00202 for (chan = 1; chan <= maxChan; chan++)
00203 w->fb.cf[chan] = ((float)chan / maxChan)*ms + mlo;
00204
00205
00206 w->fb.loChan = (short *)mymalloc((nv2 + 1) * sizeof(short));
00207 for(k = 1, chan = 1; k <= nv2; k++){
00208 if (k < w->fb.klo || k > w->fb.khi) w->fb.loChan[k] = -1;
00209 else {
00210 melk = Mel(k, w->fb.fres);
00211 while (w->fb.cf[chan] < melk && chan <= maxChan) ++chan;
00212 w->fb.loChan[k] = chan - 1;
00213 }
00214 }
00215
00216
00217 w->fb.loWt = (float *)mymalloc((nv2 + 1) * sizeof(float));
00218 for(k = 1; k <= nv2; k++) {
00219 chan = w->fb.loChan[k];
00220 if (k < w->fb.klo || k > w->fb.khi) w->fb.loWt[k] = 0.0;
00221 else {
00222 if (chan > 0)
00223 w->fb.loWt[k] = (w->fb.cf[chan + 1] - Mel(k, w->fb.fres)) / (w->fb.cf[chan + 1] - w->fb.cf[chan]);
00224 else
00225 w->fb.loWt[k] = (w->fb.cf[1] - Mel(k, w->fb.fres)) / (w->fb.cf[1] - mlo);
00226 }
00227 }
00228
00229
00230 w->fb.Re = (float *)mymalloc((w->fb.fftN + 1) * sizeof(float));
00231 w->fb.Im = (float *)mymalloc((w->fb.fftN + 1) * sizeof(float));
00232
00233 w->sqrt2var = sqrt(2.0 / para->fbank_num);
00234
00235 }
00236
00242 void
00243 FreeFBank(FBankInfo *fb)
00244 {
00245 free(fb->cf);
00246 free(fb->loChan);
00247 free(fb->loWt);
00248 free(fb->Re);
00249 free(fb->Im);
00250 }
00251
00259 void
00260 ZMeanFrame(float *wave, int framesize)
00261 {
00262 int i;
00263 float mean;
00264
00265 mean = 0.0;
00266 for(i = 1; i <= framesize; i++) mean += wave[i];
00267 mean /= framesize;
00268 for(i = 1; i <= framesize; i++) wave[i] -= mean;
00269 }
00270
00279 float CalcLogRawE(float *wave, int framesize)
00280 {
00281 int i;
00282 double raw_E = 0.0;
00283 float energy;
00284
00285 for(i = 1; i <= framesize; i++)
00286 raw_E += wave[i] * wave[i];
00287 energy = (float)log(raw_E);
00288
00289 return(energy);
00290 }
00291
00299 void PreEmphasise (float *wave, int framesize, float preEmph)
00300 {
00301 int i;
00302
00303 for(i = framesize; i >= 2; i--)
00304 wave[i] -= wave[i - 1] * preEmph;
00305 wave[1] *= 1.0 - preEmph;
00306 }
00307
00315 void Hamming(float *wave, int framesize, MFCCWork *w)
00316 {
00317 int i;
00318 #ifdef MFCC_SINCOS_TABLE
00319 for(i = 1; i <= framesize; i++)
00320 wave[i] *= w->costbl_hamming[i-1];
00321 #else
00322 float a;
00323 a = 2 * PI / (framesize - 1);
00324 for(i = 1; i <= framesize; i++)
00325 wave[i] *= 0.54 - 0.46 * cos(a * (i - 1));
00326 #endif
00327 }
00328
00337 void FFT(float *xRe, float *xIm, int p, MFCCWork *w)
00338 {
00339 int i, ip, j, k, m, me, me1, n, nv2;
00340 double uRe, uIm, vRe, vIm, wRe, wIm, tRe, tIm;
00341
00342 n = 1<<p;
00343 nv2 = n / 2;
00344
00345 j = 0;
00346 for(i = 0; i < n-1; i++){
00347 if(j > i){
00348 tRe = xRe[j]; tIm = xIm[j];
00349 xRe[j] = xRe[i]; xIm[j] = xIm[i];
00350 xRe[i] = tRe; xIm[i] = tIm;
00351 }
00352 k = nv2;
00353 while(j >= k){
00354 j -= k; k /= 2;
00355 }
00356 j += k;
00357 }
00358
00359 for(m = 1; m <= p; m++){
00360 me = 1<<m; me1 = me / 2;
00361 uRe = 1.0; uIm = 0.0;
00362 #ifdef MFCC_SINCOS_TABLE
00363 wRe = w->costbl_fft[m-1]; wIm = w->sintbl_fft[m-1];
00364 #else
00365 wRe = cos(PI / me1); wIm = -sin(PI / me1);
00366 #endif
00367 for(j = 0; j < me1; j++){
00368 for(i = j; i < n; i += me){
00369 ip = i + me1;
00370 tRe = xRe[ip] * uRe - xIm[ip] * uIm;
00371 tIm = xRe[ip] * uIm + xIm[ip] * uRe;
00372 xRe[ip] = xRe[i] - tRe; xIm[ip] = xIm[i] - tIm;
00373 xRe[i] += tRe; xIm[i] += tIm;
00374 }
00375 vRe = uRe * wRe - uIm * wIm; vIm = uRe * wIm + uIm * wRe;
00376 uRe = vRe; uIm = vIm;
00377 }
00378 }
00379 }
00380
00381
00389 void
00390 MakeFBank(float *wave, MFCCWork *w, Value *para)
00391 {
00392 int k, bin, i;
00393 double Re, Im, A, P, NP, H, temp;
00394
00395 for(k = 1; k <= para->framesize; k++){
00396 w->fb.Re[k - 1] = wave[k]; w->fb.Im[k - 1] = 0.0;
00397 }
00398 for(k = para->framesize + 1; k <= w->fb.fftN; k++){
00399 w->fb.Re[k - 1] = 0.0; w->fb.Im[k - 1] = 0.0;
00400 }
00401
00402
00403 FFT(w->fb.Re, w->fb.Im, w->fb.n, w);
00404
00405 if (w->ssbuf != NULL) {
00406
00407 for(k = 1; k <= w->fb.fftN; k++){
00408 Re = w->fb.Re[k - 1]; Im = w->fb.Im[k - 1];
00409 P = sqrt(Re * Re + Im * Im);
00410 NP = w->ssbuf[k - 1];
00411 if((P * P - w->ss_alpha * NP * NP) < 0){
00412 H = w->ss_floor;
00413 }else{
00414 H = sqrt(P * P - w->ss_alpha * NP * NP) / P;
00415 }
00416 w->fb.Re[k - 1] = H * Re;
00417 w->fb.Im[k - 1] = H * Im;
00418 }
00419 }
00420
00421
00422 for(i = 1; i <= para->fbank_num; i++)
00423 w->fbank[i] = 0.0;
00424
00425 for(k = w->fb.klo; k <= w->fb.khi; k++){
00426 Re = w->fb.Re[k-1]; Im = w->fb.Im[k-1];
00427 A = sqrt(Re * Re + Im * Im);
00428 bin = w->fb.loChan[k];
00429 Re = w->fb.loWt[k] * A;
00430 if(bin > 0) w->fbank[bin] += Re;
00431 if(bin < para->fbank_num) w->fbank[bin + 1] += A - Re;
00432 }
00433
00434
00435 for(bin = 1; bin <= para->fbank_num; bin++){
00436 temp = w->fbank[bin];
00437 if(temp < 1.0) temp = 1.0;
00438 w->fbank[bin] = log(temp);
00439 }
00440 }
00441
00450 float CalcC0(MFCCWork *w, Value *para)
00451 {
00452 int i;
00453 float S;
00454
00455 S = 0.0;
00456 for(i = 1; i <= para->fbank_num; i++)
00457 S += w->fbank[i];
00458 return S * w->sqrt2var;
00459 }
00460
00468 void MakeMFCC(float *mfcc, Value *para, MFCCWork *w)
00469 {
00470 #ifdef MFCC_SINCOS_TABLE
00471 int i, j, k;
00472 k = 0;
00473
00474 for(i = 0; i < para->mfcc_dim; i++){
00475 mfcc[i] = 0.0;
00476 for(j = 1; j <= para->fbank_num; j++)
00477 mfcc[i] += w->fbank[j] * w->costbl_makemfcc[k++];
00478 mfcc[i] *= w->sqrt2var;
00479 }
00480 #else
00481 int i, j;
00482 float B, C;
00483
00484 B = PI / para->fbank_num;
00485
00486 for(i = 1; i <= para->mfcc_dim; i++){
00487 mfcc[i - 1] = 0.0;
00488 C = i * B;
00489 for(j = 1; j <= para->fbank_num; j++)
00490 mfcc[i - 1] += w->fbank[j] * cos(C * (j - 0.5));
00491 mfcc[i - 1] *= w->sqrt2var;
00492 }
00493 #endif
00494 }
00495
00503 void WeightCepstrum (float *mfcc, Value *para, MFCCWork *w)
00504 {
00505 #ifdef MFCC_SINCOS_TABLE
00506 int i;
00507 for(i=0;i<para->mfcc_dim;i++) {
00508 mfcc[i] *= w->sintbl_wcep[i];
00509 }
00510 #else
00511 int i;
00512 float a, b, *cepWin;
00513
00514 cepWin = (float *)mymalloc(para->mfcc_dim * sizeof(float));
00515 a = PI / para->lifter;
00516 b = para->lifter / 2.0;
00517
00518 for(i = 0; i < para->mfcc_dim; i++){
00519 cepWin[i] = 1.0 + b * sin((i + 1) * a);
00520 mfcc[i] *= cepWin[i];
00521 }
00522
00523 free(cepWin);
00524 #endif
00525 }
00526
00527
00528
00529
00530
00531
00540 MFCCWork *
00541 WMP_work_new(Value *para)
00542 {
00543 MFCCWork *w;
00544
00545
00546 w = (MFCCWork *)mymalloc(sizeof(MFCCWork));
00547 memset(w, 0, sizeof(MFCCWork));
00548
00549
00550 InitFBank(w, para);
00551
00552 #ifdef MFCC_SINCOS_TABLE
00553
00554 make_costbl_hamming(w, para->framesize);
00555 make_fft_table(w, w->fb.n);
00556 if (para->mfcc_dim >= 0) {
00557 make_costbl_makemfcc(w, para->fbank_num, para->mfcc_dim);
00558 make_sintbl_wcep(w, para->lifter, para->mfcc_dim);
00559 }
00560 #endif
00561
00562
00563 w->fbank = (double *)mymalloc((para->fbank_num+1)*sizeof(double));
00564 w->bf = (float *)mymalloc(w->fb.fftN * sizeof(float));
00565 w->bflen = w->fb.fftN;
00566
00567 return w;
00568 }
00569
00578 void
00579 WMP_calc(MFCCWork *w, float *mfcc, Value *para)
00580 {
00581 float energy = 0.0;
00582 float c0 = 0.0;
00583 int p;
00584
00585 if (para->zmeanframe) {
00586 ZMeanFrame(w->bf, para->framesize);
00587 }
00588
00589 if (para->energy && para->raw_e) {
00590
00591 energy = CalcLogRawE(w->bf, para->framesize);
00592 }
00593
00594 PreEmphasise(w->bf, para->framesize, para->preEmph);
00595
00596 Hamming(w->bf, para->framesize, w);
00597 if (para->energy && ! para->raw_e) {
00598
00599 energy = CalcLogRawE(w->bf, para->framesize);
00600 }
00601
00602 MakeFBank(w->bf, w, para);
00603
00604 if (para->c0) c0 = CalcC0(w, para);
00605
00606 MakeMFCC(mfcc, para, w);
00607
00608 WeightCepstrum(mfcc, para, w);
00609
00610 p = para->mfcc_dim;
00611 if (para->c0) mfcc[p++] = c0;
00612 if (para->energy) mfcc[p++] = energy;
00613 }
00614
00620 void
00621 WMP_free(MFCCWork *w)
00622 {
00623 if (w->fbank) {
00624 FreeFBank(&(w->fb));
00625 free(w->fbank);
00626 free(w->bf);
00627 w->fbank = NULL;
00628 w->bf = NULL;
00629 }
00630 #ifdef MFCC_SINCOS_TABLE
00631 if (w->costbl_hamming) {
00632 free(w->costbl_hamming);
00633 w->costbl_hamming = NULL;
00634 }
00635 if (w->costbl_fft) {
00636 free(w->costbl_fft);
00637 w->costbl_fft = NULL;
00638 }
00639 if (w->sintbl_fft) {
00640 free(w->sintbl_fft);
00641 w->sintbl_fft = NULL;
00642 }
00643 if (w->costbl_makemfcc) {
00644 free(w->costbl_makemfcc);
00645 w->costbl_makemfcc = NULL;
00646 }
00647 if (w->sintbl_wcep) {
00648 free(w->sintbl_wcep);
00649 w->sintbl_wcep = NULL;
00650 }
00651 #endif
00652 free(w);
00653 }
00654