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