SphinxBase  0.6
fe_sigproc.c
1 /* -*- c-basic-offset: 4; indent-tabs-mode: nil -*- */
2 /* ====================================================================
3  * Copyright (c) 1996-2004 Carnegie Mellon University. All rights
4  * reserved.
5  *
6  * Redistribution and use in source and binary forms, with or without
7  * modification, are permitted provided that the following conditions
8  * are met:
9  *
10  * 1. Redistributions of source code must retain the above copyright
11  * notice, this list of conditions and the following disclaimer.
12  *
13  * 2. Redistributions in binary form must reproduce the above copyright
14  * notice, this list of conditions and the following disclaimer in
15  * the documentation and/or other materials provided with the
16  * distribution.
17  *
18  * This work was supported in part by funding from the Defense Advanced
19  * Research Projects Agency and the National Science Foundation of the
20  * United States of America, and the CMU Sphinx Speech Consortium.
21  *
22  * THIS SOFTWARE IS PROVIDED BY CARNEGIE MELLON UNIVERSITY ``AS IS'' AND
23  * ANY EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
24  * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
25  * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL CARNEGIE MELLON UNIVERSITY
26  * NOR ITS EMPLOYEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
27  * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
28  * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
29  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
30  * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
31  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
32  * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33  *
34  * ====================================================================
35  *
36  */
37 
38 #include <stdio.h>
39 #include <math.h>
40 #include <string.h>
41 #include <stdlib.h>
42 #include <assert.h>
43 
44 #ifdef HAVE_CONFIG_H
45 #include <config.h>
46 #endif
47 
48 #ifdef _MSC_VER
49 #pragma warning (disable: 4244)
50 #endif
51 
52 #include "sphinxbase/prim_type.h"
53 #include "sphinxbase/ckd_alloc.h"
54 #include "sphinxbase/byteorder.h"
55 #include "sphinxbase/fixpoint.h"
56 #include "sphinxbase/fe.h"
57 #include "sphinxbase/genrand.h"
58 #include "sphinxbase/libutil.h"
59 #include "sphinxbase/err.h"
60 
61 #include "fe_internal.h"
62 #include "fe_warp.h"
63 
64 /* Use extra precision for cosines, Hamming window, pre-emphasis
65  * coefficient, twiddle factors. */
66 #ifdef FIXED_POINT
67 #define FLOAT2COS(x) FLOAT2FIX_ANY(x,30)
68 #define COSMUL(x,y) FIXMUL_ANY(x,y,30)
69 #else
70 #define FLOAT2COS(x) (x)
71 #define COSMUL(x,y) ((x)*(y))
72 #endif
73 
74 #ifdef FIXED_POINT
75 /* Internal log-addition table for natural log with radix point at 8
76  * bits. Each entry is 256 * log(1 + e^{-n/256}). This is used in the
77  * log-add computation:
78  *
79  * e^z = e^x + e^y
80  * e^z = e^x(1 + e^{y-x}) = e^y(1 + e^{x-y})
81  * z = x + log(1 + e^{y-x}) = y + log(1 + e^{x-y})
82  *
83  * So when y > x, z = y + logadd_table[-(x-y)]
84  * when x > y, z = x + logadd_table[-(y-x)]
85  */
86 static const unsigned char fe_logadd_table[] = {
87 177, 177, 176, 176, 175, 175, 174, 174, 173, 173,
88 172, 172, 172, 171, 171, 170, 170, 169, 169, 168,
89 168, 167, 167, 166, 166, 165, 165, 164, 164, 163,
90 163, 162, 162, 161, 161, 161, 160, 160, 159, 159,
91 158, 158, 157, 157, 156, 156, 155, 155, 155, 154,
92 154, 153, 153, 152, 152, 151, 151, 151, 150, 150,
93 149, 149, 148, 148, 147, 147, 147, 146, 146, 145,
94 145, 144, 144, 144, 143, 143, 142, 142, 141, 141,
95 141, 140, 140, 139, 139, 138, 138, 138, 137, 137,
96 136, 136, 136, 135, 135, 134, 134, 134, 133, 133,
97 132, 132, 131, 131, 131, 130, 130, 129, 129, 129,
98 128, 128, 128, 127, 127, 126, 126, 126, 125, 125,
99 124, 124, 124, 123, 123, 123, 122, 122, 121, 121,
100 121, 120, 120, 119, 119, 119, 118, 118, 118, 117,
101 117, 117, 116, 116, 115, 115, 115, 114, 114, 114,
102 113, 113, 113, 112, 112, 112, 111, 111, 110, 110,
103 110, 109, 109, 109, 108, 108, 108, 107, 107, 107,
104 106, 106, 106, 105, 105, 105, 104, 104, 104, 103,
105 103, 103, 102, 102, 102, 101, 101, 101, 100, 100,
106 100, 99, 99, 99, 98, 98, 98, 97, 97, 97,
107 96, 96, 96, 96, 95, 95, 95, 94, 94, 94,
108 93, 93, 93, 92, 92, 92, 92, 91, 91, 91,
109 90, 90, 90, 89, 89, 89, 89, 88, 88, 88,
110 87, 87, 87, 87, 86, 86, 86, 85, 85, 85,
111 85, 84, 84, 84, 83, 83, 83, 83, 82, 82,
112 82, 82, 81, 81, 81, 80, 80, 80, 80, 79,
113 79, 79, 79, 78, 78, 78, 78, 77, 77, 77,
114 77, 76, 76, 76, 75, 75, 75, 75, 74, 74,
115 74, 74, 73, 73, 73, 73, 72, 72, 72, 72,
116 71, 71, 71, 71, 71, 70, 70, 70, 70, 69,
117 69, 69, 69, 68, 68, 68, 68, 67, 67, 67,
118 67, 67, 66, 66, 66, 66, 65, 65, 65, 65,
119 64, 64, 64, 64, 64, 63, 63, 63, 63, 63,
120 62, 62, 62, 62, 61, 61, 61, 61, 61, 60,
121 60, 60, 60, 60, 59, 59, 59, 59, 59, 58,
122 58, 58, 58, 58, 57, 57, 57, 57, 57, 56,
123 56, 56, 56, 56, 55, 55, 55, 55, 55, 54,
124 54, 54, 54, 54, 53, 53, 53, 53, 53, 52,
125 52, 52, 52, 52, 52, 51, 51, 51, 51, 51,
126 50, 50, 50, 50, 50, 50, 49, 49, 49, 49,
127 49, 49, 48, 48, 48, 48, 48, 48, 47, 47,
128 47, 47, 47, 47, 46, 46, 46, 46, 46, 46,
129 45, 45, 45, 45, 45, 45, 44, 44, 44, 44,
130 44, 44, 43, 43, 43, 43, 43, 43, 43, 42,
131 42, 42, 42, 42, 42, 41, 41, 41, 41, 41,
132 41, 41, 40, 40, 40, 40, 40, 40, 40, 39,
133 39, 39, 39, 39, 39, 39, 38, 38, 38, 38,
134 38, 38, 38, 37, 37, 37, 37, 37, 37, 37,
135 37, 36, 36, 36, 36, 36, 36, 36, 35, 35,
136 35, 35, 35, 35, 35, 35, 34, 34, 34, 34,
137 34, 34, 34, 34, 33, 33, 33, 33, 33, 33,
138 33, 33, 32, 32, 32, 32, 32, 32, 32, 32,
139 32, 31, 31, 31, 31, 31, 31, 31, 31, 31,
140 30, 30, 30, 30, 30, 30, 30, 30, 30, 29,
141 29, 29, 29, 29, 29, 29, 29, 29, 28, 28,
142 28, 28, 28, 28, 28, 28, 28, 28, 27, 27,
143 27, 27, 27, 27, 27, 27, 27, 27, 26, 26,
144 26, 26, 26, 26, 26, 26, 26, 26, 25, 25,
145 25, 25, 25, 25, 25, 25, 25, 25, 25, 24,
146 24, 24, 24, 24, 24, 24, 24, 24, 24, 24,
147 23, 23, 23, 23, 23, 23, 23, 23, 23, 23,
148 23, 23, 22, 22, 22, 22, 22, 22, 22, 22,
149 22, 22, 22, 22, 21, 21, 21, 21, 21, 21,
150 21, 21, 21, 21, 21, 21, 21, 20, 20, 20,
151 20, 20, 20, 20, 20, 20, 20, 20, 20, 20,
152 19, 19, 19, 19, 19, 19, 19, 19, 19, 19,
153 19, 19, 19, 19, 18, 18, 18, 18, 18, 18,
154 18, 18, 18, 18, 18, 18, 18, 18, 18, 17,
155 17, 17, 17, 17, 17, 17, 17, 17, 17, 17,
156 17, 17, 17, 17, 16, 16, 16, 16, 16, 16,
157 16, 16, 16, 16, 16, 16, 16, 16, 16, 16,
158 16, 15, 15, 15, 15, 15, 15, 15, 15, 15,
159 15, 15, 15, 15, 15, 15, 15, 15, 14, 14,
160 14, 14, 14, 14, 14, 14, 14, 14, 14, 14,
161 14, 14, 14, 14, 14, 14, 14, 13, 13, 13,
162 13, 13, 13, 13, 13, 13, 13, 13, 13, 13,
163 13, 13, 13, 13, 13, 13, 13, 12, 12, 12,
164 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
165 12, 12, 12, 12, 12, 12, 12, 12, 12, 11,
166 11, 11, 11, 11, 11, 11, 11, 11, 11, 11,
167 11, 11, 11, 11, 11, 11, 11, 11, 11, 11,
168 11, 11, 11, 10, 10, 10, 10, 10, 10, 10,
169 10, 10, 10, 10, 10, 10, 10, 10, 10, 10,
170 10, 10, 10, 10, 10, 10, 10, 10, 10, 9,
171 9, 9, 9, 9, 9, 9, 9, 9, 9, 9,
172 9, 9, 9, 9, 9, 9, 9, 9, 9, 9,
173 9, 9, 9, 9, 9, 9, 9, 9, 8, 8,
174 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
175 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
176 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
177 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
178 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
179 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
180 7, 7, 7, 7, 7, 7, 7, 7, 6, 6,
181 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
182 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
183 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
184 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
185 6, 5, 5, 5, 5, 5, 5, 5, 5, 5,
186 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
187 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
188 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
189 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
190 5, 5, 5, 4, 4, 4, 4, 4, 4, 4,
191 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
192 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
193 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
194 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
195 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
196 4, 4, 4, 4, 4, 4, 4, 4, 3, 3,
197 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
198 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
199 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
200 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
201 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
202 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
203 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
204 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
205 3, 3, 3, 3, 2, 2, 2, 2, 2, 2,
206 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
207 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
208 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
209 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
210 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
211 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
212 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
213 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
214 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
215 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
216 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
217 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
218 2, 2, 2, 2, 2, 2, 1, 1, 1, 1,
219 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
220 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
221 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
222 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
223 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
224 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
225 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
226 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
227 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
228 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
229 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
230 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
231 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
232 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
233 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
234 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
235 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
236 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
237 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
238 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
239 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
240 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
241 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
242 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
243 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
244 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
245 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
246 1, 1, 1, 1, 1, 1, 1, 0
247 };
248 static const int fe_logadd_table_size = sizeof(fe_logadd_table) / sizeof(fe_logadd_table[0]);
249 
250 static fixed32
251 fe_log_add(fixed32 x, fixed32 y)
252 {
253  fixed32 d, r;
254 
255  if (x > y) {
256  d = (x - y) >> (DEFAULT_RADIX - 8);
257  r = x;
258  }
259  else {
260  d = (y - x) >> (DEFAULT_RADIX - 8);
261  r = y;
262  }
263  if (d > fe_logadd_table_size)
264  return r;
265  else {
266  r += ((fixed32)fe_logadd_table[d] << (DEFAULT_RADIX - 8));
267 /*
268  printf("%d + %d = %d | %f + %f = %f | %f + %f = %f\n",
269  x, y, r, FIX2FLOAT(x), FIX2FLOAT(y), FIX2FLOAT(r),
270  exp(FIX2FLOAT(x)), exp(FIX2FLOAT(y)), exp(FIX2FLOAT(r)));
271 */
272  return r;
273  }
274 }
275 
276 static fixed32
277 fe_log(float32 x)
278 {
279  if (x <= 0) {
280  return MIN_FIXLOG;
281  }
282  else {
283  return FLOAT2FIX(log(x));
284  }
285 }
286 #endif
287 
288 static float32
289 fe_mel(melfb_t *mel, float32 x)
290 {
291  float32 warped = fe_warp_unwarped_to_warped(mel, x);
292 
293  return (float32) (2595.0 * log10(1.0 + warped / 700.0));
294 }
295 
296 static float32
297 fe_melinv(melfb_t *mel, float32 x)
298 {
299  float32 warped = (float32) (700.0 * (pow(10.0, x / 2595.0) - 1.0));
300  return fe_warp_warped_to_unwarped(mel, warped);
301 }
302 
303 int32
304 fe_build_melfilters(melfb_t *mel_fb)
305 {
306  float32 melmin, melmax, melbw, fftfreq;
307  int n_coeffs, i, j;
308 
309  /* Filter coefficient matrix, in flattened form. */
310  mel_fb->spec_start = ckd_malloc(mel_fb->num_filters * sizeof(*mel_fb->spec_start));
311  mel_fb->filt_start = ckd_malloc(mel_fb->num_filters * sizeof(*mel_fb->filt_start));
312  mel_fb->filt_width = ckd_malloc(mel_fb->num_filters * sizeof(*mel_fb->filt_width));
313 
314  /* First calculate the widths of each filter. */
315  /* Minimum and maximum frequencies in mel scale. */
316  melmin = fe_mel(mel_fb, mel_fb->lower_filt_freq);
317  melmax = fe_mel(mel_fb, mel_fb->upper_filt_freq);
318 
319  /* Width of filters in mel scale */
320  melbw = (melmax - melmin) / (mel_fb->num_filters + 1);
321  if (mel_fb->doublewide) {
322  melmin -= melbw;
323  melmax += melbw;
324  if ((fe_melinv(mel_fb, melmin) < 0) ||
325  (fe_melinv(mel_fb, melmax) > mel_fb->sampling_rate / 2)) {
326  E_WARN
327  ("Out of Range: low filter edge = %f (%f)\n",
328  fe_melinv(mel_fb, melmin), 0.0);
329  E_WARN
330  (" high filter edge = %f (%f)\n",
331  fe_melinv(mel_fb, melmax), mel_fb->sampling_rate / 2);
332  return FE_INVALID_PARAM_ERROR;
333  }
334  }
335 
336  /* DFT point spacing */
337  fftfreq = mel_fb->sampling_rate / (float32) mel_fb->fft_size;
338 
339  /* Count and place filter coefficients. */
340  n_coeffs = 0;
341  for (i = 0; i < mel_fb->num_filters; ++i) {
342  float32 freqs[3];
343 
344  /* Left, center, right frequencies in Hertz */
345  for (j = 0; j < 3; ++j) {
346  if (mel_fb->doublewide)
347  freqs[j] = fe_melinv(mel_fb, (i + j * 2) * melbw + melmin);
348  else
349  freqs[j] = fe_melinv(mel_fb, (i + j) * melbw + melmin);
350  /* Round them to DFT points if requested */
351  if (mel_fb->round_filters)
352  freqs[j] = ((int)(freqs[j] / fftfreq + 0.5)) * fftfreq;
353  }
354 
355  /* spec_start is the start of this filter in the power spectrum. */
356  mel_fb->spec_start[i] = -1;
357  /* There must be a better way... */
358  for (j = 0; j < mel_fb->fft_size/2+1; ++j) {
359  float32 hz = j * fftfreq;
360  if (hz < freqs[0])
361  continue;
362  else if (hz > freqs[2] || j == mel_fb->fft_size/2) {
363  /* filt_width is the width in DFT points of this filter. */
364  mel_fb->filt_width[i] = j - mel_fb->spec_start[i];
365  /* filt_start is the start of this filter in the filt_coeffs array. */
366  mel_fb->filt_start[i] = n_coeffs;
367  n_coeffs += mel_fb->filt_width[i];
368  break;
369  }
370  if (mel_fb->spec_start[i] == -1)
371  mel_fb->spec_start[i] = j;
372  }
373  }
374 
375  /* Now go back and allocate the coefficient array. */
376  mel_fb->filt_coeffs = ckd_malloc(n_coeffs * sizeof(*mel_fb->filt_coeffs));
377 
378  /* And now generate the coefficients. */
379  n_coeffs = 0;
380  for (i = 0; i < mel_fb->num_filters; ++i) {
381  float32 freqs[3];
382 
383  /* Left, center, right frequencies in Hertz */
384  for (j = 0; j < 3; ++j) {
385  if (mel_fb->doublewide)
386  freqs[j] = fe_melinv(mel_fb, (i + j * 2) * melbw + melmin);
387  else
388  freqs[j] = fe_melinv(mel_fb, (i + j) * melbw + melmin);
389  /* Round them to DFT points if requested */
390  if (mel_fb->round_filters)
391  freqs[j] = ((int)(freqs[j] / fftfreq + 0.5)) * fftfreq;
392  }
393 
394  for (j = 0; j < mel_fb->filt_width[i]; ++j) {
395  float32 hz, loslope, hislope;
396 
397  hz = (mel_fb->spec_start[i] + j) * fftfreq;
398  if (hz < freqs[0] || hz > freqs[2]) {
399  E_FATAL("WTF, %f < %f > %f\n", freqs[0], hz, freqs[2]);
400  }
401  loslope = (hz - freqs[0]) / (freqs[1] - freqs[0]);
402  hislope = (freqs[2] - hz) / (freqs[2] - freqs[1]);
403  if (mel_fb->unit_area) {
404  loslope *= 2 / (freqs[2] - freqs[0]);
405  hislope *= 2 / (freqs[2] - freqs[0]);
406  }
407  if (loslope < hislope) {
408 #ifdef FIXED_POINT
409  mel_fb->filt_coeffs[n_coeffs] = fe_log(loslope);
410 #else
411  mel_fb->filt_coeffs[n_coeffs] = loslope;
412 #endif
413  }
414  else {
415 #ifdef FIXED_POINT
416  mel_fb->filt_coeffs[n_coeffs] = fe_log(hislope);
417 #else
418  mel_fb->filt_coeffs[n_coeffs] = hislope;
419 #endif
420  }
421  ++n_coeffs;
422  }
423  }
424 
425 
426  return FE_SUCCESS;
427 }
428 
429 int32
430 fe_compute_melcosine(melfb_t * mel_fb)
431 {
432 
433  float64 freqstep;
434  int32 i, j;
435 
436  mel_fb->mel_cosine =
437  (mfcc_t **) ckd_calloc_2d(mel_fb->num_cepstra,
438  mel_fb->num_filters,
439  sizeof(mfcc_t));
440 
441  freqstep = M_PI / mel_fb->num_filters;
442  /* NOTE: The first row vector is actually unnecessary but we leave
443  * it in to avoid confusion. */
444  for (i = 0; i < mel_fb->num_cepstra; i++) {
445  for (j = 0; j < mel_fb->num_filters; j++) {
446  float64 cosine;
447 
448  cosine = cos(freqstep * i * (j + 0.5));
449  mel_fb->mel_cosine[i][j] = FLOAT2COS(cosine);
450  }
451  }
452 
453  /* Also precompute normalization constants for unitary DCT. */
454  mel_fb->sqrt_inv_n = FLOAT2COS(sqrt(1.0 / mel_fb->num_filters));
455  mel_fb->sqrt_inv_2n = FLOAT2COS(sqrt(2.0 / mel_fb->num_filters));
456 
457  /* And liftering weights */
458  if (mel_fb->lifter_val) {
459  mel_fb->lifter = calloc(mel_fb->num_cepstra, sizeof(*mel_fb->lifter));
460  for (i = 0; i < mel_fb->num_cepstra; ++i) {
461  mel_fb->lifter[i] = FLOAT2MFCC(1 + mel_fb->lifter_val / 2
462  * sin(i * M_PI / mel_fb->lifter_val));
463  }
464  }
465 
466  return (0);
467 }
468 
469 static void
470 fe_pre_emphasis(int16 const *in, frame_t * out, int32 len,
471  float32 factor, int16 prior)
472 {
473  int i;
474 
475 #if defined(FIXED16)
476  int16 fxd_alpha = (int16)(factor * 0x8000);
477  int32 tmp1, tmp2;
478 
479  tmp1 = (int32)in[0] << 15;
480  tmp2 = (int32)prior * fxd_alpha;
481  out[0] = (int16)((tmp1 - tmp2) >> 15);
482  for (i = 1; i < len; ++i) {
483  tmp1 = (int32)in[i] << 15;
484  tmp2 = (int32)in[i-1] * fxd_alpha;
485  out[i] = (int16)((tmp1 - tmp2) >> 15);
486  }
487 #elif defined(FIXED_POINT)
488  fixed32 fxd_alpha = FLOAT2FIX(factor);
489  out[0] = ((fixed32)in[0] << DEFAULT_RADIX) - (prior * fxd_alpha);
490  for (i = 1; i < len; ++i)
491  out[i] = ((fixed32)in[i] << DEFAULT_RADIX)
492  - (fixed32)in[i-1] * fxd_alpha;
493 #else
494  out[0] = (frame_t) in[0] - (frame_t) prior * factor;
495  for (i = 1; i < len; i++)
496  out[i] = (frame_t) in[i] - (frame_t) in[i-1] * factor;
497 #endif
498 }
499 
500 static void
501 fe_short_to_frame(int16 const *in, frame_t * out, int32 len)
502 {
503  int i;
504 
505 #if defined(FIXED16)
506  memcpy(out, in, len * sizeof(*out));
507 #elif defined(FIXED_POINT)
508  for (i = 0; i < len; i++)
509  out[i] = (int32) in[i] << DEFAULT_RADIX;
510 #else /* FIXED_POINT */
511  for (i = 0; i < len; i++)
512  out[i] = (frame_t) in[i];
513 #endif /* FIXED_POINT */
514 }
515 
516 void
517 fe_create_hamming(window_t * in, int32 in_len)
518 {
519  int i;
520 
521  /* Symmetric, so we only create the first half of it. */
522  for (i = 0; i < in_len / 2; i++) {
523  float64 hamm;
524  hamm = (0.54 - 0.46 * cos(2 * M_PI * i /
525  ((float64) in_len - 1.0)));
526 #ifdef FIXED16
527  in[i] = (int16)(hamm * 0x8000);
528 #else
529  in[i] = FLOAT2COS(hamm);
530 #endif
531  }
532 }
533 
534 static void
535 fe_hamming_window(frame_t * in, window_t * window, int32 in_len, int32 remove_dc)
536 {
537  int i;
538 
539  if (remove_dc) {
540 #ifdef FIXED16
541  int32 mean = 0; /* Use int32 to avoid possibility of overflow */
542 #else
543  frame_t mean = 0;
544 #endif
545 
546  for (i = 0; i < in_len; i++)
547  mean += in[i];
548  mean /= in_len;
549  for (i = 0; i < in_len; i++)
550  in[i] -= (frame_t)mean;
551  }
552 
553 #ifdef FIXED16
554  for (i = 0; i < in_len/2; i++) {
555  int32 tmp1, tmp2;
556 
557  tmp1 = (int32)in[i] * window[i];
558  tmp2 = (int32)in[in_len-1-i] * window[i];
559  in[i] = (int16)(tmp1 >> 15);
560  in[in_len-1-i] = (int16)(tmp2 >> 15);
561  }
562 #else
563  for (i = 0; i < in_len/2; i++) {
564  in[i] = COSMUL(in[i], window[i]);
565  in[in_len-1-i] = COSMUL(in[in_len-1-i], window[i]);
566  }
567 #endif
568 }
569 
570 static int
571 fe_spch_to_frame(fe_t *fe, int len)
572 {
573  /* Copy to the frame buffer. */
574  if (fe->pre_emphasis_alpha != 0.0) {
575  fe_pre_emphasis(fe->spch, fe->frame, len,
576  fe->pre_emphasis_alpha, fe->prior);
577  if (len >= fe->frame_shift)
578  fe->prior = fe->spch[fe->frame_shift - 1];
579  else
580  fe->prior = fe->spch[len - 1];
581  }
582  else
583  fe_short_to_frame(fe->spch, fe->frame, len);
584 
585  /* Zero pad up to FFT size. */
586  memset(fe->frame + len, 0,
587  (fe->fft_size - len) * sizeof(*fe->frame));
588 
589  /* Window. */
590  fe_hamming_window(fe->frame, fe->hamming_window, fe->frame_size,
591  fe->remove_dc);
592 
593  return len;
594 }
595 
596 int
597 fe_read_frame(fe_t *fe, int16 const *in, int32 len)
598 {
599  int i;
600 
601  if (len > fe->frame_size)
602  len = fe->frame_size;
603 
604  /* Read it into the raw speech buffer. */
605  memcpy(fe->spch, in, len * sizeof(*in));
606  /* Swap and dither if necessary. */
607  if (fe->swap)
608  for (i = 0; i < len; ++i)
609  SWAP_INT16(&fe->spch[i]);
610  if (fe->dither)
611  for (i = 0; i < len; ++i)
612  fe->spch[i] += (int16) ((!(s3_rand_int31() % 4)) ? 1 : 0);
613 
614  return fe_spch_to_frame(fe, len);
615 }
616 
617 int
618 fe_shift_frame(fe_t *fe, int16 const *in, int32 len)
619 {
620  int offset, i;
621 
622  if (len > fe->frame_shift)
623  len = fe->frame_shift;
624  offset = fe->frame_size - fe->frame_shift;
625 
626  /* Shift data into the raw speech buffer. */
627  memmove(fe->spch, fe->spch + fe->frame_shift,
628  offset * sizeof(*fe->spch));
629  memcpy(fe->spch + offset, in, len * sizeof(*fe->spch));
630  /* Swap and dither if necessary. */
631  if (fe->swap)
632  for (i = 0; i < len; ++i)
633  SWAP_INT16(&fe->spch[offset + i]);
634  if (fe->dither)
635  for (i = 0; i < len; ++i)
636  fe->spch[offset + i]
637  += (int16) ((!(s3_rand_int31() % 4)) ? 1 : 0);
638 
639  return fe_spch_to_frame(fe, offset + len);
640 }
641 
645 void
646 fe_create_twiddle(fe_t *fe)
647 {
648  int i;
649 
650  for (i = 0; i < fe->fft_size / 4; ++i) {
651  float64 a = 2 * M_PI * i / fe->fft_size;
652 #ifdef FIXED16
653  fe->ccc[i] = (int16)(cos(a) * 0x8000);
654  fe->sss[i] = (int16)(sin(a) * 0x8000);
655 #elif defined(FIXED_POINT)
656  fe->ccc[i] = FLOAT2COS(cos(a));
657  fe->sss[i] = FLOAT2COS(sin(a));
658 #else
659  fe->ccc[i] = cos(a);
660  fe->sss[i] = sin(a);
661 #endif
662  }
663 }
664 
665 /* Translated from the FORTRAN (obviously) from "Real-Valued Fast
666  * Fourier Transform Algorithms" by Henrik V. Sorensen et al., IEEE
667  * Transactions on Acoustics, Speech, and Signal Processing, vol. 35,
668  * no.6. The 16-bit version does a version of "block floating
669  * point" in order to avoid rounding errors.
670  */
671 #if defined(FIXED16)
672 static int
673 fe_fft_real(fe_t *fe)
674 {
675  int i, j, k, m, n, lz;
676  frame_t *x, xt, max;
677 
678  x = fe->frame;
679  m = fe->fft_order;
680  n = fe->fft_size;
681 
682  /* Bit-reverse the input. */
683  j = 0;
684  for (i = 0; i < n - 1; ++i) {
685  if (i < j) {
686  xt = x[j];
687  x[j] = x[i];
688  x[i] = xt;
689  }
690  k = n / 2;
691  while (k <= j) {
692  j -= k;
693  k /= 2;
694  }
695  j += k;
696  }
697  /* Determine how many bits of dynamic range are in the input. */
698  max = 0;
699  for (i = 0; i < n; ++i)
700  if (abs(x[i]) > max)
701  max = abs(x[i]);
702  /* The FFT has a gain of M bits, so we need to attenuate the input
703  * by M bits minus the number of leading zeroes in the input's
704  * range in order to avoid overflows. */
705  for (lz = 0; lz < m; ++lz)
706  if (max & (1 << (15-lz)))
707  break;
708 
709  /* Basic butterflies (2-point FFT, real twiddle factors):
710  * x[i] = x[i] + 1 * x[i+1]
711  * x[i+1] = x[i] + -1 * x[i+1]
712  */
713  /* The quantization error introduced by attenuating the input at
714  * any given stage of the FFT has a cascading effect, so we hold
715  * off on it until it's absolutely necessary. */
716  for (i = 0; i < n; i += 2) {
717  int atten = (lz == 0);
718  xt = x[i] >> atten;
719  x[i] = xt + (x[i + 1] >> atten);
720  x[i + 1] = xt - (x[i + 1] >> atten);
721  }
722 
723  /* The rest of the butterflies, in stages from 1..m */
724  for (k = 1; k < m; ++k) {
725  int n1, n2, n4;
726  /* Start attenuating once we hit the number of leading zeros. */
727  int atten = (k >= lz);
728 
729  n4 = k - 1;
730  n2 = k;
731  n1 = k + 1;
732  /* Stride over each (1 << (k+1)) points */
733  for (i = 0; i < n; i += (1 << n1)) {
734  /* Basic butterfly with real twiddle factors:
735  * x[i] = x[i] + 1 * x[i + (1<<k)]
736  * x[i + (1<<k)] = x[i] + -1 * x[i + (1<<k)]
737  */
738  xt = x[i] >> atten;
739  x[i] = xt + (x[i + (1 << n2)] >> atten);
740  x[i + (1 << n2)] = xt - (x[i + (1 << n2)] >> atten);
741 
742  /* The other ones with real twiddle factors:
743  * x[i + (1<<k) + (1<<(k-1))]
744  * = 0 * x[i + (1<<k-1)] + -1 * x[i + (1<<k) + (1<<k-1)]
745  * x[i + (1<<(k-1))]
746  * = 1 * x[i + (1<<k-1)] + 0 * x[i + (1<<k) + (1<<k-1)]
747  */
748  x[i + (1 << n2) + (1 << n4)] = -x[i + (1 << n2) + (1 << n4)] >> atten;
749  x[i + (1 << n4)] = x[i + (1 << n4)] >> atten;
750 
751  /* Butterflies with complex twiddle factors.
752  * There are (1<<k-1) of them.
753  */
754  for (j = 1; j < (1 << n4); ++j) {
755  frame_t cc, ss, t1, t2;
756  int i1, i2, i3, i4;
757 
758  i1 = i + j;
759  i2 = i + (1 << n2) - j;
760  i3 = i + (1 << n2) + j;
761  i4 = i + (1 << n2) + (1 << n2) - j;
762 
763  /*
764  * cc = real(W[j * n / (1<<(k+1))])
765  * ss = imag(W[j * n / (1<<(k+1))])
766  */
767  cc = fe->ccc[j << (m - n1)];
768  ss = fe->sss[j << (m - n1)];
769 
770  /* There are some symmetry properties which allow us
771  * to get away with only four multiplications here. */
772  {
773  int32 tmp1, tmp2;
774  tmp1 = (int32)x[i3] * cc + (int32)x[i4] * ss;
775  tmp2 = (int32)x[i3] * ss - (int32)x[i4] * cc;
776  t1 = (int16)(tmp1 >> 15) >> atten;
777  t2 = (int16)(tmp2 >> 15) >> atten;
778  }
779 
780  x[i4] = (x[i2] >> atten) - t2;
781  x[i3] = (-x[i2] >> atten) - t2;
782  x[i2] = (x[i1] >> atten) - t1;
783  x[i1] = (x[i1] >> atten) + t1;
784  }
785  }
786  }
787 
788  /* Return the residual scaling factor. */
789  return lz;
790 }
791 #else /* !FIXED16 */
792 static int
793 fe_fft_real(fe_t *fe)
794 {
795  int i, j, k, m, n;
796  frame_t *x, xt;
797 
798  x = fe->frame;
799  m = fe->fft_order;
800  n = fe->fft_size;
801 
802  /* Bit-reverse the input. */
803  j = 0;
804  for (i = 0; i < n - 1; ++i) {
805  if (i < j) {
806  xt = x[j];
807  x[j] = x[i];
808  x[i] = xt;
809  }
810  k = n / 2;
811  while (k <= j) {
812  j -= k;
813  k /= 2;
814  }
815  j += k;
816  }
817 
818  /* Basic butterflies (2-point FFT, real twiddle factors):
819  * x[i] = x[i] + 1 * x[i+1]
820  * x[i+1] = x[i] + -1 * x[i+1]
821  */
822  for (i = 0; i < n; i += 2) {
823  xt = x[i];
824  x[i] = (xt + x[i + 1]);
825  x[i + 1] = (xt - x[i + 1]);
826  }
827 
828  /* The rest of the butterflies, in stages from 1..m */
829  for (k = 1; k < m; ++k) {
830  int n1, n2, n4;
831 
832  n4 = k - 1;
833  n2 = k;
834  n1 = k + 1;
835  /* Stride over each (1 << (k+1)) points */
836  for (i = 0; i < n; i += (1 << n1)) {
837  /* Basic butterfly with real twiddle factors:
838  * x[i] = x[i] + 1 * x[i + (1<<k)]
839  * x[i + (1<<k)] = x[i] + -1 * x[i + (1<<k)]
840  */
841  xt = x[i];
842  x[i] = (xt + x[i + (1 << n2)]);
843  x[i + (1 << n2)] = (xt - x[i + (1 << n2)]);
844 
845  /* The other ones with real twiddle factors:
846  * x[i + (1<<k) + (1<<(k-1))]
847  * = 0 * x[i + (1<<k-1)] + -1 * x[i + (1<<k) + (1<<k-1)]
848  * x[i + (1<<(k-1))]
849  * = 1 * x[i + (1<<k-1)] + 0 * x[i + (1<<k) + (1<<k-1)]
850  */
851  x[i + (1 << n2) + (1 << n4)] = -x[i + (1 << n2) + (1 << n4)];
852  x[i + (1 << n4)] = x[i + (1 << n4)];
853 
854  /* Butterflies with complex twiddle factors.
855  * There are (1<<k-1) of them.
856  */
857  for (j = 1; j < (1 << n4); ++j) {
858  frame_t cc, ss, t1, t2;
859  int i1, i2, i3, i4;
860 
861  i1 = i + j;
862  i2 = i + (1 << n2) - j;
863  i3 = i + (1 << n2) + j;
864  i4 = i + (1 << n2) + (1 << n2) - j;
865 
866  /*
867  * cc = real(W[j * n / (1<<(k+1))])
868  * ss = imag(W[j * n / (1<<(k+1))])
869  */
870  cc = fe->ccc[j << (m - n1)];
871  ss = fe->sss[j << (m - n1)];
872 
873  /* There are some symmetry properties which allow us
874  * to get away with only four multiplications here. */
875  t1 = COSMUL(x[i3], cc) + COSMUL(x[i4], ss);
876  t2 = COSMUL(x[i3], ss) - COSMUL(x[i4], cc);
877 
878  x[i4] = (x[i2] - t2);
879  x[i3] = (-x[i2] - t2);
880  x[i2] = (x[i1] - t1);
881  x[i1] = (x[i1] + t1);
882  }
883  }
884  }
885 
886  /* This isn't used, but return it for completeness. */
887  return m;
888 }
889 #endif /* !FIXED16 */
890 
891 static void
892 fe_spec_magnitude(fe_t *fe)
893 {
894  frame_t *fft;
895  powspec_t *spec;
896  int32 j, scale, fftsize;
897 
898  /* Do FFT and get the scaling factor back (only actually used in
899  * fixed-point). Note the scaling factor is expressed in bits. */
900  scale = fe_fft_real(fe);
901 
902  /* Convenience pointers to make things less awkward below. */
903  fft = fe->frame;
904  spec = fe->spec;
905  fftsize = fe->fft_size;
906 
907  /* We need to scale things up the rest of the way to N. */
908  scale = fe->fft_order - scale;
909 
910  /* The first point (DC coefficient) has no imaginary part */
911  {
912 #ifdef FIXED16
913  spec[0] = fixlog(abs(fft[0]) << scale) * 2;
914 #elif defined(FIXED_POINT)
915  spec[0] = FIXLN(abs(fft[0]) << scale) * 2;
916 #else
917  spec[0] = fft[0] * fft[0];
918 #endif
919  }
920 
921  for (j = 1; j <= fftsize / 2; j++) {
922 #ifdef FIXED16
923  int32 rr = fixlog(abs(fft[j]) << scale) * 2;
924  int32 ii = fixlog(abs(fft[fftsize - j]) << scale) * 2;
925  spec[j] = fe_log_add(rr, ii);
926 #elif defined(FIXED_POINT)
927  int32 rr = FIXLN(abs(fft[j]) << scale) * 2;
928  int32 ii = FIXLN(abs(fft[fftsize - j]) << scale) * 2;
929  spec[j] = fe_log_add(rr, ii);
930 #else
931  spec[j] = fft[j] * fft[j] + fft[fftsize - j] * fft[fftsize - j];
932 #endif
933  }
934 }
935 
936 static void
937 fe_mel_spec(fe_t * fe)
938 {
939  int whichfilt;
940  powspec_t *spec, *mfspec;
941 
942  /* Convenience poitners. */
943  spec = fe->spec;
944  mfspec = fe->mfspec;
945 
946  for (whichfilt = 0; whichfilt < fe->mel_fb->num_filters; whichfilt++) {
947  int spec_start, filt_start, i;
948 
949  spec_start = fe->mel_fb->spec_start[whichfilt];
950  filt_start = fe->mel_fb->filt_start[whichfilt];
951 
952 #ifdef FIXED_POINT
953  mfspec[whichfilt] = spec[spec_start] + fe->mel_fb->filt_coeffs[filt_start];
954  for (i = 1; i < fe->mel_fb->filt_width[whichfilt]; i++) {
955  mfspec[whichfilt] = fe_log_add(mfspec[whichfilt],
956  spec[spec_start + i] +
957  fe->mel_fb->filt_coeffs[filt_start + i]);
958  }
959 #else /* !FIXED_POINT */
960  mfspec[whichfilt] = 0;
961  for (i = 0; i < fe->mel_fb->filt_width[whichfilt]; i++)
962  mfspec[whichfilt] +=
963  spec[spec_start + i] * fe->mel_fb->filt_coeffs[filt_start + i];
964 #endif /* !FIXED_POINT */
965  }
966 }
967 
968 static void
969 fe_mel_cep(fe_t * fe, mfcc_t *mfcep)
970 {
971  int32 i;
972  powspec_t *mfspec;
973 
974  /* Convenience pointer. */
975  mfspec = fe->mfspec;
976 
977  for (i = 0; i < fe->mel_fb->num_filters; ++i) {
978 #ifndef FIXED_POINT /* It's already in log domain for fixed point */
979  if (mfspec[i] > 0)
980  mfspec[i] = log(mfspec[i]);
981  else /* This number should be smaller than anything
982  * else, but not too small, so as to avoid
983  * infinities in the inverse transform (this is
984  * the frequency-domain equivalent of
985  * dithering) */
986  mfspec[i] = -10.0;
987 #endif /* !FIXED_POINT */
988  }
989 
990  /* If we are doing LOG_SPEC, then do nothing. */
991  if (fe->log_spec == RAW_LOG_SPEC) {
992  for (i = 0; i < fe->feature_dimension; i++) {
993  mfcep[i] = (mfcc_t) mfspec[i];
994  }
995  }
996  /* For smoothed spectrum, do DCT-II followed by (its inverse) DCT-III */
997  else if (fe->log_spec == SMOOTH_LOG_SPEC) {
998  /* FIXME: This is probably broken for fixed-point. */
999  fe_dct2(fe, mfspec, mfcep, 0);
1000  fe_dct3(fe, mfcep, mfspec);
1001  for (i = 0; i < fe->feature_dimension; i++) {
1002  mfcep[i] = (mfcc_t) mfspec[i];
1003  }
1004  }
1005  else if (fe->transform == DCT_II)
1006  fe_dct2(fe, mfspec, mfcep, FALSE);
1007  else if (fe->transform == DCT_HTK)
1008  fe_dct2(fe, mfspec, mfcep, TRUE);
1009  else
1010  fe_spec2cep(fe, mfspec, mfcep);
1011 
1012  return;
1013 }
1014 
1015 void
1016 fe_spec2cep(fe_t * fe, const powspec_t * mflogspec, mfcc_t * mfcep)
1017 {
1018  int32 i, j, beta;
1019 
1020  /* Compute C0 separately (its basis vector is 1) to avoid
1021  * costly multiplications. */
1022  mfcep[0] = mflogspec[0] / 2; /* beta = 0.5 */
1023  for (j = 1; j < fe->mel_fb->num_filters; j++)
1024  mfcep[0] += mflogspec[j]; /* beta = 1.0 */
1025  mfcep[0] /= (frame_t) fe->mel_fb->num_filters;
1026 
1027  for (i = 1; i < fe->num_cepstra; ++i) {
1028  mfcep[i] = 0;
1029  for (j = 0; j < fe->mel_fb->num_filters; j++) {
1030  if (j == 0)
1031  beta = 1; /* 0.5 */
1032  else
1033  beta = 2; /* 1.0 */
1034  mfcep[i] += COSMUL(mflogspec[j],
1035  fe->mel_fb->mel_cosine[i][j]) * beta;
1036  }
1037  /* Note that this actually normalizes by num_filters, like the
1038  * original Sphinx front-end, due to the doubled 'beta' factor
1039  * above. */
1040  mfcep[i] /= (frame_t) fe->mel_fb->num_filters * 2;
1041  }
1042 }
1043 
1044 void
1045 fe_dct2(fe_t * fe, const powspec_t * mflogspec, mfcc_t * mfcep, int htk)
1046 {
1047  int32 i, j;
1048 
1049  /* Compute C0 separately (its basis vector is 1) to avoid
1050  * costly multiplications. */
1051  mfcep[0] = mflogspec[0];
1052  for (j = 1; j < fe->mel_fb->num_filters; j++)
1053  mfcep[0] += mflogspec[j];
1054  if (htk)
1055  mfcep[0] = COSMUL(mfcep[0], fe->mel_fb->sqrt_inv_2n);
1056  else /* sqrt(1/N) = sqrt(2/N) * 1/sqrt(2) */
1057  mfcep[0] = COSMUL(mfcep[0], fe->mel_fb->sqrt_inv_n);
1058 
1059  for (i = 1; i < fe->num_cepstra; ++i) {
1060  mfcep[i] = 0;
1061  for (j = 0; j < fe->mel_fb->num_filters; j++) {
1062  mfcep[i] += COSMUL(mflogspec[j],
1063  fe->mel_fb->mel_cosine[i][j]);
1064  }
1065  mfcep[i] = COSMUL(mfcep[i], fe->mel_fb->sqrt_inv_2n);
1066  }
1067 }
1068 
1069 void
1070 fe_lifter(fe_t *fe, mfcc_t *mfcep)
1071 {
1072  int32 i;
1073 
1074  if (fe->mel_fb->lifter_val == 0)
1075  return;
1076 
1077  for (i = 0; i < fe->num_cepstra; ++i) {
1078  mfcep[i] = MFCCMUL(mfcep[i], fe->mel_fb->lifter[i]);
1079  }
1080 }
1081 
1082 void
1083 fe_dct3(fe_t * fe, const mfcc_t * mfcep, powspec_t * mflogspec)
1084 {
1085  int32 i, j;
1086 
1087  for (i = 0; i < fe->mel_fb->num_filters; ++i) {
1088  mflogspec[i] = COSMUL(mfcep[0], SQRT_HALF);
1089  for (j = 1; j < fe->num_cepstra; j++) {
1090  mflogspec[i] += COSMUL(mfcep[j],
1091  fe->mel_fb->mel_cosine[j][i]);
1092  }
1093  mflogspec[i] = COSMUL(mflogspec[i], fe->mel_fb->sqrt_inv_2n);
1094  }
1095 }
1096 
1097 int32
1098 fe_write_frame(fe_t * fe, mfcc_t * fea)
1099 {
1100  fe_spec_magnitude(fe);
1101  fe_mel_spec(fe);
1102  fe_mel_cep(fe, fea);
1103  fe_lifter(fe, fea);
1104 
1105  return 1;
1106 }
1107 
1108 void *
1109 fe_create_2d(int32 d1, int32 d2, int32 elem_size)
1110 {
1111  return (void *)ckd_calloc_2d(d1, d2, elem_size);
1112 }
1113 
1114 void
1115 fe_free_2d(void *arr)
1116 {
1117  ckd_free_2d((void **)arr);
1118 }