diff options
Diffstat (limited to 'media/sphinxbase/src/libsphinxbase/fe/fe_sigproc.c')
| -rw-r--r-- | media/sphinxbase/src/libsphinxbase/fe/fe_sigproc.c | 1377 |
1 files changed, 1377 insertions, 0 deletions
diff --git a/media/sphinxbase/src/libsphinxbase/fe/fe_sigproc.c b/media/sphinxbase/src/libsphinxbase/fe/fe_sigproc.c new file mode 100644 index 000000000..577640f62 --- /dev/null +++ b/media/sphinxbase/src/libsphinxbase/fe/fe_sigproc.c @@ -0,0 +1,1377 @@ +/* -*- c-basic-offset: 4; indent-tabs-mode: nil -*- */ +/* ==================================================================== + * Copyright (c) 1996-2004 Carnegie Mellon University. All rights + * reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in + * the documentation and/or other materials provided with the + * distribution. + * + * This work was supported in part by funding from the Defense Advanced + * Research Projects Agency and the National Science Foundation of the + * United States of America, and the CMU Sphinx Speech Consortium. + * + * THIS SOFTWARE IS PROVIDED BY CARNEGIE MELLON UNIVERSITY ``AS IS'' AND + * ANY EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, + * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR + * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL CARNEGIE MELLON UNIVERSITY + * NOR ITS EMPLOYEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT + * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + * + * ==================================================================== + * + */ + +#include <stdio.h> +#include <math.h> +#include <string.h> +#include <stdlib.h> +#include <assert.h> + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#ifdef _MSC_VER +#pragma warning (disable: 4244) +#endif + +/** + * Windows math.h does not contain M_PI + */ +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +#include "sphinxbase/prim_type.h" +#include "sphinxbase/ckd_alloc.h" +#include "sphinxbase/byteorder.h" +#include "sphinxbase/fixpoint.h" +#include "sphinxbase/fe.h" +#include "sphinxbase/genrand.h" +#include "sphinxbase/err.h" + +#include "fe_internal.h" +#include "fe_warp.h" + +/* Use extra precision for cosines, Hamming window, pre-emphasis + * coefficient, twiddle factors. */ +#ifdef FIXED_POINT +#define FLOAT2COS(x) FLOAT2FIX_ANY(x,30) +#define COSMUL(x,y) FIXMUL_ANY(x,y,30) +#else +#define FLOAT2COS(x) (x) +#define COSMUL(x,y) ((x)*(y)) +#endif + +#ifdef FIXED_POINT + +/* Internal log-addition table for natural log with radix point at 8 + * bits. Each entry is 256 * log(1 + e^{-n/256}). This is used in the + * log-add computation: + * + * e^z = e^x + e^y + * e^z = e^x(1 + e^{y-x}) = e^y(1 + e^{x-y}) + * z = x + log(1 + e^{y-x}) = y + log(1 + e^{x-y}) + * + * So when y > x, z = y + logadd_table[-(x-y)] + * when x > y, z = x + logadd_table[-(y-x)] + */ +static const unsigned char fe_logadd_table[] = { + 177, 177, 176, 176, 175, 175, 174, 174, 173, 173, + 172, 172, 172, 171, 171, 170, 170, 169, 169, 168, + 168, 167, 167, 166, 166, 165, 165, 164, 164, 163, + 163, 162, 162, 161, 161, 161, 160, 160, 159, 159, + 158, 158, 157, 157, 156, 156, 155, 155, 155, 154, + 154, 153, 153, 152, 152, 151, 151, 151, 150, 150, + 149, 149, 148, 148, 147, 147, 147, 146, 146, 145, + 145, 144, 144, 144, 143, 143, 142, 142, 141, 141, + 141, 140, 140, 139, 139, 138, 138, 138, 137, 137, + 136, 136, 136, 135, 135, 134, 134, 134, 133, 133, + 132, 132, 131, 131, 131, 130, 130, 129, 129, 129, + 128, 128, 128, 127, 127, 126, 126, 126, 125, 125, + 124, 124, 124, 123, 123, 123, 122, 122, 121, 121, + 121, 120, 120, 119, 119, 119, 118, 118, 118, 117, + 117, 117, 116, 116, 115, 115, 115, 114, 114, 114, + 113, 113, 113, 112, 112, 112, 111, 111, 110, 110, + 110, 109, 109, 109, 108, 108, 108, 107, 107, 107, + 106, 106, 106, 105, 105, 105, 104, 104, 104, 103, + 103, 103, 102, 102, 102, 101, 101, 101, 100, 100, + 100, 99, 99, 99, 98, 98, 98, 97, 97, 97, + 96, 96, 96, 96, 95, 95, 95, 94, 94, 94, + 93, 93, 93, 92, 92, 92, 92, 91, 91, 91, + 90, 90, 90, 89, 89, 89, 89, 88, 88, 88, + 87, 87, 87, 87, 86, 86, 86, 85, 85, 85, + 85, 84, 84, 84, 83, 83, 83, 83, 82, 82, + 82, 82, 81, 81, 81, 80, 80, 80, 80, 79, + 79, 79, 79, 78, 78, 78, 78, 77, 77, 77, + 77, 76, 76, 76, 75, 75, 75, 75, 74, 74, + 74, 74, 73, 73, 73, 73, 72, 72, 72, 72, + 71, 71, 71, 71, 71, 70, 70, 70, 70, 69, + 69, 69, 69, 68, 68, 68, 68, 67, 67, 67, + 67, 67, 66, 66, 66, 66, 65, 65, 65, 65, + 64, 64, 64, 64, 64, 63, 63, 63, 63, 63, + 62, 62, 62, 62, 61, 61, 61, 61, 61, 60, + 60, 60, 60, 60, 59, 59, 59, 59, 59, 58, + 58, 58, 58, 58, 57, 57, 57, 57, 57, 56, + 56, 56, 56, 56, 55, 55, 55, 55, 55, 54, + 54, 54, 54, 54, 53, 53, 53, 53, 53, 52, + 52, 52, 52, 52, 52, 51, 51, 51, 51, 51, + 50, 50, 50, 50, 50, 50, 49, 49, 49, 49, + 49, 49, 48, 48, 48, 48, 48, 48, 47, 47, + 47, 47, 47, 47, 46, 46, 46, 46, 46, 46, + 45, 45, 45, 45, 45, 45, 44, 44, 44, 44, + 44, 44, 43, 43, 43, 43, 43, 43, 43, 42, + 42, 42, 42, 42, 42, 41, 41, 41, 41, 41, + 41, 41, 40, 40, 40, 40, 40, 40, 40, 39, + 39, 39, 39, 39, 39, 39, 38, 38, 38, 38, + 38, 38, 38, 37, 37, 37, 37, 37, 37, 37, + 37, 36, 36, 36, 36, 36, 36, 36, 35, 35, + 35, 35, 35, 35, 35, 35, 34, 34, 34, 34, + 34, 34, 34, 34, 33, 33, 33, 33, 33, 33, + 33, 33, 32, 32, 32, 32, 32, 32, 32, 32, + 32, 31, 31, 31, 31, 31, 31, 31, 31, 31, + 30, 30, 30, 30, 30, 30, 30, 30, 30, 29, + 29, 29, 29, 29, 29, 29, 29, 29, 28, 28, + 28, 28, 28, 28, 28, 28, 28, 28, 27, 27, + 27, 27, 27, 27, 27, 27, 27, 27, 26, 26, + 26, 26, 26, 26, 26, 26, 26, 26, 25, 25, + 25, 25, 25, 25, 25, 25, 25, 25, 25, 24, + 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, + 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, + 23, 23, 22, 22, 22, 22, 22, 22, 22, 22, + 22, 22, 22, 22, 21, 21, 21, 21, 21, 21, + 21, 21, 21, 21, 21, 21, 21, 20, 20, 20, + 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, + 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, + 19, 19, 19, 19, 18, 18, 18, 18, 18, 18, + 18, 18, 18, 18, 18, 18, 18, 18, 18, 17, + 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, + 17, 17, 17, 17, 16, 16, 16, 16, 16, 16, + 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, + 16, 15, 15, 15, 15, 15, 15, 15, 15, 15, + 15, 15, 15, 15, 15, 15, 15, 15, 14, 14, + 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, + 14, 14, 14, 14, 14, 14, 14, 13, 13, 13, + 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, + 13, 13, 13, 13, 13, 13, 13, 12, 12, 12, + 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, + 12, 12, 12, 12, 12, 12, 12, 12, 12, 11, + 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, + 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, + 11, 11, 11, 10, 10, 10, 10, 10, 10, 10, + 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, + 10, 10, 10, 10, 10, 10, 10, 10, 10, 9, + 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, + 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, + 9, 9, 9, 9, 9, 9, 9, 9, 8, 8, + 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, + 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, + 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, + 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, + 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, + 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, + 7, 7, 7, 7, 7, 7, 7, 7, 6, 6, + 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, + 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, + 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, + 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, + 6, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, + 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, + 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, + 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, + 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, + 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, + 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, + 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, + 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, + 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, + 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, + 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 0 +}; + +static const int fe_logadd_table_size = + sizeof(fe_logadd_table) / sizeof(fe_logadd_table[0]); + +fixed32 +fe_log_add(fixed32 x, fixed32 y) +{ + fixed32 d, r; + + if (x > y) { + d = (x - y) >> (DEFAULT_RADIX - 8); + r = x; + } + else { + d = (y - x) >> (DEFAULT_RADIX - 8); + r = y; + } + + if (r <= MIN_FIXLOG) + return MIN_FIXLOG; + else if (d > fe_logadd_table_size - 1) + return r; + else { + r += ((fixed32) fe_logadd_table[d] << (DEFAULT_RADIX - 8)); +/* printf("%d - %d = %d | %f - %f = %f | %f - %f = %f\n", + x, y, r, FIX2FLOAT(x), FIX2FLOAT(y), FIX2FLOAT(r), + exp(FIX2FLOAT(x)), exp(FIX2FLOAT(y)), exp(FIX2FLOAT(r))); +*/ + return r; + } +} + +/* + * log_sub for spectral subtraction, similar to logadd but we had + * to smooth function around zero with fixlog in order to improve + * table interpolation properties + * + * The table is created with the file included into distribution + * + * e^z = e^x - e^y + * e^z = e^x (1 - e^(-(x - y))) + * z = x + log(1 - e^(-(x - y))) + * z = x + fixlog(a) + (log(1 - e^(- a)) - log(a)) + * + * Input radix is 8 output radix is 10 + */ +static const uint16 fe_logsub_table[] = { +1, 3, 5, 7, 9, 11, 13, 15, 17, 19, +21, 23, 25, 27, 29, 31, 33, 35, 37, 39, +41, 43, 45, 47, 49, 51, 53, 55, 56, 58, +60, 62, 64, 66, 68, 70, 72, 74, 76, 78, +80, 82, 84, 86, 88, 90, 92, 94, 95, 97, +99, 101, 103, 105, 107, 109, 111, 113, 115, 117, +119, 121, 122, 124, 126, 128, 130, 132, 134, 136, +138, 140, 142, 143, 145, 147, 149, 151, 153, 155, +157, 159, 161, 162, 164, 166, 168, 170, 172, 174, +176, 178, 179, 181, 183, 185, 187, 189, 191, 193, +194, 196, 198, 200, 202, 204, 206, 207, 209, 211, +213, 215, 217, 219, 220, 222, 224, 226, 228, 230, +232, 233, 235, 237, 239, 241, 243, 244, 246, 248, +250, 252, 254, 255, 257, 259, 261, 263, 265, 266, +268, 270, 272, 274, 275, 277, 279, 281, 283, 284, +286, 288, 290, 292, 294, 295, 297, 299, 301, 302, +304, 306, 308, 310, 311, 313, 315, 317, 319, 320, +322, 324, 326, 327, 329, 331, 333, 335, 336, 338, +340, 342, 343, 345, 347, 349, 350, 352, 354, 356, +357, 359, 361, 363, 364, 366, 368, 370, 371, 373, +375, 377, 378, 380, 382, 384, 385, 387, 389, 391, +392, 394, 396, 397, 399, 401, 403, 404, 406, 408, +410, 411, 413, 415, 416, 418, 420, 422, 423, 425, +427, 428, 430, 432, 433, 435, 437, 439, 440, 442, +444, 445, 447, 449, 450, 452, 454, 455, 457, 459, +460, 462, 464, 465, 467, 469, 471, 472, 474, 476, +477, 479, 481, 482, 484, 486, 487, 489, 490, 492, +494, 495, 497, 499, 500, 502, 504, 505, 507, 509, +510, 512, 514, 515, 517, 518, 520, 522, 523, 525, +527, 528, 530, 532, 533, 535, 536, 538, 540, 541, +543, 544, 546, 548, 549, 551, 553, 554, 556, 557, +559, 561, 562, 564, 565, 567, 569, 570, 572, 573, +575, 577, 578, 580, 581, 583, 585, 586, 588, 589, +591, 592, 594, 596, 597, 599, 600, 602, 603, 605, +607, 608, 610, 611, 613, 614, 616, 618, 619, 621, +622, 624, 625, 627, 628, 630, 632, 633, 635, 636, +638, 639, 641, 642, 644, 645, 647, 649, 650, 652, +653, 655, 656, 658, 659, 661, 662, 664, 665, 667, +668, 670, 671, 673, 674, 676, 678, 679, 681, 682, +684, 685, 687, 688, 690, 691, 693, 694, 696, 697, +699, 700, 702, 703, 705, 706, 708, 709, 711, 712, +714, 715, 717, 718, 719, 721, 722, 724, 725, 727, +728, 730, 731, 733, 734, 736, 737, 739, 740, 742, +743, 745, 746, 747, 749, 750, 752, 753, 755, 756, +758, 759, 761, 762, 763, 765, 766, 768, 769, 771, +772, 774, 775, 776, 778, 779, 781, 782, 784, 785, +786, 788, 789, 791, 792, 794, 795, 796, 798, 799, +801, 802, 804, 805, 806, 808, 809, 811, 812, 813, +815, 816, 818, 819, 820, 822, 823, 825, 826, 827, +829, 830, 832, 833, 834, 836, 837, 839, 840, 841, +843, 844, 846, 847, 848, 850, 851, 852, 854, 855, +857, 858, 859, 861, 862, 863, 865, 866, 868, 869, +870, 872, 873, 874, 876, 877, 878, 880, 881, 883, +884, 885, 887, 888, 889, 891, 892, 893, 895, 896, +897, 899, 900, 901, 903, 904, 905, 907, 908, 909, +911, 912, 913, 915, 916, 917, 919, 920, 921, 923, +924, 925, 927, 928, 929, 931, 932, 933, 935, 936, +937, 939, 940, 941, 942, 944, 945, 946, 948, 949, +950, 952, 953, 954, 956, 957, 958, 959, 961, 962, +963, 965, 966, 967, 968, 970, 971, 972, 974, 975, +976, 977, 979, 980, 981, 983, 984, 985, 986, 988, +989, 990, 991, 993, 994, 995, 997, 998, 999, 1000, +1002, 1003, 1004, 1005, 1007, 1008, 1009, 1010, 1012, 1013, +1014, 1015, 1017, 1018, 1019, 1020, 1022, 1023, 1024, 1025, +1027, 1028, 1029, 1030, 1032, 1033, 1034, 1035, 1037, 1038, +1039, 1040, 1041, 1043, 1044, 1045, 1046, 1048, 1049, 1050, +1051, 1052, 1054, 1055, 1056, 1057, 1059, 1060, 1061, 1062, +1063, 1065, 1066, 1067, 1068, 1069, 1071, 1072, 1073, 1074, +1076, 1077, 1078, 1079, 1080, 1082, 1083, 1084, 1085, 1086, +1087, 1089, 1090, 1091, 1092, 1093, 1095, 1096, 1097, 1098, +1099, 1101, 1102, 1103, 1104, 1105, 1106, 1108, 1109, 1110, +1111, 1112, 1114, 1115, 1116, 1117, 1118, 1119, 1121, 1122, +1123, 1124, 1125, 1126, 1128, 1129, 1130, 1131, 1132, 1133, +1135, 1136, 1137, 1138, 1139, 1140, 1141, 1143, 1144, 1145, +1146, 1147, 1148, 1149, 1151, 1152, 1153, 1154, 1155, 1156, +1157, 1159, 1160, 1161, 1162, 1163, 1164, 1165, 1167, 1168, +1169, 1170, 1171, 1172, 1173, 1174, 1176, 1177, 1178, 1179, +1180, 1181, 1182, 1183, 1185, 1186, 1187, 1188, 1189, 1190, +1191, 1192, 1193, 1195, 1196, 1197, 1198, 1199, 1200, 1201, +1202, 1203, 1205, 1206, 1207, 1208, 1209, 1210, 1211, 1212, +1213, 1214, 1216, 1217, 1218, 1219, 1220, 1221, 1222, 1223, +1224, 1225, 1226, 1228, 1229, 1230, 1231, 1232, 1233, 1234, +1235, 1236, 1237, 1238, 1239, 1240, 1242, 1243, 1244, 1245, +1246, 1247, 1248, 1249, 1250, 1251, 1252, 1253, 1254, 1255, +1256, 1258, 1259, 1260, 1261, 1262, 1263, 1264, 1265, 1266, +1267, 1268, 1269, 1270, 1271, 1272, 1273, 1274, 1275, 1277, +1278, 1279, 1280, 1281, 1282, 1283, 1284, 1285, 1286, 1287, +1288, 1289, 1290, 1291, 1292, 1293, 1294, 1295, 1296, 1297, +1298, 1299, 1300, 1301, 1302, 1303, 1305, 1306, 1307, 1308, +1309, 1310, 1311, 1312, 1313, 1314, 1315, 1316, 1317, 1318, +1319, 1320, 1321, 1322, 1323, 1324, 1325, 1326, 1327, 1328, +1329, 1330, 1331, 1332, 1333, 1334, 1335, 1336, 1337, 1338, +1339, 1340, 1341, 1342, 1343, 1344, 1345, 1346, 1347, 1348, +1349, 1350, 1351, 1352, 1353, 1354, 1355, 1356, 1357, 1358, +1359, 1360, 1361, 1362, 1363, 1364, 1365, 1366, 1367, 1368, +1369, 1370, 1371, 1372, 1372, 1373, 1374, 1375, 1376, 1377, +1378, 1379, 1380, 1381, 1382, 1383, 1384, 1385, 1386, 1387, +1388, 1389, 1390, 1391, 1392, 1393, 1394, 1395, 1396, 1397, +1398, 1399, 1399, 1400, 1401, 1402, 1403, 1404, 1405, 1406, +1407, 1408, 1409, 1410, 1411, 1412, 1413, 1414, 1415, 1416, +1417, 1418, 1418, 1419, 1420, 1421, 1422, 1423, 1424, 1425, +1426, 1427, 1428, 1429, 1430, 1431, 1432, 1432, 1433, 1434, +1435, 1436, 1437, 1438, 1439, 1440, 1441, 1442, 1443, 1444, +1444, 1445, 1446, 1447, 1448, 1449, 1450, 1451, 1452, 1453, +1454, 1455, 1455, 1456, 1457, 1458, 1459, 1460, 1461, 1462, +1463, 1464, 1465, 1466, 1466, 1467, 1468, 1469, 1470, 1471, +1472, 1473, 1474, 1475, 1475, 1476, 1477, 1478, 1479, 1480, +1481, 1482, 1483, 1483, 1484, 1485, 1486, 1487, 1488, 1489, +1490, 1491, 1491, 1492, 1493, 1494, 1495, 1496, 1497, 1498, +1499, 1499, 1500, 1501, 1502, 1503, 1504, 1505, 1506, 1506, +1507, 1508, 1509, 1510, 1511, 1512, 1513, 1513, 1514, 1515, +1516, 1517, 1518, 1519, 1520, 1520, 1521, 1522, 1523, 1524, +1525, 1526, 1526, 1527, 1528, 1529, 1530, 1531, 1532, 1532, +1533, 1534, 1535, 1536, 1537, 1538, 1538, 1539, 1540, 1541, +1542, 1543, 1544, 1544, 1545, 1546, 1547, 1548, 1549, 1550, +1550, 1551, 1552, 1553, 1554, 1555, 1555, 1556, 1557, 1558, +1559, 1560, 1560, 1561, 1562, 1563, 1564, 1565, 1565, 1566, +1567, 1568, 1569, 1570, 1570, 1571, 1572, 1573, 1574, 1575, +1575, 1576, 1577, 1578, 1579, 1580, 1580, 1581, 1582, 1583, +1584, 1584, 1585, 1586, 1587, 1588, 1589, 1589, 1590, 1591, +1592, 1593, 1593, 1594, 1595, 1596, 1597, 1598, 1598, 1599, +1600, 1601, 1602, 1602, 1603, 1604, 1605, 1606, 1606, 1607, +1608, 1609, 1610, 1610, 1611, 1612, 1613, 1614, 1614, 1615, +1616, 1617, 1618, 1618, 1619, 1620, 1621, 1622, 1622, 1623, +1624, 1625, 1626, 1626, 1627, 1628, 1629, 1630, 1630, 1631, +1632, 1633, 1634, 1634, 1635, 1636, 1637, 1637, 1638, 1639, +1640, 1641, 1641, 1642, 1643, 1644, 1645, 1645, 1646, 1647, +1648, 1648, 1649, 1650, 1651, 1652, 1652, 1653, 1654, 1655, +1655, 1656, 1657, 1658, 1658, 1659, 1660, 1661, 1662, 1662, +1663, 1664, 1665, 1665, 1666, 1667, 1668, 1668, 1669, 1670, +1671, 1671, 1672, 1673, 1674, 1675, 1675, 1676, 1677, 1678, +1678, 1679, 1680, 1681, 1681, 1682, 1683, 1684, 1684, 1685, +1686, 1687, 1687, 1688, 1689, 1690, 1690, 1691, 1692, 1693, +1693, 1694, 1695, 1696, 1696, 1697, 1698, 1699, 1699, 1700, +1701, 1702, 1702, 1703, 1704, 1705, 1705, 1706, 1707, 1707, +1708, 1709, 1710, 1710, 1711, 1712, 1713, 1713, 1714, 1715, +1716, 1716, 1717, 1718, 1718, 1719, 1720, 1721, 1721, 1722, +1723, 1724, 1724, 1725, 1726, 1727, 1727, 1728, 1729, 1729, +1730, 1731, 1732, 1732, 1733, 1734, 1734, 1735, 1736, 1737, +1737, 1738, 1739, 1740, 1740, 1741, 1742, 1742, 1743, 1744, +1745, 1745, 1746, 1747, 1747, 1748, 1749, 1749, 1750, 1751, +1752, 1752, 1753, 1754, 1754, 1755, 1756, 1757, 1757, 1758, +1759, 1759, 1760, 1761, 1762, 1762, 1763, 1764, 1764, 1765, +1766, 1766, 1767, 1768, 1769, 1769, 1770, 1771, 1771, 1772, +1773, 1773, 1774, 1775, 1776, 1776, 1777, 1778, 1778, 1779, +1780, 1780, 1781, 1782, 1782, 1783, 1784, 1784, 1785, 1786, +1787, 1787, 1788, 1789, 1789, 1790, 1791, 1791, 1792, 1793, +1793, 1794, 1795, 1795, 1796, 1797, 1798, 1798, 1799, 1800, +1800, 1801, 1802, 1802, 1803, 1804, 1804, 1805, 1806, 1806, +1807, 1808, 1808, 1809, 1810, 1810, 1811, 1812, 1812, 1813, +1814, 1814, 1815, 1816, 1816, 1817, 1818, 1818, 1819, 1820, +1820, 1821, 1822, 1822, 1823, 1824, 1824, 1825, 1826, 1826, +1827, 1828, 1828, 1829, 1830, 1830, 1831, 1832, 1832, 1833, +1834, 1834, 1835, 1836, 1836, 1837, 1838, 1838, 1839, 1840, +1840, 1841, 1842, 1842, 1843, 1844, 1844, 1845, 1845, 1846, +1847, 1847, 1848, 1849, 1849, 1850, 1851, 1851, 1852, 1853, +1853, 1854, 1855, 1855, 1856, 1857, 1857, 1858, 1858, 1859, +1860, 1860, 1861, 1862, 1862, 1863, 1864, 1864, 1865, 1866, +1866, 1867, 1867, 1868, 1869, 1869, 1870, 1871, 1871, 1872, +1873, 1873, 1874, 1874, 1875, 1876, 1876, 1877, 1878, 1878, +1879, 1879, 1880, 1881, 1881, 1882, 1883, 1883, 1884, 1885, +1885, 1886, 1886, 1887, 1888, 1888, 1889, 1890, 1890, 1891, +1891, 1892, 1893, 1893, 1894, 1895, 1895, 1896, 1896, 1897, +1898, 1898, 1899, 1900, 1900, 1901, 1901, 1902, 1903, 1903, +1904, 1904, 1905, 1906, 1906, 1907, 1908, 1908, 1909, 1909, +1910, 1911, 1911, 1912, 1912, 1913, 1914, 1914, 1915, 1916, +1916, 1917, 1917, 1918, 1919, 1919, 1920, 1920, 1921, 1922, +1922, 1923, 1923, 1924, 1925, 1925, 1926, 1926, 1927, 1928, +1928, 1929, 1929, 1930, 1931, 1931, 1932, 1932, 1933, 1934, +1934, 1935, 1935, 1936, 1937, 1937, 1938, 1938, 1939, 1940, +1940, 1941, 1941, 1942, 1943, 1943, 1944, 1944, 1945, 1946, +1946, 1947, 1947, 1948, 1949, 1949, 1950, 1950, 1951, 1952, +1952, 1953, 1953, 1954, 1955, 1955, 1956, 1956, 1957, 1957, +1958, 1959, 1959, 1960, 1960, 1961, 1962, 1962, 1963, 1963, +1964, 1964, 1965, 1966, 1966, 1967, 1967, 1968, 1969, 1969, +1970, 1970, 1971, 1971, 1972, 1973, 1973, 1974, 1974, 1975, +1976, 1976, 1977, 1977, 1978, 1978, 1979, 1980, 1980, 1981, +1981, 1982, 1982, 1983, 1984, 1984, 1985, 1985, 1986, 1986, +1987, 1988, 1988, 1989, 1989, 1990, 1990, 1991, 1992, 1992, +1993, 1993, 1994, 1994, 1995, 1996, 1996, 1997, 1997, 1998, +1998, 1999, 1999, 2000, 2001, 2001, 2002, 2002, 2003, 2003, +2004, 2005, 2005, 2006, 2006, 2007, 2007, 2008, 2008, 2009, +2010, 2010, 2011, 2011, 2012, 2012, 2013, 2014, 2014, 2015, +2015, 2016, 2016, 2017, 2017, 2018, 2019, 2019, 2020, 2020, +2021, 2021, 2022, 2022, 2023, 2023, 2024, 2025, 2025, 2026, +2026, 2027, 2027, 2028, 2028, 2029, 2030, 2030, 2031, 2031, +2032, 2032, 2033, 2033, 2034, 2034, 2035, 2036, 2036, 2037, +2037, 2038, 2038, 2039, 2039, 2040, 2040, 2041, 2042, 2042, +2043, 2043, 2044, 2044, 2045, 2045, 2046, 2046, 2047, 2048, +2048, 2049, 2049, 2050, 2050, 2051, 2051, 2052, 2052, 2053, +2053, 2054, 2054, 2055, 2056, 2056, 2057, 2057, 2058, 2058, +2059, 2059, 2060, 2060, 2061, 2061, 2062, 2062, 2063, 2064, +2064, 2065, 2065, 2066, 2066, 2067, 2067, 2068, 2068, 2069, +2069, 2070, 2070, 2071, 2072, 2072, 2073, 2073, 2074, 2074, +2075, 2075, 2076, 2076, 2077, 2077, 2078, 2078, 2079, 2079, +2080, 2080, 2081 +}; + +static const int fe_logsub_table_size = + sizeof(fe_logsub_table) / sizeof(fe_logsub_table[0]); + +fixed32 +fe_log_sub(fixed32 x, fixed32 y) +{ + fixed32 d, r; + + if (x < MIN_FIXLOG || y >= x) + return MIN_FIXLOG; + + d = (x - y) >> (DEFAULT_RADIX - 8); + + if (d > fe_logsub_table_size - 1) + return x; + + r = fe_logsub_table[d] << (DEFAULT_RADIX - 10); +/* + printf("diff=%d\n", + x + FIXLN(x-y) - r - + (x + FLOAT2FIX(logf(-expm1f(FIX2FLOAT(y - x)))))); +*/ + return x + FIXLN(x-y) - r; +} + +static fixed32 +fe_log(float32 x) +{ + if (x <= 0) { + return MIN_FIXLOG; + } + else { + return FLOAT2FIX(log(x)); + } +} +#endif + +static float32 +fe_mel(melfb_t * mel, float32 x) +{ + float32 warped = fe_warp_unwarped_to_warped(mel, x); + + return (float32) (2595.0 * log10(1.0 + warped / 700.0)); +} + +static float32 +fe_melinv(melfb_t * mel, float32 x) +{ + float32 warped = (float32) (700.0 * (pow(10.0, x / 2595.0) - 1.0)); + return fe_warp_warped_to_unwarped(mel, warped); +} + +int32 +fe_build_melfilters(melfb_t * mel_fb) +{ + float32 melmin, melmax, melbw, fftfreq; + int n_coeffs, i, j; + + + /* Filter coefficient matrix, in flattened form. */ + mel_fb->spec_start = + ckd_calloc(mel_fb->num_filters, sizeof(*mel_fb->spec_start)); + mel_fb->filt_start = + ckd_calloc(mel_fb->num_filters, sizeof(*mel_fb->filt_start)); + mel_fb->filt_width = + ckd_calloc(mel_fb->num_filters, sizeof(*mel_fb->filt_width)); + + /* First calculate the widths of each filter. */ + /* Minimum and maximum frequencies in mel scale. */ + melmin = fe_mel(mel_fb, mel_fb->lower_filt_freq); + melmax = fe_mel(mel_fb, mel_fb->upper_filt_freq); + + /* Width of filters in mel scale */ + melbw = (melmax - melmin) / (mel_fb->num_filters + 1); + if (mel_fb->doublewide) { + melmin -= melbw; + melmax += melbw; + if ((fe_melinv(mel_fb, melmin) < 0) || + (fe_melinv(mel_fb, melmax) > mel_fb->sampling_rate / 2)) { + E_WARN + ("Out of Range: low filter edge = %f (%f)\n", + fe_melinv(mel_fb, melmin), 0.0); + E_WARN + (" high filter edge = %f (%f)\n", + fe_melinv(mel_fb, melmax), mel_fb->sampling_rate / 2); + return FE_INVALID_PARAM_ERROR; + } + } + + /* DFT point spacing */ + fftfreq = mel_fb->sampling_rate / (float32) mel_fb->fft_size; + + /* Count and place filter coefficients. */ + n_coeffs = 0; + for (i = 0; i < mel_fb->num_filters; ++i) { + float32 freqs[3]; + + /* Left, center, right frequencies in Hertz */ + for (j = 0; j < 3; ++j) { + if (mel_fb->doublewide) + freqs[j] = fe_melinv(mel_fb, (i + j * 2) * melbw + melmin); + else + freqs[j] = fe_melinv(mel_fb, (i + j) * melbw + melmin); + /* Round them to DFT points if requested */ + if (mel_fb->round_filters) + freqs[j] = ((int) (freqs[j] / fftfreq + 0.5)) * fftfreq; + } + + /* spec_start is the start of this filter in the power spectrum. */ + mel_fb->spec_start[i] = -1; + /* There must be a better way... */ + for (j = 0; j < mel_fb->fft_size / 2 + 1; ++j) { + float32 hz = j * fftfreq; + if (hz < freqs[0]) + continue; + else if (hz > freqs[2] || j == mel_fb->fft_size / 2) { + /* filt_width is the width in DFT points of this filter. */ + mel_fb->filt_width[i] = j - mel_fb->spec_start[i]; + /* filt_start is the start of this filter in the filt_coeffs array. */ + mel_fb->filt_start[i] = n_coeffs; + n_coeffs += mel_fb->filt_width[i]; + break; + } + if (mel_fb->spec_start[i] == -1) + mel_fb->spec_start[i] = j; + } + } + + /* Now go back and allocate the coefficient array. */ + mel_fb->filt_coeffs = + ckd_malloc(n_coeffs * sizeof(*mel_fb->filt_coeffs)); + + /* And now generate the coefficients. */ + n_coeffs = 0; + for (i = 0; i < mel_fb->num_filters; ++i) { + float32 freqs[3]; + + /* Left, center, right frequencies in Hertz */ + for (j = 0; j < 3; ++j) { + if (mel_fb->doublewide) + freqs[j] = fe_melinv(mel_fb, (i + j * 2) * melbw + melmin); + else + freqs[j] = fe_melinv(mel_fb, (i + j) * melbw + melmin); + /* Round them to DFT points if requested */ + if (mel_fb->round_filters) + freqs[j] = ((int) (freqs[j] / fftfreq + 0.5)) * fftfreq; + } + + for (j = 0; j < mel_fb->filt_width[i]; ++j) { + float32 hz, loslope, hislope; + + hz = (mel_fb->spec_start[i] + j) * fftfreq; + if (hz < freqs[0] || hz > freqs[2]) { + E_FATAL + ("Failed to create filterbank, frequency range does not match. " + "Sample rate %f, FFT size %d, lowerf %f < freq %f > upperf %f.\n", + mel_fb->sampling_rate, mel_fb->fft_size, freqs[0], hz, + freqs[2]); + } + loslope = (hz - freqs[0]) / (freqs[1] - freqs[0]); + hislope = (freqs[2] - hz) / (freqs[2] - freqs[1]); + if (mel_fb->unit_area) { + loslope *= 2 / (freqs[2] - freqs[0]); + hislope *= 2 / (freqs[2] - freqs[0]); + } + if (loslope < hislope) { +#ifdef FIXED_POINT + mel_fb->filt_coeffs[n_coeffs] = fe_log(loslope); +#else + mel_fb->filt_coeffs[n_coeffs] = loslope; +#endif + } + else { +#ifdef FIXED_POINT + mel_fb->filt_coeffs[n_coeffs] = fe_log(hislope); +#else + mel_fb->filt_coeffs[n_coeffs] = hislope; +#endif + } + ++n_coeffs; + } + } + + return FE_SUCCESS; +} + +int32 +fe_compute_melcosine(melfb_t * mel_fb) +{ + + float64 freqstep; + int32 i, j; + + mel_fb->mel_cosine = + (mfcc_t **) ckd_calloc_2d(mel_fb->num_cepstra, + mel_fb->num_filters, sizeof(mfcc_t)); + + freqstep = M_PI / mel_fb->num_filters; + /* NOTE: The first row vector is actually unnecessary but we leave + * it in to avoid confusion. */ + for (i = 0; i < mel_fb->num_cepstra; i++) { + for (j = 0; j < mel_fb->num_filters; j++) { + float64 cosine; + + cosine = cos(freqstep * i * (j + 0.5)); + mel_fb->mel_cosine[i][j] = FLOAT2COS(cosine); + } + } + + /* Also precompute normalization constants for unitary DCT. */ + mel_fb->sqrt_inv_n = FLOAT2COS(sqrt(1.0 / mel_fb->num_filters)); + mel_fb->sqrt_inv_2n = FLOAT2COS(sqrt(2.0 / mel_fb->num_filters)); + + /* And liftering weights */ + if (mel_fb->lifter_val) { + mel_fb->lifter = + calloc(mel_fb->num_cepstra, sizeof(*mel_fb->lifter)); + for (i = 0; i < mel_fb->num_cepstra; ++i) { + mel_fb->lifter[i] = FLOAT2MFCC(1 + mel_fb->lifter_val / 2 + * sin(i * M_PI / + mel_fb->lifter_val)); + } + } + + return (0); +} + +static void +fe_pre_emphasis(int16 const *in, frame_t * out, int32 len, + float32 factor, int16 prior) +{ + int i; + +#if defined(FIXED16) + int16 fxd_alpha = (int16) (factor * 0x8000); + int32 tmp1, tmp2; + + tmp1 = (int32) in[0] << 15; + tmp2 = (int32) prior *fxd_alpha; + out[0] = (int16) ((tmp1 - tmp2) >> 15); + for (i = 1; i < len; ++i) { + tmp1 = (int32) in[i] << 15; + tmp2 = (int32) in[i - 1] * fxd_alpha; + out[i] = (int16) ((tmp1 - tmp2) >> 15); + } +#elif defined(FIXED_POINT) + fixed32 fxd_alpha = FLOAT2FIX(factor); + out[0] = ((fixed32) in[0] << DEFAULT_RADIX) - (prior * fxd_alpha); + for (i = 1; i < len; ++i) + out[i] = ((fixed32) in[i] << DEFAULT_RADIX) + - (fixed32) in[i - 1] * fxd_alpha; +#else + out[0] = (frame_t) in[0] - (frame_t) prior *factor; + for (i = 1; i < len; i++) + out[i] = (frame_t) in[i] - (frame_t) in[i - 1] * factor; +#endif +} + +static void +fe_short_to_frame(int16 const *in, frame_t * out, int32 len) +{ + int i; + +#if defined(FIXED16) + memcpy(out, in, len * sizeof(*out)); +#elif defined(FIXED_POINT) + for (i = 0; i < len; i++) + out[i] = (int32) in[i] << DEFAULT_RADIX; +#else /* FIXED_POINT */ + for (i = 0; i < len; i++) + out[i] = (frame_t) in[i]; +#endif /* FIXED_POINT */ +} + +void +fe_create_hamming(window_t * in, int32 in_len) +{ + int i; + + /* Symmetric, so we only create the first half of it. */ + for (i = 0; i < in_len / 2; i++) { + float64 hamm; + hamm = (0.54 - 0.46 * cos(2 * M_PI * i / + ((float64) in_len - 1.0))); +#ifdef FIXED16 + in[i] = (int16) (hamm * 0x8000); +#else + in[i] = FLOAT2COS(hamm); +#endif + } +} + +static void +fe_hamming_window(frame_t * in, window_t * window, int32 in_len, + int32 remove_dc) +{ + int i; + + if (remove_dc) { +#ifdef FIXED16 + int32 mean = 0; /* Use int32 to avoid possibility of overflow */ +#else + frame_t mean = 0; +#endif + + for (i = 0; i < in_len; i++) + mean += in[i]; + mean /= in_len; + for (i = 0; i < in_len; i++) + in[i] -= (frame_t) mean; + } + +#ifdef FIXED16 + for (i = 0; i < in_len / 2; i++) { + int32 tmp1, tmp2; + + tmp1 = (int32) in[i] * window[i]; + tmp2 = (int32) in[in_len - 1 - i] * window[i]; + in[i] = (int16) (tmp1 >> 15); + in[in_len - 1 - i] = (int16) (tmp2 >> 15); + } +#else + for (i = 0; i < in_len / 2; i++) { + in[i] = COSMUL(in[i], window[i]); + in[in_len - 1 - i] = COSMUL(in[in_len - 1 - i], window[i]); + } +#endif +} + +static int +fe_spch_to_frame(fe_t * fe, int len) +{ + /* Copy to the frame buffer. */ + if (fe->pre_emphasis_alpha != 0.0) { + fe_pre_emphasis(fe->spch, fe->frame, len, + fe->pre_emphasis_alpha, fe->prior); + if (len >= fe->frame_shift) + fe->prior = fe->spch[fe->frame_shift - 1]; + else + fe->prior = fe->spch[len - 1]; + } + else + fe_short_to_frame(fe->spch, fe->frame, len); + + /* Zero pad up to FFT size. */ + memset(fe->frame + len, 0, (fe->fft_size - len) * sizeof(*fe->frame)); + + /* Window. */ + fe_hamming_window(fe->frame, fe->hamming_window, fe->frame_size, + fe->remove_dc); + + return len; +} + +int +fe_read_frame(fe_t * fe, int16 const *in, int32 len) +{ + int i; + + if (len > fe->frame_size) + len = fe->frame_size; + + /* Read it into the raw speech buffer. */ + memcpy(fe->spch, in, len * sizeof(*in)); + /* Swap and dither if necessary. */ + if (fe->swap) + for (i = 0; i < len; ++i) + SWAP_INT16(&fe->spch[i]); + if (fe->dither) + for (i = 0; i < len; ++i) + fe->spch[i] += (int16) ((!(s3_rand_int31() % 4)) ? 1 : 0); + + return fe_spch_to_frame(fe, len); +} + +int +fe_shift_frame(fe_t * fe, int16 const *in, int32 len) +{ + int offset, i; + + if (len > fe->frame_shift) + len = fe->frame_shift; + offset = fe->frame_size - fe->frame_shift; + + /* Shift data into the raw speech buffer. */ + memmove(fe->spch, fe->spch + fe->frame_shift, + offset * sizeof(*fe->spch)); + memcpy(fe->spch + offset, in, len * sizeof(*fe->spch)); + /* Swap and dither if necessary. */ + if (fe->swap) + for (i = 0; i < len; ++i) + SWAP_INT16(&fe->spch[offset + i]); + if (fe->dither) + for (i = 0; i < len; ++i) + fe->spch[offset + i] + += (int16) ((!(s3_rand_int31() % 4)) ? 1 : 0); + + return fe_spch_to_frame(fe, offset + len); +} + +/** + * Create arrays of twiddle factors. + */ +void +fe_create_twiddle(fe_t * fe) +{ + int i; + + for (i = 0; i < fe->fft_size / 4; ++i) { + float64 a = 2 * M_PI * i / fe->fft_size; +#ifdef FIXED16 + fe->ccc[i] = (int16) (cos(a) * 0x8000); + fe->sss[i] = (int16) (sin(a) * 0x8000); +#elif defined(FIXED_POINT) + fe->ccc[i] = FLOAT2COS(cos(a)); + fe->sss[i] = FLOAT2COS(sin(a)); +#else + fe->ccc[i] = cos(a); + fe->sss[i] = sin(a); +#endif + } +} + + +/* Translated from the FORTRAN (obviously) from "Real-Valued Fast + * Fourier Transform Algorithms" by Henrik V. Sorensen et al., IEEE + * Transactions on Acoustics, Speech, and Signal Processing, vol. 35, + * no.6. The 16-bit version does a version of "block floating + * point" in order to avoid rounding errors. + */ +#if defined(FIXED16) +static int +fe_fft_real(fe_t * fe) +{ + int i, j, k, m, n, lz; + frame_t *x, xt, max; + + x = fe->frame; + m = fe->fft_order; + n = fe->fft_size; + + /* Bit-reverse the input. */ + j = 0; + for (i = 0; i < n - 1; ++i) { + if (i < j) { + xt = x[j]; + x[j] = x[i]; + x[i] = xt; + } + k = n / 2; + while (k <= j) { + j -= k; + k /= 2; + } + j += k; + } + /* Determine how many bits of dynamic range are in the input. */ + max = 0; + for (i = 0; i < n; ++i) + if (abs(x[i]) > max) + max = abs(x[i]); + /* The FFT has a gain of M bits, so we need to attenuate the input + * by M bits minus the number of leading zeroes in the input's + * range in order to avoid overflows. */ + for (lz = 0; lz < m; ++lz) + if (max & (1 << (15 - lz))) + break; + + /* Basic butterflies (2-point FFT, real twiddle factors): + * x[i] = x[i] + 1 * x[i+1] + * x[i+1] = x[i] + -1 * x[i+1] + */ + /* The quantization error introduced by attenuating the input at + * any given stage of the FFT has a cascading effect, so we hold + * off on it until it's absolutely necessary. */ + for (i = 0; i < n; i += 2) { + int atten = (lz == 0); + xt = x[i] >> atten; + x[i] = xt + (x[i + 1] >> atten); + x[i + 1] = xt - (x[i + 1] >> atten); + } + + /* The rest of the butterflies, in stages from 1..m */ + for (k = 1; k < m; ++k) { + int n1, n2, n4; + /* Start attenuating once we hit the number of leading zeros. */ + int atten = (k >= lz); + + n4 = k - 1; + n2 = k; + n1 = k + 1; + /* Stride over each (1 << (k+1)) points */ + for (i = 0; i < n; i += (1 << n1)) { + /* Basic butterfly with real twiddle factors: + * x[i] = x[i] + 1 * x[i + (1<<k)] + * x[i + (1<<k)] = x[i] + -1 * x[i + (1<<k)] + */ + xt = x[i] >> atten; + x[i] = xt + (x[i + (1 << n2)] >> atten); + x[i + (1 << n2)] = xt - (x[i + (1 << n2)] >> atten); + + /* The other ones with real twiddle factors: + * x[i + (1<<k) + (1<<(k-1))] + * = 0 * x[i + (1<<k-1)] + -1 * x[i + (1<<k) + (1<<k-1)] + * x[i + (1<<(k-1))] + * = 1 * x[i + (1<<k-1)] + 0 * x[i + (1<<k) + (1<<k-1)] + */ + x[i + (1 << n2) + (1 << n4)] = + -x[i + (1 << n2) + (1 << n4)] >> atten; + x[i + (1 << n4)] = x[i + (1 << n4)] >> atten; + + /* Butterflies with complex twiddle factors. + * There are (1<<k-1) of them. + */ + for (j = 1; j < (1 << n4); ++j) { + frame_t cc, ss, t1, t2; + int i1, i2, i3, i4; + + i1 = i + j; + i2 = i + (1 << n2) - j; + i3 = i + (1 << n2) + j; + i4 = i + (1 << n2) + (1 << n2) - j; + + /* + * cc = real(W[j * n / (1<<(k+1))]) + * ss = imag(W[j * n / (1<<(k+1))]) + */ + cc = fe->ccc[j << (m - n1)]; + ss = fe->sss[j << (m - n1)]; + + /* There are some symmetry properties which allow us + * to get away with only four multiplications here. */ + { + int32 tmp1, tmp2; + tmp1 = (int32) x[i3] * cc + (int32) x[i4] * ss; + tmp2 = (int32) x[i3] * ss - (int32) x[i4] * cc; + t1 = (int16) (tmp1 >> 15) >> atten; + t2 = (int16) (tmp2 >> 15) >> atten; + } + + x[i4] = (x[i2] >> atten) - t2; + x[i3] = (-x[i2] >> atten) - t2; + x[i2] = (x[i1] >> atten) - t1; + x[i1] = (x[i1] >> atten) + t1; + } + } + } + + /* Return the residual scaling factor. */ + return lz; +} +#else /* !FIXED16 */ +static int +fe_fft_real(fe_t * fe) +{ + int i, j, k, m, n; + frame_t *x, xt; + + x = fe->frame; + m = fe->fft_order; + n = fe->fft_size; + + /* Bit-reverse the input. */ + j = 0; + for (i = 0; i < n - 1; ++i) { + if (i < j) { + xt = x[j]; + x[j] = x[i]; + x[i] = xt; + } + k = n / 2; + while (k <= j) { + j -= k; + k /= 2; + } + j += k; + } + + /* Basic butterflies (2-point FFT, real twiddle factors): + * x[i] = x[i] + 1 * x[i+1] + * x[i+1] = x[i] + -1 * x[i+1] + */ + for (i = 0; i < n; i += 2) { + xt = x[i]; + x[i] = (xt + x[i + 1]); + x[i + 1] = (xt - x[i + 1]); + } + + /* The rest of the butterflies, in stages from 1..m */ + for (k = 1; k < m; ++k) { + int n1, n2, n4; + + n4 = k - 1; + n2 = k; + n1 = k + 1; + /* Stride over each (1 << (k+1)) points */ + for (i = 0; i < n; i += (1 << n1)) { + /* Basic butterfly with real twiddle factors: + * x[i] = x[i] + 1 * x[i + (1<<k)] + * x[i + (1<<k)] = x[i] + -1 * x[i + (1<<k)] + */ + xt = x[i]; + x[i] = (xt + x[i + (1 << n2)]); + x[i + (1 << n2)] = (xt - x[i + (1 << n2)]); + + /* The other ones with real twiddle factors: + * x[i + (1<<k) + (1<<(k-1))] + * = 0 * x[i + (1<<k-1)] + -1 * x[i + (1<<k) + (1<<k-1)] + * x[i + (1<<(k-1))] + * = 1 * x[i + (1<<k-1)] + 0 * x[i + (1<<k) + (1<<k-1)] + */ + x[i + (1 << n2) + (1 << n4)] = -x[i + (1 << n2) + (1 << n4)]; + x[i + (1 << n4)] = x[i + (1 << n4)]; + + /* Butterflies with complex twiddle factors. + * There are (1<<k-1) of them. + */ + for (j = 1; j < (1 << n4); ++j) { + frame_t cc, ss, t1, t2; + int i1, i2, i3, i4; + + i1 = i + j; + i2 = i + (1 << n2) - j; + i3 = i + (1 << n2) + j; + i4 = i + (1 << n2) + (1 << n2) - j; + + /* + * cc = real(W[j * n / (1<<(k+1))]) + * ss = imag(W[j * n / (1<<(k+1))]) + */ + cc = fe->ccc[j << (m - n1)]; + ss = fe->sss[j << (m - n1)]; + + /* There are some symmetry properties which allow us + * to get away with only four multiplications here. */ + t1 = COSMUL(x[i3], cc) + COSMUL(x[i4], ss); + t2 = COSMUL(x[i3], ss) - COSMUL(x[i4], cc); + + x[i4] = (x[i2] - t2); + x[i3] = (-x[i2] - t2); + x[i2] = (x[i1] - t1); + x[i1] = (x[i1] + t1); + } + } + } + + /* This isn't used, but return it for completeness. */ + return m; +} +#endif /* !FIXED16 */ + +static void +fe_spec_magnitude(fe_t * fe) +{ + frame_t *fft; + powspec_t *spec; + int32 j, scale, fftsize; + + /* Do FFT and get the scaling factor back (only actually used in + * fixed-point). Note the scaling factor is expressed in bits. */ + scale = fe_fft_real(fe); + + /* Convenience pointers to make things less awkward below. */ + fft = fe->frame; + spec = fe->spec; + fftsize = fe->fft_size; + + /* We need to scale things up the rest of the way to N. */ + scale = fe->fft_order - scale; + + /* The first point (DC coefficient) has no imaginary part */ + { +#ifdef FIXED16 + spec[0] = fixlog(abs(fft[0]) << scale) * 2; +#elif defined(FIXED_POINT) + spec[0] = FIXLN(abs(fft[0]) << scale) * 2; +#else + spec[0] = fft[0] * fft[0]; +#endif + } + + for (j = 1; j <= fftsize / 2; j++) { +#ifdef FIXED16 + int32 rr = fixlog(abs(fft[j]) << scale) * 2; + int32 ii = fixlog(abs(fft[fftsize - j]) << scale) * 2; + spec[j] = fe_log_add(rr, ii); +#elif defined(FIXED_POINT) + int32 rr = FIXLN(abs(fft[j]) << scale) * 2; + int32 ii = FIXLN(abs(fft[fftsize - j]) << scale) * 2; + spec[j] = fe_log_add(rr, ii); +#else + spec[j] = fft[j] * fft[j] + fft[fftsize - j] * fft[fftsize - j]; +#endif + } +} + +static void +fe_mel_spec(fe_t * fe) +{ + int whichfilt; + powspec_t *spec, *mfspec; + + /* Convenience poitners. */ + spec = fe->spec; + mfspec = fe->mfspec; + for (whichfilt = 0; whichfilt < fe->mel_fb->num_filters; whichfilt++) { + int spec_start, filt_start, i; + + spec_start = fe->mel_fb->spec_start[whichfilt]; + filt_start = fe->mel_fb->filt_start[whichfilt]; + +#ifdef FIXED_POINT + mfspec[whichfilt] = + spec[spec_start] + fe->mel_fb->filt_coeffs[filt_start]; + for (i = 1; i < fe->mel_fb->filt_width[whichfilt]; i++) { + mfspec[whichfilt] = fe_log_add(mfspec[whichfilt], + spec[spec_start + i] + + fe->mel_fb-> + filt_coeffs[filt_start + i]); + } +#else /* !FIXED_POINT */ + mfspec[whichfilt] = 0; + for (i = 0; i < fe->mel_fb->filt_width[whichfilt]; i++) + mfspec[whichfilt] += + spec[spec_start + i] * fe->mel_fb->filt_coeffs[filt_start + + i]; +#endif /* !FIXED_POINT */ + } + +} + +#define LOG_FLOOR 1e-4 + +static void +fe_mel_cep(fe_t * fe, mfcc_t * mfcep) +{ + int32 i; + powspec_t *mfspec; + + /* Convenience pointer. */ + mfspec = fe->mfspec; + + for (i = 0; i < fe->mel_fb->num_filters; ++i) { +#ifndef FIXED_POINT /* It's already in log domain for fixed point */ + mfspec[i] = log(mfspec[i] + LOG_FLOOR); +#endif /* !FIXED_POINT */ + } + + /* If we are doing LOG_SPEC, then do nothing. */ + if (fe->log_spec == RAW_LOG_SPEC) { + for (i = 0; i < fe->feature_dimension; i++) { + mfcep[i] = (mfcc_t) mfspec[i]; + } + } + /* For smoothed spectrum, do DCT-II followed by (its inverse) DCT-III */ + else if (fe->log_spec == SMOOTH_LOG_SPEC) { + /* FIXME: This is probably broken for fixed-point. */ + fe_dct2(fe, mfspec, mfcep, 0); + fe_dct3(fe, mfcep, mfspec); + for (i = 0; i < fe->feature_dimension; i++) { + mfcep[i] = (mfcc_t) mfspec[i]; + } + } + else if (fe->transform == DCT_II) + fe_dct2(fe, mfspec, mfcep, FALSE); + else if (fe->transform == DCT_HTK) + fe_dct2(fe, mfspec, mfcep, TRUE); + else + fe_spec2cep(fe, mfspec, mfcep); + + return; +} + +void +fe_spec2cep(fe_t * fe, const powspec_t * mflogspec, mfcc_t * mfcep) +{ + int32 i, j, beta; + + /* Compute C0 separately (its basis vector is 1) to avoid + * costly multiplications. */ + mfcep[0] = mflogspec[0] / 2; /* beta = 0.5 */ + for (j = 1; j < fe->mel_fb->num_filters; j++) + mfcep[0] += mflogspec[j]; /* beta = 1.0 */ + mfcep[0] /= (frame_t) fe->mel_fb->num_filters; + + for (i = 1; i < fe->num_cepstra; ++i) { + mfcep[i] = 0; + for (j = 0; j < fe->mel_fb->num_filters; j++) { + if (j == 0) + beta = 1; /* 0.5 */ + else + beta = 2; /* 1.0 */ + mfcep[i] += COSMUL(mflogspec[j], + fe->mel_fb->mel_cosine[i][j]) * beta; + } + /* Note that this actually normalizes by num_filters, like the + * original Sphinx front-end, due to the doubled 'beta' factor + * above. */ + mfcep[i] /= (frame_t) fe->mel_fb->num_filters * 2; + } +} + +void +fe_dct2(fe_t * fe, const powspec_t * mflogspec, mfcc_t * mfcep, int htk) +{ + int32 i, j; + + /* Compute C0 separately (its basis vector is 1) to avoid + * costly multiplications. */ + mfcep[0] = mflogspec[0]; + for (j = 1; j < fe->mel_fb->num_filters; j++) + mfcep[0] += mflogspec[j]; + if (htk) + mfcep[0] = COSMUL(mfcep[0], fe->mel_fb->sqrt_inv_2n); + else /* sqrt(1/N) = sqrt(2/N) * 1/sqrt(2) */ + mfcep[0] = COSMUL(mfcep[0], fe->mel_fb->sqrt_inv_n); + + for (i = 1; i < fe->num_cepstra; ++i) { + mfcep[i] = 0; + for (j = 0; j < fe->mel_fb->num_filters; j++) { + mfcep[i] += COSMUL(mflogspec[j], fe->mel_fb->mel_cosine[i][j]); + } + mfcep[i] = COSMUL(mfcep[i], fe->mel_fb->sqrt_inv_2n); + } +} + +void +fe_lifter(fe_t * fe, mfcc_t * mfcep) +{ + int32 i; + + if (fe->mel_fb->lifter_val == 0) + return; + + for (i = 0; i < fe->num_cepstra; ++i) { + mfcep[i] = MFCCMUL(mfcep[i], fe->mel_fb->lifter[i]); + } +} + +void +fe_dct3(fe_t * fe, const mfcc_t * mfcep, powspec_t * mflogspec) +{ + int32 i, j; + + for (i = 0; i < fe->mel_fb->num_filters; ++i) { + mflogspec[i] = COSMUL(mfcep[0], SQRT_HALF); + for (j = 1; j < fe->num_cepstra; j++) { + mflogspec[i] += COSMUL(mfcep[j], fe->mel_fb->mel_cosine[j][i]); + } + mflogspec[i] = COSMUL(mflogspec[i], fe->mel_fb->sqrt_inv_2n); + } +} + +void +fe_write_frame(fe_t * fe, mfcc_t * fea) +{ + int32 is_speech; + + fe_spec_magnitude(fe); + fe_mel_spec(fe); + fe_track_snr(fe, &is_speech); + fe_mel_cep(fe, fea); + fe_lifter(fe, fea); + fe_vad_hangover(fe, fea, is_speech); +} + + +void * +fe_create_2d(int32 d1, int32 d2, int32 elem_size) +{ + return (void *) ckd_calloc_2d(d1, d2, elem_size); +} + +void +fe_free_2d(void *arr) +{ + ckd_free_2d((void **) arr); +} |