SphinxBase  5prealpha
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
57 #endif
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 
82 /* Internal log-addition table for natural log with radix point at 8
83  * bits. Each entry is 256 * log(1 + e^{-n/256}). This is used in the
84  * log-add computation:
85  *
86  * e^z = e^x + e^y
87  * e^z = e^x(1 + e^{y-x}) = e^y(1 + e^{x-y})
88  * z = x + log(1 + e^{y-x}) = y + log(1 + e^{x-y})
89  *
90  * So when y > x, z = y + logadd_table[-(x-y)]
91  * when x > y, z = x + logadd_table[-(y-x)]
92  */
93 static const unsigned char fe_logadd_table[] = {
94  177, 177, 176, 176, 175, 175, 174, 174, 173, 173,
95  172, 172, 172, 171, 171, 170, 170, 169, 169, 168,
96  168, 167, 167, 166, 166, 165, 165, 164, 164, 163,
97  163, 162, 162, 161, 161, 161, 160, 160, 159, 159,
98  158, 158, 157, 157, 156, 156, 155, 155, 155, 154,
99  154, 153, 153, 152, 152, 151, 151, 151, 150, 150,
100  149, 149, 148, 148, 147, 147, 147, 146, 146, 145,
101  145, 144, 144, 144, 143, 143, 142, 142, 141, 141,
102  141, 140, 140, 139, 139, 138, 138, 138, 137, 137,
103  136, 136, 136, 135, 135, 134, 134, 134, 133, 133,
104  132, 132, 131, 131, 131, 130, 130, 129, 129, 129,
105  128, 128, 128, 127, 127, 126, 126, 126, 125, 125,
106  124, 124, 124, 123, 123, 123, 122, 122, 121, 121,
107  121, 120, 120, 119, 119, 119, 118, 118, 118, 117,
108  117, 117, 116, 116, 115, 115, 115, 114, 114, 114,
109  113, 113, 113, 112, 112, 112, 111, 111, 110, 110,
110  110, 109, 109, 109, 108, 108, 108, 107, 107, 107,
111  106, 106, 106, 105, 105, 105, 104, 104, 104, 103,
112  103, 103, 102, 102, 102, 101, 101, 101, 100, 100,
113  100, 99, 99, 99, 98, 98, 98, 97, 97, 97,
114  96, 96, 96, 96, 95, 95, 95, 94, 94, 94,
115  93, 93, 93, 92, 92, 92, 92, 91, 91, 91,
116  90, 90, 90, 89, 89, 89, 89, 88, 88, 88,
117  87, 87, 87, 87, 86, 86, 86, 85, 85, 85,
118  85, 84, 84, 84, 83, 83, 83, 83, 82, 82,
119  82, 82, 81, 81, 81, 80, 80, 80, 80, 79,
120  79, 79, 79, 78, 78, 78, 78, 77, 77, 77,
121  77, 76, 76, 76, 75, 75, 75, 75, 74, 74,
122  74, 74, 73, 73, 73, 73, 72, 72, 72, 72,
123  71, 71, 71, 71, 71, 70, 70, 70, 70, 69,
124  69, 69, 69, 68, 68, 68, 68, 67, 67, 67,
125  67, 67, 66, 66, 66, 66, 65, 65, 65, 65,
126  64, 64, 64, 64, 64, 63, 63, 63, 63, 63,
127  62, 62, 62, 62, 61, 61, 61, 61, 61, 60,
128  60, 60, 60, 60, 59, 59, 59, 59, 59, 58,
129  58, 58, 58, 58, 57, 57, 57, 57, 57, 56,
130  56, 56, 56, 56, 55, 55, 55, 55, 55, 54,
131  54, 54, 54, 54, 53, 53, 53, 53, 53, 52,
132  52, 52, 52, 52, 52, 51, 51, 51, 51, 51,
133  50, 50, 50, 50, 50, 50, 49, 49, 49, 49,
134  49, 49, 48, 48, 48, 48, 48, 48, 47, 47,
135  47, 47, 47, 47, 46, 46, 46, 46, 46, 46,
136  45, 45, 45, 45, 45, 45, 44, 44, 44, 44,
137  44, 44, 43, 43, 43, 43, 43, 43, 43, 42,
138  42, 42, 42, 42, 42, 41, 41, 41, 41, 41,
139  41, 41, 40, 40, 40, 40, 40, 40, 40, 39,
140  39, 39, 39, 39, 39, 39, 38, 38, 38, 38,
141  38, 38, 38, 37, 37, 37, 37, 37, 37, 37,
142  37, 36, 36, 36, 36, 36, 36, 36, 35, 35,
143  35, 35, 35, 35, 35, 35, 34, 34, 34, 34,
144  34, 34, 34, 34, 33, 33, 33, 33, 33, 33,
145  33, 33, 32, 32, 32, 32, 32, 32, 32, 32,
146  32, 31, 31, 31, 31, 31, 31, 31, 31, 31,
147  30, 30, 30, 30, 30, 30, 30, 30, 30, 29,
148  29, 29, 29, 29, 29, 29, 29, 29, 28, 28,
149  28, 28, 28, 28, 28, 28, 28, 28, 27, 27,
150  27, 27, 27, 27, 27, 27, 27, 27, 26, 26,
151  26, 26, 26, 26, 26, 26, 26, 26, 25, 25,
152  25, 25, 25, 25, 25, 25, 25, 25, 25, 24,
153  24, 24, 24, 24, 24, 24, 24, 24, 24, 24,
154  23, 23, 23, 23, 23, 23, 23, 23, 23, 23,
155  23, 23, 22, 22, 22, 22, 22, 22, 22, 22,
156  22, 22, 22, 22, 21, 21, 21, 21, 21, 21,
157  21, 21, 21, 21, 21, 21, 21, 20, 20, 20,
158  20, 20, 20, 20, 20, 20, 20, 20, 20, 20,
159  19, 19, 19, 19, 19, 19, 19, 19, 19, 19,
160  19, 19, 19, 19, 18, 18, 18, 18, 18, 18,
161  18, 18, 18, 18, 18, 18, 18, 18, 18, 17,
162  17, 17, 17, 17, 17, 17, 17, 17, 17, 17,
163  17, 17, 17, 17, 16, 16, 16, 16, 16, 16,
164  16, 16, 16, 16, 16, 16, 16, 16, 16, 16,
165  16, 15, 15, 15, 15, 15, 15, 15, 15, 15,
166  15, 15, 15, 15, 15, 15, 15, 15, 14, 14,
167  14, 14, 14, 14, 14, 14, 14, 14, 14, 14,
168  14, 14, 14, 14, 14, 14, 14, 13, 13, 13,
169  13, 13, 13, 13, 13, 13, 13, 13, 13, 13,
170  13, 13, 13, 13, 13, 13, 13, 12, 12, 12,
171  12, 12, 12, 12, 12, 12, 12, 12, 12, 12,
172  12, 12, 12, 12, 12, 12, 12, 12, 12, 11,
173  11, 11, 11, 11, 11, 11, 11, 11, 11, 11,
174  11, 11, 11, 11, 11, 11, 11, 11, 11, 11,
175  11, 11, 11, 10, 10, 10, 10, 10, 10, 10,
176  10, 10, 10, 10, 10, 10, 10, 10, 10, 10,
177  10, 10, 10, 10, 10, 10, 10, 10, 10, 9,
178  9, 9, 9, 9, 9, 9, 9, 9, 9, 9,
179  9, 9, 9, 9, 9, 9, 9, 9, 9, 9,
180  9, 9, 9, 9, 9, 9, 9, 9, 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  8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
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, 7, 7,
187  7, 7, 7, 7, 7, 7, 7, 7, 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, 6, 6, 6, 6, 6, 6, 6, 6, 6,
192  6, 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, 5, 5, 5, 5, 5, 5, 5,
197  5, 5, 5, 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, 4, 4,
203  4, 4, 4, 4, 4, 4, 4, 4, 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, 3, 3, 3, 3, 3, 3,
212  3, 3, 3, 3, 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, 2, 2, 2, 2,
225  2, 2, 2, 2, 2, 2, 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, 1, 1, 1,
253  1, 1, 1, 1, 1, 1, 1, 0
254 };
255 
256 static const int fe_logadd_table_size =
257  sizeof(fe_logadd_table) / sizeof(fe_logadd_table[0]);
258 
259 fixed32
260 fe_log_add(fixed32 x, fixed32 y)
261 {
262  fixed32 d, r;
263 
264  if (x > y) {
265  d = (x - y) >> (DEFAULT_RADIX - 8);
266  r = x;
267  }
268  else {
269  d = (y - x) >> (DEFAULT_RADIX - 8);
270  r = y;
271  }
272 
273  if (r <= MIN_FIXLOG)
274  return MIN_FIXLOG;
275  else if (d > fe_logadd_table_size - 1)
276  return r;
277  else {
278  r += ((fixed32) fe_logadd_table[d] << (DEFAULT_RADIX - 8));
279 /* printf("%d - %d = %d | %f - %f = %f | %f - %f = %f\n",
280  x, y, r, FIX2FLOAT(x), FIX2FLOAT(y), FIX2FLOAT(r),
281  exp(FIX2FLOAT(x)), exp(FIX2FLOAT(y)), exp(FIX2FLOAT(r)));
282 */
283  return r;
284  }
285 }
286 
287 /*
288  * log_sub for spectral subtraction, similar to logadd but we had
289  * to smooth function around zero with fixlog in order to improve
290  * table interpolation properties
291  *
292  * The table is created with the file included into distribution
293  *
294  * e^z = e^x - e^y
295  * e^z = e^x (1 - e^(-(x - y)))
296  * z = x + log(1 - e^(-(x - y)))
297  * z = x + fixlog(a) + (log(1 - e^(- a)) - log(a))
298  *
299  * Input radix is 8 output radix is 10
300  */
301 static const uint16 fe_logsub_table[] = {
302 1, 3, 5, 7, 9, 11, 13, 15, 17, 19,
303 21, 23, 25, 27, 29, 31, 33, 35, 37, 39,
304 41, 43, 45, 47, 49, 51, 53, 55, 56, 58,
305 60, 62, 64, 66, 68, 70, 72, 74, 76, 78,
306 80, 82, 84, 86, 88, 90, 92, 94, 95, 97,
307 99, 101, 103, 105, 107, 109, 111, 113, 115, 117,
308 119, 121, 122, 124, 126, 128, 130, 132, 134, 136,
309 138, 140, 142, 143, 145, 147, 149, 151, 153, 155,
310 157, 159, 161, 162, 164, 166, 168, 170, 172, 174,
311 176, 178, 179, 181, 183, 185, 187, 189, 191, 193,
312 194, 196, 198, 200, 202, 204, 206, 207, 209, 211,
313 213, 215, 217, 219, 220, 222, 224, 226, 228, 230,
314 232, 233, 235, 237, 239, 241, 243, 244, 246, 248,
315 250, 252, 254, 255, 257, 259, 261, 263, 265, 266,
316 268, 270, 272, 274, 275, 277, 279, 281, 283, 284,
317 286, 288, 290, 292, 294, 295, 297, 299, 301, 302,
318 304, 306, 308, 310, 311, 313, 315, 317, 319, 320,
319 322, 324, 326, 327, 329, 331, 333, 335, 336, 338,
320 340, 342, 343, 345, 347, 349, 350, 352, 354, 356,
321 357, 359, 361, 363, 364, 366, 368, 370, 371, 373,
322 375, 377, 378, 380, 382, 384, 385, 387, 389, 391,
323 392, 394, 396, 397, 399, 401, 403, 404, 406, 408,
324 410, 411, 413, 415, 416, 418, 420, 422, 423, 425,
325 427, 428, 430, 432, 433, 435, 437, 439, 440, 442,
326 444, 445, 447, 449, 450, 452, 454, 455, 457, 459,
327 460, 462, 464, 465, 467, 469, 471, 472, 474, 476,
328 477, 479, 481, 482, 484, 486, 487, 489, 490, 492,
329 494, 495, 497, 499, 500, 502, 504, 505, 507, 509,
330 510, 512, 514, 515, 517, 518, 520, 522, 523, 525,
331 527, 528, 530, 532, 533, 535, 536, 538, 540, 541,
332 543, 544, 546, 548, 549, 551, 553, 554, 556, 557,
333 559, 561, 562, 564, 565, 567, 569, 570, 572, 573,
334 575, 577, 578, 580, 581, 583, 585, 586, 588, 589,
335 591, 592, 594, 596, 597, 599, 600, 602, 603, 605,
336 607, 608, 610, 611, 613, 614, 616, 618, 619, 621,
337 622, 624, 625, 627, 628, 630, 632, 633, 635, 636,
338 638, 639, 641, 642, 644, 645, 647, 649, 650, 652,
339 653, 655, 656, 658, 659, 661, 662, 664, 665, 667,
340 668, 670, 671, 673, 674, 676, 678, 679, 681, 682,
341 684, 685, 687, 688, 690, 691, 693, 694, 696, 697,
342 699, 700, 702, 703, 705, 706, 708, 709, 711, 712,
343 714, 715, 717, 718, 719, 721, 722, 724, 725, 727,
344 728, 730, 731, 733, 734, 736, 737, 739, 740, 742,
345 743, 745, 746, 747, 749, 750, 752, 753, 755, 756,
346 758, 759, 761, 762, 763, 765, 766, 768, 769, 771,
347 772, 774, 775, 776, 778, 779, 781, 782, 784, 785,
348 786, 788, 789, 791, 792, 794, 795, 796, 798, 799,
349 801, 802, 804, 805, 806, 808, 809, 811, 812, 813,
350 815, 816, 818, 819, 820, 822, 823, 825, 826, 827,
351 829, 830, 832, 833, 834, 836, 837, 839, 840, 841,
352 843, 844, 846, 847, 848, 850, 851, 852, 854, 855,
353 857, 858, 859, 861, 862, 863, 865, 866, 868, 869,
354 870, 872, 873, 874, 876, 877, 878, 880, 881, 883,
355 884, 885, 887, 888, 889, 891, 892, 893, 895, 896,
356 897, 899, 900, 901, 903, 904, 905, 907, 908, 909,
357 911, 912, 913, 915, 916, 917, 919, 920, 921, 923,
358 924, 925, 927, 928, 929, 931, 932, 933, 935, 936,
359 937, 939, 940, 941, 942, 944, 945, 946, 948, 949,
360 950, 952, 953, 954, 956, 957, 958, 959, 961, 962,
361 963, 965, 966, 967, 968, 970, 971, 972, 974, 975,
362 976, 977, 979, 980, 981, 983, 984, 985, 986, 988,
363 989, 990, 991, 993, 994, 995, 997, 998, 999, 1000,
364 1002, 1003, 1004, 1005, 1007, 1008, 1009, 1010, 1012, 1013,
365 1014, 1015, 1017, 1018, 1019, 1020, 1022, 1023, 1024, 1025,
366 1027, 1028, 1029, 1030, 1032, 1033, 1034, 1035, 1037, 1038,
367 1039, 1040, 1041, 1043, 1044, 1045, 1046, 1048, 1049, 1050,
368 1051, 1052, 1054, 1055, 1056, 1057, 1059, 1060, 1061, 1062,
369 1063, 1065, 1066, 1067, 1068, 1069, 1071, 1072, 1073, 1074,
370 1076, 1077, 1078, 1079, 1080, 1082, 1083, 1084, 1085, 1086,
371 1087, 1089, 1090, 1091, 1092, 1093, 1095, 1096, 1097, 1098,
372 1099, 1101, 1102, 1103, 1104, 1105, 1106, 1108, 1109, 1110,
373 1111, 1112, 1114, 1115, 1116, 1117, 1118, 1119, 1121, 1122,
374 1123, 1124, 1125, 1126, 1128, 1129, 1130, 1131, 1132, 1133,
375 1135, 1136, 1137, 1138, 1139, 1140, 1141, 1143, 1144, 1145,
376 1146, 1147, 1148, 1149, 1151, 1152, 1153, 1154, 1155, 1156,
377 1157, 1159, 1160, 1161, 1162, 1163, 1164, 1165, 1167, 1168,
378 1169, 1170, 1171, 1172, 1173, 1174, 1176, 1177, 1178, 1179,
379 1180, 1181, 1182, 1183, 1185, 1186, 1187, 1188, 1189, 1190,
380 1191, 1192, 1193, 1195, 1196, 1197, 1198, 1199, 1200, 1201,
381 1202, 1203, 1205, 1206, 1207, 1208, 1209, 1210, 1211, 1212,
382 1213, 1214, 1216, 1217, 1218, 1219, 1220, 1221, 1222, 1223,
383 1224, 1225, 1226, 1228, 1229, 1230, 1231, 1232, 1233, 1234,
384 1235, 1236, 1237, 1238, 1239, 1240, 1242, 1243, 1244, 1245,
385 1246, 1247, 1248, 1249, 1250, 1251, 1252, 1253, 1254, 1255,
386 1256, 1258, 1259, 1260, 1261, 1262, 1263, 1264, 1265, 1266,
387 1267, 1268, 1269, 1270, 1271, 1272, 1273, 1274, 1275, 1277,
388 1278, 1279, 1280, 1281, 1282, 1283, 1284, 1285, 1286, 1287,
389 1288, 1289, 1290, 1291, 1292, 1293, 1294, 1295, 1296, 1297,
390 1298, 1299, 1300, 1301, 1302, 1303, 1305, 1306, 1307, 1308,
391 1309, 1310, 1311, 1312, 1313, 1314, 1315, 1316, 1317, 1318,
392 1319, 1320, 1321, 1322, 1323, 1324, 1325, 1326, 1327, 1328,
393 1329, 1330, 1331, 1332, 1333, 1334, 1335, 1336, 1337, 1338,
394 1339, 1340, 1341, 1342, 1343, 1344, 1345, 1346, 1347, 1348,
395 1349, 1350, 1351, 1352, 1353, 1354, 1355, 1356, 1357, 1358,
396 1359, 1360, 1361, 1362, 1363, 1364, 1365, 1366, 1367, 1368,
397 1369, 1370, 1371, 1372, 1372, 1373, 1374, 1375, 1376, 1377,
398 1378, 1379, 1380, 1381, 1382, 1383, 1384, 1385, 1386, 1387,
399 1388, 1389, 1390, 1391, 1392, 1393, 1394, 1395, 1396, 1397,
400 1398, 1399, 1399, 1400, 1401, 1402, 1403, 1404, 1405, 1406,
401 1407, 1408, 1409, 1410, 1411, 1412, 1413, 1414, 1415, 1416,
402 1417, 1418, 1418, 1419, 1420, 1421, 1422, 1423, 1424, 1425,
403 1426, 1427, 1428, 1429, 1430, 1431, 1432, 1432, 1433, 1434,
404 1435, 1436, 1437, 1438, 1439, 1440, 1441, 1442, 1443, 1444,
405 1444, 1445, 1446, 1447, 1448, 1449, 1450, 1451, 1452, 1453,
406 1454, 1455, 1455, 1456, 1457, 1458, 1459, 1460, 1461, 1462,
407 1463, 1464, 1465, 1466, 1466, 1467, 1468, 1469, 1470, 1471,
408 1472, 1473, 1474, 1475, 1475, 1476, 1477, 1478, 1479, 1480,
409 1481, 1482, 1483, 1483, 1484, 1485, 1486, 1487, 1488, 1489,
410 1490, 1491, 1491, 1492, 1493, 1494, 1495, 1496, 1497, 1498,
411 1499, 1499, 1500, 1501, 1502, 1503, 1504, 1505, 1506, 1506,
412 1507, 1508, 1509, 1510, 1511, 1512, 1513, 1513, 1514, 1515,
413 1516, 1517, 1518, 1519, 1520, 1520, 1521, 1522, 1523, 1524,
414 1525, 1526, 1526, 1527, 1528, 1529, 1530, 1531, 1532, 1532,
415 1533, 1534, 1535, 1536, 1537, 1538, 1538, 1539, 1540, 1541,
416 1542, 1543, 1544, 1544, 1545, 1546, 1547, 1548, 1549, 1550,
417 1550, 1551, 1552, 1553, 1554, 1555, 1555, 1556, 1557, 1558,
418 1559, 1560, 1560, 1561, 1562, 1563, 1564, 1565, 1565, 1566,
419 1567, 1568, 1569, 1570, 1570, 1571, 1572, 1573, 1574, 1575,
420 1575, 1576, 1577, 1578, 1579, 1580, 1580, 1581, 1582, 1583,
421 1584, 1584, 1585, 1586, 1587, 1588, 1589, 1589, 1590, 1591,
422 1592, 1593, 1593, 1594, 1595, 1596, 1597, 1598, 1598, 1599,
423 1600, 1601, 1602, 1602, 1603, 1604, 1605, 1606, 1606, 1607,
424 1608, 1609, 1610, 1610, 1611, 1612, 1613, 1614, 1614, 1615,
425 1616, 1617, 1618, 1618, 1619, 1620, 1621, 1622, 1622, 1623,
426 1624, 1625, 1626, 1626, 1627, 1628, 1629, 1630, 1630, 1631,
427 1632, 1633, 1634, 1634, 1635, 1636, 1637, 1637, 1638, 1639,
428 1640, 1641, 1641, 1642, 1643, 1644, 1645, 1645, 1646, 1647,
429 1648, 1648, 1649, 1650, 1651, 1652, 1652, 1653, 1654, 1655,
430 1655, 1656, 1657, 1658, 1658, 1659, 1660, 1661, 1662, 1662,
431 1663, 1664, 1665, 1665, 1666, 1667, 1668, 1668, 1669, 1670,
432 1671, 1671, 1672, 1673, 1674, 1675, 1675, 1676, 1677, 1678,
433 1678, 1679, 1680, 1681, 1681, 1682, 1683, 1684, 1684, 1685,
434 1686, 1687, 1687, 1688, 1689, 1690, 1690, 1691, 1692, 1693,
435 1693, 1694, 1695, 1696, 1696, 1697, 1698, 1699, 1699, 1700,
436 1701, 1702, 1702, 1703, 1704, 1705, 1705, 1706, 1707, 1707,
437 1708, 1709, 1710, 1710, 1711, 1712, 1713, 1713, 1714, 1715,
438 1716, 1716, 1717, 1718, 1718, 1719, 1720, 1721, 1721, 1722,
439 1723, 1724, 1724, 1725, 1726, 1727, 1727, 1728, 1729, 1729,
440 1730, 1731, 1732, 1732, 1733, 1734, 1734, 1735, 1736, 1737,
441 1737, 1738, 1739, 1740, 1740, 1741, 1742, 1742, 1743, 1744,
442 1745, 1745, 1746, 1747, 1747, 1748, 1749, 1749, 1750, 1751,
443 1752, 1752, 1753, 1754, 1754, 1755, 1756, 1757, 1757, 1758,
444 1759, 1759, 1760, 1761, 1762, 1762, 1763, 1764, 1764, 1765,
445 1766, 1766, 1767, 1768, 1769, 1769, 1770, 1771, 1771, 1772,
446 1773, 1773, 1774, 1775, 1776, 1776, 1777, 1778, 1778, 1779,
447 1780, 1780, 1781, 1782, 1782, 1783, 1784, 1784, 1785, 1786,
448 1787, 1787, 1788, 1789, 1789, 1790, 1791, 1791, 1792, 1793,
449 1793, 1794, 1795, 1795, 1796, 1797, 1798, 1798, 1799, 1800,
450 1800, 1801, 1802, 1802, 1803, 1804, 1804, 1805, 1806, 1806,
451 1807, 1808, 1808, 1809, 1810, 1810, 1811, 1812, 1812, 1813,
452 1814, 1814, 1815, 1816, 1816, 1817, 1818, 1818, 1819, 1820,
453 1820, 1821, 1822, 1822, 1823, 1824, 1824, 1825, 1826, 1826,
454 1827, 1828, 1828, 1829, 1830, 1830, 1831, 1832, 1832, 1833,
455 1834, 1834, 1835, 1836, 1836, 1837, 1838, 1838, 1839, 1840,
456 1840, 1841, 1842, 1842, 1843, 1844, 1844, 1845, 1845, 1846,
457 1847, 1847, 1848, 1849, 1849, 1850, 1851, 1851, 1852, 1853,
458 1853, 1854, 1855, 1855, 1856, 1857, 1857, 1858, 1858, 1859,
459 1860, 1860, 1861, 1862, 1862, 1863, 1864, 1864, 1865, 1866,
460 1866, 1867, 1867, 1868, 1869, 1869, 1870, 1871, 1871, 1872,
461 1873, 1873, 1874, 1874, 1875, 1876, 1876, 1877, 1878, 1878,
462 1879, 1879, 1880, 1881, 1881, 1882, 1883, 1883, 1884, 1885,
463 1885, 1886, 1886, 1887, 1888, 1888, 1889, 1890, 1890, 1891,
464 1891, 1892, 1893, 1893, 1894, 1895, 1895, 1896, 1896, 1897,
465 1898, 1898, 1899, 1900, 1900, 1901, 1901, 1902, 1903, 1903,
466 1904, 1904, 1905, 1906, 1906, 1907, 1908, 1908, 1909, 1909,
467 1910, 1911, 1911, 1912, 1912, 1913, 1914, 1914, 1915, 1916,
468 1916, 1917, 1917, 1918, 1919, 1919, 1920, 1920, 1921, 1922,
469 1922, 1923, 1923, 1924, 1925, 1925, 1926, 1926, 1927, 1928,
470 1928, 1929, 1929, 1930, 1931, 1931, 1932, 1932, 1933, 1934,
471 1934, 1935, 1935, 1936, 1937, 1937, 1938, 1938, 1939, 1940,
472 1940, 1941, 1941, 1942, 1943, 1943, 1944, 1944, 1945, 1946,
473 1946, 1947, 1947, 1948, 1949, 1949, 1950, 1950, 1951, 1952,
474 1952, 1953, 1953, 1954, 1955, 1955, 1956, 1956, 1957, 1957,
475 1958, 1959, 1959, 1960, 1960, 1961, 1962, 1962, 1963, 1963,
476 1964, 1964, 1965, 1966, 1966, 1967, 1967, 1968, 1969, 1969,
477 1970, 1970, 1971, 1971, 1972, 1973, 1973, 1974, 1974, 1975,
478 1976, 1976, 1977, 1977, 1978, 1978, 1979, 1980, 1980, 1981,
479 1981, 1982, 1982, 1983, 1984, 1984, 1985, 1985, 1986, 1986,
480 1987, 1988, 1988, 1989, 1989, 1990, 1990, 1991, 1992, 1992,
481 1993, 1993, 1994, 1994, 1995, 1996, 1996, 1997, 1997, 1998,
482 1998, 1999, 1999, 2000, 2001, 2001, 2002, 2002, 2003, 2003,
483 2004, 2005, 2005, 2006, 2006, 2007, 2007, 2008, 2008, 2009,
484 2010, 2010, 2011, 2011, 2012, 2012, 2013, 2014, 2014, 2015,
485 2015, 2016, 2016, 2017, 2017, 2018, 2019, 2019, 2020, 2020,
486 2021, 2021, 2022, 2022, 2023, 2023, 2024, 2025, 2025, 2026,
487 2026, 2027, 2027, 2028, 2028, 2029, 2030, 2030, 2031, 2031,
488 2032, 2032, 2033, 2033, 2034, 2034, 2035, 2036, 2036, 2037,
489 2037, 2038, 2038, 2039, 2039, 2040, 2040, 2041, 2042, 2042,
490 2043, 2043, 2044, 2044, 2045, 2045, 2046, 2046, 2047, 2048,
491 2048, 2049, 2049, 2050, 2050, 2051, 2051, 2052, 2052, 2053,
492 2053, 2054, 2054, 2055, 2056, 2056, 2057, 2057, 2058, 2058,
493 2059, 2059, 2060, 2060, 2061, 2061, 2062, 2062, 2063, 2064,
494 2064, 2065, 2065, 2066, 2066, 2067, 2067, 2068, 2068, 2069,
495 2069, 2070, 2070, 2071, 2072, 2072, 2073, 2073, 2074, 2074,
496 2075, 2075, 2076, 2076, 2077, 2077, 2078, 2078, 2079, 2079,
497 2080, 2080, 2081
498 };
499 
500 static const int fe_logsub_table_size =
501  sizeof(fe_logsub_table) / sizeof(fe_logsub_table[0]);
502 
503 fixed32
504 fe_log_sub(fixed32 x, fixed32 y)
505 {
506  fixed32 d, r;
507 
508  if (x < MIN_FIXLOG || y >= x)
509  return MIN_FIXLOG;
510 
511  d = (x - y) >> (DEFAULT_RADIX - 8);
512 
513  if (d > fe_logsub_table_size - 1)
514  return x;
515 
516  r = fe_logsub_table[d] << (DEFAULT_RADIX - 10);
517 /*
518  printf("diff=%d\n",
519  x + FIXLN(x-y) - r -
520  (x + FLOAT2FIX(logf(-expm1f(FIX2FLOAT(y - x))))));
521 */
522  return x + FIXLN(x-y) - r;
523 }
524 
525 static fixed32
526 fe_log(float32 x)
527 {
528  if (x <= 0) {
529  return MIN_FIXLOG;
530  }
531  else {
532  return FLOAT2FIX(log(x));
533  }
534 }
535 #endif
536 
537 static float32
538 fe_mel(melfb_t * mel, float32 x)
539 {
540  float32 warped = fe_warp_unwarped_to_warped(mel, x);
541 
542  return (float32) (2595.0 * log10(1.0 + warped / 700.0));
543 }
544 
545 static float32
546 fe_melinv(melfb_t * mel, float32 x)
547 {
548  float32 warped = (float32) (700.0 * (pow(10.0, x / 2595.0) - 1.0));
549  return fe_warp_warped_to_unwarped(mel, warped);
550 }
551 
552 int32
553 fe_build_melfilters(melfb_t * mel_fb)
554 {
555  float32 melmin, melmax, melbw, fftfreq;
556  int n_coeffs, i, j;
557 
558 
559  /* Filter coefficient matrix, in flattened form. */
560  mel_fb->spec_start =
561  ckd_calloc(mel_fb->num_filters, sizeof(*mel_fb->spec_start));
562  mel_fb->filt_start =
563  ckd_calloc(mel_fb->num_filters, sizeof(*mel_fb->filt_start));
564  mel_fb->filt_width =
565  ckd_calloc(mel_fb->num_filters, sizeof(*mel_fb->filt_width));
566 
567  /* First calculate the widths of each filter. */
568  /* Minimum and maximum frequencies in mel scale. */
569  melmin = fe_mel(mel_fb, mel_fb->lower_filt_freq);
570  melmax = fe_mel(mel_fb, mel_fb->upper_filt_freq);
571 
572  /* Width of filters in mel scale */
573  melbw = (melmax - melmin) / (mel_fb->num_filters + 1);
574  if (mel_fb->doublewide) {
575  melmin -= melbw;
576  melmax += melbw;
577  if ((fe_melinv(mel_fb, melmin) < 0) ||
578  (fe_melinv(mel_fb, melmax) > mel_fb->sampling_rate / 2)) {
579  E_WARN
580  ("Out of Range: low filter edge = %f (%f)\n",
581  fe_melinv(mel_fb, melmin), 0.0);
582  E_WARN
583  (" high filter edge = %f (%f)\n",
584  fe_melinv(mel_fb, melmax), mel_fb->sampling_rate / 2);
585  return FE_INVALID_PARAM_ERROR;
586  }
587  }
588 
589  /* DFT point spacing */
590  fftfreq = mel_fb->sampling_rate / (float32) mel_fb->fft_size;
591 
592  /* Count and place filter coefficients. */
593  n_coeffs = 0;
594  for (i = 0; i < mel_fb->num_filters; ++i) {
595  float32 freqs[3];
596 
597  /* Left, center, right frequencies in Hertz */
598  for (j = 0; j < 3; ++j) {
599  if (mel_fb->doublewide)
600  freqs[j] = fe_melinv(mel_fb, (i + j * 2) * melbw + melmin);
601  else
602  freqs[j] = fe_melinv(mel_fb, (i + j) * melbw + melmin);
603  /* Round them to DFT points if requested */
604  if (mel_fb->round_filters)
605  freqs[j] = ((int) (freqs[j] / fftfreq + 0.5)) * fftfreq;
606  }
607 
608  /* spec_start is the start of this filter in the power spectrum. */
609  mel_fb->spec_start[i] = -1;
610  /* There must be a better way... */
611  for (j = 0; j < mel_fb->fft_size / 2 + 1; ++j) {
612  float32 hz = j * fftfreq;
613  if (hz < freqs[0])
614  continue;
615  else if (hz > freqs[2] || j == mel_fb->fft_size / 2) {
616  /* filt_width is the width in DFT points of this filter. */
617  mel_fb->filt_width[i] = j - mel_fb->spec_start[i];
618  /* filt_start is the start of this filter in the filt_coeffs array. */
619  mel_fb->filt_start[i] = n_coeffs;
620  n_coeffs += mel_fb->filt_width[i];
621  break;
622  }
623  if (mel_fb->spec_start[i] == -1)
624  mel_fb->spec_start[i] = j;
625  }
626  }
627 
628  /* Now go back and allocate the coefficient array. */
629  mel_fb->filt_coeffs =
630  ckd_malloc(n_coeffs * sizeof(*mel_fb->filt_coeffs));
631 
632  /* And now generate the coefficients. */
633  n_coeffs = 0;
634  for (i = 0; i < mel_fb->num_filters; ++i) {
635  float32 freqs[3];
636 
637  /* Left, center, right frequencies in Hertz */
638  for (j = 0; j < 3; ++j) {
639  if (mel_fb->doublewide)
640  freqs[j] = fe_melinv(mel_fb, (i + j * 2) * melbw + melmin);
641  else
642  freqs[j] = fe_melinv(mel_fb, (i + j) * melbw + melmin);
643  /* Round them to DFT points if requested */
644  if (mel_fb->round_filters)
645  freqs[j] = ((int) (freqs[j] / fftfreq + 0.5)) * fftfreq;
646  }
647 
648  for (j = 0; j < mel_fb->filt_width[i]; ++j) {
649  float32 hz, loslope, hislope;
650 
651  hz = (mel_fb->spec_start[i] + j) * fftfreq;
652  if (hz < freqs[0] || hz > freqs[2]) {
653  E_FATAL
654  ("Failed to create filterbank, frequency range does not match. "
655  "Sample rate %f, FFT size %d, lowerf %f < freq %f > upperf %f.\n",
656  mel_fb->sampling_rate, mel_fb->fft_size, freqs[0], hz,
657  freqs[2]);
658  }
659  loslope = (hz - freqs[0]) / (freqs[1] - freqs[0]);
660  hislope = (freqs[2] - hz) / (freqs[2] - freqs[1]);
661  if (mel_fb->unit_area) {
662  loslope *= 2 / (freqs[2] - freqs[0]);
663  hislope *= 2 / (freqs[2] - freqs[0]);
664  }
665  if (loslope < hislope) {
666 #ifdef FIXED_POINT
667  mel_fb->filt_coeffs[n_coeffs] = fe_log(loslope);
668 #else
669  mel_fb->filt_coeffs[n_coeffs] = loslope;
670 #endif
671  }
672  else {
673 #ifdef FIXED_POINT
674  mel_fb->filt_coeffs[n_coeffs] = fe_log(hislope);
675 #else
676  mel_fb->filt_coeffs[n_coeffs] = hislope;
677 #endif
678  }
679  ++n_coeffs;
680  }
681  }
682 
683  return FE_SUCCESS;
684 }
685 
686 int32
687 fe_compute_melcosine(melfb_t * mel_fb)
688 {
689 
690  float64 freqstep;
691  int32 i, j;
692 
693  mel_fb->mel_cosine =
694  (mfcc_t **) ckd_calloc_2d(mel_fb->num_cepstra,
695  mel_fb->num_filters, sizeof(mfcc_t));
696 
697  freqstep = M_PI / mel_fb->num_filters;
698  /* NOTE: The first row vector is actually unnecessary but we leave
699  * it in to avoid confusion. */
700  for (i = 0; i < mel_fb->num_cepstra; i++) {
701  for (j = 0; j < mel_fb->num_filters; j++) {
702  float64 cosine;
703 
704  cosine = cos(freqstep * i * (j + 0.5));
705  mel_fb->mel_cosine[i][j] = FLOAT2COS(cosine);
706  }
707  }
708 
709  /* Also precompute normalization constants for unitary DCT. */
710  mel_fb->sqrt_inv_n = FLOAT2COS(sqrt(1.0 / mel_fb->num_filters));
711  mel_fb->sqrt_inv_2n = FLOAT2COS(sqrt(2.0 / mel_fb->num_filters));
712 
713  /* And liftering weights */
714  if (mel_fb->lifter_val) {
715  mel_fb->lifter =
716  calloc(mel_fb->num_cepstra, sizeof(*mel_fb->lifter));
717  for (i = 0; i < mel_fb->num_cepstra; ++i) {
718  mel_fb->lifter[i] = FLOAT2MFCC(1 + mel_fb->lifter_val / 2
719  * sin(i * M_PI /
720  mel_fb->lifter_val));
721  }
722  }
723 
724  return (0);
725 }
726 
727 static void
728 fe_pre_emphasis(int16 const *in, frame_t * out, int32 len,
729  float32 factor, int16 prior)
730 {
731  int i;
732 
733 #if defined(FIXED_POINT)
734  fixed32 fxd_alpha = FLOAT2FIX(factor);
735  out[0] = ((fixed32) in[0] << DEFAULT_RADIX) - (prior * fxd_alpha);
736  for (i = 1; i < len; ++i)
737  out[i] = ((fixed32) in[i] << DEFAULT_RADIX)
738  - (fixed32) in[i - 1] * fxd_alpha;
739 #else
740  out[0] = (frame_t) in[0] - (frame_t) prior *factor;
741  for (i = 1; i < len; i++)
742  out[i] = (frame_t) in[i] - (frame_t) in[i - 1] * factor;
743 #endif
744 }
745 
746 static void
747 fe_short_to_frame(int16 const *in, frame_t * out, int32 len)
748 {
749  int i;
750 
751 #if defined(FIXED_POINT)
752  for (i = 0; i < len; i++)
753  out[i] = (int32) in[i] << DEFAULT_RADIX;
754 #else /* FIXED_POINT */
755  for (i = 0; i < len; i++)
756  out[i] = (frame_t) in[i];
757 #endif /* FIXED_POINT */
758 }
759 
760 void
761 fe_create_hamming(window_t * in, int32 in_len)
762 {
763  int i;
764 
765  /* Symmetric, so we only create the first half of it. */
766  for (i = 0; i < in_len / 2; i++) {
767  float64 hamm;
768  hamm = (0.54 - 0.46 * cos(2 * M_PI * i /
769  ((float64) in_len - 1.0)));
770  in[i] = FLOAT2COS(hamm);
771  }
772 }
773 
774 static void
775 fe_hamming_window(frame_t * in, window_t * window, int32 in_len,
776  int32 remove_dc)
777 {
778  int i;
779 
780  if (remove_dc) {
781  frame_t mean = 0;
782 
783  for (i = 0; i < in_len; i++)
784  mean += in[i];
785  mean /= in_len;
786  for (i = 0; i < in_len; i++)
787  in[i] -= (frame_t) mean;
788  }
789 
790  for (i = 0; i < in_len / 2; i++) {
791  in[i] = COSMUL(in[i], window[i]);
792  in[in_len - 1 - i] = COSMUL(in[in_len - 1 - i], window[i]);
793  }
794 }
795 
796 static int
797 fe_spch_to_frame(fe_t * fe, int len)
798 {
799  /* Copy to the frame buffer. */
800  if (fe->pre_emphasis_alpha != 0.0) {
801  fe_pre_emphasis(fe->spch, fe->frame, len,
802  fe->pre_emphasis_alpha, fe->pre_emphasis_prior);
803  if (len >= fe->frame_shift)
804  fe->pre_emphasis_prior = fe->spch[fe->frame_shift - 1];
805  else
806  fe->pre_emphasis_prior = fe->spch[len - 1];
807  }
808  else
809  fe_short_to_frame(fe->spch, fe->frame, len);
810 
811  /* Zero pad up to FFT size. */
812  memset(fe->frame + len, 0, (fe->fft_size - len) * sizeof(*fe->frame));
813 
814  /* Window. */
815  fe_hamming_window(fe->frame, fe->hamming_window, fe->frame_size,
816  fe->remove_dc);
817 
818  return len;
819 }
820 
821 int
822 fe_read_frame(fe_t * fe, int16 const *in, int32 len)
823 {
824  int i;
825 
826  if (len > fe->frame_size)
827  len = fe->frame_size;
828 
829  /* Read it into the raw speech buffer. */
830  memcpy(fe->spch, in, len * sizeof(*in));
831  /* Swap and dither if necessary. */
832  if (fe->swap)
833  for (i = 0; i < len; ++i)
834  SWAP_INT16(&fe->spch[i]);
835  if (fe->dither)
836  for (i = 0; i < len; ++i)
837  fe->spch[i] += (int16) ((!(s3_rand_int31() % 4)) ? 1 : 0);
838 
839  return fe_spch_to_frame(fe, len);
840 }
841 
842 int
843 fe_shift_frame(fe_t * fe, int16 const *in, int32 len)
844 {
845  int offset, i;
846 
847  if (len > fe->frame_shift)
848  len = fe->frame_shift;
849  offset = fe->frame_size - fe->frame_shift;
850 
851  /* Shift data into the raw speech buffer. */
852  memmove(fe->spch, fe->spch + fe->frame_shift,
853  offset * sizeof(*fe->spch));
854  memcpy(fe->spch + offset, in, len * sizeof(*fe->spch));
855  /* Swap and dither if necessary. */
856  if (fe->swap)
857  for (i = 0; i < len; ++i)
858  SWAP_INT16(&fe->spch[offset + i]);
859  if (fe->dither)
860  for (i = 0; i < len; ++i)
861  fe->spch[offset + i]
862  += (int16) ((!(s3_rand_int31() % 4)) ? 1 : 0);
863 
864  return fe_spch_to_frame(fe, offset + len);
865 }
866 
870 void
871 fe_create_twiddle(fe_t * fe)
872 {
873  int i;
874 
875  for (i = 0; i < fe->fft_size / 4; ++i) {
876  float64 a = 2 * M_PI * i / fe->fft_size;
877 #if defined(FIXED_POINT)
878  fe->ccc[i] = FLOAT2COS(cos(a));
879  fe->sss[i] = FLOAT2COS(sin(a));
880 #else
881  fe->ccc[i] = cos(a);
882  fe->sss[i] = sin(a);
883 #endif
884  }
885 }
886 
887 
888 static int
889 fe_fft_real(fe_t * fe)
890 {
891  int i, j, k, m, n;
892  frame_t *x, xt;
893 
894  x = fe->frame;
895  m = fe->fft_order;
896  n = fe->fft_size;
897 
898  /* Bit-reverse the input. */
899  j = 0;
900  for (i = 0; i < n - 1; ++i) {
901  if (i < j) {
902  xt = x[j];
903  x[j] = x[i];
904  x[i] = xt;
905  }
906  k = n / 2;
907  while (k <= j) {
908  j -= k;
909  k /= 2;
910  }
911  j += k;
912  }
913 
914  /* Basic butterflies (2-point FFT, real twiddle factors):
915  * x[i] = x[i] + 1 * x[i+1]
916  * x[i+1] = x[i] + -1 * x[i+1]
917  */
918  for (i = 0; i < n; i += 2) {
919  xt = x[i];
920  x[i] = (xt + x[i + 1]);
921  x[i + 1] = (xt - x[i + 1]);
922  }
923 
924  /* The rest of the butterflies, in stages from 1..m */
925  for (k = 1; k < m; ++k) {
926  int n1, n2, n4;
927 
928  n4 = k - 1;
929  n2 = k;
930  n1 = k + 1;
931  /* Stride over each (1 << (k+1)) points */
932  for (i = 0; i < n; i += (1 << n1)) {
933  /* Basic butterfly with real twiddle factors:
934  * x[i] = x[i] + 1 * x[i + (1<<k)]
935  * x[i + (1<<k)] = x[i] + -1 * x[i + (1<<k)]
936  */
937  xt = x[i];
938  x[i] = (xt + x[i + (1 << n2)]);
939  x[i + (1 << n2)] = (xt - x[i + (1 << n2)]);
940 
941  /* The other ones with real twiddle factors:
942  * x[i + (1<<k) + (1<<(k-1))]
943  * = 0 * x[i + (1<<k-1)] + -1 * x[i + (1<<k) + (1<<k-1)]
944  * x[i + (1<<(k-1))]
945  * = 1 * x[i + (1<<k-1)] + 0 * x[i + (1<<k) + (1<<k-1)]
946  */
947  x[i + (1 << n2) + (1 << n4)] = -x[i + (1 << n2) + (1 << n4)];
948  x[i + (1 << n4)] = x[i + (1 << n4)];
949 
950  /* Butterflies with complex twiddle factors.
951  * There are (1<<k-1) of them.
952  */
953  for (j = 1; j < (1 << n4); ++j) {
954  frame_t cc, ss, t1, t2;
955  int i1, i2, i3, i4;
956 
957  i1 = i + j;
958  i2 = i + (1 << n2) - j;
959  i3 = i + (1 << n2) + j;
960  i4 = i + (1 << n2) + (1 << n2) - j;
961 
962  /*
963  * cc = real(W[j * n / (1<<(k+1))])
964  * ss = imag(W[j * n / (1<<(k+1))])
965  */
966  cc = fe->ccc[j << (m - n1)];
967  ss = fe->sss[j << (m - n1)];
968 
969  /* There are some symmetry properties which allow us
970  * to get away with only four multiplications here. */
971  t1 = COSMUL(x[i3], cc) + COSMUL(x[i4], ss);
972  t2 = COSMUL(x[i3], ss) - COSMUL(x[i4], cc);
973 
974  x[i4] = (x[i2] - t2);
975  x[i3] = (-x[i2] - t2);
976  x[i2] = (x[i1] - t1);
977  x[i1] = (x[i1] + t1);
978  }
979  }
980  }
981 
982  /* This isn't used, but return it for completeness. */
983  return m;
984 }
985 
986 static void
987 fe_spec_magnitude(fe_t * fe)
988 {
989  frame_t *fft;
990  powspec_t *spec;
991  int32 j, scale, fftsize;
992 
993  /* Do FFT and get the scaling factor back (only actually used in
994  * fixed-point). Note the scaling factor is expressed in bits. */
995  scale = fe_fft_real(fe);
996 
997  /* Convenience pointers to make things less awkward below. */
998  fft = fe->frame;
999  spec = fe->spec;
1000  fftsize = fe->fft_size;
1001 
1002  /* We need to scale things up the rest of the way to N. */
1003  scale = fe->fft_order - scale;
1004 
1005  /* The first point (DC coefficient) has no imaginary part */
1006  {
1007 #if defined(FIXED_POINT)
1008  spec[0] = FIXLN(abs(fft[0]) << scale) * 2;
1009 #else
1010  spec[0] = fft[0] * fft[0];
1011 #endif
1012  }
1013 
1014  for (j = 1; j <= fftsize / 2; j++) {
1015 #if defined(FIXED_POINT)
1016  int32 rr = FIXLN(abs(fft[j]) << scale) * 2;
1017  int32 ii = FIXLN(abs(fft[fftsize - j]) << scale) * 2;
1018  spec[j] = fe_log_add(rr, ii);
1019 #else
1020  spec[j] = fft[j] * fft[j] + fft[fftsize - j] * fft[fftsize - j];
1021 #endif
1022  }
1023 }
1024 
1025 static void
1026 fe_mel_spec(fe_t * fe)
1027 {
1028  int whichfilt;
1029  powspec_t *spec, *mfspec;
1030 
1031  /* Convenience poitners. */
1032  spec = fe->spec;
1033  mfspec = fe->mfspec;
1034  for (whichfilt = 0; whichfilt < fe->mel_fb->num_filters; whichfilt++) {
1035  int spec_start, filt_start, i;
1036 
1037  spec_start = fe->mel_fb->spec_start[whichfilt];
1038  filt_start = fe->mel_fb->filt_start[whichfilt];
1039 
1040 #ifdef FIXED_POINT
1041  mfspec[whichfilt] =
1042  spec[spec_start] + fe->mel_fb->filt_coeffs[filt_start];
1043  for (i = 1; i < fe->mel_fb->filt_width[whichfilt]; i++) {
1044  mfspec[whichfilt] = fe_log_add(mfspec[whichfilt],
1045  spec[spec_start + i] +
1046  fe->mel_fb->
1047  filt_coeffs[filt_start + i]);
1048  }
1049 #else /* !FIXED_POINT */
1050  mfspec[whichfilt] = 0;
1051  for (i = 0; i < fe->mel_fb->filt_width[whichfilt]; i++)
1052  mfspec[whichfilt] +=
1053  spec[spec_start + i] * fe->mel_fb->filt_coeffs[filt_start +
1054  i];
1055 #endif /* !FIXED_POINT */
1056  }
1057 
1058 }
1059 
1060 #define LOG_FLOOR 1e-4
1061 
1062 static void
1063 fe_mel_cep(fe_t * fe, mfcc_t * mfcep)
1064 {
1065  int32 i;
1066  powspec_t *mfspec;
1067 
1068  /* Convenience pointer. */
1069  mfspec = fe->mfspec;
1070 
1071  for (i = 0; i < fe->mel_fb->num_filters; ++i) {
1072 #ifndef FIXED_POINT /* It's already in log domain for fixed point */
1073  mfspec[i] = log(mfspec[i] + LOG_FLOOR);
1074 #endif /* !FIXED_POINT */
1075  }
1076 
1077  /* If we are doing LOG_SPEC, then do nothing. */
1078  if (fe->log_spec == RAW_LOG_SPEC) {
1079  for (i = 0; i < fe->feature_dimension; i++) {
1080  mfcep[i] = (mfcc_t) mfspec[i];
1081  }
1082  }
1083  /* For smoothed spectrum, do DCT-II followed by (its inverse) DCT-III */
1084  else if (fe->log_spec == SMOOTH_LOG_SPEC) {
1085  /* FIXME: This is probably broken for fixed-point. */
1086  fe_dct2(fe, mfspec, mfcep, 0);
1087  fe_dct3(fe, mfcep, mfspec);
1088  for (i = 0; i < fe->feature_dimension; i++) {
1089  mfcep[i] = (mfcc_t) mfspec[i];
1090  }
1091  }
1092  else if (fe->transform == DCT_II)
1093  fe_dct2(fe, mfspec, mfcep, FALSE);
1094  else if (fe->transform == DCT_HTK)
1095  fe_dct2(fe, mfspec, mfcep, TRUE);
1096  else
1097  fe_spec2cep(fe, mfspec, mfcep);
1098 
1099  return;
1100 }
1101 
1102 void
1103 fe_spec2cep(fe_t * fe, const powspec_t * mflogspec, mfcc_t * mfcep)
1104 {
1105  int32 i, j, beta;
1106 
1107  /* Compute C0 separately (its basis vector is 1) to avoid
1108  * costly multiplications. */
1109  mfcep[0] = mflogspec[0] / 2; /* beta = 0.5 */
1110  for (j = 1; j < fe->mel_fb->num_filters; j++)
1111  mfcep[0] += mflogspec[j]; /* beta = 1.0 */
1112  mfcep[0] /= (frame_t) fe->mel_fb->num_filters;
1113 
1114  for (i = 1; i < fe->num_cepstra; ++i) {
1115  mfcep[i] = 0;
1116  for (j = 0; j < fe->mel_fb->num_filters; j++) {
1117  if (j == 0)
1118  beta = 1; /* 0.5 */
1119  else
1120  beta = 2; /* 1.0 */
1121  mfcep[i] += COSMUL(mflogspec[j],
1122  fe->mel_fb->mel_cosine[i][j]) * beta;
1123  }
1124  /* Note that this actually normalizes by num_filters, like the
1125  * original Sphinx front-end, due to the doubled 'beta' factor
1126  * above. */
1127  mfcep[i] /= (frame_t) fe->mel_fb->num_filters * 2;
1128  }
1129 }
1130 
1131 void
1132 fe_dct2(fe_t * fe, const powspec_t * mflogspec, mfcc_t * mfcep, int htk)
1133 {
1134  int32 i, j;
1135 
1136  /* Compute C0 separately (its basis vector is 1) to avoid
1137  * costly multiplications. */
1138  mfcep[0] = mflogspec[0];
1139  for (j = 1; j < fe->mel_fb->num_filters; j++)
1140  mfcep[0] += mflogspec[j];
1141  if (htk)
1142  mfcep[0] = COSMUL(mfcep[0], fe->mel_fb->sqrt_inv_2n);
1143  else /* sqrt(1/N) = sqrt(2/N) * 1/sqrt(2) */
1144  mfcep[0] = COSMUL(mfcep[0], fe->mel_fb->sqrt_inv_n);
1145 
1146  for (i = 1; i < fe->num_cepstra; ++i) {
1147  mfcep[i] = 0;
1148  for (j = 0; j < fe->mel_fb->num_filters; j++) {
1149  mfcep[i] += COSMUL(mflogspec[j], fe->mel_fb->mel_cosine[i][j]);
1150  }
1151  mfcep[i] = COSMUL(mfcep[i], fe->mel_fb->sqrt_inv_2n);
1152  }
1153 }
1154 
1155 void
1156 fe_lifter(fe_t * fe, mfcc_t * mfcep)
1157 {
1158  int32 i;
1159 
1160  if (fe->mel_fb->lifter_val == 0)
1161  return;
1162 
1163  for (i = 0; i < fe->num_cepstra; ++i) {
1164  mfcep[i] = MFCCMUL(mfcep[i], fe->mel_fb->lifter[i]);
1165  }
1166 }
1167 
1168 void
1169 fe_dct3(fe_t * fe, const mfcc_t * mfcep, powspec_t * mflogspec)
1170 {
1171  int32 i, j;
1172 
1173  for (i = 0; i < fe->mel_fb->num_filters; ++i) {
1174  mflogspec[i] = COSMUL(mfcep[0], SQRT_HALF);
1175  for (j = 1; j < fe->num_cepstra; j++) {
1176  mflogspec[i] += COSMUL(mfcep[j], fe->mel_fb->mel_cosine[j][i]);
1177  }
1178  mflogspec[i] = COSMUL(mflogspec[i], fe->mel_fb->sqrt_inv_2n);
1179  }
1180 }
1181 
1182 void
1183 fe_write_frame(fe_t * fe, mfcc_t * feat, int32 store_pcm)
1184 {
1185  int32 is_speech;
1186 
1187  fe_spec_magnitude(fe);
1188  fe_mel_spec(fe);
1189  fe_track_snr(fe, &is_speech);
1190  fe_mel_cep(fe, feat);
1191  fe_lifter(fe, feat);
1192  fe_vad_hangover(fe, feat, is_speech, store_pcm);
1193 }
1194 
1195 
1196 void *
1197 fe_create_2d(int32 d1, int32 d2, int32 elem_size)
1198 {
1199  return (void *) ckd_calloc_2d(d1, d2, elem_size);
1200 }
1201 
1202 void
1203 fe_free_2d(void *arr)
1204 {
1205  ckd_free_2d((void **) arr);
1206 }
ckd_calloc_2d
#define ckd_calloc_2d(d1, d2, sz)
Macro for ckd_calloc_2d
Definition: ckd_alloc.h:270
ckd_calloc
#define ckd_calloc(n, sz)
Macros to simplify the use of above functions.
Definition: ckd_alloc.h:248
melfb_s
Base Struct to hold all structure for MFCC computation.
Definition: fe_internal.h:75
ckd_alloc.h
Sphinx&#39;s memory allocation/deallocation routines.
prim_type.h
Basic type definitions used in Sphinx.
err.h
Implementation of logging routines.
E_WARN
#define E_WARN(...)
Print warning message to error log.
Definition: err.h:109
genrand.h
High performance prortable random generator created by Takuji Nishimura and Makoto Matsumoto...
ckd_malloc
#define ckd_malloc(sz)
Macro for ckd_malloc
Definition: ckd_alloc.h:253
ckd_free_2d
SPHINXBASE_EXPORT void ckd_free_2d(void *ptr)
Free a 2-D array (ptr) previously allocated by ckd_calloc_2d.
Definition: ckd_alloc.c:255
E_FATAL
#define E_FATAL(...)
Exit with non-zero status after error message.
Definition: err.h:81
fe_s
Structure for the front-end computation.
Definition: fe_internal.h:117