00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038 #include <stdio.h>
00039 #include <math.h>
00040 #include <string.h>
00041 #include <stdlib.h>
00042 #include <assert.h>
00043
00044 #ifdef HAVE_CONFIG_H
00045 #include <config.h>
00046 #endif
00047
00048 #ifdef _MSC_VER
00049 #pragma warning (disable: 4244)
00050 #endif
00051
00052 #include "prim_type.h"
00053 #include "ckd_alloc.h"
00054 #include "byteorder.h"
00055 #include "fixpoint.h"
00056 #include "fe.h"
00057 #include "fe_internal.h"
00058 #include "fe_warp.h"
00059 #include "genrand.h"
00060 #include "err.h"
00061
00062
00063
00064 #ifdef FIXED_POINT
00065 #define FLOAT2COS(x) FLOAT2FIX_ANY(x,30)
00066 #define COSMUL(x,y) FIXMUL_ANY(x,y,30)
00067 #else
00068 #define FLOAT2COS(x) (x)
00069 #define COSMUL(x,y) ((x)*(y))
00070 #endif
00071
00072 #ifdef FIXED_POINT
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084 static const unsigned char fe_logadd_table[] = {
00085 177, 177, 176, 176, 175, 175, 174, 174, 173, 173,
00086 172, 172, 172, 171, 171, 170, 170, 169, 169, 168,
00087 168, 167, 167, 166, 166, 165, 165, 164, 164, 163,
00088 163, 162, 162, 161, 161, 161, 160, 160, 159, 159,
00089 158, 158, 157, 157, 156, 156, 155, 155, 155, 154,
00090 154, 153, 153, 152, 152, 151, 151, 151, 150, 150,
00091 149, 149, 148, 148, 147, 147, 147, 146, 146, 145,
00092 145, 144, 144, 144, 143, 143, 142, 142, 141, 141,
00093 141, 140, 140, 139, 139, 138, 138, 138, 137, 137,
00094 136, 136, 136, 135, 135, 134, 134, 134, 133, 133,
00095 132, 132, 131, 131, 131, 130, 130, 129, 129, 129,
00096 128, 128, 128, 127, 127, 126, 126, 126, 125, 125,
00097 124, 124, 124, 123, 123, 123, 122, 122, 121, 121,
00098 121, 120, 120, 119, 119, 119, 118, 118, 118, 117,
00099 117, 117, 116, 116, 115, 115, 115, 114, 114, 114,
00100 113, 113, 113, 112, 112, 112, 111, 111, 110, 110,
00101 110, 109, 109, 109, 108, 108, 108, 107, 107, 107,
00102 106, 106, 106, 105, 105, 105, 104, 104, 104, 103,
00103 103, 103, 102, 102, 102, 101, 101, 101, 100, 100,
00104 100, 99, 99, 99, 98, 98, 98, 97, 97, 97,
00105 96, 96, 96, 96, 95, 95, 95, 94, 94, 94,
00106 93, 93, 93, 92, 92, 92, 92, 91, 91, 91,
00107 90, 90, 90, 89, 89, 89, 89, 88, 88, 88,
00108 87, 87, 87, 87, 86, 86, 86, 85, 85, 85,
00109 85, 84, 84, 84, 83, 83, 83, 83, 82, 82,
00110 82, 82, 81, 81, 81, 80, 80, 80, 80, 79,
00111 79, 79, 79, 78, 78, 78, 78, 77, 77, 77,
00112 77, 76, 76, 76, 75, 75, 75, 75, 74, 74,
00113 74, 74, 73, 73, 73, 73, 72, 72, 72, 72,
00114 71, 71, 71, 71, 71, 70, 70, 70, 70, 69,
00115 69, 69, 69, 68, 68, 68, 68, 67, 67, 67,
00116 67, 67, 66, 66, 66, 66, 65, 65, 65, 65,
00117 64, 64, 64, 64, 64, 63, 63, 63, 63, 63,
00118 62, 62, 62, 62, 61, 61, 61, 61, 61, 60,
00119 60, 60, 60, 60, 59, 59, 59, 59, 59, 58,
00120 58, 58, 58, 58, 57, 57, 57, 57, 57, 56,
00121 56, 56, 56, 56, 55, 55, 55, 55, 55, 54,
00122 54, 54, 54, 54, 53, 53, 53, 53, 53, 52,
00123 52, 52, 52, 52, 52, 51, 51, 51, 51, 51,
00124 50, 50, 50, 50, 50, 50, 49, 49, 49, 49,
00125 49, 49, 48, 48, 48, 48, 48, 48, 47, 47,
00126 47, 47, 47, 47, 46, 46, 46, 46, 46, 46,
00127 45, 45, 45, 45, 45, 45, 44, 44, 44, 44,
00128 44, 44, 43, 43, 43, 43, 43, 43, 43, 42,
00129 42, 42, 42, 42, 42, 41, 41, 41, 41, 41,
00130 41, 41, 40, 40, 40, 40, 40, 40, 40, 39,
00131 39, 39, 39, 39, 39, 39, 38, 38, 38, 38,
00132 38, 38, 38, 37, 37, 37, 37, 37, 37, 37,
00133 37, 36, 36, 36, 36, 36, 36, 36, 35, 35,
00134 35, 35, 35, 35, 35, 35, 34, 34, 34, 34,
00135 34, 34, 34, 34, 33, 33, 33, 33, 33, 33,
00136 33, 33, 32, 32, 32, 32, 32, 32, 32, 32,
00137 32, 31, 31, 31, 31, 31, 31, 31, 31, 31,
00138 30, 30, 30, 30, 30, 30, 30, 30, 30, 29,
00139 29, 29, 29, 29, 29, 29, 29, 29, 28, 28,
00140 28, 28, 28, 28, 28, 28, 28, 28, 27, 27,
00141 27, 27, 27, 27, 27, 27, 27, 27, 26, 26,
00142 26, 26, 26, 26, 26, 26, 26, 26, 25, 25,
00143 25, 25, 25, 25, 25, 25, 25, 25, 25, 24,
00144 24, 24, 24, 24, 24, 24, 24, 24, 24, 24,
00145 23, 23, 23, 23, 23, 23, 23, 23, 23, 23,
00146 23, 23, 22, 22, 22, 22, 22, 22, 22, 22,
00147 22, 22, 22, 22, 21, 21, 21, 21, 21, 21,
00148 21, 21, 21, 21, 21, 21, 21, 20, 20, 20,
00149 20, 20, 20, 20, 20, 20, 20, 20, 20, 20,
00150 19, 19, 19, 19, 19, 19, 19, 19, 19, 19,
00151 19, 19, 19, 19, 18, 18, 18, 18, 18, 18,
00152 18, 18, 18, 18, 18, 18, 18, 18, 18, 17,
00153 17, 17, 17, 17, 17, 17, 17, 17, 17, 17,
00154 17, 17, 17, 17, 16, 16, 16, 16, 16, 16,
00155 16, 16, 16, 16, 16, 16, 16, 16, 16, 16,
00156 16, 15, 15, 15, 15, 15, 15, 15, 15, 15,
00157 15, 15, 15, 15, 15, 15, 15, 15, 14, 14,
00158 14, 14, 14, 14, 14, 14, 14, 14, 14, 14,
00159 14, 14, 14, 14, 14, 14, 14, 13, 13, 13,
00160 13, 13, 13, 13, 13, 13, 13, 13, 13, 13,
00161 13, 13, 13, 13, 13, 13, 13, 12, 12, 12,
00162 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
00163 12, 12, 12, 12, 12, 12, 12, 12, 12, 11,
00164 11, 11, 11, 11, 11, 11, 11, 11, 11, 11,
00165 11, 11, 11, 11, 11, 11, 11, 11, 11, 11,
00166 11, 11, 11, 10, 10, 10, 10, 10, 10, 10,
00167 10, 10, 10, 10, 10, 10, 10, 10, 10, 10,
00168 10, 10, 10, 10, 10, 10, 10, 10, 10, 9,
00169 9, 9, 9, 9, 9, 9, 9, 9, 9, 9,
00170 9, 9, 9, 9, 9, 9, 9, 9, 9, 9,
00171 9, 9, 9, 9, 9, 9, 9, 9, 8, 8,
00172 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
00173 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
00174 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
00175 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
00176 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
00177 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
00178 7, 7, 7, 7, 7, 7, 7, 7, 6, 6,
00179 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
00180 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
00181 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
00182 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
00183 6, 5, 5, 5, 5, 5, 5, 5, 5, 5,
00184 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
00185 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
00186 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
00187 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
00188 5, 5, 5, 4, 4, 4, 4, 4, 4, 4,
00189 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
00190 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
00191 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
00192 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
00193 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
00194 4, 4, 4, 4, 4, 4, 4, 4, 3, 3,
00195 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
00196 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
00197 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
00198 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
00199 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
00200 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
00201 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
00202 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
00203 3, 3, 3, 3, 2, 2, 2, 2, 2, 2,
00204 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
00205 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
00206 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
00207 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
00208 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
00209 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
00210 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
00211 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
00212 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
00213 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
00214 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
00215 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
00216 2, 2, 2, 2, 2, 2, 1, 1, 1, 1,
00217 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
00218 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
00219 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
00220 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
00221 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
00222 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
00223 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
00224 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
00225 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
00226 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
00227 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
00228 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
00229 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
00230 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
00231 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
00232 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
00233 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
00234 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
00235 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
00236 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
00237 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
00238 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
00239 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
00240 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
00241 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
00242 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
00243 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
00244 1, 1, 1, 1, 1, 1, 1, 0
00245 };
00246 static const int fe_logadd_table_size = sizeof(fe_logadd_table) / sizeof(fe_logadd_table[0]);
00247
00248 static fixed32
00249 fe_log_add(fixed32 x, fixed32 y)
00250 {
00251 fixed32 d, r;
00252
00253 if (x > y) {
00254 d = (x - y) >> (DEFAULT_RADIX - 8);
00255 r = x;
00256 }
00257 else {
00258 d = (y - x) >> (DEFAULT_RADIX - 8);
00259 r = y;
00260 }
00261 if (d > fe_logadd_table_size)
00262 return r;
00263 else {
00264 r += ((fixed32)fe_logadd_table[d] << (DEFAULT_RADIX - 8));
00265
00266
00267
00268
00269
00270 return r;
00271 }
00272 }
00273
00274 static fixed32
00275 fe_log(float32 x)
00276 {
00277 if (x <= 0) {
00278 return MIN_FIXLOG;
00279 }
00280 else {
00281 return FLOAT2FIX(log(x));
00282 }
00283 }
00284 #endif
00285
00286 static float32
00287 fe_mel(melfb_t *mel, float32 x)
00288 {
00289 float32 warped = fe_warp_unwarped_to_warped(mel, x);
00290
00291 return (float32) (2595.0 * log10(1.0 + warped / 700.0));
00292 }
00293
00294 static float32
00295 fe_melinv(melfb_t *mel, float32 x)
00296 {
00297 float32 warped = (float32) (700.0 * (pow(10.0, x / 2595.0) - 1.0));
00298 return fe_warp_warped_to_unwarped(mel, warped);
00299 }
00300
00301 int32
00302 fe_build_melfilters(melfb_t *mel_fb)
00303 {
00304 float32 melmin, melmax, melbw, fftfreq;
00305 int n_coeffs, i, j;
00306
00307
00308 mel_fb->spec_start = ckd_malloc(mel_fb->num_filters * sizeof(*mel_fb->spec_start));
00309 mel_fb->filt_start = ckd_malloc(mel_fb->num_filters * sizeof(*mel_fb->filt_start));
00310 mel_fb->filt_width = ckd_malloc(mel_fb->num_filters * sizeof(*mel_fb->filt_width));
00311
00312
00313
00314 melmin = fe_mel(mel_fb, mel_fb->lower_filt_freq);
00315 melmax = fe_mel(mel_fb, mel_fb->upper_filt_freq);
00316
00317
00318 melbw = (melmax - melmin) / (mel_fb->num_filters + 1);
00319 if (mel_fb->doublewide) {
00320 melmin -= melbw;
00321 melmax += melbw;
00322 if ((fe_melinv(mel_fb, melmin) < 0) ||
00323 (fe_melinv(mel_fb, melmax) > mel_fb->sampling_rate / 2)) {
00324 E_WARN
00325 ("Out of Range: low filter edge = %f (%f)\n",
00326 fe_melinv(mel_fb, melmin), 0.0);
00327 E_WARN
00328 (" high filter edge = %f (%f)\n",
00329 fe_melinv(mel_fb, melmax), mel_fb->sampling_rate / 2);
00330 return FE_INVALID_PARAM_ERROR;
00331 }
00332 }
00333
00334
00335 fftfreq = mel_fb->sampling_rate / (float32) mel_fb->fft_size;
00336
00337
00338 n_coeffs = 0;
00339 for (i = 0; i < mel_fb->num_filters; ++i) {
00340 float32 freqs[3];
00341
00342
00343 for (j = 0; j < 3; ++j) {
00344 if (mel_fb->doublewide)
00345 freqs[j] = fe_melinv(mel_fb, (i + j * 2) * melbw + melmin);
00346 else
00347 freqs[j] = fe_melinv(mel_fb, (i + j) * melbw + melmin);
00348
00349 if (mel_fb->round_filters)
00350 freqs[j] = ((int)(freqs[j] / fftfreq + 0.5)) * fftfreq;
00351 }
00352
00353
00354 mel_fb->spec_start[i] = -1;
00355
00356 for (j = 0; j < mel_fb->fft_size/2+1; ++j) {
00357 float32 hz = j * fftfreq;
00358 if (hz < freqs[0])
00359 continue;
00360 else if (hz > freqs[2] || j == mel_fb->fft_size/2) {
00361
00362 mel_fb->filt_width[i] = j - mel_fb->spec_start[i];
00363
00364 mel_fb->filt_start[i] = n_coeffs;
00365 n_coeffs += mel_fb->filt_width[i];
00366 break;
00367 }
00368 if (mel_fb->spec_start[i] == -1)
00369 mel_fb->spec_start[i] = j;
00370 }
00371 }
00372
00373
00374 mel_fb->filt_coeffs = ckd_malloc(n_coeffs * sizeof(*mel_fb->filt_coeffs));
00375
00376
00377 n_coeffs = 0;
00378 for (i = 0; i < mel_fb->num_filters; ++i) {
00379 float32 freqs[3];
00380
00381
00382 for (j = 0; j < 3; ++j) {
00383 if (mel_fb->doublewide)
00384 freqs[j] = fe_melinv(mel_fb, (i + j * 2) * melbw + melmin);
00385 else
00386 freqs[j] = fe_melinv(mel_fb, (i + j) * melbw + melmin);
00387
00388 if (mel_fb->round_filters)
00389 freqs[j] = ((int)(freqs[j] / fftfreq + 0.5)) * fftfreq;
00390 }
00391
00392 for (j = 0; j < mel_fb->filt_width[i]; ++j) {
00393 float32 hz, loslope, hislope;
00394
00395 hz = (mel_fb->spec_start[i] + j) * fftfreq;
00396 if (hz < freqs[0] || hz > freqs[2]) {
00397 E_FATAL("WTF, %f < %f > %f\n", freqs[0], hz, freqs[2]);
00398 }
00399 loslope = (hz - freqs[0]) / (freqs[1] - freqs[0]);
00400 hislope = (freqs[2] - hz) / (freqs[2] - freqs[1]);
00401 if (mel_fb->unit_area) {
00402 loslope *= 2 / (freqs[2] - freqs[0]);
00403 hislope *= 2 / (freqs[2] - freqs[0]);
00404 }
00405 if (loslope < hislope) {
00406 #ifdef FIXED_POINT
00407 mel_fb->filt_coeffs[n_coeffs] = fe_log(loslope);
00408 #else
00409 mel_fb->filt_coeffs[n_coeffs] = loslope;
00410 #endif
00411 }
00412 else {
00413 #ifdef FIXED_POINT
00414 mel_fb->filt_coeffs[n_coeffs] = fe_log(hislope);
00415 #else
00416 mel_fb->filt_coeffs[n_coeffs] = hislope;
00417 #endif
00418 }
00419 ++n_coeffs;
00420 }
00421 }
00422
00423
00424 return FE_SUCCESS;
00425 }
00426
00427 int32
00428 fe_compute_melcosine(melfb_t * mel_fb)
00429 {
00430
00431 float64 freqstep;
00432 int32 i, j;
00433
00434 mel_fb->mel_cosine =
00435 (mfcc_t **) ckd_calloc_2d(mel_fb->num_cepstra,
00436 mel_fb->num_filters,
00437 sizeof(mfcc_t));
00438
00439 freqstep = M_PI / mel_fb->num_filters;
00440
00441
00442 for (i = 0; i < mel_fb->num_cepstra; i++) {
00443 for (j = 0; j < mel_fb->num_filters; j++) {
00444 float64 cosine;
00445
00446 cosine = cos(freqstep * i * (j + 0.5));
00447 mel_fb->mel_cosine[i][j] = FLOAT2COS(cosine);
00448 }
00449 }
00450
00451
00452 mel_fb->sqrt_inv_n = FLOAT2COS(sqrt(1.0 / mel_fb->num_filters));
00453 mel_fb->sqrt_inv_2n = FLOAT2COS(sqrt(2.0 / mel_fb->num_filters));
00454
00455
00456 if (mel_fb->lifter_val) {
00457 mel_fb->lifter = calloc(mel_fb->num_cepstra, sizeof(*mel_fb->lifter));
00458 for (i = 0; i < mel_fb->num_cepstra; ++i) {
00459 mel_fb->lifter[i] = FLOAT2MFCC(1 + mel_fb->lifter_val / 2
00460 * sin(i * M_PI / mel_fb->lifter_val));
00461 }
00462 }
00463
00464 return (0);
00465 }
00466
00467 static void
00468 fe_pre_emphasis(int16 const *in, frame_t * out, int32 len,
00469 float32 factor, int16 prior)
00470 {
00471 int i;
00472
00473 #if defined(FIXED16)
00474 int16 fxd_alpha = (int16)(factor * 0x8000);
00475 int32 tmp1, tmp2;
00476
00477 tmp1 = (int32)in[0] << 15;
00478 tmp2 = (int32)prior * fxd_alpha;
00479 out[0] = (int16)((tmp1 - tmp2) >> 15);
00480 for (i = 1; i < len; ++i) {
00481 tmp1 = (int32)in[i] << 15;
00482 tmp2 = (int32)in[i-1] * fxd_alpha;
00483 out[i] = (int16)((tmp1 - tmp2) >> 15);
00484 }
00485 #elif defined(FIXED_POINT)
00486 fixed32 fxd_alpha = FLOAT2FIX(factor);
00487 out[0] = ((fixed32)in[0] << DEFAULT_RADIX) - (prior * fxd_alpha);
00488 for (i = 1; i < len; ++i)
00489 out[i] = ((fixed32)in[i] << DEFAULT_RADIX)
00490 - (fixed32)in[i-1] * fxd_alpha;
00491 #else
00492 out[0] = (frame_t) in[0] - (frame_t) prior * factor;
00493 for (i = 1; i < len; i++)
00494 out[i] = (frame_t) in[i] - (frame_t) in[i-1] * factor;
00495 #endif
00496 }
00497
00498 static void
00499 fe_short_to_frame(int16 const *in, frame_t * out, int32 len)
00500 {
00501 int i;
00502
00503 #if defined(FIXED16)
00504 memcpy(out, in, len * sizeof(*out));
00505 #elif defined(FIXED_POINT)
00506 for (i = 0; i < len; i++)
00507 out[i] = (int32) in[i] << DEFAULT_RADIX;
00508 #else
00509 for (i = 0; i < len; i++)
00510 out[i] = (frame_t) in[i];
00511 #endif
00512 }
00513
00514 void
00515 fe_create_hamming(window_t * in, int32 in_len)
00516 {
00517 int i;
00518
00519
00520 for (i = 0; i < in_len / 2; i++) {
00521 float64 hamm;
00522 hamm = (0.54 - 0.46 * cos(2 * M_PI * i /
00523 ((float64) in_len - 1.0)));
00524 #ifdef FIXED16
00525 in[i] = (int16)(hamm * 0x8000);
00526 #else
00527 in[i] = FLOAT2COS(hamm);
00528 #endif
00529 }
00530 }
00531
00532 static void
00533 fe_hamming_window(frame_t * in, window_t * window, int32 in_len, int32 remove_dc)
00534 {
00535 int i;
00536
00537 if (remove_dc) {
00538 #ifdef FIXED16
00539 int32 mean = 0;
00540 #else
00541 frame_t mean = 0;
00542 #endif
00543
00544 for (i = 0; i < in_len; i++)
00545 mean += in[i];
00546 mean /= in_len;
00547 for (i = 0; i < in_len; i++)
00548 in[i] -= (frame_t)mean;
00549 }
00550
00551 #ifdef FIXED16
00552 for (i = 0; i < in_len/2; i++) {
00553 int32 tmp1, tmp2;
00554
00555 tmp1 = (int32)in[i] * window[i];
00556 tmp2 = (int32)in[in_len-1-i] * window[i];
00557 in[i] = (int16)(tmp1 >> 15);
00558 in[in_len-1-i] = (int16)(tmp2 >> 15);
00559 }
00560 #else
00561 for (i = 0; i < in_len/2; i++) {
00562 in[i] = COSMUL(in[i], window[i]);
00563 in[in_len-1-i] = COSMUL(in[in_len-1-i], window[i]);
00564 }
00565 #endif
00566 }
00567
00568 static int
00569 fe_spch_to_frame(fe_t *fe, int len)
00570 {
00571
00572 if (fe->pre_emphasis_alpha != 0.0) {
00573 fe_pre_emphasis(fe->spch, fe->frame, len,
00574 fe->pre_emphasis_alpha, fe->prior);
00575 if (len >= fe->frame_shift)
00576 fe->prior = fe->spch[fe->frame_shift - 1];
00577 else
00578 fe->prior = fe->spch[len - 1];
00579 }
00580 else
00581 fe_short_to_frame(fe->spch, fe->frame, len);
00582
00583
00584 memset(fe->frame + len, 0,
00585 (fe->fft_size - len) * sizeof(*fe->frame));
00586
00587
00588 fe_hamming_window(fe->frame, fe->hamming_window, fe->frame_size,
00589 fe->remove_dc);
00590
00591 return len;
00592 }
00593
00594 int
00595 fe_read_frame(fe_t *fe, int16 const *in, int32 len)
00596 {
00597 int i;
00598
00599 if (len > fe->frame_size)
00600 len = fe->frame_size;
00601
00602
00603 memcpy(fe->spch, in, len * sizeof(*in));
00604
00605 if (fe->swap)
00606 for (i = 0; i < len; ++i)
00607 SWAP_INT16(&fe->spch[i]);
00608 if (fe->dither)
00609 for (i = 0; i < len; ++i)
00610 fe->spch[i] += (int16) ((!(s3_rand_int31() % 4)) ? 1 : 0);
00611
00612 return fe_spch_to_frame(fe, len);
00613 }
00614
00615 int
00616 fe_shift_frame(fe_t *fe, int16 const *in, int32 len)
00617 {
00618 int offset, i;
00619
00620 if (len > fe->frame_shift)
00621 len = fe->frame_shift;
00622 offset = fe->frame_size - fe->frame_shift;
00623
00624
00625 memmove(fe->spch, fe->spch + fe->frame_shift,
00626 offset * sizeof(*fe->spch));
00627 memcpy(fe->spch + offset, in, len * sizeof(*fe->spch));
00628
00629 if (fe->swap)
00630 for (i = 0; i < len; ++i)
00631 SWAP_INT16(&fe->spch[offset + i]);
00632 if (fe->dither)
00633 for (i = 0; i < len; ++i)
00634 fe->spch[offset + i]
00635 += (int16) ((!(s3_rand_int31() % 4)) ? 1 : 0);
00636
00637 return fe_spch_to_frame(fe, offset + len);
00638 }
00639
00643 void
00644 fe_create_twiddle(fe_t *fe)
00645 {
00646 int i;
00647
00648 for (i = 0; i < fe->fft_size / 4; ++i) {
00649 float64 a = 2 * M_PI * i / fe->fft_size;
00650 #ifdef FIXED16
00651 fe->ccc[i] = (int16)(cos(a) * 0x8000);
00652 fe->sss[i] = (int16)(sin(a) * 0x8000);
00653 #elif defined(FIXED_POINT)
00654 fe->ccc[i] = FLOAT2COS(cos(a));
00655 fe->sss[i] = FLOAT2COS(sin(a));
00656 #else
00657 fe->ccc[i] = cos(a);
00658 fe->sss[i] = sin(a);
00659 #endif
00660 }
00661 }
00662
00663
00664
00665
00666
00667
00668
00669 #if defined(FIXED16)
00670 static int
00671 fe_fft_real(fe_t *fe)
00672 {
00673 int i, j, k, m, n, lz;
00674 frame_t *x, xt, max;
00675
00676 x = fe->frame;
00677 m = fe->fft_order;
00678 n = fe->fft_size;
00679
00680
00681 j = 0;
00682 for (i = 0; i < n - 1; ++i) {
00683 if (i < j) {
00684 xt = x[j];
00685 x[j] = x[i];
00686 x[i] = xt;
00687 }
00688 k = n / 2;
00689 while (k <= j) {
00690 j -= k;
00691 k /= 2;
00692 }
00693 j += k;
00694 }
00695
00696 max = 0;
00697 for (i = 0; i < n; ++i)
00698 if (abs(x[i]) > max)
00699 max = abs(x[i]);
00700
00701
00702
00703 for (lz = 0; lz < m; ++lz)
00704 if (max & (1 << (15-lz)))
00705 break;
00706
00707
00708
00709
00710
00711
00712
00713
00714 for (i = 0; i < n; i += 2) {
00715 int atten = (lz == 0);
00716 xt = x[i] >> atten;
00717 x[i] = xt + (x[i + 1] >> atten);
00718 x[i + 1] = xt - (x[i + 1] >> atten);
00719 }
00720
00721
00722 for (k = 1; k < m; ++k) {
00723 int n1, n2, n4;
00724
00725 int atten = (k >= lz);
00726
00727 n4 = k - 1;
00728 n2 = k;
00729 n1 = k + 1;
00730
00731 for (i = 0; i < n; i += (1 << n1)) {
00732
00733
00734
00735
00736 xt = x[i] >> atten;
00737 x[i] = xt + (x[i + (1 << n2)] >> atten);
00738 x[i + (1 << n2)] = xt - (x[i + (1 << n2)] >> atten);
00739
00740
00741
00742
00743
00744
00745
00746 x[i + (1 << n2) + (1 << n4)] = -x[i + (1 << n2) + (1 << n4)] >> atten;
00747 x[i + (1 << n4)] = x[i + (1 << n4)] >> atten;
00748
00749
00750
00751
00752 for (j = 1; j < (1 << n4); ++j) {
00753 frame_t cc, ss, t1, t2;
00754 int i1, i2, i3, i4;
00755
00756 i1 = i + j;
00757 i2 = i + (1 << n2) - j;
00758 i3 = i + (1 << n2) + j;
00759 i4 = i + (1 << n2) + (1 << n2) - j;
00760
00761
00762
00763
00764
00765 cc = fe->ccc[j << (m - n1)];
00766 ss = fe->sss[j << (m - n1)];
00767
00768
00769
00770 {
00771 int32 tmp1, tmp2;
00772 tmp1 = (int32)x[i3] * cc + (int32)x[i4] * ss;
00773 tmp2 = (int32)x[i3] * ss - (int32)x[i4] * cc;
00774 t1 = (int16)(tmp1 >> 15) >> atten;
00775 t2 = (int16)(tmp2 >> 15) >> atten;
00776 }
00777
00778 x[i4] = (x[i2] >> atten) - t2;
00779 x[i3] = (-x[i2] >> atten) - t2;
00780 x[i2] = (x[i1] >> atten) - t1;
00781 x[i1] = (x[i1] >> atten) + t1;
00782 }
00783 }
00784 }
00785
00786
00787 return lz;
00788 }
00789 #else
00790 static int
00791 fe_fft_real(fe_t *fe)
00792 {
00793 int i, j, k, m, n;
00794 frame_t *x, xt;
00795
00796 x = fe->frame;
00797 m = fe->fft_order;
00798 n = fe->fft_size;
00799
00800
00801 j = 0;
00802 for (i = 0; i < n - 1; ++i) {
00803 if (i < j) {
00804 xt = x[j];
00805 x[j] = x[i];
00806 x[i] = xt;
00807 }
00808 k = n / 2;
00809 while (k <= j) {
00810 j -= k;
00811 k /= 2;
00812 }
00813 j += k;
00814 }
00815
00816
00817
00818
00819
00820 for (i = 0; i < n; i += 2) {
00821 xt = x[i];
00822 x[i] = (xt + x[i + 1]);
00823 x[i + 1] = (xt - x[i + 1]);
00824 }
00825
00826
00827 for (k = 1; k < m; ++k) {
00828 int n1, n2, n4;
00829
00830 n4 = k - 1;
00831 n2 = k;
00832 n1 = k + 1;
00833
00834 for (i = 0; i < n; i += (1 << n1)) {
00835
00836
00837
00838
00839 xt = x[i];
00840 x[i] = (xt + x[i + (1 << n2)]);
00841 x[i + (1 << n2)] = (xt - x[i + (1 << n2)]);
00842
00843
00844
00845
00846
00847
00848
00849 x[i + (1 << n2) + (1 << n4)] = -x[i + (1 << n2) + (1 << n4)];
00850 x[i + (1 << n4)] = x[i + (1 << n4)];
00851
00852
00853
00854
00855 for (j = 1; j < (1 << n4); ++j) {
00856 frame_t cc, ss, t1, t2;
00857 int i1, i2, i3, i4;
00858
00859 i1 = i + j;
00860 i2 = i + (1 << n2) - j;
00861 i3 = i + (1 << n2) + j;
00862 i4 = i + (1 << n2) + (1 << n2) - j;
00863
00864
00865
00866
00867
00868 cc = fe->ccc[j << (m - n1)];
00869 ss = fe->sss[j << (m - n1)];
00870
00871
00872
00873 t1 = COSMUL(x[i3], cc) + COSMUL(x[i4], ss);
00874 t2 = COSMUL(x[i3], ss) - COSMUL(x[i4], cc);
00875
00876 x[i4] = (x[i2] - t2);
00877 x[i3] = (-x[i2] - t2);
00878 x[i2] = (x[i1] - t1);
00879 x[i1] = (x[i1] + t1);
00880 }
00881 }
00882 }
00883
00884
00885 return m;
00886 }
00887 #endif
00888
00889 static void
00890 fe_spec_magnitude(fe_t *fe)
00891 {
00892 frame_t *fft;
00893 powspec_t *spec;
00894 int32 j, scale, fftsize;
00895
00896
00897
00898 scale = fe_fft_real(fe);
00899
00900
00901 fft = fe->frame;
00902 spec = fe->spec;
00903 fftsize = fe->fft_size;
00904
00905
00906 scale = fe->fft_order - scale;
00907
00908
00909 {
00910 #ifdef FIXED16
00911 spec[0] = fixlog(abs(fft[0]) << scale) * 2;
00912 #elif defined(FIXED_POINT)
00913 spec[0] = FIXLN(abs(fft[0]) << scale) * 2;
00914 #else
00915 spec[0] = fft[0] * fft[0];
00916 #endif
00917 }
00918
00919 for (j = 1; j <= fftsize / 2; j++) {
00920 #ifdef FIXED16
00921 int32 rr = fixlog(abs(fft[j]) << scale) * 2;
00922 int32 ii = fixlog(abs(fft[fftsize - j]) << scale) * 2;
00923 spec[j] = fe_log_add(rr, ii);
00924 #elif defined(FIXED_POINT)
00925 int32 rr = FIXLN(abs(fft[j]) << scale) * 2;
00926 int32 ii = FIXLN(abs(fft[fftsize - j]) << scale) * 2;
00927 spec[j] = fe_log_add(rr, ii);
00928 #else
00929 spec[j] = fft[j] * fft[j] + fft[fftsize - j] * fft[fftsize - j];
00930 #endif
00931 }
00932 }
00933
00934 static void
00935 fe_mel_spec(fe_t * fe)
00936 {
00937 int whichfilt;
00938 powspec_t *spec, *mfspec;
00939
00940
00941 spec = fe->spec;
00942 mfspec = fe->mfspec;
00943
00944 for (whichfilt = 0; whichfilt < fe->mel_fb->num_filters; whichfilt++) {
00945 int spec_start, filt_start, i;
00946
00947 spec_start = fe->mel_fb->spec_start[whichfilt];
00948 filt_start = fe->mel_fb->filt_start[whichfilt];
00949
00950 #ifdef FIXED_POINT
00951 mfspec[whichfilt] = spec[spec_start] + fe->mel_fb->filt_coeffs[filt_start];
00952 for (i = 1; i < fe->mel_fb->filt_width[whichfilt]; i++) {
00953 mfspec[whichfilt] = fe_log_add(mfspec[whichfilt],
00954 spec[spec_start + i] +
00955 fe->mel_fb->filt_coeffs[filt_start + i]);
00956 }
00957 #else
00958 mfspec[whichfilt] = 0;
00959 for (i = 0; i < fe->mel_fb->filt_width[whichfilt]; i++)
00960 mfspec[whichfilt] +=
00961 spec[spec_start + i] * fe->mel_fb->filt_coeffs[filt_start + i];
00962 #endif
00963 }
00964 }
00965
00966 static void
00967 fe_mel_cep(fe_t * fe, mfcc_t *mfcep)
00968 {
00969 int32 i;
00970 powspec_t *mfspec;
00971
00972
00973 mfspec = fe->mfspec;
00974
00975 for (i = 0; i < fe->mel_fb->num_filters; ++i) {
00976 #ifndef FIXED_POINT
00977 if (mfspec[i] > 0)
00978 mfspec[i] = log(mfspec[i]);
00979 else
00980
00981
00982
00983
00984 mfspec[i] = -10.0;
00985 #endif
00986 }
00987
00988
00989 if (fe->log_spec == RAW_LOG_SPEC) {
00990 for (i = 0; i < fe->feature_dimension; i++) {
00991 mfcep[i] = (mfcc_t) mfspec[i];
00992 }
00993 }
00994
00995 else if (fe->log_spec == SMOOTH_LOG_SPEC) {
00996
00997 fe_dct2(fe, mfspec, mfcep, 0);
00998 fe_dct3(fe, mfcep, mfspec);
00999 for (i = 0; i < fe->feature_dimension; i++) {
01000 mfcep[i] = (mfcc_t) mfspec[i];
01001 }
01002 }
01003 else if (fe->transform == DCT_II)
01004 fe_dct2(fe, mfspec, mfcep, 0);
01005 else if (fe->transform == DCT_HTK)
01006 fe_dct2(fe, mfspec, mfcep, 1);
01007 else
01008 fe_spec2cep(fe, mfspec, mfcep);
01009
01010 return;
01011 }
01012
01013 void
01014 fe_spec2cep(fe_t * fe, const powspec_t * mflogspec, mfcc_t * mfcep)
01015 {
01016 int32 i, j, beta;
01017
01018
01019
01020 mfcep[0] = mflogspec[0] / 2;
01021 for (j = 1; j < fe->mel_fb->num_filters; j++)
01022 mfcep[0] += mflogspec[j];
01023 mfcep[0] /= (frame_t) fe->mel_fb->num_filters;
01024
01025 for (i = 1; i < fe->num_cepstra; ++i) {
01026 mfcep[i] = 0;
01027 for (j = 0; j < fe->mel_fb->num_filters; j++) {
01028 if (j == 0)
01029 beta = 1;
01030 else
01031 beta = 2;
01032 mfcep[i] += COSMUL(mflogspec[j],
01033 fe->mel_fb->mel_cosine[i][j]) * beta;
01034 }
01035
01036
01037
01038 mfcep[i] /= (frame_t) fe->mel_fb->num_filters * 2;
01039 }
01040 }
01041
01042 void
01043 fe_dct2(fe_t * fe, const powspec_t * mflogspec, mfcc_t * mfcep, int htk)
01044 {
01045 int32 i, j;
01046
01047
01048
01049 mfcep[0] = mflogspec[0];
01050 for (j = 1; j < fe->mel_fb->num_filters; j++)
01051 mfcep[0] += mflogspec[j];
01052 if (htk)
01053 mfcep[0] = COSMUL(mfcep[0], fe->mel_fb->sqrt_inv_2n);
01054 else
01055 mfcep[0] = COSMUL(mfcep[0], fe->mel_fb->sqrt_inv_n);
01056
01057 for (i = 1; i < fe->num_cepstra; ++i) {
01058 mfcep[i] = 0;
01059 for (j = 0; j < fe->mel_fb->num_filters; j++) {
01060 mfcep[i] += COSMUL(mflogspec[j],
01061 fe->mel_fb->mel_cosine[i][j]);
01062 }
01063 mfcep[i] = COSMUL(mfcep[i], fe->mel_fb->sqrt_inv_2n);
01064 }
01065 }
01066
01067 void
01068 fe_lifter(fe_t *fe, mfcc_t *mfcep)
01069 {
01070 int32 i;
01071
01072 if (fe->mel_fb->lifter_val == 0)
01073 return;
01074
01075 for (i = 0; i < fe->num_cepstra; ++i) {
01076 mfcep[i] = MFCCMUL(mfcep[i], fe->mel_fb->lifter[i]);
01077 }
01078 }
01079
01080 void
01081 fe_dct3(fe_t * fe, const mfcc_t * mfcep, powspec_t * mflogspec)
01082 {
01083 int32 i, j;
01084
01085 for (i = 0; i < fe->mel_fb->num_filters; ++i) {
01086 mflogspec[i] = COSMUL(mfcep[0], SQRT_HALF);
01087 for (j = 1; j < fe->num_cepstra; j++) {
01088 mflogspec[i] += COSMUL(mfcep[j],
01089 fe->mel_fb->mel_cosine[j][i]);
01090 }
01091 mflogspec[i] = COSMUL(mflogspec[i], fe->mel_fb->sqrt_inv_2n);
01092 }
01093 }
01094
01095 int32
01096 fe_write_frame(fe_t * fe, mfcc_t * fea)
01097 {
01098 fe_spec_magnitude(fe);
01099 fe_mel_spec(fe);
01100 fe_mel_cep(fe, fea);
01101 fe_lifter(fe, fea);
01102
01103 return 0;
01104 }
01105
01106 void *
01107 fe_create_2d(int32 d1, int32 d2, int32 elem_size)
01108 {
01109 return (void *)ckd_calloc_2d(d1, d2, elem_size);
01110 }
01111
01112 void
01113 fe_free_2d(void *arr)
01114 {
01115 ckd_free_2d((void **)arr);
01116 }