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, 0 insertions, 1377 deletions
diff --git a/media/sphinxbase/src/libsphinxbase/fe/fe_sigproc.c b/media/sphinxbase/src/libsphinxbase/fe/fe_sigproc.c deleted file mode 100644 index 577640f62..000000000 --- a/media/sphinxbase/src/libsphinxbase/fe/fe_sigproc.c +++ /dev/null @@ -1,1377 +0,0 @@ -/* -*- 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); -} |