• Main Page
  • Related Pages
  • Data Structures
  • Files
  • File List
  • Globals

src/libsphinxbase/fe/fe_sigproc.c

00001 /* -*- c-basic-offset: 4; indent-tabs-mode: nil -*- */
00002 /* ====================================================================
00003  * Copyright (c) 1996-2004 Carnegie Mellon University.  All rights 
00004  * reserved.
00005  *
00006  * Redistribution and use in source and binary forms, with or without
00007  * modification, are permitted provided that the following conditions
00008  * are met:
00009  *
00010  * 1. Redistributions of source code must retain the above copyright
00011  *    notice, this list of conditions and the following disclaimer. 
00012  *
00013  * 2. Redistributions in binary form must reproduce the above copyright
00014  *    notice, this list of conditions and the following disclaimer in
00015  *    the documentation and/or other materials provided with the
00016  *    distribution.
00017  *
00018  * This work was supported in part by funding from the Defense Advanced 
00019  * Research Projects Agency and the National Science Foundation of the 
00020  * United States of America, and the CMU Sphinx Speech Consortium.
00021  *
00022  * THIS SOFTWARE IS PROVIDED BY CARNEGIE MELLON UNIVERSITY ``AS IS'' AND 
00023  * ANY EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 
00024  * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00025  * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL CARNEGIE MELLON UNIVERSITY
00026  * NOR ITS EMPLOYEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
00027  * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 
00028  * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 
00029  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 
00030  * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 
00031  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 
00032  * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
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 /* Use extra precision for cosines, Hamming window, pre-emphasis
00063  * coefficient, twiddle factors. */
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 /* Internal log-addition table for natural log with radix point at 8
00074  * bits.  Each entry is 256 * log(1 + e^{-n/256}).  This is used in the
00075  * log-add computation:
00076  *
00077  * e^z = e^x + e^y
00078  * e^z = e^x(1 + e^{y-x})     = e^y(1 + e^{x-y})
00079  * z   = x + log(1 + e^{y-x}) = y + log(1 + e^{x-y})
00080  *
00081  * So when y > x, z = y + logadd_table[-(x-y)]
00082  *    when x > y, z = x + logadd_table[-(y-x)]
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         printf("%d + %d = %d | %f + %f = %f | %f + %f = %f\n",
00267                x, y, r, FIX2FLOAT(x), FIX2FLOAT(y), FIX2FLOAT(r),
00268                exp(FIX2FLOAT(x)), exp(FIX2FLOAT(y)), exp(FIX2FLOAT(r)));
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     /* Filter coefficient matrix, in flattened form. */
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     /* First calculate the widths of each filter. */
00313     /* Minimum and maximum frequencies in mel scale. */
00314     melmin = fe_mel(mel_fb, mel_fb->lower_filt_freq);
00315     melmax = fe_mel(mel_fb, mel_fb->upper_filt_freq);
00316 
00317     /* Width of filters in mel scale */
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     /* DFT point spacing */
00335     fftfreq = mel_fb->sampling_rate / (float32) mel_fb->fft_size;
00336 
00337     /* Count and place filter coefficients. */
00338     n_coeffs = 0;
00339     for (i = 0; i < mel_fb->num_filters; ++i) {
00340         float32 freqs[3];
00341 
00342         /* Left, center, right frequencies in Hertz */
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             /* Round them to DFT points if requested */
00349             if (mel_fb->round_filters)
00350                 freqs[j] = ((int)(freqs[j] / fftfreq + 0.5)) * fftfreq;
00351         }
00352 
00353         /* spec_start is the start of this filter in the power spectrum. */
00354         mel_fb->spec_start[i] = -1;
00355         /* There must be a better way... */
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                 /* filt_width is the width in DFT points of this filter. */
00362                 mel_fb->filt_width[i] = j - mel_fb->spec_start[i];
00363                 /* filt_start is the start of this filter in the filt_coeffs array. */
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     /* Now go back and allocate the coefficient array. */
00374     mel_fb->filt_coeffs = ckd_malloc(n_coeffs * sizeof(*mel_fb->filt_coeffs));
00375 
00376     /* And now generate the coefficients. */
00377     n_coeffs = 0;
00378     for (i = 0; i < mel_fb->num_filters; ++i) {
00379         float32 freqs[3];
00380 
00381         /* Left, center, right frequencies in Hertz */
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             /* Round them to DFT points if requested */
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     /* NOTE: The first row vector is actually unnecessary but we leave
00441      * it in to avoid confusion. */
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     /* Also precompute normalization constants for unitary DCT. */
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     /* And liftering weights */
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                           /* FIXED_POINT */
00509     for (i = 0; i < len; i++)
00510         out[i] = (frame_t) in[i];
00511 #endif                          /* FIXED_POINT */
00512 }
00513 
00514 void
00515 fe_create_hamming(window_t * in, int32 in_len)
00516 {
00517     int i;
00518 
00519     /* Symmetric, so we only create the first half of it. */
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; /* Use int32 to avoid possibility of overflow */
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     /* Copy to the frame buffer. */
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     /* Zero pad up to FFT size. */
00584     memset(fe->frame + len, 0,
00585            (fe->fft_size - len) * sizeof(*fe->frame));
00586 
00587     /* Window. */
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     /* Read it into the raw speech buffer. */
00603     memcpy(fe->spch, in, len * sizeof(*in));
00604     /* Swap and dither if necessary. */
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     /* Shift data into the raw speech buffer. */
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     /* Swap and dither if necessary. */
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 /* Translated from the FORTRAN (obviously) from "Real-Valued Fast
00664  * Fourier Transform Algorithms" by Henrik V. Sorensen et al., IEEE
00665  * Transactions on Acoustics, Speech, and Signal Processing, vol. 35,
00666  * no.6.  The 16-bit version does a version of "block floating
00667  * point" in order to avoid rounding errors.
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     /* Bit-reverse the input. */
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     /* Determine how many bits of dynamic range are in the input. */
00696     max = 0;
00697     for (i = 0; i < n; ++i)
00698         if (abs(x[i]) > max)
00699             max = abs(x[i]);
00700     /* The FFT has a gain of M bits, so we need to attenuate the input
00701      * by M bits minus the number of leading zeroes in the input's
00702      * range in order to avoid overflows.  */
00703     for (lz = 0; lz < m; ++lz)
00704         if (max & (1 << (15-lz)))
00705             break;
00706 
00707     /* Basic butterflies (2-point FFT, real twiddle factors):
00708      * x[i]   = x[i] +  1 * x[i+1]
00709      * x[i+1] = x[i] + -1 * x[i+1]
00710      */
00711     /* The quantization error introduced by attenuating the input at
00712      * any given stage of the FFT has a cascading effect, so we hold
00713      * off on it until it's absolutely necessary. */
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     /* The rest of the butterflies, in stages from 1..m */
00722     for (k = 1; k < m; ++k) {
00723         int n1, n2, n4;
00724         /* Start attenuating once we hit the number of leading zeros. */
00725         int atten = (k >= lz);
00726 
00727         n4 = k - 1;
00728         n2 = k;
00729         n1 = k + 1;
00730         /* Stride over each (1 << (k+1)) points */
00731         for (i = 0; i < n; i += (1 << n1)) {
00732             /* Basic butterfly with real twiddle factors:
00733              * x[i]          = x[i] +  1 * x[i + (1<<k)]
00734              * x[i + (1<<k)] = x[i] + -1 * x[i + (1<<k)]
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             /* The other ones with real twiddle factors:
00741              * x[i + (1<<k) + (1<<(k-1))]
00742              *   = 0 * x[i + (1<<k-1)] + -1 * x[i + (1<<k) + (1<<k-1)]
00743              * x[i + (1<<(k-1))]
00744              *   = 1 * x[i + (1<<k-1)] +  0 * x[i + (1<<k) + (1<<k-1)]
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             /* Butterflies with complex twiddle factors.
00750              * There are (1<<k-1) of them.
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                  * cc = real(W[j * n / (1<<(k+1))])
00763                  * ss = imag(W[j * n / (1<<(k+1))])
00764                  */
00765                 cc = fe->ccc[j << (m - n1)];
00766                 ss = fe->sss[j << (m - n1)];
00767 
00768                 /* There are some symmetry properties which allow us
00769                  * to get away with only four multiplications here. */
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     /* Return the residual scaling factor. */
00787     return lz;
00788 }
00789 #else /* !FIXED16 */
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     /* Bit-reverse the input. */
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     /* Basic butterflies (2-point FFT, real twiddle factors):
00817      * x[i]   = x[i] +  1 * x[i+1]
00818      * x[i+1] = x[i] + -1 * x[i+1]
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     /* The rest of the butterflies, in stages from 1..m */
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         /* Stride over each (1 << (k+1)) points */
00834         for (i = 0; i < n; i += (1 << n1)) {
00835             /* Basic butterfly with real twiddle factors:
00836              * x[i]          = x[i] +  1 * x[i + (1<<k)]
00837              * x[i + (1<<k)] = x[i] + -1 * x[i + (1<<k)]
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             /* The other ones with real twiddle factors:
00844              * x[i + (1<<k) + (1<<(k-1))]
00845              *   = 0 * x[i + (1<<k-1)] + -1 * x[i + (1<<k) + (1<<k-1)]
00846              * x[i + (1<<(k-1))]
00847              *   = 1 * x[i + (1<<k-1)] +  0 * x[i + (1<<k) + (1<<k-1)]
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             /* Butterflies with complex twiddle factors.
00853              * There are (1<<k-1) of them.
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                  * cc = real(W[j * n / (1<<(k+1))])
00866                  * ss = imag(W[j * n / (1<<(k+1))])
00867                  */
00868                 cc = fe->ccc[j << (m - n1)];
00869                 ss = fe->sss[j << (m - n1)];
00870 
00871                 /* There are some symmetry properties which allow us
00872                  * to get away with only four multiplications here. */
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     /* This isn't used, but return it for completeness. */
00885     return m;
00886 }
00887 #endif /* !FIXED16 */
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     /* Do FFT and get the scaling factor back (only actually used in
00897      * fixed-point).  Note the scaling factor is expressed in bits. */
00898     scale = fe_fft_real(fe);
00899 
00900     /* Convenience pointers to make things less awkward below. */
00901     fft = fe->frame;
00902     spec = fe->spec;
00903     fftsize = fe->fft_size;
00904 
00905     /* We need to scale things up the rest of the way to N. */
00906     scale = fe->fft_order - scale;
00907 
00908     /* The first point (DC coefficient) has no imaginary part */
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     /* Convenience poitners. */
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                           /* !FIXED_POINT */
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                          /* !FIXED_POINT */
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     /* Convenience pointer. */
00973     mfspec = fe->mfspec;
00974 
00975     for (i = 0; i < fe->mel_fb->num_filters; ++i) {
00976 #ifndef FIXED_POINT /* It's already in log domain for fixed point */
00977         if (mfspec[i] > 0)
00978             mfspec[i] = log(mfspec[i]);
00979         else                    /* This number should be smaller than anything
00980                                  * else, but not too small, so as to avoid
00981                                  * infinities in the inverse transform (this is
00982                                  * the frequency-domain equivalent of
00983                                  * dithering) */
00984             mfspec[i] = -10.0;
00985 #endif                          /* !FIXED_POINT */
00986     }
00987 
00988     /* If we are doing LOG_SPEC, then do nothing. */
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     /* For smoothed spectrum, do DCT-II followed by (its inverse) DCT-III */
00995     else if (fe->log_spec == SMOOTH_LOG_SPEC) {
00996         /* FIXME: This is probably broken for fixed-point. */
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, FALSE);
01005     else if (fe->transform == DCT_HTK)
01006         fe_dct2(fe, mfspec, mfcep, TRUE);
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     /* Compute C0 separately (its basis vector is 1) to avoid
01019      * costly multiplications. */
01020     mfcep[0] = mflogspec[0] / 2; /* beta = 0.5 */
01021     for (j = 1; j < fe->mel_fb->num_filters; j++)
01022         mfcep[0] += mflogspec[j]; /* beta = 1.0 */
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;       /* 0.5 */
01030             else
01031                 beta = 2;       /* 1.0 */
01032             mfcep[i] += COSMUL(mflogspec[j],
01033                                fe->mel_fb->mel_cosine[i][j]) * beta;
01034         }
01035         /* Note that this actually normalizes by num_filters, like the
01036          * original Sphinx front-end, due to the doubled 'beta' factor
01037          * above.  */
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     /* Compute C0 separately (its basis vector is 1) to avoid
01048      * costly multiplications. */
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 /* sqrt(1/N) = sqrt(2/N) * 1/sqrt(2) */
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 1;
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 }

Generated on Mon Aug 29 2011 for SphinxBase by  doxygen 1.7.1