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
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 */
93static 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
256static const int fe_logadd_table_size =
257 sizeof(fe_logadd_table) / sizeof(fe_logadd_table[0]);
258
259fixed32
260fe_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 */
301static const uint16 fe_logsub_table[] = {
3021, 3, 5, 7, 9, 11, 13, 15, 17, 19,
30321, 23, 25, 27, 29, 31, 33, 35, 37, 39,
30441, 43, 45, 47, 49, 51, 53, 55, 56, 58,
30560, 62, 64, 66, 68, 70, 72, 74, 76, 78,
30680, 82, 84, 86, 88, 90, 92, 94, 95, 97,
30799, 101, 103, 105, 107, 109, 111, 113, 115, 117,
308119, 121, 122, 124, 126, 128, 130, 132, 134, 136,
309138, 140, 142, 143, 145, 147, 149, 151, 153, 155,
310157, 159, 161, 162, 164, 166, 168, 170, 172, 174,
311176, 178, 179, 181, 183, 185, 187, 189, 191, 193,
312194, 196, 198, 200, 202, 204, 206, 207, 209, 211,
313213, 215, 217, 219, 220, 222, 224, 226, 228, 230,
314232, 233, 235, 237, 239, 241, 243, 244, 246, 248,
315250, 252, 254, 255, 257, 259, 261, 263, 265, 266,
316268, 270, 272, 274, 275, 277, 279, 281, 283, 284,
317286, 288, 290, 292, 294, 295, 297, 299, 301, 302,
318304, 306, 308, 310, 311, 313, 315, 317, 319, 320,
319322, 324, 326, 327, 329, 331, 333, 335, 336, 338,
320340, 342, 343, 345, 347, 349, 350, 352, 354, 356,
321357, 359, 361, 363, 364, 366, 368, 370, 371, 373,
322375, 377, 378, 380, 382, 384, 385, 387, 389, 391,
323392, 394, 396, 397, 399, 401, 403, 404, 406, 408,
324410, 411, 413, 415, 416, 418, 420, 422, 423, 425,
325427, 428, 430, 432, 433, 435, 437, 439, 440, 442,
326444, 445, 447, 449, 450, 452, 454, 455, 457, 459,
327460, 462, 464, 465, 467, 469, 471, 472, 474, 476,
328477, 479, 481, 482, 484, 486, 487, 489, 490, 492,
329494, 495, 497, 499, 500, 502, 504, 505, 507, 509,
330510, 512, 514, 515, 517, 518, 520, 522, 523, 525,
331527, 528, 530, 532, 533, 535, 536, 538, 540, 541,
332543, 544, 546, 548, 549, 551, 553, 554, 556, 557,
333559, 561, 562, 564, 565, 567, 569, 570, 572, 573,
334575, 577, 578, 580, 581, 583, 585, 586, 588, 589,
335591, 592, 594, 596, 597, 599, 600, 602, 603, 605,
336607, 608, 610, 611, 613, 614, 616, 618, 619, 621,
337622, 624, 625, 627, 628, 630, 632, 633, 635, 636,
338638, 639, 641, 642, 644, 645, 647, 649, 650, 652,
339653, 655, 656, 658, 659, 661, 662, 664, 665, 667,
340668, 670, 671, 673, 674, 676, 678, 679, 681, 682,
341684, 685, 687, 688, 690, 691, 693, 694, 696, 697,
342699, 700, 702, 703, 705, 706, 708, 709, 711, 712,
343714, 715, 717, 718, 719, 721, 722, 724, 725, 727,
344728, 730, 731, 733, 734, 736, 737, 739, 740, 742,
345743, 745, 746, 747, 749, 750, 752, 753, 755, 756,
346758, 759, 761, 762, 763, 765, 766, 768, 769, 771,
347772, 774, 775, 776, 778, 779, 781, 782, 784, 785,
348786, 788, 789, 791, 792, 794, 795, 796, 798, 799,
349801, 802, 804, 805, 806, 808, 809, 811, 812, 813,
350815, 816, 818, 819, 820, 822, 823, 825, 826, 827,
351829, 830, 832, 833, 834, 836, 837, 839, 840, 841,
352843, 844, 846, 847, 848, 850, 851, 852, 854, 855,
353857, 858, 859, 861, 862, 863, 865, 866, 868, 869,
354870, 872, 873, 874, 876, 877, 878, 880, 881, 883,
355884, 885, 887, 888, 889, 891, 892, 893, 895, 896,
356897, 899, 900, 901, 903, 904, 905, 907, 908, 909,
357911, 912, 913, 915, 916, 917, 919, 920, 921, 923,
358924, 925, 927, 928, 929, 931, 932, 933, 935, 936,
359937, 939, 940, 941, 942, 944, 945, 946, 948, 949,
360950, 952, 953, 954, 956, 957, 958, 959, 961, 962,
361963, 965, 966, 967, 968, 970, 971, 972, 974, 975,
362976, 977, 979, 980, 981, 983, 984, 985, 986, 988,
363989, 990, 991, 993, 994, 995, 997, 998, 999, 1000,
3641002, 1003, 1004, 1005, 1007, 1008, 1009, 1010, 1012, 1013,
3651014, 1015, 1017, 1018, 1019, 1020, 1022, 1023, 1024, 1025,
3661027, 1028, 1029, 1030, 1032, 1033, 1034, 1035, 1037, 1038,
3671039, 1040, 1041, 1043, 1044, 1045, 1046, 1048, 1049, 1050,
3681051, 1052, 1054, 1055, 1056, 1057, 1059, 1060, 1061, 1062,
3691063, 1065, 1066, 1067, 1068, 1069, 1071, 1072, 1073, 1074,
3701076, 1077, 1078, 1079, 1080, 1082, 1083, 1084, 1085, 1086,
3711087, 1089, 1090, 1091, 1092, 1093, 1095, 1096, 1097, 1098,
3721099, 1101, 1102, 1103, 1104, 1105, 1106, 1108, 1109, 1110,
3731111, 1112, 1114, 1115, 1116, 1117, 1118, 1119, 1121, 1122,
3741123, 1124, 1125, 1126, 1128, 1129, 1130, 1131, 1132, 1133,
3751135, 1136, 1137, 1138, 1139, 1140, 1141, 1143, 1144, 1145,
3761146, 1147, 1148, 1149, 1151, 1152, 1153, 1154, 1155, 1156,
3771157, 1159, 1160, 1161, 1162, 1163, 1164, 1165, 1167, 1168,
3781169, 1170, 1171, 1172, 1173, 1174, 1176, 1177, 1178, 1179,
3791180, 1181, 1182, 1183, 1185, 1186, 1187, 1188, 1189, 1190,
3801191, 1192, 1193, 1195, 1196, 1197, 1198, 1199, 1200, 1201,
3811202, 1203, 1205, 1206, 1207, 1208, 1209, 1210, 1211, 1212,
3821213, 1214, 1216, 1217, 1218, 1219, 1220, 1221, 1222, 1223,
3831224, 1225, 1226, 1228, 1229, 1230, 1231, 1232, 1233, 1234,
3841235, 1236, 1237, 1238, 1239, 1240, 1242, 1243, 1244, 1245,
3851246, 1247, 1248, 1249, 1250, 1251, 1252, 1253, 1254, 1255,
3861256, 1258, 1259, 1260, 1261, 1262, 1263, 1264, 1265, 1266,
3871267, 1268, 1269, 1270, 1271, 1272, 1273, 1274, 1275, 1277,
3881278, 1279, 1280, 1281, 1282, 1283, 1284, 1285, 1286, 1287,
3891288, 1289, 1290, 1291, 1292, 1293, 1294, 1295, 1296, 1297,
3901298, 1299, 1300, 1301, 1302, 1303, 1305, 1306, 1307, 1308,
3911309, 1310, 1311, 1312, 1313, 1314, 1315, 1316, 1317, 1318,
3921319, 1320, 1321, 1322, 1323, 1324, 1325, 1326, 1327, 1328,
3931329, 1330, 1331, 1332, 1333, 1334, 1335, 1336, 1337, 1338,
3941339, 1340, 1341, 1342, 1343, 1344, 1345, 1346, 1347, 1348,
3951349, 1350, 1351, 1352, 1353, 1354, 1355, 1356, 1357, 1358,
3961359, 1360, 1361, 1362, 1363, 1364, 1365, 1366, 1367, 1368,
3971369, 1370, 1371, 1372, 1372, 1373, 1374, 1375, 1376, 1377,
3981378, 1379, 1380, 1381, 1382, 1383, 1384, 1385, 1386, 1387,
3991388, 1389, 1390, 1391, 1392, 1393, 1394, 1395, 1396, 1397,
4001398, 1399, 1399, 1400, 1401, 1402, 1403, 1404, 1405, 1406,
4011407, 1408, 1409, 1410, 1411, 1412, 1413, 1414, 1415, 1416,
4021417, 1418, 1418, 1419, 1420, 1421, 1422, 1423, 1424, 1425,
4031426, 1427, 1428, 1429, 1430, 1431, 1432, 1432, 1433, 1434,
4041435, 1436, 1437, 1438, 1439, 1440, 1441, 1442, 1443, 1444,
4051444, 1445, 1446, 1447, 1448, 1449, 1450, 1451, 1452, 1453,
4061454, 1455, 1455, 1456, 1457, 1458, 1459, 1460, 1461, 1462,
4071463, 1464, 1465, 1466, 1466, 1467, 1468, 1469, 1470, 1471,
4081472, 1473, 1474, 1475, 1475, 1476, 1477, 1478, 1479, 1480,
4091481, 1482, 1483, 1483, 1484, 1485, 1486, 1487, 1488, 1489,
4101490, 1491, 1491, 1492, 1493, 1494, 1495, 1496, 1497, 1498,
4111499, 1499, 1500, 1501, 1502, 1503, 1504, 1505, 1506, 1506,
4121507, 1508, 1509, 1510, 1511, 1512, 1513, 1513, 1514, 1515,
4131516, 1517, 1518, 1519, 1520, 1520, 1521, 1522, 1523, 1524,
4141525, 1526, 1526, 1527, 1528, 1529, 1530, 1531, 1532, 1532,
4151533, 1534, 1535, 1536, 1537, 1538, 1538, 1539, 1540, 1541,
4161542, 1543, 1544, 1544, 1545, 1546, 1547, 1548, 1549, 1550,
4171550, 1551, 1552, 1553, 1554, 1555, 1555, 1556, 1557, 1558,
4181559, 1560, 1560, 1561, 1562, 1563, 1564, 1565, 1565, 1566,
4191567, 1568, 1569, 1570, 1570, 1571, 1572, 1573, 1574, 1575,
4201575, 1576, 1577, 1578, 1579, 1580, 1580, 1581, 1582, 1583,
4211584, 1584, 1585, 1586, 1587, 1588, 1589, 1589, 1590, 1591,
4221592, 1593, 1593, 1594, 1595, 1596, 1597, 1598, 1598, 1599,
4231600, 1601, 1602, 1602, 1603, 1604, 1605, 1606, 1606, 1607,
4241608, 1609, 1610, 1610, 1611, 1612, 1613, 1614, 1614, 1615,
4251616, 1617, 1618, 1618, 1619, 1620, 1621, 1622, 1622, 1623,
4261624, 1625, 1626, 1626, 1627, 1628, 1629, 1630, 1630, 1631,
4271632, 1633, 1634, 1634, 1635, 1636, 1637, 1637, 1638, 1639,
4281640, 1641, 1641, 1642, 1643, 1644, 1645, 1645, 1646, 1647,
4291648, 1648, 1649, 1650, 1651, 1652, 1652, 1653, 1654, 1655,
4301655, 1656, 1657, 1658, 1658, 1659, 1660, 1661, 1662, 1662,
4311663, 1664, 1665, 1665, 1666, 1667, 1668, 1668, 1669, 1670,
4321671, 1671, 1672, 1673, 1674, 1675, 1675, 1676, 1677, 1678,
4331678, 1679, 1680, 1681, 1681, 1682, 1683, 1684, 1684, 1685,
4341686, 1687, 1687, 1688, 1689, 1690, 1690, 1691, 1692, 1693,
4351693, 1694, 1695, 1696, 1696, 1697, 1698, 1699, 1699, 1700,
4361701, 1702, 1702, 1703, 1704, 1705, 1705, 1706, 1707, 1707,
4371708, 1709, 1710, 1710, 1711, 1712, 1713, 1713, 1714, 1715,
4381716, 1716, 1717, 1718, 1718, 1719, 1720, 1721, 1721, 1722,
4391723, 1724, 1724, 1725, 1726, 1727, 1727, 1728, 1729, 1729,
4401730, 1731, 1732, 1732, 1733, 1734, 1734, 1735, 1736, 1737,
4411737, 1738, 1739, 1740, 1740, 1741, 1742, 1742, 1743, 1744,
4421745, 1745, 1746, 1747, 1747, 1748, 1749, 1749, 1750, 1751,
4431752, 1752, 1753, 1754, 1754, 1755, 1756, 1757, 1757, 1758,
4441759, 1759, 1760, 1761, 1762, 1762, 1763, 1764, 1764, 1765,
4451766, 1766, 1767, 1768, 1769, 1769, 1770, 1771, 1771, 1772,
4461773, 1773, 1774, 1775, 1776, 1776, 1777, 1778, 1778, 1779,
4471780, 1780, 1781, 1782, 1782, 1783, 1784, 1784, 1785, 1786,
4481787, 1787, 1788, 1789, 1789, 1790, 1791, 1791, 1792, 1793,
4491793, 1794, 1795, 1795, 1796, 1797, 1798, 1798, 1799, 1800,
4501800, 1801, 1802, 1802, 1803, 1804, 1804, 1805, 1806, 1806,
4511807, 1808, 1808, 1809, 1810, 1810, 1811, 1812, 1812, 1813,
4521814, 1814, 1815, 1816, 1816, 1817, 1818, 1818, 1819, 1820,
4531820, 1821, 1822, 1822, 1823, 1824, 1824, 1825, 1826, 1826,
4541827, 1828, 1828, 1829, 1830, 1830, 1831, 1832, 1832, 1833,
4551834, 1834, 1835, 1836, 1836, 1837, 1838, 1838, 1839, 1840,
4561840, 1841, 1842, 1842, 1843, 1844, 1844, 1845, 1845, 1846,
4571847, 1847, 1848, 1849, 1849, 1850, 1851, 1851, 1852, 1853,
4581853, 1854, 1855, 1855, 1856, 1857, 1857, 1858, 1858, 1859,
4591860, 1860, 1861, 1862, 1862, 1863, 1864, 1864, 1865, 1866,
4601866, 1867, 1867, 1868, 1869, 1869, 1870, 1871, 1871, 1872,
4611873, 1873, 1874, 1874, 1875, 1876, 1876, 1877, 1878, 1878,
4621879, 1879, 1880, 1881, 1881, 1882, 1883, 1883, 1884, 1885,
4631885, 1886, 1886, 1887, 1888, 1888, 1889, 1890, 1890, 1891,
4641891, 1892, 1893, 1893, 1894, 1895, 1895, 1896, 1896, 1897,
4651898, 1898, 1899, 1900, 1900, 1901, 1901, 1902, 1903, 1903,
4661904, 1904, 1905, 1906, 1906, 1907, 1908, 1908, 1909, 1909,
4671910, 1911, 1911, 1912, 1912, 1913, 1914, 1914, 1915, 1916,
4681916, 1917, 1917, 1918, 1919, 1919, 1920, 1920, 1921, 1922,
4691922, 1923, 1923, 1924, 1925, 1925, 1926, 1926, 1927, 1928,
4701928, 1929, 1929, 1930, 1931, 1931, 1932, 1932, 1933, 1934,
4711934, 1935, 1935, 1936, 1937, 1937, 1938, 1938, 1939, 1940,
4721940, 1941, 1941, 1942, 1943, 1943, 1944, 1944, 1945, 1946,
4731946, 1947, 1947, 1948, 1949, 1949, 1950, 1950, 1951, 1952,
4741952, 1953, 1953, 1954, 1955, 1955, 1956, 1956, 1957, 1957,
4751958, 1959, 1959, 1960, 1960, 1961, 1962, 1962, 1963, 1963,
4761964, 1964, 1965, 1966, 1966, 1967, 1967, 1968, 1969, 1969,
4771970, 1970, 1971, 1971, 1972, 1973, 1973, 1974, 1974, 1975,
4781976, 1976, 1977, 1977, 1978, 1978, 1979, 1980, 1980, 1981,
4791981, 1982, 1982, 1983, 1984, 1984, 1985, 1985, 1986, 1986,
4801987, 1988, 1988, 1989, 1989, 1990, 1990, 1991, 1992, 1992,
4811993, 1993, 1994, 1994, 1995, 1996, 1996, 1997, 1997, 1998,
4821998, 1999, 1999, 2000, 2001, 2001, 2002, 2002, 2003, 2003,
4832004, 2005, 2005, 2006, 2006, 2007, 2007, 2008, 2008, 2009,
4842010, 2010, 2011, 2011, 2012, 2012, 2013, 2014, 2014, 2015,
4852015, 2016, 2016, 2017, 2017, 2018, 2019, 2019, 2020, 2020,
4862021, 2021, 2022, 2022, 2023, 2023, 2024, 2025, 2025, 2026,
4872026, 2027, 2027, 2028, 2028, 2029, 2030, 2030, 2031, 2031,
4882032, 2032, 2033, 2033, 2034, 2034, 2035, 2036, 2036, 2037,
4892037, 2038, 2038, 2039, 2039, 2040, 2040, 2041, 2042, 2042,
4902043, 2043, 2044, 2044, 2045, 2045, 2046, 2046, 2047, 2048,
4912048, 2049, 2049, 2050, 2050, 2051, 2051, 2052, 2052, 2053,
4922053, 2054, 2054, 2055, 2056, 2056, 2057, 2057, 2058, 2058,
4932059, 2059, 2060, 2060, 2061, 2061, 2062, 2062, 2063, 2064,
4942064, 2065, 2065, 2066, 2066, 2067, 2067, 2068, 2068, 2069,
4952069, 2070, 2070, 2071, 2072, 2072, 2073, 2073, 2074, 2074,
4962075, 2075, 2076, 2076, 2077, 2077, 2078, 2078, 2079, 2079,
4972080, 2080, 2081
498};
499
500static const int fe_logsub_table_size =
501 sizeof(fe_logsub_table) / sizeof(fe_logsub_table[0]);
502
503fixed32
504fe_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
525static fixed32
526fe_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
537static float32
538fe_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
545static float32
546fe_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
552int32
553fe_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
686int32
687fe_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
727static void
728fe_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
746static void
747fe_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
760void
761fe_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
774static void
775fe_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
796static int
797fe_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
821int
822fe_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
842int
843fe_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
870void
871fe_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
888static int
889fe_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
986static void
987fe_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
1025static void
1026fe_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
1062static void
1063fe_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
1102void
1103fe_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
1131void
1132fe_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
1155void
1156fe_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
1168void
1169fe_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
1182void
1183fe_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
1196void *
1197fe_create_2d(int32 d1, int32 d2, int32 elem_size)
1198{
1199 return (void *) ckd_calloc_2d(d1, d2, elem_size);
1200}
1201
1202void
1203fe_free_2d(void *arr)
1204{
1205 ckd_free_2d((void **) arr);
1206}
Sphinx's memory allocation/deallocation routines.
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
#define ckd_malloc(sz)
Macro for ckd_malloc
Definition: ckd_alloc.h:253
#define ckd_calloc_2d(d1, d2, sz)
Macro for ckd_calloc_2d
Definition: ckd_alloc.h:270
#define ckd_calloc(n, sz)
Macros to simplify the use of above functions.
Definition: ckd_alloc.h:248
Implementation of logging routines.
#define E_FATAL(...)
Exit with non-zero status after error message.
Definition: err.h:81
#define E_WARN(...)
Print warning message to error log.
Definition: err.h:109
High performance prortable random generator created by Takuji Nishimura and Makoto Matsumoto.
Basic type definitions used in Sphinx.
Structure for the front-end computation.
Definition: fe_internal.h:117
Base Struct to hold all structure for MFCC computation.
Definition: fe_internal.h:75