00001
00022
00023
00024
00025
00026
00027
00028
00029 #include <sent/stddefs.h>
00030 #include <sent/mfcc.h>
00031
00032
00033 #ifdef MFCC_SINCOS_TABLE
00034
00035
00036
00037 static double *costbl_hamming;
00038 static int costbl_hamming_len = 0;
00039
00040 static double *costbl_fft;
00041 static double *sintbl_fft;
00042 static int tbllen = 0;
00043
00044 static double *costbl_makemfcc;
00045 static int costbl_makemfcc_len = 0;
00046
00047 static double *sintbl_wcep;
00048 static int sintbl_wcep_len = 0;
00049
00050 #endif
00051
00052 static float sqrt2var;
00053
00054 #ifdef MFCC_SINCOS_TABLE
00055
00056
00062 void
00063 make_costbl_hamming(int framesize)
00064 {
00065 int i;
00066 float a;
00067
00068 if (costbl_hamming_len == framesize) return;
00069 if (costbl_hamming_len != 0) {
00070 j_printerr("WARNING!! frame size changed!! (%d -> %d)\n",
00071 costbl_hamming_len, framesize);
00072 free(costbl_hamming);
00073 }
00074 costbl_hamming = (double *)mymalloc(sizeof(double) * framesize);
00075
00076 a = 2.0 * PI / (framesize - 1);
00077 for(i=1;i<=framesize;i++) {
00078
00079 costbl_hamming[i-1] = 0.54 - 0.46 * cos(a * (i - 1));
00080 }
00081 costbl_hamming_len = framesize;
00082 #ifdef MFCC_TABLE_DEBUG
00083 j_printerr("generated Hamming cos table (%d bytes)\n",
00084 costbl_hamming_len * sizeof(double));
00085 #endif
00086 }
00087
00093 void
00094 make_fft_table(int n)
00095 {
00096 int m;
00097 int me, me1;
00098
00099 if (tbllen == n) return;
00100
00101 if (tbllen != 0) {
00102 j_printerr("WARNING!! fft size changed!! (%d -> %d)\n", tbllen, n);
00103 free(costbl_fft);
00104 free(sintbl_fft);
00105 }
00106 costbl_fft = (double *)mymalloc(sizeof(double) * n);
00107 sintbl_fft = (double *)mymalloc(sizeof(double) * n);
00108 for (m = 1; m <= n; m++) {
00109 me = 1 << m;
00110 me1 = me / 2;
00111 costbl_fft[m-1] = cos(PI / me1);
00112 sintbl_fft[m-1] = -sin(PI / me1);
00113 }
00114 tbllen = n;
00115 #ifdef MFCC_TABLE_DEBUG
00116 j_printerr("generated FFT sin/cos table (%d bytes)\n", tbllen * sizeof(double));
00117 #endif
00118 }
00119
00126 void
00127 make_costbl_makemfcc(int fbank_num, int mfcc_dim)
00128 {
00129 int size;
00130 int i, j, k;
00131 float B, C;
00132
00133 size = fbank_num * mfcc_dim;
00134 if (costbl_makemfcc_len == size) return;
00135 if (costbl_makemfcc_len != 0) {
00136 j_printerr("WARNING!! fbank num or mfcc dim changed!! (%d -> %d)\n",
00137 costbl_makemfcc_len, size);
00138 free(costbl_makemfcc);
00139 }
00140 costbl_makemfcc = (double *)mymalloc(sizeof(double) * size);
00141
00142 B = PI / fbank_num;
00143 k = 0;
00144 for(i=1;i<=mfcc_dim;i++) {
00145 C = i * B;
00146 for(j=1;j<=fbank_num;j++) {
00147 costbl_makemfcc[k] = cos(C * (j - 0.5));
00148 k++;
00149 }
00150 }
00151 costbl_makemfcc_len = size;
00152 #ifdef MFCC_TABLE_DEBUG
00153 j_printerr("generated MakeMFCC cos table (%d bytes)\n",
00154 costbl_makemfcc_len * sizeof(double));
00155 #endif
00156 }
00157
00164 void
00165 make_sintbl_wcep(int lifter, int mfcc_dim)
00166 {
00167 int i;
00168 float a, b;
00169
00170 if (sintbl_wcep_len == mfcc_dim) return;
00171 if (sintbl_wcep_len != 0) {
00172 j_printerr("WARNING!! mfcc dim changed!! (%d -> %d)\n",
00173 sintbl_wcep_len, mfcc_dim);
00174 free(sintbl_wcep);
00175 }
00176 sintbl_wcep = (double *)mymalloc(sizeof(double) * mfcc_dim);
00177 a = PI / lifter;
00178 b = lifter / 2.0;
00179 for(i=0;i<mfcc_dim;i++) {
00180 sintbl_wcep[i] = 1.0 + b * sin((i+1) * a);
00181 }
00182 sintbl_wcep_len = mfcc_dim;
00183 #ifdef MFCC_TABLE_DEBUG
00184 j_printerr("generated WeightCepstrum sin table (%d bytes)\n",
00185 sintbl_wcep_len * sizeof(double));
00186 #endif
00187 }
00188
00189 #endif
00190
00199 float Mel(int k, float fres)
00200 {
00201 return(1127 * log(1 + (k-1) * fres));
00202 }
00203
00211 FBankInfo InitFBank(Value para)
00212 {
00213 FBankInfo fb;
00214 float mlo, mhi, ms, melk;
00215 int k, chan, maxChan, nv2;
00216
00217
00218 fb.fftN = 2; fb.n = 1;
00219 while(para.framesize > fb.fftN){
00220 fb.fftN *= 2; fb.n++;
00221 }
00222
00223 nv2 = fb.fftN / 2;
00224 fb.fres = 1.0E7 / (para.smp_period * fb.fftN * 700.0);
00225 maxChan = para.fbank_num + 1;
00226 fb.klo = 2; fb.khi = nv2;
00227 mlo = 0; mhi = Mel(nv2 + 1, fb.fres);
00228
00229
00230 if (para.lopass >= 0) {
00231 mlo = 1127*log(1+(float)para.lopass/700.0);
00232 fb.klo = ((float)para.lopass * para.smp_period * 1.0e-7 * fb.fftN) + 2.5;
00233 if (fb.klo<2) fb.klo = 2;
00234 }
00235
00236 if (para.hipass >= 0) {
00237 mhi = 1127*log(1+(float)para.hipass/700.0);
00238 fb.khi = ((float)para.hipass * para.smp_period * 1.0e-7 * fb.fftN) + 0.5;
00239 if (fb.khi>nv2) fb.khi = nv2;
00240 }
00241
00242
00243 if((fb.cf = (float *)mymalloc((maxChan + 1) * sizeof(float))) == NULL){
00244 j_error("InitFBank: failed to malloc\n");
00245 }
00246 ms = mhi - mlo;
00247 for (chan = 1; chan <= maxChan; chan++)
00248 fb.cf[chan] = ((float)chan / maxChan)*ms + mlo;
00249
00250
00251 if((fb.loChan = (short *)mymalloc((nv2 + 1) * sizeof(short))) == NULL){
00252 j_error("InitFBank: failed to malloc\n");
00253 }
00254 for(k = 1, chan = 1; k <= nv2; k++){
00255 if (k < fb.klo || k > fb.khi) fb.loChan[k] = -1;
00256 else {
00257 melk = Mel(k, fb.fres);
00258 while (fb.cf[chan] < melk && chan <= maxChan) ++chan;
00259 fb.loChan[k] = chan - 1;
00260 }
00261 }
00262
00263
00264 if((fb.loWt = (float *)mymalloc((nv2 + 1) * sizeof(float))) == NULL){
00265 j_error("InitFBank: failed to malloc\n");
00266 }
00267 for(k = 1; k <= nv2; k++) {
00268 chan = fb.loChan[k];
00269 if (k < fb.klo || k > fb.khi) fb.loWt[k] = 0.0;
00270 else {
00271 if (chan > 0)
00272 fb.loWt[k] = (fb.cf[chan + 1] - Mel(k, fb.fres)) / (fb.cf[chan + 1] - fb.cf[chan]);
00273 else
00274 fb.loWt[k] = (fb.cf[1] - Mel(k, fb.fres)) / (fb.cf[1] - mlo);
00275 }
00276 }
00277
00278
00279 if((fb.Re = (float *)mymalloc((fb.fftN + 1) * sizeof(float))) == NULL){
00280 j_error("InitFBank: failed to malloc\n");
00281 }
00282 if((fb.Im = (float *)mymalloc((fb.fftN + 1) * sizeof(float))) == NULL){
00283 j_error("InitFBank: failed to malloc\n");
00284 }
00285
00286 #ifdef MFCC_SINCOS_TABLE
00287
00288 make_costbl_hamming(para.framesize);
00289 make_fft_table(fb.n);
00290 make_costbl_makemfcc(para.fbank_num, para.mfcc_dim);
00291 make_sintbl_wcep(para.lifter, para.mfcc_dim);
00292 #endif
00293
00294 sqrt2var = sqrt(2.0 / para.fbank_num);
00295
00296 return(fb);
00297 }
00298
00304 void
00305 FreeFBank(FBankInfo fb)
00306 {
00307 free(fb.cf);
00308 free(fb.loChan);
00309 free(fb.loWt);
00310 free(fb.Re);
00311 free(fb.Im);
00312 }
00313
00321 void
00322 ZMeanFrame(float *wave, int framesize)
00323 {
00324 int i;
00325 float mean;
00326
00327 mean = 0.0;
00328 for(i = 1; i <= framesize; i++) mean += wave[i];
00329 mean /= framesize;
00330 for(i = 1; i <= framesize; i++) wave[i] -= mean;
00331 }
00332
00341 float CalcLogRawE(float *wave, int framesize)
00342 {
00343 int i;
00344 double raw_E = 0.0;
00345 float energy;
00346
00347 for(i = 1; i <= framesize; i++)
00348 raw_E += wave[i] * wave[i];
00349 energy = (float)log(raw_E);
00350
00351 return(energy);
00352 }
00353
00360 void PreEmphasise (float *wave, Value para)
00361 {
00362 int i;
00363
00364 for(i = para.framesize; i >= 2; i--)
00365 wave[i] -= wave[i - 1] * para.preEmph;
00366 wave[1] *= 1.0 - para.preEmph;
00367 }
00368
00375 void Hamming(float *wave, int framesize)
00376 {
00377 int i;
00378 #ifdef MFCC_SINCOS_TABLE
00379 for(i = 1; i <= framesize; i++)
00380 wave[i] *= costbl_hamming[i-1];
00381 #else
00382 float a;
00383 a = 2 * PI / (framesize - 1);
00384 for(i = 1; i <= framesize; i++)
00385 wave[i] *= 0.54 - 0.46 * cos(a * (i - 1));
00386 #endif
00387 }
00388
00396 void FFT(float *xRe, float *xIm, int p)
00397 {
00398 int i, ip, j, k, m, me, me1, n, nv2;
00399 double uRe, uIm, vRe, vIm, wRe, wIm, tRe, tIm;
00400
00401 n = 1<<p;
00402 nv2 = n / 2;
00403
00404 j = 0;
00405 for(i = 0; i < n-1; i++){
00406 if(j > i){
00407 tRe = xRe[j]; tIm = xIm[j];
00408 xRe[j] = xRe[i]; xIm[j] = xIm[i];
00409 xRe[i] = tRe; xIm[i] = tIm;
00410 }
00411 k = nv2;
00412 while(j >= k){
00413 j -= k; k /= 2;
00414 }
00415 j += k;
00416 }
00417
00418 for(m = 1; m <= p; m++){
00419 me = 1<<m; me1 = me / 2;
00420 uRe = 1.0; uIm = 0.0;
00421 #ifdef MFCC_SINCOS_TABLE
00422 wRe = costbl_fft[m-1]; wIm = sintbl_fft[m-1];
00423 #else
00424 wRe = cos(PI / me1); wIm = -sin(PI / me1);
00425 #endif
00426 for(j = 0; j < me1; j++){
00427 for(i = j; i < n; i += me){
00428 ip = i + me1;
00429 tRe = xRe[ip] * uRe - xIm[ip] * uIm;
00430 tIm = xRe[ip] * uIm + xIm[ip] * uRe;
00431 xRe[ip] = xRe[i] - tRe; xIm[ip] = xIm[i] - tIm;
00432 xRe[i] += tRe; xIm[i] += tIm;
00433 }
00434 vRe = uRe * wRe - uIm * wIm; vIm = uRe * wIm + uIm * wRe;
00435 uRe = vRe; uIm = vIm;
00436 }
00437 }
00438 }
00439
00440
00450 void MakeFBank(float *wave, double *fbank, FBankInfo fb, Value para, float *ssbuf)
00451 {
00452 int k, bin, i;
00453 double Re, Im, A, P, NP, H, temp;
00454
00455 for(k = 1; k <= para.framesize; k++){
00456 fb.Re[k - 1] = wave[k]; fb.Im[k - 1] = 0.0;
00457 }
00458 for(k = para.framesize + 1; k <= fb.fftN; k++){
00459 fb.Re[k - 1] = 0.0; fb.Im[k - 1] = 0.0;
00460 }
00461
00462
00463 FFT(fb.Re, fb.Im, fb.n);
00464
00465 if (ssbuf != NULL) {
00466
00467 for(k = 1; k <= fb.fftN; k++){
00468 Re = fb.Re[k - 1]; Im = fb.Im[k - 1];
00469 P = sqrt(Re * Re + Im * Im);
00470 NP = ssbuf[k - 1];
00471 if((P * P - para.ss_alpha * NP * NP) < 0){
00472 H = para.ss_floor;
00473 }else{
00474 H = sqrt(P * P - para.ss_alpha * NP * NP) / P;
00475 }
00476 fb.Re[k - 1] = H * Re;
00477 fb.Im[k - 1] = H * Im;
00478 }
00479 }
00480
00481
00482 for(i = 1; i <= para.fbank_num; i++)
00483 fbank[i] = 0.0;
00484
00485 for(k = fb.klo; k <= fb.khi; k++){
00486 Re = fb.Re[k-1]; Im = fb.Im[k-1];
00487 A = sqrt(Re * Re + Im * Im);
00488 bin = fb.loChan[k];
00489 Re = fb.loWt[k] * A;
00490 if(bin > 0) fbank[bin] += Re;
00491 if(bin < para.fbank_num) fbank[bin + 1] += A - Re;
00492 }
00493
00494
00495 for(bin = 1; bin <= para.fbank_num; bin++){
00496 temp = fbank[bin];
00497 if(temp < 1.0) temp = 1.0;
00498 fbank[bin] = log(temp);
00499 }
00500 }
00501
00510 float CalcC0(double *fbank, Value para)
00511 {
00512 int i;
00513 float S;
00514
00515 S = 0.0;
00516 for(i = 1; i <= para.fbank_num; i++)
00517 S += fbank[i];
00518 return S * sqrt2var;
00519 }
00520
00528 void MakeMFCC(double *fbank, float *mfcc, Value para)
00529 {
00530 #ifdef MFCC_SINCOS_TABLE
00531 int i, j, k;
00532 k = 0;
00533
00534 for(i = 0; i < para.mfcc_dim; i++){
00535 mfcc[i] = 0.0;
00536 for(j = 1; j <= para.fbank_num; j++)
00537 mfcc[i] += fbank[j] * costbl_makemfcc[k++];
00538 mfcc[i] *= sqrt2var;
00539 }
00540 #else
00541 int i, j;
00542 float B, C;
00543
00544 B = PI / para.fbank_num;
00545
00546 for(i = 1; i <= para.mfcc_dim; i++){
00547 mfcc[i - 1] = 0.0;
00548 C = i * B;
00549 for(j = 1; j <= para.fbank_num; j++)
00550 mfcc[i - 1] += fbank[j] * cos(C * (j - 0.5));
00551 mfcc[i - 1] *= sqrt2var;
00552 }
00553 #endif
00554 }
00555
00562 void WeightCepstrum (float *mfcc, Value para)
00563 {
00564 #ifdef MFCC_SINCOS_TABLE
00565 int i;
00566 for(i=0;i<para.mfcc_dim;i++) {
00567 mfcc[i] *= sintbl_wcep[i];
00568 }
00569 #else
00570 int i;
00571 float a, b, *cepWin;
00572
00573 if((cepWin = (float *)mymalloc(para.mfcc_dim * sizeof(float))) == NULL){
00574 j_error("WeightCepstrum: failed to malloc\n");
00575 }
00576 a = PI / para.lifter;
00577 b = para.lifter / 2.0;
00578
00579 for(i = 0; i < para.mfcc_dim; i++){
00580 cepWin[i] = 1.0 + b * sin((i + 1) * a);
00581 mfcc[i] *= cepWin[i];
00582 }
00583
00584 free(cepWin);
00585 #endif
00586 }
00587
00588
00589
00590
00591
00592
00593
00594 static double *fbank;
00595 static FBankInfo fb;
00596
00604 void
00605 WMP_calc_init(Value para, float **bf, int *bflen)
00606 {
00607
00608 fb = InitFBank(para);
00609
00610 if((fbank = (double *)mymalloc((para.fbank_num+1)*sizeof(double))) == NULL){
00611 j_error("Error: WMP_calc_init failed\n");
00612 }
00613
00614 if((*bf = (float *)mymalloc(fb.fftN * sizeof(float))) == NULL){
00615 j_error("Error: WMP_calc_init failed\n");
00616 }
00617
00618 *bflen = fb.fftN;
00619
00620 }
00621
00631 void
00632 WMP_calc(float *mfcc, float *bf, Value para, float *ssbuf)
00633 {
00634 float energy = 0.0;
00635 float c0 = 0.0;
00636 int p;
00637
00638 if (para.zmeanframe) {
00639 ZMeanFrame(bf, para.framesize);
00640 }
00641
00642 if (para.energy && para.raw_e) {
00643
00644 energy = CalcLogRawE(bf, para.framesize);
00645 }
00646
00647 PreEmphasise(bf, para);
00648
00649 Hamming(bf, para.framesize);
00650 if (para.energy && ! para.raw_e) {
00651
00652 energy = CalcLogRawE(bf, para.framesize);
00653 }
00654
00655 MakeFBank(bf, fbank, fb, para, ssbuf);
00656
00657 if (para.c0) c0 = CalcC0(fbank, para);
00658
00659 MakeMFCC(fbank, mfcc, para);
00660
00661 WeightCepstrum(mfcc, para);
00662
00663 p = para.mfcc_dim;
00664 if (para.c0) mfcc[p++] = c0;
00665 if (para.energy) mfcc[p++] = energy;
00666 }
00667
00673 void
00674 WMP_calc_fin(float *bf)
00675 {
00676 FreeFBank(fb);
00677 free(fbank);
00678 free(bf);
00679 }