Home | History | Annotate | Download | only in ns
      1 /*
      2  *  Copyright (c) 2012 The WebRTC project authors. All Rights Reserved.
      3  *
      4  *  Use of this source code is governed by a BSD-style license
      5  *  that can be found in the LICENSE file in the root of the source
      6  *  tree. An additional intellectual property rights grant can be found
      7  *  in the file PATENTS.  All contributing project authors may
      8  *  be found in the AUTHORS file in the root of the source tree.
      9  */
     10 
     11 #include "webrtc/modules/audio_processing/ns/noise_suppression_x.h"
     12 
     13 #include <assert.h>
     14 #include <math.h>
     15 #include <stdlib.h>
     16 #include <string.h>
     17 
     18 #include "webrtc/common_audio/signal_processing/include/real_fft.h"
     19 #include "webrtc/modules/audio_processing/ns/nsx_core.h"
     20 #include "webrtc/system_wrappers/include/cpu_features_wrapper.h"
     21 
     22 #if (defined WEBRTC_DETECT_NEON || defined WEBRTC_HAS_NEON)
     23 /* Tables are defined in ARM assembly files. */
     24 extern const int16_t WebRtcNsx_kLogTable[9];
     25 extern const int16_t WebRtcNsx_kCounterDiv[201];
     26 extern const int16_t WebRtcNsx_kLogTableFrac[256];
     27 #else
     28 static const int16_t WebRtcNsx_kLogTable[9] = {
     29   0, 177, 355, 532, 710, 887, 1065, 1242, 1420
     30 };
     31 
     32 static const int16_t WebRtcNsx_kCounterDiv[201] = {
     33   32767, 16384, 10923, 8192, 6554, 5461, 4681, 4096, 3641, 3277, 2979, 2731,
     34   2521, 2341, 2185, 2048, 1928, 1820, 1725, 1638, 1560, 1489, 1425, 1365, 1311,
     35   1260, 1214, 1170, 1130, 1092, 1057, 1024, 993, 964, 936, 910, 886, 862, 840,
     36   819, 799, 780, 762, 745, 728, 712, 697, 683, 669, 655, 643, 630, 618, 607,
     37   596, 585, 575, 565, 555, 546, 537, 529, 520, 512, 504, 496, 489, 482, 475,
     38   468, 462, 455, 449, 443, 437, 431, 426, 420, 415, 410, 405, 400, 395, 390,
     39   386, 381, 377, 372, 368, 364, 360, 356, 352, 349, 345, 341, 338, 334, 331,
     40   328, 324, 321, 318, 315, 312, 309, 306, 303, 301, 298, 295, 293, 290, 287,
     41   285, 282, 280, 278, 275, 273, 271, 269, 266, 264, 262, 260, 258, 256, 254,
     42   252, 250, 248, 246, 245, 243, 241, 239, 237, 236, 234, 232, 231, 229, 228,
     43   226, 224, 223, 221, 220, 218, 217, 216, 214, 213, 211, 210, 209, 207, 206,
     44   205, 204, 202, 201, 200, 199, 197, 196, 195, 194, 193, 192, 191, 189, 188,
     45   187, 186, 185, 184, 183, 182, 181, 180, 179, 178, 177, 176, 175, 174, 173,
     46   172, 172, 171, 170, 169, 168, 167, 166, 165, 165, 164, 163
     47 };
     48 
     49 static const int16_t WebRtcNsx_kLogTableFrac[256] = {
     50   0,   1,   3,   4,   6,   7,   9,  10,  11,  13,  14,  16,  17,  18,  20,  21,
     51   22,  24,  25,  26,  28,  29,  30,  32,  33,  34,  36,  37,  38,  40,  41,  42,
     52   44,  45,  46,  47,  49,  50,  51,  52,  54,  55,  56,  57,  59,  60,  61,  62,
     53   63,  65,  66,  67,  68,  69,  71,  72,  73,  74,  75,  77,  78,  79,  80,  81,
     54   82,  84,  85,  86,  87,  88,  89,  90,  92,  93,  94,  95,  96,  97,  98,  99,
     55   100, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 116,
     56   117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131,
     57   132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146,
     58   147, 148, 149, 150, 151, 152, 153, 154, 155, 155, 156, 157, 158, 159, 160,
     59   161, 162, 163, 164, 165, 166, 167, 168, 169, 169, 170, 171, 172, 173, 174,
     60   175, 176, 177, 178, 178, 179, 180, 181, 182, 183, 184, 185, 185, 186, 187,
     61   188, 189, 190, 191, 192, 192, 193, 194, 195, 196, 197, 198, 198, 199, 200,
     62   201, 202, 203, 203, 204, 205, 206, 207, 208, 208, 209, 210, 211, 212, 212,
     63   213, 214, 215, 216, 216, 217, 218, 219, 220, 220, 221, 222, 223, 224, 224,
     64   225, 226, 227, 228, 228, 229, 230, 231, 231, 232, 233, 234, 234, 235, 236,
     65   237, 238, 238, 239, 240, 241, 241, 242, 243, 244, 244, 245, 246, 247, 247,
     66   248, 249, 249, 250, 251, 252, 252, 253, 254, 255, 255
     67 };
     68 #endif  // WEBRTC_DETECT_NEON || WEBRTC_HAS_NEON
     69 
     70 // Skip first frequency bins during estimation. (0 <= value < 64)
     71 static const size_t kStartBand = 5;
     72 
     73 // hybrib Hanning & flat window
     74 static const int16_t kBlocks80w128x[128] = {
     75   0,    536,   1072,   1606,   2139,   2669,   3196,   3720,   4240,   4756,   5266,
     76   5771,   6270,   6762,   7246,   7723,   8192,   8652,   9102,   9543,   9974,  10394,
     77   10803,  11200,  11585,  11958,  12318,  12665,  12998,  13318,  13623,  13913,  14189,
     78   14449,  14694,  14924,  15137,  15334,  15515,  15679,  15826,  15956,  16069,  16165,
     79   16244,  16305,  16349,  16375,  16384,  16384,  16384,  16384,  16384,  16384,  16384,
     80   16384,  16384,  16384,  16384,  16384,  16384,  16384,  16384,  16384,  16384,  16384,
     81   16384,  16384,  16384,  16384,  16384,  16384,  16384,  16384,  16384,  16384,  16384,
     82   16384,  16384,  16384,  16384,  16375,  16349,  16305,  16244,  16165,  16069,  15956,
     83   15826,  15679,  15515,  15334,  15137,  14924,  14694,  14449,  14189,  13913,  13623,
     84   13318,  12998,  12665,  12318,  11958,  11585,  11200,  10803,  10394,   9974,   9543,
     85   9102,   8652,   8192,   7723,   7246,   6762,   6270,   5771,   5266,   4756,   4240,
     86   3720,   3196,   2669,   2139,   1606,   1072,    536
     87 };
     88 
     89 // hybrib Hanning & flat window
     90 static const int16_t kBlocks160w256x[256] = {
     91   0,   268,   536,   804,  1072,  1339,  1606,  1872,
     92   2139,  2404,  2669,  2933,  3196,  3459,  3720,  3981,
     93   4240,  4499,  4756,  5012,  5266,  5520,  5771,  6021,
     94   6270,  6517,  6762,  7005,  7246,  7486,  7723,  7959,
     95   8192,  8423,  8652,  8878,  9102,  9324,  9543,  9760,
     96   9974, 10185, 10394, 10600, 10803, 11003, 11200, 11394,
     97   11585, 11773, 11958, 12140, 12318, 12493, 12665, 12833,
     98   12998, 13160, 13318, 13472, 13623, 13770, 13913, 14053,
     99   14189, 14321, 14449, 14574, 14694, 14811, 14924, 15032,
    100   15137, 15237, 15334, 15426, 15515, 15599, 15679, 15754,
    101   15826, 15893, 15956, 16015, 16069, 16119, 16165, 16207,
    102   16244, 16277, 16305, 16329, 16349, 16364, 16375, 16382,
    103   16384, 16384, 16384, 16384, 16384, 16384, 16384, 16384,
    104   16384, 16384, 16384, 16384, 16384, 16384, 16384, 16384,
    105   16384, 16384, 16384, 16384, 16384, 16384, 16384, 16384,
    106   16384, 16384, 16384, 16384, 16384, 16384, 16384, 16384,
    107   16384, 16384, 16384, 16384, 16384, 16384, 16384, 16384,
    108   16384, 16384, 16384, 16384, 16384, 16384, 16384, 16384,
    109   16384, 16384, 16384, 16384, 16384, 16384, 16384, 16384,
    110   16384, 16384, 16384, 16384, 16384, 16384, 16384, 16384,
    111   16384, 16382, 16375, 16364, 16349, 16329, 16305, 16277,
    112   16244, 16207, 16165, 16119, 16069, 16015, 15956, 15893,
    113   15826, 15754, 15679, 15599, 15515, 15426, 15334, 15237,
    114   15137, 15032, 14924, 14811, 14694, 14574, 14449, 14321,
    115   14189, 14053, 13913, 13770, 13623, 13472, 13318, 13160,
    116   12998, 12833, 12665, 12493, 12318, 12140, 11958, 11773,
    117   11585, 11394, 11200, 11003, 10803, 10600, 10394, 10185,
    118   9974,  9760,  9543,  9324,  9102,  8878,  8652,  8423,
    119   8192,  7959,  7723,  7486,  7246,  7005,  6762,  6517,
    120   6270,  6021,  5771,  5520,  5266,  5012,  4756,  4499,
    121   4240,  3981,  3720,  3459,  3196,  2933,  2669,  2404,
    122   2139,  1872,  1606,  1339,  1072,   804,   536,   268
    123 };
    124 
    125 // Gain factor1 table: Input value in Q8 and output value in Q13
    126 // original floating point code
    127 //  if (gain > blim) {
    128 //    factor1 = 1.0 + 1.3 * (gain - blim);
    129 //    if (gain * factor1 > 1.0) {
    130 //      factor1 = 1.0 / gain;
    131 //    }
    132 //  } else {
    133 //    factor1 = 1.0;
    134 //  }
    135 static const int16_t kFactor1Table[257] = {
    136   8192, 8192, 8192, 8192, 8192, 8192, 8192,
    137   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    138   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    139   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    140   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    141   8192, 8192, 8233, 8274, 8315, 8355, 8396, 8436, 8475, 8515, 8554, 8592, 8631, 8669,
    142   8707, 8745, 8783, 8820, 8857, 8894, 8931, 8967, 9003, 9039, 9075, 9111, 9146, 9181,
    143   9216, 9251, 9286, 9320, 9354, 9388, 9422, 9456, 9489, 9523, 9556, 9589, 9622, 9655,
    144   9687, 9719, 9752, 9784, 9816, 9848, 9879, 9911, 9942, 9973, 10004, 10035, 10066,
    145   10097, 10128, 10158, 10188, 10218, 10249, 10279, 10308, 10338, 10368, 10397, 10426,
    146   10456, 10485, 10514, 10543, 10572, 10600, 10629, 10657, 10686, 10714, 10742, 10770,
    147   10798, 10826, 10854, 10882, 10847, 10810, 10774, 10737, 10701, 10666, 10631, 10596,
    148   10562, 10527, 10494, 10460, 10427, 10394, 10362, 10329, 10297, 10266, 10235, 10203,
    149   10173, 10142, 10112, 10082, 10052, 10023, 9994, 9965, 9936, 9908, 9879, 9851, 9824,
    150   9796, 9769, 9742, 9715, 9689, 9662, 9636, 9610, 9584, 9559, 9534, 9508, 9484, 9459,
    151   9434, 9410, 9386, 9362, 9338, 9314, 9291, 9268, 9245, 9222, 9199, 9176, 9154, 9132,
    152   9110, 9088, 9066, 9044, 9023, 9002, 8980, 8959, 8939, 8918, 8897, 8877, 8857, 8836,
    153   8816, 8796, 8777, 8757, 8738, 8718, 8699, 8680, 8661, 8642, 8623, 8605, 8586, 8568,
    154   8550, 8532, 8514, 8496, 8478, 8460, 8443, 8425, 8408, 8391, 8373, 8356, 8339, 8323,
    155   8306, 8289, 8273, 8256, 8240, 8224, 8208, 8192
    156 };
    157 
    158 // For Factor2 tables
    159 // original floating point code
    160 // if (gain > blim) {
    161 //   factor2 = 1.0;
    162 // } else {
    163 //   factor2 = 1.0 - 0.3 * (blim - gain);
    164 //   if (gain <= inst->denoiseBound) {
    165 //     factor2 = 1.0 - 0.3 * (blim - inst->denoiseBound);
    166 //   }
    167 // }
    168 //
    169 // Gain factor table: Input value in Q8 and output value in Q13
    170 static const int16_t kFactor2Aggressiveness1[257] = {
    171   7577, 7577, 7577, 7577, 7577, 7577,
    172   7577, 7577, 7577, 7577, 7577, 7577, 7577, 7577, 7577, 7577, 7577, 7596, 7614, 7632,
    173   7650, 7667, 7683, 7699, 7715, 7731, 7746, 7761, 7775, 7790, 7804, 7818, 7832, 7845,
    174   7858, 7871, 7884, 7897, 7910, 7922, 7934, 7946, 7958, 7970, 7982, 7993, 8004, 8016,
    175   8027, 8038, 8049, 8060, 8070, 8081, 8091, 8102, 8112, 8122, 8132, 8143, 8152, 8162,
    176   8172, 8182, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    177   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    178   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    179   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    180   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    181   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    182   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    183   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    184   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    185   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    186   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    187   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    188   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    189   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192
    190 };
    191 
    192 // Gain factor table: Input value in Q8 and output value in Q13
    193 static const int16_t kFactor2Aggressiveness2[257] = {
    194   7270, 7270, 7270, 7270, 7270, 7306,
    195   7339, 7369, 7397, 7424, 7448, 7472, 7495, 7517, 7537, 7558, 7577, 7596, 7614, 7632,
    196   7650, 7667, 7683, 7699, 7715, 7731, 7746, 7761, 7775, 7790, 7804, 7818, 7832, 7845,
    197   7858, 7871, 7884, 7897, 7910, 7922, 7934, 7946, 7958, 7970, 7982, 7993, 8004, 8016,
    198   8027, 8038, 8049, 8060, 8070, 8081, 8091, 8102, 8112, 8122, 8132, 8143, 8152, 8162,
    199   8172, 8182, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    200   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    201   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    202   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    203   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    204   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    205   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    206   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    207   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    208   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    209   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    210   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    211   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    212   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192
    213 };
    214 
    215 // Gain factor table: Input value in Q8 and output value in Q13
    216 static const int16_t kFactor2Aggressiveness3[257] = {
    217   7184, 7184, 7184, 7229, 7270, 7306,
    218   7339, 7369, 7397, 7424, 7448, 7472, 7495, 7517, 7537, 7558, 7577, 7596, 7614, 7632,
    219   7650, 7667, 7683, 7699, 7715, 7731, 7746, 7761, 7775, 7790, 7804, 7818, 7832, 7845,
    220   7858, 7871, 7884, 7897, 7910, 7922, 7934, 7946, 7958, 7970, 7982, 7993, 8004, 8016,
    221   8027, 8038, 8049, 8060, 8070, 8081, 8091, 8102, 8112, 8122, 8132, 8143, 8152, 8162,
    222   8172, 8182, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    223   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    224   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    225   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    226   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    227   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    228   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    229   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    230   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    231   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    232   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    233   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    234   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
    235   8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192
    236 };
    237 
    238 // sum of log2(i) from table index to inst->anaLen2 in Q5
    239 // Note that the first table value is invalid, since log2(0) = -infinity
    240 static const int16_t kSumLogIndex[66] = {
    241   0,  22917,  22917,  22885,  22834,  22770,  22696,  22613,
    242   22524,  22428,  22326,  22220,  22109,  21994,  21876,  21754,
    243   21629,  21501,  21370,  21237,  21101,  20963,  20822,  20679,
    244   20535,  20388,  20239,  20089,  19937,  19783,  19628,  19470,
    245   19312,  19152,  18991,  18828,  18664,  18498,  18331,  18164,
    246   17994,  17824,  17653,  17480,  17306,  17132,  16956,  16779,
    247   16602,  16423,  16243,  16063,  15881,  15699,  15515,  15331,
    248   15146,  14960,  14774,  14586,  14398,  14209,  14019,  13829,
    249   13637,  13445
    250 };
    251 
    252 // sum of log2(i)^2 from table index to inst->anaLen2 in Q2
    253 // Note that the first table value is invalid, since log2(0) = -infinity
    254 static const int16_t kSumSquareLogIndex[66] = {
    255   0,  16959,  16959,  16955,  16945,  16929,  16908,  16881,
    256   16850,  16814,  16773,  16729,  16681,  16630,  16575,  16517,
    257   16456,  16392,  16325,  16256,  16184,  16109,  16032,  15952,
    258   15870,  15786,  15700,  15612,  15521,  15429,  15334,  15238,
    259   15140,  15040,  14938,  14834,  14729,  14622,  14514,  14404,
    260   14292,  14179,  14064,  13947,  13830,  13710,  13590,  13468,
    261   13344,  13220,  13094,  12966,  12837,  12707,  12576,  12444,
    262   12310,  12175,  12039,  11902,  11763,  11624,  11483,  11341,
    263   11198,  11054
    264 };
    265 
    266 // log2(table index) in Q12
    267 // Note that the first table value is invalid, since log2(0) = -infinity
    268 static const int16_t kLogIndex[129] = {
    269   0,      0,   4096,   6492,   8192,   9511,  10588,  11499,
    270   12288,  12984,  13607,  14170,  14684,  15157,  15595,  16003,
    271   16384,  16742,  17080,  17400,  17703,  17991,  18266,  18529,
    272   18780,  19021,  19253,  19476,  19691,  19898,  20099,  20292,
    273   20480,  20662,  20838,  21010,  21176,  21338,  21496,  21649,
    274   21799,  21945,  22087,  22226,  22362,  22495,  22625,  22752,
    275   22876,  22998,  23117,  23234,  23349,  23462,  23572,  23680,
    276   23787,  23892,  23994,  24095,  24195,  24292,  24388,  24483,
    277   24576,  24668,  24758,  24847,  24934,  25021,  25106,  25189,
    278   25272,  25354,  25434,  25513,  25592,  25669,  25745,  25820,
    279   25895,  25968,  26041,  26112,  26183,  26253,  26322,  26390,
    280   26458,  26525,  26591,  26656,  26721,  26784,  26848,  26910,
    281   26972,  27033,  27094,  27154,  27213,  27272,  27330,  27388,
    282   27445,  27502,  27558,  27613,  27668,  27722,  27776,  27830,
    283   27883,  27935,  27988,  28039,  28090,  28141,  28191,  28241,
    284   28291,  28340,  28388,  28437,  28484,  28532,  28579,  28626,
    285   28672
    286 };
    287 
    288 // determinant of estimation matrix in Q0 corresponding to the log2 tables above
    289 // Note that the first table value is invalid, since log2(0) = -infinity
    290 static const int16_t kDeterminantEstMatrix[66] = {
    291   0,  29814,  25574,  22640,  20351,  18469,  16873,  15491,
    292   14277,  13199,  12233,  11362,  10571,   9851,   9192,   8587,
    293   8030,   7515,   7038,   6596,   6186,   5804,   5448,   5115,
    294   4805,   4514,   4242,   3988,   3749,   3524,   3314,   3116,
    295   2930,   2755,   2590,   2435,   2289,   2152,   2022,   1900,
    296   1785,   1677,   1575,   1478,   1388,   1302,   1221,   1145,
    297   1073,   1005,    942,    881,    825,    771,    721,    674,
    298   629,    587,    547,    510,    475,    442,    411,    382,
    299   355,    330
    300 };
    301 
    302 // Update the noise estimation information.
    303 static void UpdateNoiseEstimate(NoiseSuppressionFixedC* inst, int offset) {
    304   int32_t tmp32no1 = 0;
    305   int32_t tmp32no2 = 0;
    306   int16_t tmp16 = 0;
    307   const int16_t kExp2Const = 11819; // Q13
    308 
    309   size_t i = 0;
    310 
    311   tmp16 = WebRtcSpl_MaxValueW16(inst->noiseEstLogQuantile + offset,
    312                                    inst->magnLen);
    313   // Guarantee a Q-domain as high as possible and still fit in int16
    314   inst->qNoise = 14 - (int) WEBRTC_SPL_MUL_16_16_RSFT_WITH_ROUND(
    315                    kExp2Const, tmp16, 21);
    316   for (i = 0; i < inst->magnLen; i++) {
    317     // inst->quantile[i]=exp(inst->lquantile[offset+i]);
    318     // in Q21
    319     tmp32no2 = kExp2Const * inst->noiseEstLogQuantile[offset + i];
    320     tmp32no1 = (0x00200000 | (tmp32no2 & 0x001FFFFF)); // 2^21 + frac
    321     tmp16 = (int16_t)(tmp32no2 >> 21);
    322     tmp16 -= 21;// shift 21 to get result in Q0
    323     tmp16 += (int16_t) inst->qNoise; //shift to get result in Q(qNoise)
    324     if (tmp16 < 0) {
    325       tmp32no1 >>= -tmp16;
    326     } else {
    327       tmp32no1 <<= tmp16;
    328     }
    329     inst->noiseEstQuantile[i] = WebRtcSpl_SatW32ToW16(tmp32no1);
    330   }
    331 }
    332 
    333 // Noise Estimation
    334 static void NoiseEstimationC(NoiseSuppressionFixedC* inst,
    335                              uint16_t* magn,
    336                              uint32_t* noise,
    337                              int16_t* q_noise) {
    338   int16_t lmagn[HALF_ANAL_BLOCKL], counter, countDiv;
    339   int16_t countProd, delta, zeros, frac;
    340   int16_t log2, tabind, logval, tmp16, tmp16no1, tmp16no2;
    341   const int16_t log2_const = 22713; // Q15
    342   const int16_t width_factor = 21845;
    343 
    344   size_t i, s, offset;
    345 
    346   tabind = inst->stages - inst->normData;
    347   assert(tabind < 9);
    348   assert(tabind > -9);
    349   if (tabind < 0) {
    350     logval = -WebRtcNsx_kLogTable[-tabind];
    351   } else {
    352     logval = WebRtcNsx_kLogTable[tabind];
    353   }
    354 
    355   // lmagn(i)=log(magn(i))=log(2)*log2(magn(i))
    356   // magn is in Q(-stages), and the real lmagn values are:
    357   // real_lmagn(i)=log(magn(i)*2^stages)=log(magn(i))+log(2^stages)
    358   // lmagn in Q8
    359   for (i = 0; i < inst->magnLen; i++) {
    360     if (magn[i]) {
    361       zeros = WebRtcSpl_NormU32((uint32_t)magn[i]);
    362       frac = (int16_t)((((uint32_t)magn[i] << zeros)
    363                               & 0x7FFFFFFF) >> 23);
    364       // log2(magn(i))
    365       assert(frac < 256);
    366       log2 = (int16_t)(((31 - zeros) << 8)
    367                              + WebRtcNsx_kLogTableFrac[frac]);
    368       // log2(magn(i))*log(2)
    369       lmagn[i] = (int16_t)((log2 * log2_const) >> 15);
    370       // + log(2^stages)
    371       lmagn[i] += logval;
    372     } else {
    373       lmagn[i] = logval;//0;
    374     }
    375   }
    376 
    377   // loop over simultaneous estimates
    378   for (s = 0; s < SIMULT; s++) {
    379     offset = s * inst->magnLen;
    380 
    381     // Get counter values from state
    382     counter = inst->noiseEstCounter[s];
    383     assert(counter < 201);
    384     countDiv = WebRtcNsx_kCounterDiv[counter];
    385     countProd = (int16_t)(counter * countDiv);
    386 
    387     // quant_est(...)
    388     for (i = 0; i < inst->magnLen; i++) {
    389       // compute delta
    390       if (inst->noiseEstDensity[offset + i] > 512) {
    391         // Get the value for delta by shifting intead of dividing.
    392         int factor = WebRtcSpl_NormW16(inst->noiseEstDensity[offset + i]);
    393         delta = (int16_t)(FACTOR_Q16 >> (14 - factor));
    394       } else {
    395         delta = FACTOR_Q7;
    396         if (inst->blockIndex < END_STARTUP_LONG) {
    397           // Smaller step size during startup. This prevents from using
    398           // unrealistic values causing overflow.
    399           delta = FACTOR_Q7_STARTUP;
    400         }
    401       }
    402 
    403       // update log quantile estimate
    404       tmp16 = (int16_t)((delta * countDiv) >> 14);
    405       if (lmagn[i] > inst->noiseEstLogQuantile[offset + i]) {
    406         // +=QUANTILE*delta/(inst->counter[s]+1) QUANTILE=0.25, =1 in Q2
    407         // CounterDiv=1/(inst->counter[s]+1) in Q15
    408         tmp16 += 2;
    409         inst->noiseEstLogQuantile[offset + i] += tmp16 / 4;
    410       } else {
    411         tmp16 += 1;
    412         // *(1-QUANTILE), in Q2 QUANTILE=0.25, 1-0.25=0.75=3 in Q2
    413         // TODO(bjornv): investigate why we need to truncate twice.
    414         tmp16no2 = (int16_t)((tmp16 / 2) * 3 / 2);
    415         inst->noiseEstLogQuantile[offset + i] -= tmp16no2;
    416         if (inst->noiseEstLogQuantile[offset + i] < logval) {
    417           // This is the smallest fixed point representation we can
    418           // have, hence we limit the output.
    419           inst->noiseEstLogQuantile[offset + i] = logval;
    420         }
    421       }
    422 
    423       // update density estimate
    424       if (WEBRTC_SPL_ABS_W16(lmagn[i] - inst->noiseEstLogQuantile[offset + i])
    425           < WIDTH_Q8) {
    426         tmp16no1 = (int16_t)WEBRTC_SPL_MUL_16_16_RSFT_WITH_ROUND(
    427                      inst->noiseEstDensity[offset + i], countProd, 15);
    428         tmp16no2 = (int16_t)WEBRTC_SPL_MUL_16_16_RSFT_WITH_ROUND(
    429                      width_factor, countDiv, 15);
    430         inst->noiseEstDensity[offset + i] = tmp16no1 + tmp16no2;
    431       }
    432     }  // end loop over magnitude spectrum
    433 
    434     if (counter >= END_STARTUP_LONG) {
    435       inst->noiseEstCounter[s] = 0;
    436       if (inst->blockIndex >= END_STARTUP_LONG) {
    437         UpdateNoiseEstimate(inst, offset);
    438       }
    439     }
    440     inst->noiseEstCounter[s]++;
    441 
    442   }  // end loop over simultaneous estimates
    443 
    444   // Sequentially update the noise during startup
    445   if (inst->blockIndex < END_STARTUP_LONG) {
    446     UpdateNoiseEstimate(inst, offset);
    447   }
    448 
    449   for (i = 0; i < inst->magnLen; i++) {
    450     noise[i] = (uint32_t)(inst->noiseEstQuantile[i]); // Q(qNoise)
    451   }
    452   (*q_noise) = (int16_t)inst->qNoise;
    453 }
    454 
    455 // Filter the data in the frequency domain, and create spectrum.
    456 static void PrepareSpectrumC(NoiseSuppressionFixedC* inst, int16_t* freq_buf) {
    457   size_t i = 0, j = 0;
    458 
    459   for (i = 0; i < inst->magnLen; i++) {
    460     inst->real[i] = (int16_t)((inst->real[i] *
    461         (int16_t)(inst->noiseSupFilter[i])) >> 14);  // Q(normData-stages)
    462     inst->imag[i] = (int16_t)((inst->imag[i] *
    463         (int16_t)(inst->noiseSupFilter[i])) >> 14);  // Q(normData-stages)
    464   }
    465 
    466   freq_buf[0] = inst->real[0];
    467   freq_buf[1] = -inst->imag[0];
    468   for (i = 1, j = 2; i < inst->anaLen2; i += 1, j += 2) {
    469     freq_buf[j] = inst->real[i];
    470     freq_buf[j + 1] = -inst->imag[i];
    471   }
    472   freq_buf[inst->anaLen] = inst->real[inst->anaLen2];
    473   freq_buf[inst->anaLen + 1] = -inst->imag[inst->anaLen2];
    474 }
    475 
    476 // Denormalize the real-valued signal |in|, the output from inverse FFT.
    477 static void DenormalizeC(NoiseSuppressionFixedC* inst,
    478                          int16_t* in,
    479                          int factor) {
    480   size_t i = 0;
    481   int32_t tmp32 = 0;
    482   for (i = 0; i < inst->anaLen; i += 1) {
    483     tmp32 = WEBRTC_SPL_SHIFT_W32((int32_t)in[i],
    484                                  factor - inst->normData);
    485     inst->real[i] = WebRtcSpl_SatW32ToW16(tmp32); // Q0
    486   }
    487 }
    488 
    489 // For the noise supression process, synthesis, read out fully processed
    490 // segment, and update synthesis buffer.
    491 static void SynthesisUpdateC(NoiseSuppressionFixedC* inst,
    492                              int16_t* out_frame,
    493                              int16_t gain_factor) {
    494   size_t i = 0;
    495   int16_t tmp16a = 0;
    496   int16_t tmp16b = 0;
    497   int32_t tmp32 = 0;
    498 
    499   // synthesis
    500   for (i = 0; i < inst->anaLen; i++) {
    501     tmp16a = (int16_t)WEBRTC_SPL_MUL_16_16_RSFT_WITH_ROUND(
    502                  inst->window[i], inst->real[i], 14); // Q0, window in Q14
    503     tmp32 = WEBRTC_SPL_MUL_16_16_RSFT_WITH_ROUND(tmp16a, gain_factor, 13); // Q0
    504     // Down shift with rounding
    505     tmp16b = WebRtcSpl_SatW32ToW16(tmp32); // Q0
    506     inst->synthesisBuffer[i] = WebRtcSpl_AddSatW16(inst->synthesisBuffer[i],
    507                                                    tmp16b); // Q0
    508   }
    509 
    510   // read out fully processed segment
    511   for (i = 0; i < inst->blockLen10ms; i++) {
    512     out_frame[i] = inst->synthesisBuffer[i]; // Q0
    513   }
    514 
    515   // update synthesis buffer
    516   memcpy(inst->synthesisBuffer, inst->synthesisBuffer + inst->blockLen10ms,
    517       (inst->anaLen - inst->blockLen10ms) * sizeof(*inst->synthesisBuffer));
    518   WebRtcSpl_ZerosArrayW16(inst->synthesisBuffer
    519       + inst->anaLen - inst->blockLen10ms, inst->blockLen10ms);
    520 }
    521 
    522 // Update analysis buffer for lower band, and window data before FFT.
    523 static void AnalysisUpdateC(NoiseSuppressionFixedC* inst,
    524                             int16_t* out,
    525                             int16_t* new_speech) {
    526   size_t i = 0;
    527 
    528   // For lower band update analysis buffer.
    529   memcpy(inst->analysisBuffer, inst->analysisBuffer + inst->blockLen10ms,
    530       (inst->anaLen - inst->blockLen10ms) * sizeof(*inst->analysisBuffer));
    531   memcpy(inst->analysisBuffer + inst->anaLen - inst->blockLen10ms, new_speech,
    532       inst->blockLen10ms * sizeof(*inst->analysisBuffer));
    533 
    534   // Window data before FFT.
    535   for (i = 0; i < inst->anaLen; i++) {
    536     out[i] = (int16_t)WEBRTC_SPL_MUL_16_16_RSFT_WITH_ROUND(
    537                inst->window[i], inst->analysisBuffer[i], 14); // Q0
    538   }
    539 }
    540 
    541 // Normalize the real-valued signal |in|, the input to forward FFT.
    542 static void NormalizeRealBufferC(NoiseSuppressionFixedC* inst,
    543                                  const int16_t* in,
    544                                  int16_t* out) {
    545   size_t i = 0;
    546   assert(inst->normData >= 0);
    547   for (i = 0; i < inst->anaLen; ++i) {
    548     out[i] = in[i] << inst->normData;  // Q(normData)
    549   }
    550 }
    551 
    552 // Declare function pointers.
    553 NoiseEstimation WebRtcNsx_NoiseEstimation;
    554 PrepareSpectrum WebRtcNsx_PrepareSpectrum;
    555 SynthesisUpdate WebRtcNsx_SynthesisUpdate;
    556 AnalysisUpdate WebRtcNsx_AnalysisUpdate;
    557 Denormalize WebRtcNsx_Denormalize;
    558 NormalizeRealBuffer WebRtcNsx_NormalizeRealBuffer;
    559 
    560 #if (defined WEBRTC_DETECT_NEON || defined WEBRTC_HAS_NEON)
    561 // Initialize function pointers for ARM Neon platform.
    562 static void WebRtcNsx_InitNeon(void) {
    563   WebRtcNsx_NoiseEstimation = WebRtcNsx_NoiseEstimationNeon;
    564   WebRtcNsx_PrepareSpectrum = WebRtcNsx_PrepareSpectrumNeon;
    565   WebRtcNsx_SynthesisUpdate = WebRtcNsx_SynthesisUpdateNeon;
    566   WebRtcNsx_AnalysisUpdate = WebRtcNsx_AnalysisUpdateNeon;
    567 }
    568 #endif
    569 
    570 #if defined(MIPS32_LE)
    571 // Initialize function pointers for MIPS platform.
    572 static void WebRtcNsx_InitMips(void) {
    573   WebRtcNsx_PrepareSpectrum = WebRtcNsx_PrepareSpectrum_mips;
    574   WebRtcNsx_SynthesisUpdate = WebRtcNsx_SynthesisUpdate_mips;
    575   WebRtcNsx_AnalysisUpdate = WebRtcNsx_AnalysisUpdate_mips;
    576   WebRtcNsx_NormalizeRealBuffer = WebRtcNsx_NormalizeRealBuffer_mips;
    577 #if defined(MIPS_DSP_R1_LE)
    578   WebRtcNsx_Denormalize = WebRtcNsx_Denormalize_mips;
    579 #endif
    580 }
    581 #endif
    582 
    583 void WebRtcNsx_CalcParametricNoiseEstimate(NoiseSuppressionFixedC* inst,
    584                                            int16_t pink_noise_exp_avg,
    585                                            int32_t pink_noise_num_avg,
    586                                            int freq_index,
    587                                            uint32_t* noise_estimate,
    588                                            uint32_t* noise_estimate_avg) {
    589   int32_t tmp32no1 = 0;
    590   int32_t tmp32no2 = 0;
    591 
    592   int16_t int_part = 0;
    593   int16_t frac_part = 0;
    594 
    595   // Use pink noise estimate
    596   // noise_estimate = 2^(pinkNoiseNumerator + pinkNoiseExp * log2(j))
    597   assert(freq_index >= 0);
    598   assert(freq_index < 129);
    599   tmp32no2 = (pink_noise_exp_avg * kLogIndex[freq_index]) >> 15;  // Q11
    600   tmp32no1 = pink_noise_num_avg - tmp32no2; // Q11
    601 
    602   // Calculate output: 2^tmp32no1
    603   // Output in Q(minNorm-stages)
    604   tmp32no1 += (inst->minNorm - inst->stages) << 11;
    605   if (tmp32no1 > 0) {
    606     int_part = (int16_t)(tmp32no1 >> 11);
    607     frac_part = (int16_t)(tmp32no1 & 0x000007ff); // Q11
    608     // Piecewise linear approximation of 'b' in
    609     // 2^(int_part+frac_part) = 2^int_part * (1 + b)
    610     // 'b' is given in Q11 and below stored in frac_part.
    611     if (frac_part >> 10) {
    612       // Upper fractional part
    613       tmp32no2 = (2048 - frac_part) * 1244;  // Q21
    614       tmp32no2 = 2048 - (tmp32no2 >> 10);
    615     } else {
    616       // Lower fractional part
    617       tmp32no2 = (frac_part * 804) >> 10;
    618     }
    619     // Shift fractional part to Q(minNorm-stages)
    620     tmp32no2 = WEBRTC_SPL_SHIFT_W32(tmp32no2, int_part - 11);
    621     *noise_estimate_avg = (1 << int_part) + (uint32_t)tmp32no2;
    622     // Scale up to initMagnEst, which is not block averaged
    623     *noise_estimate = (*noise_estimate_avg) * (uint32_t)(inst->blockIndex + 1);
    624   }
    625 }
    626 
    627 // Initialize state
    628 int32_t WebRtcNsx_InitCore(NoiseSuppressionFixedC* inst, uint32_t fs) {
    629   int i;
    630 
    631   //check for valid pointer
    632   if (inst == NULL) {
    633     return -1;
    634   }
    635   //
    636 
    637   // Initialization of struct
    638   if (fs == 8000 || fs == 16000 || fs == 32000 || fs == 48000) {
    639     inst->fs = fs;
    640   } else {
    641     return -1;
    642   }
    643 
    644   if (fs == 8000) {
    645     inst->blockLen10ms = 80;
    646     inst->anaLen = 128;
    647     inst->stages = 7;
    648     inst->window = kBlocks80w128x;
    649     inst->thresholdLogLrt = 131072; //default threshold for LRT feature
    650     inst->maxLrt = 0x0040000;
    651     inst->minLrt = 52429;
    652   } else {
    653     inst->blockLen10ms = 160;
    654     inst->anaLen = 256;
    655     inst->stages = 8;
    656     inst->window = kBlocks160w256x;
    657     inst->thresholdLogLrt = 212644; //default threshold for LRT feature
    658     inst->maxLrt = 0x0080000;
    659     inst->minLrt = 104858;
    660   }
    661   inst->anaLen2 = inst->anaLen / 2;
    662   inst->magnLen = inst->anaLen2 + 1;
    663 
    664   if (inst->real_fft != NULL) {
    665     WebRtcSpl_FreeRealFFT(inst->real_fft);
    666   }
    667   inst->real_fft = WebRtcSpl_CreateRealFFT(inst->stages);
    668   if (inst->real_fft == NULL) {
    669     return -1;
    670   }
    671 
    672   WebRtcSpl_ZerosArrayW16(inst->analysisBuffer, ANAL_BLOCKL_MAX);
    673   WebRtcSpl_ZerosArrayW16(inst->synthesisBuffer, ANAL_BLOCKL_MAX);
    674 
    675   // for HB processing
    676   WebRtcSpl_ZerosArrayW16(inst->dataBufHBFX[0],
    677                           NUM_HIGH_BANDS_MAX * ANAL_BLOCKL_MAX);
    678   // for quantile noise estimation
    679   WebRtcSpl_ZerosArrayW16(inst->noiseEstQuantile, HALF_ANAL_BLOCKL);
    680   for (i = 0; i < SIMULT * HALF_ANAL_BLOCKL; i++) {
    681     inst->noiseEstLogQuantile[i] = 2048; // Q8
    682     inst->noiseEstDensity[i] = 153; // Q9
    683   }
    684   for (i = 0; i < SIMULT; i++) {
    685     inst->noiseEstCounter[i] = (int16_t)(END_STARTUP_LONG * (i + 1)) / SIMULT;
    686   }
    687 
    688   // Initialize suppression filter with ones
    689   WebRtcSpl_MemSetW16((int16_t*)inst->noiseSupFilter, 16384, HALF_ANAL_BLOCKL);
    690 
    691   // Set the aggressiveness: default
    692   inst->aggrMode = 0;
    693 
    694   //initialize variables for new method
    695   inst->priorNonSpeechProb = 8192; // Q14(0.5) prior probability for speech/noise
    696   for (i = 0; i < HALF_ANAL_BLOCKL; i++) {
    697     inst->prevMagnU16[i] = 0;
    698     inst->prevNoiseU32[i] = 0; //previous noise-spectrum
    699     inst->logLrtTimeAvgW32[i] = 0; //smooth LR ratio
    700     inst->avgMagnPause[i] = 0; //conservative noise spectrum estimate
    701     inst->initMagnEst[i] = 0; //initial average magnitude spectrum
    702   }
    703 
    704   //feature quantities
    705   inst->thresholdSpecDiff = 50; //threshold for difference feature: determined on-line
    706   inst->thresholdSpecFlat = 20480; //threshold for flatness: determined on-line
    707   inst->featureLogLrt = inst->thresholdLogLrt; //average LRT factor (= threshold)
    708   inst->featureSpecFlat = inst->thresholdSpecFlat; //spectral flatness (= threshold)
    709   inst->featureSpecDiff = inst->thresholdSpecDiff; //spectral difference (= threshold)
    710   inst->weightLogLrt = 6; //default weighting par for LRT feature
    711   inst->weightSpecFlat = 0; //default weighting par for spectral flatness feature
    712   inst->weightSpecDiff = 0; //default weighting par for spectral difference feature
    713 
    714   inst->curAvgMagnEnergy = 0; //window time-average of input magnitude spectrum
    715   inst->timeAvgMagnEnergy = 0; //normalization for spectral difference
    716   inst->timeAvgMagnEnergyTmp = 0; //normalization for spectral difference
    717 
    718   //histogram quantities: used to estimate/update thresholds for features
    719   WebRtcSpl_ZerosArrayW16(inst->histLrt, HIST_PAR_EST);
    720   WebRtcSpl_ZerosArrayW16(inst->histSpecDiff, HIST_PAR_EST);
    721   WebRtcSpl_ZerosArrayW16(inst->histSpecFlat, HIST_PAR_EST);
    722 
    723   inst->blockIndex = -1; //frame counter
    724 
    725   //inst->modelUpdate    = 500;   //window for update
    726   inst->modelUpdate = (1 << STAT_UPDATES); //window for update
    727   inst->cntThresUpdate = 0; //counter feature thresholds updates
    728 
    729   inst->sumMagn = 0;
    730   inst->magnEnergy = 0;
    731   inst->prevQMagn = 0;
    732   inst->qNoise = 0;
    733   inst->prevQNoise = 0;
    734 
    735   inst->energyIn = 0;
    736   inst->scaleEnergyIn = 0;
    737 
    738   inst->whiteNoiseLevel = 0;
    739   inst->pinkNoiseNumerator = 0;
    740   inst->pinkNoiseExp = 0;
    741   inst->minNorm = 15; // Start with full scale
    742   inst->zeroInputSignal = 0;
    743 
    744   //default mode
    745   WebRtcNsx_set_policy_core(inst, 0);
    746 
    747 #ifdef NS_FILEDEBUG
    748   inst->infile = fopen("indebug.pcm", "wb");
    749   inst->outfile = fopen("outdebug.pcm", "wb");
    750   inst->file1 = fopen("file1.pcm", "wb");
    751   inst->file2 = fopen("file2.pcm", "wb");
    752   inst->file3 = fopen("file3.pcm", "wb");
    753   inst->file4 = fopen("file4.pcm", "wb");
    754   inst->file5 = fopen("file5.pcm", "wb");
    755 #endif
    756 
    757   // Initialize function pointers.
    758   WebRtcNsx_NoiseEstimation = NoiseEstimationC;
    759   WebRtcNsx_PrepareSpectrum = PrepareSpectrumC;
    760   WebRtcNsx_SynthesisUpdate = SynthesisUpdateC;
    761   WebRtcNsx_AnalysisUpdate = AnalysisUpdateC;
    762   WebRtcNsx_Denormalize = DenormalizeC;
    763   WebRtcNsx_NormalizeRealBuffer = NormalizeRealBufferC;
    764 
    765 #ifdef WEBRTC_DETECT_NEON
    766   uint64_t features = WebRtc_GetCPUFeaturesARM();
    767   if ((features & kCPUFeatureNEON) != 0) {
    768       WebRtcNsx_InitNeon();
    769   }
    770 #elif defined(WEBRTC_HAS_NEON)
    771   WebRtcNsx_InitNeon();
    772 #endif
    773 
    774 #if defined(MIPS32_LE)
    775   WebRtcNsx_InitMips();
    776 #endif
    777 
    778   inst->initFlag = 1;
    779 
    780   return 0;
    781 }
    782 
    783 int WebRtcNsx_set_policy_core(NoiseSuppressionFixedC* inst, int mode) {
    784   // allow for modes:0,1,2,3
    785   if (mode < 0 || mode > 3) {
    786     return -1;
    787   }
    788 
    789   inst->aggrMode = mode;
    790   if (mode == 0) {
    791     inst->overdrive = 256; // Q8(1.0)
    792     inst->denoiseBound = 8192; // Q14(0.5)
    793     inst->gainMap = 0; // No gain compensation
    794   } else if (mode == 1) {
    795     inst->overdrive = 256; // Q8(1.0)
    796     inst->denoiseBound = 4096; // Q14(0.25)
    797     inst->factor2Table = kFactor2Aggressiveness1;
    798     inst->gainMap = 1;
    799   } else if (mode == 2) {
    800     inst->overdrive = 282; // ~= Q8(1.1)
    801     inst->denoiseBound = 2048; // Q14(0.125)
    802     inst->factor2Table = kFactor2Aggressiveness2;
    803     inst->gainMap = 1;
    804   } else if (mode == 3) {
    805     inst->overdrive = 320; // Q8(1.25)
    806     inst->denoiseBound = 1475; // ~= Q14(0.09)
    807     inst->factor2Table = kFactor2Aggressiveness3;
    808     inst->gainMap = 1;
    809   }
    810   return 0;
    811 }
    812 
    813 // Extract thresholds for feature parameters
    814 // histograms are computed over some window_size (given by window_pars)
    815 // thresholds and weights are extracted every window
    816 // flag 0 means update histogram only, flag 1 means compute the thresholds/weights
    817 // threshold and weights are returned in: inst->priorModelPars
    818 void WebRtcNsx_FeatureParameterExtraction(NoiseSuppressionFixedC* inst,
    819                                           int flag) {
    820   uint32_t tmpU32;
    821   uint32_t histIndex;
    822   uint32_t posPeak1SpecFlatFX, posPeak2SpecFlatFX;
    823   uint32_t posPeak1SpecDiffFX, posPeak2SpecDiffFX;
    824 
    825   int32_t tmp32;
    826   int32_t fluctLrtFX, thresFluctLrtFX;
    827   int32_t avgHistLrtFX, avgSquareHistLrtFX, avgHistLrtComplFX;
    828 
    829   int16_t j;
    830   int16_t numHistLrt;
    831 
    832   int i;
    833   int useFeatureSpecFlat, useFeatureSpecDiff, featureSum;
    834   int maxPeak1, maxPeak2;
    835   int weightPeak1SpecFlat, weightPeak2SpecFlat;
    836   int weightPeak1SpecDiff, weightPeak2SpecDiff;
    837 
    838   //update histograms
    839   if (!flag) {
    840     // LRT
    841     // Type casting to UWord32 is safe since negative values will not be wrapped to larger
    842     // values than HIST_PAR_EST
    843     histIndex = (uint32_t)(inst->featureLogLrt);
    844     if (histIndex < HIST_PAR_EST) {
    845       inst->histLrt[histIndex]++;
    846     }
    847     // Spectral flatness
    848     // (inst->featureSpecFlat*20)>>10 = (inst->featureSpecFlat*5)>>8
    849     histIndex = (inst->featureSpecFlat * 5) >> 8;
    850     if (histIndex < HIST_PAR_EST) {
    851       inst->histSpecFlat[histIndex]++;
    852     }
    853     // Spectral difference
    854     histIndex = HIST_PAR_EST;
    855     if (inst->timeAvgMagnEnergy > 0) {
    856       // Guard against division by zero
    857       // If timeAvgMagnEnergy == 0 we have no normalizing statistics and
    858       // therefore can't update the histogram
    859       histIndex = ((inst->featureSpecDiff * 5) >> inst->stages) /
    860           inst->timeAvgMagnEnergy;
    861     }
    862     if (histIndex < HIST_PAR_EST) {
    863       inst->histSpecDiff[histIndex]++;
    864     }
    865   }
    866 
    867   // extract parameters for speech/noise probability
    868   if (flag) {
    869     useFeatureSpecDiff = 1;
    870     //for LRT feature:
    871     // compute the average over inst->featureExtractionParams.rangeAvgHistLrt
    872     avgHistLrtFX = 0;
    873     avgSquareHistLrtFX = 0;
    874     numHistLrt = 0;
    875     for (i = 0; i < BIN_SIZE_LRT; i++) {
    876       j = (2 * i + 1);
    877       tmp32 = inst->histLrt[i] * j;
    878       avgHistLrtFX += tmp32;
    879       numHistLrt += inst->histLrt[i];
    880       avgSquareHistLrtFX += tmp32 * j;
    881     }
    882     avgHistLrtComplFX = avgHistLrtFX;
    883     for (; i < HIST_PAR_EST; i++) {
    884       j = (2 * i + 1);
    885       tmp32 = inst->histLrt[i] * j;
    886       avgHistLrtComplFX += tmp32;
    887       avgSquareHistLrtFX += tmp32 * j;
    888     }
    889     fluctLrtFX = avgSquareHistLrtFX * numHistLrt -
    890         avgHistLrtFX * avgHistLrtComplFX;
    891     thresFluctLrtFX = THRES_FLUCT_LRT * numHistLrt;
    892     // get threshold for LRT feature:
    893     tmpU32 = (FACTOR_1_LRT_DIFF * (uint32_t)avgHistLrtFX);
    894     if ((fluctLrtFX < thresFluctLrtFX) || (numHistLrt == 0) ||
    895         (tmpU32 > (uint32_t)(100 * numHistLrt))) {
    896       //very low fluctuation, so likely noise
    897       inst->thresholdLogLrt = inst->maxLrt;
    898     } else {
    899       tmp32 = (int32_t)((tmpU32 << (9 + inst->stages)) / numHistLrt /
    900                               25);
    901       // check if value is within min/max range
    902       inst->thresholdLogLrt = WEBRTC_SPL_SAT(inst->maxLrt,
    903                                              tmp32,
    904                                              inst->minLrt);
    905     }
    906     if (fluctLrtFX < thresFluctLrtFX) {
    907       // Do not use difference feature if fluctuation of LRT feature is very low:
    908       // most likely just noise state
    909       useFeatureSpecDiff = 0;
    910     }
    911 
    912     // for spectral flatness and spectral difference: compute the main peaks of histogram
    913     maxPeak1 = 0;
    914     maxPeak2 = 0;
    915     posPeak1SpecFlatFX = 0;
    916     posPeak2SpecFlatFX = 0;
    917     weightPeak1SpecFlat = 0;
    918     weightPeak2SpecFlat = 0;
    919 
    920     // peaks for flatness
    921     for (i = 0; i < HIST_PAR_EST; i++) {
    922       if (inst->histSpecFlat[i] > maxPeak1) {
    923         // Found new "first" peak
    924         maxPeak2 = maxPeak1;
    925         weightPeak2SpecFlat = weightPeak1SpecFlat;
    926         posPeak2SpecFlatFX = posPeak1SpecFlatFX;
    927 
    928         maxPeak1 = inst->histSpecFlat[i];
    929         weightPeak1SpecFlat = inst->histSpecFlat[i];
    930         posPeak1SpecFlatFX = (uint32_t)(2 * i + 1);
    931       } else if (inst->histSpecFlat[i] > maxPeak2) {
    932         // Found new "second" peak
    933         maxPeak2 = inst->histSpecFlat[i];
    934         weightPeak2SpecFlat = inst->histSpecFlat[i];
    935         posPeak2SpecFlatFX = (uint32_t)(2 * i + 1);
    936       }
    937     }
    938 
    939     // for spectral flatness feature
    940     useFeatureSpecFlat = 1;
    941     // merge the two peaks if they are close
    942     if ((posPeak1SpecFlatFX - posPeak2SpecFlatFX < LIM_PEAK_SPACE_FLAT_DIFF)
    943         && (weightPeak2SpecFlat * LIM_PEAK_WEIGHT_FLAT_DIFF > weightPeak1SpecFlat)) {
    944       weightPeak1SpecFlat += weightPeak2SpecFlat;
    945       posPeak1SpecFlatFX = (posPeak1SpecFlatFX + posPeak2SpecFlatFX) >> 1;
    946     }
    947     //reject if weight of peaks is not large enough, or peak value too small
    948     if (weightPeak1SpecFlat < THRES_WEIGHT_FLAT_DIFF || posPeak1SpecFlatFX
    949         < THRES_PEAK_FLAT) {
    950       useFeatureSpecFlat = 0;
    951     } else { // if selected, get the threshold
    952       // compute the threshold and check if value is within min/max range
    953       inst->thresholdSpecFlat = WEBRTC_SPL_SAT(MAX_FLAT_Q10, FACTOR_2_FLAT_Q10
    954                                                * posPeak1SpecFlatFX, MIN_FLAT_Q10); //Q10
    955     }
    956     // done with flatness feature
    957 
    958     if (useFeatureSpecDiff) {
    959       //compute two peaks for spectral difference
    960       maxPeak1 = 0;
    961       maxPeak2 = 0;
    962       posPeak1SpecDiffFX = 0;
    963       posPeak2SpecDiffFX = 0;
    964       weightPeak1SpecDiff = 0;
    965       weightPeak2SpecDiff = 0;
    966       // peaks for spectral difference
    967       for (i = 0; i < HIST_PAR_EST; i++) {
    968         if (inst->histSpecDiff[i] > maxPeak1) {
    969           // Found new "first" peak
    970           maxPeak2 = maxPeak1;
    971           weightPeak2SpecDiff = weightPeak1SpecDiff;
    972           posPeak2SpecDiffFX = posPeak1SpecDiffFX;
    973 
    974           maxPeak1 = inst->histSpecDiff[i];
    975           weightPeak1SpecDiff = inst->histSpecDiff[i];
    976           posPeak1SpecDiffFX = (uint32_t)(2 * i + 1);
    977         } else if (inst->histSpecDiff[i] > maxPeak2) {
    978           // Found new "second" peak
    979           maxPeak2 = inst->histSpecDiff[i];
    980           weightPeak2SpecDiff = inst->histSpecDiff[i];
    981           posPeak2SpecDiffFX = (uint32_t)(2 * i + 1);
    982         }
    983       }
    984 
    985       // merge the two peaks if they are close
    986       if ((posPeak1SpecDiffFX - posPeak2SpecDiffFX < LIM_PEAK_SPACE_FLAT_DIFF)
    987           && (weightPeak2SpecDiff * LIM_PEAK_WEIGHT_FLAT_DIFF > weightPeak1SpecDiff)) {
    988         weightPeak1SpecDiff += weightPeak2SpecDiff;
    989         posPeak1SpecDiffFX = (posPeak1SpecDiffFX + posPeak2SpecDiffFX) >> 1;
    990       }
    991       // get the threshold value and check if value is within min/max range
    992       inst->thresholdSpecDiff = WEBRTC_SPL_SAT(MAX_DIFF, FACTOR_1_LRT_DIFF
    993                                                * posPeak1SpecDiffFX, MIN_DIFF); //5x bigger
    994       //reject if weight of peaks is not large enough
    995       if (weightPeak1SpecDiff < THRES_WEIGHT_FLAT_DIFF) {
    996         useFeatureSpecDiff = 0;
    997       }
    998       // done with spectral difference feature
    999     }
   1000 
   1001     // select the weights between the features
   1002     // inst->priorModelPars[4] is weight for LRT: always selected
   1003     featureSum = 6 / (1 + useFeatureSpecFlat + useFeatureSpecDiff);
   1004     inst->weightLogLrt = featureSum;
   1005     inst->weightSpecFlat = useFeatureSpecFlat * featureSum;
   1006     inst->weightSpecDiff = useFeatureSpecDiff * featureSum;
   1007 
   1008     // set histograms to zero for next update
   1009     WebRtcSpl_ZerosArrayW16(inst->histLrt, HIST_PAR_EST);
   1010     WebRtcSpl_ZerosArrayW16(inst->histSpecDiff, HIST_PAR_EST);
   1011     WebRtcSpl_ZerosArrayW16(inst->histSpecFlat, HIST_PAR_EST);
   1012   }  // end of flag == 1
   1013 }
   1014 
   1015 
   1016 // Compute spectral flatness on input spectrum
   1017 // magn is the magnitude spectrum
   1018 // spectral flatness is returned in inst->featureSpecFlat
   1019 void WebRtcNsx_ComputeSpectralFlatness(NoiseSuppressionFixedC* inst,
   1020                                        uint16_t* magn) {
   1021   uint32_t tmpU32;
   1022   uint32_t avgSpectralFlatnessNum, avgSpectralFlatnessDen;
   1023 
   1024   int32_t tmp32;
   1025   int32_t currentSpectralFlatness, logCurSpectralFlatness;
   1026 
   1027   int16_t zeros, frac, intPart;
   1028 
   1029   size_t i;
   1030 
   1031   // for flatness
   1032   avgSpectralFlatnessNum = 0;
   1033   avgSpectralFlatnessDen = inst->sumMagn - (uint32_t)magn[0]; // Q(normData-stages)
   1034 
   1035   // compute log of ratio of the geometric to arithmetic mean: check for log(0) case
   1036   // flatness = exp( sum(log(magn[i]))/N - log(sum(magn[i])/N) )
   1037   //          = exp( sum(log(magn[i]))/N ) * N / sum(magn[i])
   1038   //          = 2^( sum(log2(magn[i]))/N - (log2(sum(magn[i])) - log2(N)) ) [This is used]
   1039   for (i = 1; i < inst->magnLen; i++) {
   1040     // First bin is excluded from spectrum measures. Number of bins is now a power of 2
   1041     if (magn[i]) {
   1042       zeros = WebRtcSpl_NormU32((uint32_t)magn[i]);
   1043       frac = (int16_t)(((uint32_t)((uint32_t)(magn[i]) << zeros)
   1044                               & 0x7FFFFFFF) >> 23);
   1045       // log2(magn(i))
   1046       assert(frac < 256);
   1047       tmpU32 = (uint32_t)(((31 - zeros) << 8)
   1048                                 + WebRtcNsx_kLogTableFrac[frac]); // Q8
   1049       avgSpectralFlatnessNum += tmpU32; // Q8
   1050     } else {
   1051       //if at least one frequency component is zero, treat separately
   1052       tmpU32 = WEBRTC_SPL_UMUL_32_16(inst->featureSpecFlat, SPECT_FLAT_TAVG_Q14); // Q24
   1053       inst->featureSpecFlat -= tmpU32 >> 14;  // Q10
   1054       return;
   1055     }
   1056   }
   1057   //ratio and inverse log: check for case of log(0)
   1058   zeros = WebRtcSpl_NormU32(avgSpectralFlatnessDen);
   1059   frac = (int16_t)(((avgSpectralFlatnessDen << zeros) & 0x7FFFFFFF) >> 23);
   1060   // log2(avgSpectralFlatnessDen)
   1061   assert(frac < 256);
   1062   tmp32 = (int32_t)(((31 - zeros) << 8) + WebRtcNsx_kLogTableFrac[frac]); // Q8
   1063   logCurSpectralFlatness = (int32_t)avgSpectralFlatnessNum;
   1064   logCurSpectralFlatness += ((int32_t)(inst->stages - 1) << (inst->stages + 7)); // Q(8+stages-1)
   1065   logCurSpectralFlatness -= (tmp32 << (inst->stages - 1));
   1066   logCurSpectralFlatness <<= (10 - inst->stages);  // Q17
   1067   tmp32 = (int32_t)(0x00020000 | (WEBRTC_SPL_ABS_W32(logCurSpectralFlatness)
   1068                                         & 0x0001FFFF)); //Q17
   1069   intPart = 7 - (logCurSpectralFlatness >> 17);  // Add 7 for output in Q10.
   1070   if (intPart > 0) {
   1071     currentSpectralFlatness = tmp32 >> intPart;
   1072   } else {
   1073     currentSpectralFlatness = tmp32 << -intPart;
   1074   }
   1075 
   1076   //time average update of spectral flatness feature
   1077   tmp32 = currentSpectralFlatness - (int32_t)inst->featureSpecFlat; // Q10
   1078   tmp32 *= SPECT_FLAT_TAVG_Q14;  // Q24
   1079   inst->featureSpecFlat += tmp32 >> 14;  // Q10
   1080   // done with flatness feature
   1081 }
   1082 
   1083 
   1084 // Compute the difference measure between input spectrum and a template/learned noise spectrum
   1085 // magn_tmp is the input spectrum
   1086 // the reference/template spectrum is  inst->magn_avg_pause[i]
   1087 // returns (normalized) spectral difference in inst->featureSpecDiff
   1088 void WebRtcNsx_ComputeSpectralDifference(NoiseSuppressionFixedC* inst,
   1089                                          uint16_t* magnIn) {
   1090   // This is to be calculated:
   1091   // avgDiffNormMagn = var(magnIn) - cov(magnIn, magnAvgPause)^2 / var(magnAvgPause)
   1092 
   1093   uint32_t tmpU32no1, tmpU32no2;
   1094   uint32_t varMagnUFX, varPauseUFX, avgDiffNormMagnUFX;
   1095 
   1096   int32_t tmp32no1, tmp32no2;
   1097   int32_t avgPauseFX, avgMagnFX, covMagnPauseFX;
   1098   int32_t maxPause, minPause;
   1099 
   1100   int16_t tmp16no1;
   1101 
   1102   size_t i;
   1103   int norm32, nShifts;
   1104 
   1105   avgPauseFX = 0;
   1106   maxPause = 0;
   1107   minPause = inst->avgMagnPause[0]; // Q(prevQMagn)
   1108   // compute average quantities
   1109   for (i = 0; i < inst->magnLen; i++) {
   1110     // Compute mean of magn_pause
   1111     avgPauseFX += inst->avgMagnPause[i]; // in Q(prevQMagn)
   1112     maxPause = WEBRTC_SPL_MAX(maxPause, inst->avgMagnPause[i]);
   1113     minPause = WEBRTC_SPL_MIN(minPause, inst->avgMagnPause[i]);
   1114   }
   1115   // normalize by replacing div of "inst->magnLen" with "inst->stages-1" shifts
   1116   avgPauseFX >>= inst->stages - 1;
   1117   avgMagnFX = inst->sumMagn >> (inst->stages - 1);
   1118   // Largest possible deviation in magnPause for (co)var calculations
   1119   tmp32no1 = WEBRTC_SPL_MAX(maxPause - avgPauseFX, avgPauseFX - minPause);
   1120   // Get number of shifts to make sure we don't get wrap around in varPause
   1121   nShifts = WEBRTC_SPL_MAX(0, 10 + inst->stages - WebRtcSpl_NormW32(tmp32no1));
   1122 
   1123   varMagnUFX = 0;
   1124   varPauseUFX = 0;
   1125   covMagnPauseFX = 0;
   1126   for (i = 0; i < inst->magnLen; i++) {
   1127     // Compute var and cov of magn and magn_pause
   1128     tmp16no1 = (int16_t)((int32_t)magnIn[i] - avgMagnFX);
   1129     tmp32no2 = inst->avgMagnPause[i] - avgPauseFX;
   1130     varMagnUFX += (uint32_t)(tmp16no1 * tmp16no1);  // Q(2*qMagn)
   1131     tmp32no1 = tmp32no2 * tmp16no1;  // Q(prevQMagn+qMagn)
   1132     covMagnPauseFX += tmp32no1; // Q(prevQMagn+qMagn)
   1133     tmp32no1 = tmp32no2 >> nShifts;  // Q(prevQMagn-minPause).
   1134     varPauseUFX += tmp32no1 * tmp32no1;  // Q(2*(prevQMagn-minPause))
   1135   }
   1136   //update of average magnitude spectrum: Q(-2*stages) and averaging replaced by shifts
   1137   inst->curAvgMagnEnergy +=
   1138       inst->magnEnergy >> (2 * inst->normData + inst->stages - 1);
   1139 
   1140   avgDiffNormMagnUFX = varMagnUFX; // Q(2*qMagn)
   1141   if ((varPauseUFX) && (covMagnPauseFX)) {
   1142     tmpU32no1 = (uint32_t)WEBRTC_SPL_ABS_W32(covMagnPauseFX); // Q(prevQMagn+qMagn)
   1143     norm32 = WebRtcSpl_NormU32(tmpU32no1) - 16;
   1144     if (norm32 > 0) {
   1145       tmpU32no1 <<= norm32;  // Q(prevQMagn+qMagn+norm32)
   1146     } else {
   1147       tmpU32no1 >>= -norm32;  // Q(prevQMagn+qMagn+norm32)
   1148     }
   1149     tmpU32no2 = WEBRTC_SPL_UMUL(tmpU32no1, tmpU32no1); // Q(2*(prevQMagn+qMagn-norm32))
   1150 
   1151     nShifts += norm32;
   1152     nShifts <<= 1;
   1153     if (nShifts < 0) {
   1154       varPauseUFX >>= (-nShifts); // Q(2*(qMagn+norm32+minPause))
   1155       nShifts = 0;
   1156     }
   1157     if (varPauseUFX > 0) {
   1158       // Q(2*(qMagn+norm32-16+minPause))
   1159       tmpU32no1 = tmpU32no2 / varPauseUFX;
   1160       tmpU32no1 >>= nShifts;
   1161 
   1162       // Q(2*qMagn)
   1163       avgDiffNormMagnUFX -= WEBRTC_SPL_MIN(avgDiffNormMagnUFX, tmpU32no1);
   1164     } else {
   1165       avgDiffNormMagnUFX = 0;
   1166     }
   1167   }
   1168   //normalize and compute time average update of difference feature
   1169   tmpU32no1 = avgDiffNormMagnUFX >> (2 * inst->normData);
   1170   if (inst->featureSpecDiff > tmpU32no1) {
   1171     tmpU32no2 = WEBRTC_SPL_UMUL_32_16(inst->featureSpecDiff - tmpU32no1,
   1172                                       SPECT_DIFF_TAVG_Q8); // Q(8-2*stages)
   1173     inst->featureSpecDiff -= tmpU32no2 >> 8;  // Q(-2*stages)
   1174   } else {
   1175     tmpU32no2 = WEBRTC_SPL_UMUL_32_16(tmpU32no1 - inst->featureSpecDiff,
   1176                                       SPECT_DIFF_TAVG_Q8); // Q(8-2*stages)
   1177     inst->featureSpecDiff += tmpU32no2 >> 8;  // Q(-2*stages)
   1178   }
   1179 }
   1180 
   1181 // Transform input (speechFrame) to frequency domain magnitude (magnU16)
   1182 void WebRtcNsx_DataAnalysis(NoiseSuppressionFixedC* inst,
   1183                             short* speechFrame,
   1184                             uint16_t* magnU16) {
   1185   uint32_t tmpU32no1;
   1186 
   1187   int32_t   tmp_1_w32 = 0;
   1188   int32_t   tmp_2_w32 = 0;
   1189   int32_t   sum_log_magn = 0;
   1190   int32_t   sum_log_i_log_magn = 0;
   1191 
   1192   uint16_t  sum_log_magn_u16 = 0;
   1193   uint16_t  tmp_u16 = 0;
   1194 
   1195   int16_t   sum_log_i = 0;
   1196   int16_t   sum_log_i_square = 0;
   1197   int16_t   frac = 0;
   1198   int16_t   log2 = 0;
   1199   int16_t   matrix_determinant = 0;
   1200   int16_t   maxWinData;
   1201 
   1202   size_t i, j;
   1203   int zeros;
   1204   int net_norm = 0;
   1205   int right_shifts_in_magnU16 = 0;
   1206   int right_shifts_in_initMagnEst = 0;
   1207 
   1208   int16_t winData_buff[ANAL_BLOCKL_MAX * 2 + 16];
   1209   int16_t realImag_buff[ANAL_BLOCKL_MAX * 2 + 16];
   1210 
   1211   // Align the structures to 32-byte boundary for the FFT function.
   1212   int16_t* winData = (int16_t*) (((uintptr_t)winData_buff + 31) & ~31);
   1213   int16_t* realImag = (int16_t*) (((uintptr_t) realImag_buff + 31) & ~31);
   1214 
   1215   // Update analysis buffer for lower band, and window data before FFT.
   1216   WebRtcNsx_AnalysisUpdate(inst, winData, speechFrame);
   1217 
   1218   // Get input energy
   1219   inst->energyIn =
   1220       WebRtcSpl_Energy(winData, inst->anaLen, &inst->scaleEnergyIn);
   1221 
   1222   // Reset zero input flag
   1223   inst->zeroInputSignal = 0;
   1224   // Acquire norm for winData
   1225   maxWinData = WebRtcSpl_MaxAbsValueW16(winData, inst->anaLen);
   1226   inst->normData = WebRtcSpl_NormW16(maxWinData);
   1227   if (maxWinData == 0) {
   1228     // Treat zero input separately.
   1229     inst->zeroInputSignal = 1;
   1230     return;
   1231   }
   1232 
   1233   // Determine the net normalization in the frequency domain
   1234   net_norm = inst->stages - inst->normData;
   1235   // Track lowest normalization factor and use it to prevent wrap around in shifting
   1236   right_shifts_in_magnU16 = inst->normData - inst->minNorm;
   1237   right_shifts_in_initMagnEst = WEBRTC_SPL_MAX(-right_shifts_in_magnU16, 0);
   1238   inst->minNorm -= right_shifts_in_initMagnEst;
   1239   right_shifts_in_magnU16 = WEBRTC_SPL_MAX(right_shifts_in_magnU16, 0);
   1240 
   1241   // create realImag as winData interleaved with zeros (= imag. part), normalize it
   1242   WebRtcNsx_NormalizeRealBuffer(inst, winData, realImag);
   1243 
   1244   // FFT output will be in winData[].
   1245   WebRtcSpl_RealForwardFFT(inst->real_fft, realImag, winData);
   1246 
   1247   inst->imag[0] = 0; // Q(normData-stages)
   1248   inst->imag[inst->anaLen2] = 0;
   1249   inst->real[0] = winData[0]; // Q(normData-stages)
   1250   inst->real[inst->anaLen2] = winData[inst->anaLen];
   1251   // Q(2*(normData-stages))
   1252   inst->magnEnergy = (uint32_t)(inst->real[0] * inst->real[0]);
   1253   inst->magnEnergy += (uint32_t)(inst->real[inst->anaLen2] *
   1254                                  inst->real[inst->anaLen2]);
   1255   magnU16[0] = (uint16_t)WEBRTC_SPL_ABS_W16(inst->real[0]); // Q(normData-stages)
   1256   magnU16[inst->anaLen2] = (uint16_t)WEBRTC_SPL_ABS_W16(inst->real[inst->anaLen2]);
   1257   inst->sumMagn = (uint32_t)magnU16[0]; // Q(normData-stages)
   1258   inst->sumMagn += (uint32_t)magnU16[inst->anaLen2];
   1259 
   1260   if (inst->blockIndex >= END_STARTUP_SHORT) {
   1261     for (i = 1, j = 2; i < inst->anaLen2; i += 1, j += 2) {
   1262       inst->real[i] = winData[j];
   1263       inst->imag[i] = -winData[j + 1];
   1264       // magnitude spectrum
   1265       // energy in Q(2*(normData-stages))
   1266       tmpU32no1 = (uint32_t)(winData[j] * winData[j]);
   1267       tmpU32no1 += (uint32_t)(winData[j + 1] * winData[j + 1]);
   1268       inst->magnEnergy += tmpU32no1; // Q(2*(normData-stages))
   1269 
   1270       magnU16[i] = (uint16_t)WebRtcSpl_SqrtFloor(tmpU32no1); // Q(normData-stages)
   1271       inst->sumMagn += (uint32_t)magnU16[i]; // Q(normData-stages)
   1272     }
   1273   } else {
   1274     //
   1275     // Gather information during startup for noise parameter estimation
   1276     //
   1277 
   1278     // Switch initMagnEst to Q(minNorm-stages)
   1279     inst->initMagnEst[0] >>= right_shifts_in_initMagnEst;
   1280     inst->initMagnEst[inst->anaLen2] >>= right_shifts_in_initMagnEst;
   1281 
   1282     // Update initMagnEst with magnU16 in Q(minNorm-stages).
   1283     inst->initMagnEst[0] += magnU16[0] >> right_shifts_in_magnU16;
   1284     inst->initMagnEst[inst->anaLen2] +=
   1285         magnU16[inst->anaLen2] >> right_shifts_in_magnU16;
   1286 
   1287     log2 = 0;
   1288     if (magnU16[inst->anaLen2]) {
   1289       // Calculate log2(magnU16[inst->anaLen2])
   1290       zeros = WebRtcSpl_NormU32((uint32_t)magnU16[inst->anaLen2]);
   1291       frac = (int16_t)((((uint32_t)magnU16[inst->anaLen2] << zeros) &
   1292                               0x7FFFFFFF) >> 23); // Q8
   1293       // log2(magnU16(i)) in Q8
   1294       assert(frac < 256);
   1295       log2 = (int16_t)(((31 - zeros) << 8) + WebRtcNsx_kLogTableFrac[frac]);
   1296     }
   1297 
   1298     sum_log_magn = (int32_t)log2; // Q8
   1299     // sum_log_i_log_magn in Q17
   1300     sum_log_i_log_magn = (kLogIndex[inst->anaLen2] * log2) >> 3;
   1301 
   1302     for (i = 1, j = 2; i < inst->anaLen2; i += 1, j += 2) {
   1303       inst->real[i] = winData[j];
   1304       inst->imag[i] = -winData[j + 1];
   1305       // magnitude spectrum
   1306       // energy in Q(2*(normData-stages))
   1307       tmpU32no1 = (uint32_t)(winData[j] * winData[j]);
   1308       tmpU32no1 += (uint32_t)(winData[j + 1] * winData[j + 1]);
   1309       inst->magnEnergy += tmpU32no1; // Q(2*(normData-stages))
   1310 
   1311       magnU16[i] = (uint16_t)WebRtcSpl_SqrtFloor(tmpU32no1); // Q(normData-stages)
   1312       inst->sumMagn += (uint32_t)magnU16[i]; // Q(normData-stages)
   1313 
   1314       // Switch initMagnEst to Q(minNorm-stages)
   1315       inst->initMagnEst[i] >>= right_shifts_in_initMagnEst;
   1316 
   1317       // Update initMagnEst with magnU16 in Q(minNorm-stages).
   1318       inst->initMagnEst[i] += magnU16[i] >> right_shifts_in_magnU16;
   1319 
   1320       if (i >= kStartBand) {
   1321         // For pink noise estimation. Collect data neglecting lower frequency band
   1322         log2 = 0;
   1323         if (magnU16[i]) {
   1324           zeros = WebRtcSpl_NormU32((uint32_t)magnU16[i]);
   1325           frac = (int16_t)((((uint32_t)magnU16[i] << zeros) &
   1326                                   0x7FFFFFFF) >> 23);
   1327           // log2(magnU16(i)) in Q8
   1328           assert(frac < 256);
   1329           log2 = (int16_t)(((31 - zeros) << 8)
   1330                                  + WebRtcNsx_kLogTableFrac[frac]);
   1331         }
   1332         sum_log_magn += (int32_t)log2; // Q8
   1333         // sum_log_i_log_magn in Q17
   1334         sum_log_i_log_magn += (kLogIndex[i] * log2) >> 3;
   1335       }
   1336     }
   1337 
   1338     //
   1339     //compute simplified noise model during startup
   1340     //
   1341 
   1342     // Estimate White noise
   1343 
   1344     // Switch whiteNoiseLevel to Q(minNorm-stages)
   1345     inst->whiteNoiseLevel >>= right_shifts_in_initMagnEst;
   1346 
   1347     // Update the average magnitude spectrum, used as noise estimate.
   1348     tmpU32no1 = WEBRTC_SPL_UMUL_32_16(inst->sumMagn, inst->overdrive);
   1349     tmpU32no1 >>= inst->stages + 8;
   1350 
   1351     // Replacing division above with 'stages' shifts
   1352     // Shift to same Q-domain as whiteNoiseLevel
   1353     tmpU32no1 >>= right_shifts_in_magnU16;
   1354     // This operation is safe from wrap around as long as END_STARTUP_SHORT < 128
   1355     assert(END_STARTUP_SHORT < 128);
   1356     inst->whiteNoiseLevel += tmpU32no1; // Q(minNorm-stages)
   1357 
   1358     // Estimate Pink noise parameters
   1359     // Denominator used in both parameter estimates.
   1360     // The value is only dependent on the size of the frequency band (kStartBand)
   1361     // and to reduce computational complexity stored in a table (kDeterminantEstMatrix[])
   1362     assert(kStartBand < 66);
   1363     matrix_determinant = kDeterminantEstMatrix[kStartBand]; // Q0
   1364     sum_log_i = kSumLogIndex[kStartBand]; // Q5
   1365     sum_log_i_square = kSumSquareLogIndex[kStartBand]; // Q2
   1366     if (inst->fs == 8000) {
   1367       // Adjust values to shorter blocks in narrow band.
   1368       tmp_1_w32 = (int32_t)matrix_determinant;
   1369       tmp_1_w32 += (kSumLogIndex[65] * sum_log_i) >> 9;
   1370       tmp_1_w32 -= (kSumLogIndex[65] * kSumLogIndex[65]) >> 10;
   1371       tmp_1_w32 -= (int32_t)sum_log_i_square << 4;
   1372       tmp_1_w32 -= ((inst->magnLen - kStartBand) * kSumSquareLogIndex[65]) >> 2;
   1373       matrix_determinant = (int16_t)tmp_1_w32;
   1374       sum_log_i -= kSumLogIndex[65]; // Q5
   1375       sum_log_i_square -= kSumSquareLogIndex[65]; // Q2
   1376     }
   1377 
   1378     // Necessary number of shifts to fit sum_log_magn in a word16
   1379     zeros = 16 - WebRtcSpl_NormW32(sum_log_magn);
   1380     if (zeros < 0) {
   1381       zeros = 0;
   1382     }
   1383     tmp_1_w32 = sum_log_magn << 1;  // Q9
   1384     sum_log_magn_u16 = (uint16_t)(tmp_1_w32 >> zeros);  // Q(9-zeros).
   1385 
   1386     // Calculate and update pinkNoiseNumerator. Result in Q11.
   1387     tmp_2_w32 = WEBRTC_SPL_MUL_16_U16(sum_log_i_square, sum_log_magn_u16); // Q(11-zeros)
   1388     tmpU32no1 = sum_log_i_log_magn >> 12;  // Q5
   1389 
   1390     // Shift the largest value of sum_log_i and tmp32no3 before multiplication
   1391     tmp_u16 = ((uint16_t)sum_log_i << 1);  // Q6
   1392     if ((uint32_t)sum_log_i > tmpU32no1) {
   1393       tmp_u16 >>= zeros;
   1394     } else {
   1395       tmpU32no1 >>= zeros;
   1396     }
   1397     tmp_2_w32 -= (int32_t)WEBRTC_SPL_UMUL_32_16(tmpU32no1, tmp_u16); // Q(11-zeros)
   1398     matrix_determinant >>= zeros;  // Q(-zeros)
   1399     tmp_2_w32 = WebRtcSpl_DivW32W16(tmp_2_w32, matrix_determinant); // Q11
   1400     tmp_2_w32 += (int32_t)net_norm << 11;  // Q11
   1401     if (tmp_2_w32 < 0) {
   1402       tmp_2_w32 = 0;
   1403     }
   1404     inst->pinkNoiseNumerator += tmp_2_w32; // Q11
   1405 
   1406     // Calculate and update pinkNoiseExp. Result in Q14.
   1407     tmp_2_w32 = WEBRTC_SPL_MUL_16_U16(sum_log_i, sum_log_magn_u16); // Q(14-zeros)
   1408     tmp_1_w32 = sum_log_i_log_magn >> (3 + zeros);
   1409     tmp_1_w32 *= inst->magnLen - kStartBand;
   1410     tmp_2_w32 -= tmp_1_w32; // Q(14-zeros)
   1411     if (tmp_2_w32 > 0) {
   1412       // If the exponential parameter is negative force it to zero, which means a
   1413       // flat spectrum.
   1414       tmp_1_w32 = WebRtcSpl_DivW32W16(tmp_2_w32, matrix_determinant); // Q14
   1415       inst->pinkNoiseExp += WEBRTC_SPL_SAT(16384, tmp_1_w32, 0); // Q14
   1416     }
   1417   }
   1418 }
   1419 
   1420 void WebRtcNsx_DataSynthesis(NoiseSuppressionFixedC* inst, short* outFrame) {
   1421   int32_t energyOut;
   1422 
   1423   int16_t realImag_buff[ANAL_BLOCKL_MAX * 2 + 16];
   1424   int16_t rfft_out_buff[ANAL_BLOCKL_MAX * 2 + 16];
   1425 
   1426   // Align the structures to 32-byte boundary for the FFT function.
   1427   int16_t* realImag = (int16_t*) (((uintptr_t)realImag_buff + 31) & ~31);
   1428   int16_t* rfft_out = (int16_t*) (((uintptr_t) rfft_out_buff + 31) & ~31);
   1429 
   1430   int16_t tmp16no1, tmp16no2;
   1431   int16_t energyRatio;
   1432   int16_t gainFactor, gainFactor1, gainFactor2;
   1433 
   1434   size_t i;
   1435   int outCIFFT;
   1436   int scaleEnergyOut = 0;
   1437 
   1438   if (inst->zeroInputSignal) {
   1439     // synthesize the special case of zero input
   1440     // read out fully processed segment
   1441     for (i = 0; i < inst->blockLen10ms; i++) {
   1442       outFrame[i] = inst->synthesisBuffer[i]; // Q0
   1443     }
   1444     // update synthesis buffer
   1445     memcpy(inst->synthesisBuffer, inst->synthesisBuffer + inst->blockLen10ms,
   1446         (inst->anaLen - inst->blockLen10ms) * sizeof(*inst->synthesisBuffer));
   1447     WebRtcSpl_ZerosArrayW16(inst->synthesisBuffer + inst->anaLen - inst->blockLen10ms,
   1448                             inst->blockLen10ms);
   1449     return;
   1450   }
   1451 
   1452   // Filter the data in the frequency domain, and create spectrum.
   1453   WebRtcNsx_PrepareSpectrum(inst, realImag);
   1454 
   1455   // Inverse FFT output will be in rfft_out[].
   1456   outCIFFT = WebRtcSpl_RealInverseFFT(inst->real_fft, realImag, rfft_out);
   1457 
   1458   WebRtcNsx_Denormalize(inst, rfft_out, outCIFFT);
   1459 
   1460   //scale factor: only do it after END_STARTUP_LONG time
   1461   gainFactor = 8192; // 8192 = Q13(1.0)
   1462   if (inst->gainMap == 1 &&
   1463       inst->blockIndex > END_STARTUP_LONG &&
   1464       inst->energyIn > 0) {
   1465     // Q(-scaleEnergyOut)
   1466     energyOut = WebRtcSpl_Energy(inst->real, inst->anaLen, &scaleEnergyOut);
   1467     if (scaleEnergyOut == 0 && !(energyOut & 0x7f800000)) {
   1468       energyOut = WEBRTC_SPL_SHIFT_W32(energyOut, 8 + scaleEnergyOut
   1469                                        - inst->scaleEnergyIn);
   1470     } else {
   1471       // |energyIn| is currently in Q(|scaleEnergyIn|), but to later on end up
   1472       // with an |energyRatio| in Q8 we need to change the Q-domain to
   1473       // Q(-8-scaleEnergyOut).
   1474       inst->energyIn >>= 8 + scaleEnergyOut - inst->scaleEnergyIn;
   1475     }
   1476 
   1477     assert(inst->energyIn > 0);
   1478     energyRatio = (energyOut + inst->energyIn / 2) / inst->energyIn;  // Q8
   1479     // Limit the ratio to [0, 1] in Q8, i.e., [0, 256]
   1480     energyRatio = WEBRTC_SPL_SAT(256, energyRatio, 0);
   1481 
   1482     // all done in lookup tables now
   1483     assert(energyRatio < 257);
   1484     gainFactor1 = kFactor1Table[energyRatio]; // Q8
   1485     gainFactor2 = inst->factor2Table[energyRatio]; // Q8
   1486 
   1487     //combine both scales with speech/noise prob: note prior (priorSpeechProb) is not frequency dependent
   1488 
   1489     // factor = inst->priorSpeechProb*factor1 + (1.0-inst->priorSpeechProb)*factor2; // original code
   1490     tmp16no1 = (int16_t)(((16384 - inst->priorNonSpeechProb) * gainFactor1) >>
   1491         14);  // in Q13, where 16384 = Q14(1.0)
   1492     tmp16no2 = (int16_t)((inst->priorNonSpeechProb * gainFactor2) >> 14);
   1493     gainFactor = tmp16no1 + tmp16no2; // Q13
   1494   }  // out of flag_gain_map==1
   1495 
   1496   // Synthesis, read out fully processed segment, and update synthesis buffer.
   1497   WebRtcNsx_SynthesisUpdate(inst, outFrame, gainFactor);
   1498 }
   1499 
   1500 void WebRtcNsx_ProcessCore(NoiseSuppressionFixedC* inst,
   1501                            const short* const* speechFrame,
   1502                            int num_bands,
   1503                            short* const* outFrame) {
   1504   // main routine for noise suppression
   1505 
   1506   uint32_t tmpU32no1, tmpU32no2, tmpU32no3;
   1507   uint32_t satMax, maxNoiseU32;
   1508   uint32_t tmpMagnU32, tmpNoiseU32;
   1509   uint32_t nearMagnEst;
   1510   uint32_t noiseUpdateU32;
   1511   uint32_t noiseU32[HALF_ANAL_BLOCKL];
   1512   uint32_t postLocSnr[HALF_ANAL_BLOCKL];
   1513   uint32_t priorLocSnr[HALF_ANAL_BLOCKL];
   1514   uint32_t prevNearSnr[HALF_ANAL_BLOCKL];
   1515   uint32_t curNearSnr;
   1516   uint32_t priorSnr;
   1517   uint32_t noise_estimate = 0;
   1518   uint32_t noise_estimate_avg = 0;
   1519   uint32_t numerator = 0;
   1520 
   1521   int32_t tmp32no1, tmp32no2;
   1522   int32_t pink_noise_num_avg = 0;
   1523 
   1524   uint16_t tmpU16no1;
   1525   uint16_t magnU16[HALF_ANAL_BLOCKL];
   1526   uint16_t prevNoiseU16[HALF_ANAL_BLOCKL];
   1527   uint16_t nonSpeechProbFinal[HALF_ANAL_BLOCKL];
   1528   uint16_t gammaNoise, prevGammaNoise;
   1529   uint16_t noiseSupFilterTmp[HALF_ANAL_BLOCKL];
   1530 
   1531   int16_t qMagn, qNoise;
   1532   int16_t avgProbSpeechHB, gainModHB, avgFilterGainHB, gainTimeDomainHB;
   1533   int16_t pink_noise_exp_avg = 0;
   1534 
   1535   size_t i, j;
   1536   int nShifts, postShifts;
   1537   int norm32no1, norm32no2;
   1538   int flag, sign;
   1539   int q_domain_to_use = 0;
   1540 
   1541   // Code for ARMv7-Neon platform assumes the following:
   1542   assert(inst->anaLen > 0);
   1543   assert(inst->anaLen2 > 0);
   1544   assert(inst->anaLen % 16 == 0);
   1545   assert(inst->anaLen2 % 8 == 0);
   1546   assert(inst->blockLen10ms > 0);
   1547   assert(inst->blockLen10ms % 16 == 0);
   1548   assert(inst->magnLen == inst->anaLen2 + 1);
   1549 
   1550 #ifdef NS_FILEDEBUG
   1551   if (fwrite(spframe, sizeof(short),
   1552              inst->blockLen10ms, inst->infile) != inst->blockLen10ms) {
   1553     assert(false);
   1554   }
   1555 #endif
   1556 
   1557   // Check that initialization has been done
   1558   assert(inst->initFlag == 1);
   1559   assert((num_bands - 1) <= NUM_HIGH_BANDS_MAX);
   1560 
   1561   const short* const* speechFrameHB = NULL;
   1562   short* const* outFrameHB = NULL;
   1563   size_t num_high_bands = 0;
   1564   if (num_bands > 1) {
   1565     speechFrameHB = &speechFrame[1];
   1566     outFrameHB = &outFrame[1];
   1567     num_high_bands = (size_t)(num_bands - 1);
   1568   }
   1569 
   1570   // Store speechFrame and transform to frequency domain
   1571   WebRtcNsx_DataAnalysis(inst, (short*)speechFrame[0], magnU16);
   1572 
   1573   if (inst->zeroInputSignal) {
   1574     WebRtcNsx_DataSynthesis(inst, outFrame[0]);
   1575 
   1576     if (num_bands > 1) {
   1577       // update analysis buffer for H band
   1578       // append new data to buffer FX
   1579       for (i = 0; i < num_high_bands; ++i) {
   1580         int block_shift = inst->anaLen - inst->blockLen10ms;
   1581         memcpy(inst->dataBufHBFX[i], inst->dataBufHBFX[i] + inst->blockLen10ms,
   1582             block_shift * sizeof(*inst->dataBufHBFX[i]));
   1583         memcpy(inst->dataBufHBFX[i] + block_shift, speechFrameHB[i],
   1584             inst->blockLen10ms * sizeof(*inst->dataBufHBFX[i]));
   1585         for (j = 0; j < inst->blockLen10ms; j++) {
   1586           outFrameHB[i][j] = inst->dataBufHBFX[i][j]; // Q0
   1587         }
   1588       }
   1589     }  // end of H band gain computation
   1590     return;
   1591   }
   1592 
   1593   // Update block index when we have something to process
   1594   inst->blockIndex++;
   1595   //
   1596 
   1597   // Norm of magn
   1598   qMagn = inst->normData - inst->stages;
   1599 
   1600   // Compute spectral flatness on input spectrum
   1601   WebRtcNsx_ComputeSpectralFlatness(inst, magnU16);
   1602 
   1603   // quantile noise estimate
   1604   WebRtcNsx_NoiseEstimation(inst, magnU16, noiseU32, &qNoise);
   1605 
   1606   //noise estimate from previous frame
   1607   for (i = 0; i < inst->magnLen; i++) {
   1608     prevNoiseU16[i] = (uint16_t)(inst->prevNoiseU32[i] >> 11);  // Q(prevQNoise)
   1609   }
   1610 
   1611   if (inst->blockIndex < END_STARTUP_SHORT) {
   1612     // Noise Q-domain to be used later; see description at end of section.
   1613     q_domain_to_use = WEBRTC_SPL_MIN((int)qNoise, inst->minNorm - inst->stages);
   1614 
   1615     // Calculate frequency independent parts in parametric noise estimate and calculate
   1616     // the estimate for the lower frequency band (same values for all frequency bins)
   1617     if (inst->pinkNoiseExp) {
   1618       pink_noise_exp_avg = (int16_t)WebRtcSpl_DivW32W16(inst->pinkNoiseExp,
   1619                                                               (int16_t)(inst->blockIndex + 1)); // Q14
   1620       pink_noise_num_avg = WebRtcSpl_DivW32W16(inst->pinkNoiseNumerator,
   1621                                                (int16_t)(inst->blockIndex + 1)); // Q11
   1622       WebRtcNsx_CalcParametricNoiseEstimate(inst,
   1623                                             pink_noise_exp_avg,
   1624                                             pink_noise_num_avg,
   1625                                             kStartBand,
   1626                                             &noise_estimate,
   1627                                             &noise_estimate_avg);
   1628     } else {
   1629       // Use white noise estimate if we have poor pink noise parameter estimates
   1630       noise_estimate = inst->whiteNoiseLevel; // Q(minNorm-stages)
   1631       noise_estimate_avg = noise_estimate / (inst->blockIndex + 1); // Q(minNorm-stages)
   1632     }
   1633     for (i = 0; i < inst->magnLen; i++) {
   1634       // Estimate the background noise using the pink noise parameters if permitted
   1635       if ((inst->pinkNoiseExp) && (i >= kStartBand)) {
   1636         // Reset noise_estimate
   1637         noise_estimate = 0;
   1638         noise_estimate_avg = 0;
   1639         // Calculate the parametric noise estimate for current frequency bin
   1640         WebRtcNsx_CalcParametricNoiseEstimate(inst,
   1641                                               pink_noise_exp_avg,
   1642                                               pink_noise_num_avg,
   1643                                               i,
   1644                                               &noise_estimate,
   1645                                               &noise_estimate_avg);
   1646       }
   1647       // Calculate parametric Wiener filter
   1648       noiseSupFilterTmp[i] = inst->denoiseBound;
   1649       if (inst->initMagnEst[i]) {
   1650         // numerator = (initMagnEst - noise_estimate * overdrive)
   1651         // Result in Q(8+minNorm-stages)
   1652         tmpU32no1 = WEBRTC_SPL_UMUL_32_16(noise_estimate, inst->overdrive);
   1653         numerator = inst->initMagnEst[i] << 8;
   1654         if (numerator > tmpU32no1) {
   1655           // Suppression filter coefficient larger than zero, so calculate.
   1656           numerator -= tmpU32no1;
   1657 
   1658           // Determine number of left shifts in numerator for best accuracy after
   1659           // division
   1660           nShifts = WebRtcSpl_NormU32(numerator);
   1661           nShifts = WEBRTC_SPL_SAT(6, nShifts, 0);
   1662 
   1663           // Shift numerator to Q(nShifts+8+minNorm-stages)
   1664           numerator <<= nShifts;
   1665 
   1666           // Shift denominator to Q(nShifts-6+minNorm-stages)
   1667           tmpU32no1 = inst->initMagnEst[i] >> (6 - nShifts);
   1668           if (tmpU32no1 == 0) {
   1669             // This is only possible if numerator = 0, in which case
   1670             // we don't need any division.
   1671             tmpU32no1 = 1;
   1672           }
   1673           tmpU32no2 = numerator / tmpU32no1;  // Q14
   1674           noiseSupFilterTmp[i] = (uint16_t)WEBRTC_SPL_SAT(16384, tmpU32no2,
   1675               (uint32_t)(inst->denoiseBound)); // Q14
   1676         }
   1677       }
   1678       // Weight quantile noise 'noiseU32' with modeled noise 'noise_estimate_avg'
   1679       // 'noiseU32 is in Q(qNoise) and 'noise_estimate' in Q(minNorm-stages)
   1680       // To guarantee that we do not get wrap around when shifting to the same domain
   1681       // we use the lowest one. Furthermore, we need to save 6 bits for the weighting.
   1682       // 'noise_estimate_avg' can handle this operation by construction, but 'noiseU32'
   1683       // may not.
   1684 
   1685       // Shift 'noiseU32' to 'q_domain_to_use'
   1686       tmpU32no1 = noiseU32[i] >> (qNoise - q_domain_to_use);
   1687       // Shift 'noise_estimate_avg' to 'q_domain_to_use'
   1688       tmpU32no2 = noise_estimate_avg >>
   1689           (inst->minNorm - inst->stages - q_domain_to_use);
   1690       // Make a simple check to see if we have enough room for weighting 'tmpU32no1'
   1691       // without wrap around
   1692       nShifts = 0;
   1693       if (tmpU32no1 & 0xfc000000) {
   1694         tmpU32no1 >>= 6;
   1695         tmpU32no2 >>= 6;
   1696         nShifts = 6;
   1697       }
   1698       tmpU32no1 *= inst->blockIndex;
   1699       tmpU32no2 *= (END_STARTUP_SHORT - inst->blockIndex);
   1700       // Add them together and divide by startup length
   1701       noiseU32[i] = WebRtcSpl_DivU32U16(tmpU32no1 + tmpU32no2, END_STARTUP_SHORT);
   1702       // Shift back if necessary
   1703       noiseU32[i] <<= nShifts;
   1704     }
   1705     // Update new Q-domain for 'noiseU32'
   1706     qNoise = q_domain_to_use;
   1707   }
   1708   // compute average signal during END_STARTUP_LONG time:
   1709   // used to normalize spectral difference measure
   1710   if (inst->blockIndex < END_STARTUP_LONG) {
   1711     // substituting division with shift ending up in Q(-2*stages)
   1712     inst->timeAvgMagnEnergyTmp +=
   1713         inst->magnEnergy >> (2 * inst->normData + inst->stages - 1);
   1714     inst->timeAvgMagnEnergy = WebRtcSpl_DivU32U16(inst->timeAvgMagnEnergyTmp,
   1715                                                   inst->blockIndex + 1);
   1716   }
   1717 
   1718   //start processing at frames == converged+1
   1719   // STEP 1: compute prior and post SNR based on quantile noise estimates
   1720 
   1721   // compute direct decision (DD) estimate of prior SNR: needed for new method
   1722   satMax = (uint32_t)1048575;// Largest possible value without getting overflow despite shifting 12 steps
   1723   postShifts = 6 + qMagn - qNoise;
   1724   nShifts = 5 - inst->prevQMagn + inst->prevQNoise;
   1725   for (i = 0; i < inst->magnLen; i++) {
   1726     // FLOAT:
   1727     // post SNR
   1728     // postLocSnr[i] = 0.0;
   1729     // if (magn[i] > noise[i])
   1730     // {
   1731     //   postLocSnr[i] = magn[i] / (noise[i] + 0.0001);
   1732     // }
   1733     // // previous post SNR
   1734     // // previous estimate: based on previous frame with gain filter (smooth is previous filter)
   1735     //
   1736     // prevNearSnr[i] = inst->prevMagnU16[i] / (inst->noisePrev[i] + 0.0001) * (inst->smooth[i]);
   1737     //
   1738     // // DD estimate is sum of two terms: current estimate and previous estimate
   1739     // // directed decision update of priorSnr (or we actually store [2*priorSnr+1])
   1740     //
   1741     // priorLocSnr[i] = DD_PR_SNR * prevNearSnr[i] + (1.0 - DD_PR_SNR) * (postLocSnr[i] - 1.0);
   1742 
   1743     // calculate post SNR: output in Q11
   1744     postLocSnr[i] = 2048; // 1.0 in Q11
   1745     tmpU32no1 = (uint32_t)magnU16[i] << 6;  // Q(6+qMagn)
   1746     if (postShifts < 0) {
   1747       tmpU32no2 = noiseU32[i] >> -postShifts;  // Q(6+qMagn)
   1748     } else {
   1749       tmpU32no2 = noiseU32[i] << postShifts;  // Q(6+qMagn)
   1750     }
   1751     if (tmpU32no1 > tmpU32no2) {
   1752       // Current magnitude larger than noise
   1753       tmpU32no1 <<= 11;  // Q(17+qMagn)
   1754       if (tmpU32no2 > 0) {
   1755         tmpU32no1 /= tmpU32no2;  // Q11
   1756         postLocSnr[i] = WEBRTC_SPL_MIN(satMax, tmpU32no1); // Q11
   1757       } else {
   1758         postLocSnr[i] = satMax;
   1759       }
   1760     }
   1761 
   1762     // calculate prevNearSnr[i] and save for later instead of recalculating it later
   1763     // |nearMagnEst| in Q(prevQMagn + 14)
   1764     nearMagnEst = inst->prevMagnU16[i] * inst->noiseSupFilter[i];
   1765     tmpU32no1 = nearMagnEst << 3;  // Q(prevQMagn+17)
   1766     tmpU32no2 = inst->prevNoiseU32[i] >> nShifts;  // Q(prevQMagn+6)
   1767 
   1768     if (tmpU32no2 > 0) {
   1769       tmpU32no1 /= tmpU32no2;  // Q11
   1770       tmpU32no1 = WEBRTC_SPL_MIN(satMax, tmpU32no1); // Q11
   1771     } else {
   1772       tmpU32no1 = satMax; // Q11
   1773     }
   1774     prevNearSnr[i] = tmpU32no1; // Q11
   1775 
   1776     //directed decision update of priorSnr
   1777     tmpU32no1 = WEBRTC_SPL_UMUL_32_16(prevNearSnr[i], DD_PR_SNR_Q11); // Q22
   1778     tmpU32no2 = WEBRTC_SPL_UMUL_32_16(postLocSnr[i] - 2048, ONE_MINUS_DD_PR_SNR_Q11); // Q22
   1779     priorSnr = tmpU32no1 + tmpU32no2 + 512; // Q22 (added 512 for rounding)
   1780     // priorLocSnr = 1 + 2*priorSnr
   1781     priorLocSnr[i] = 2048 + (priorSnr >> 10);  // Q11
   1782   }  // end of loop over frequencies
   1783   // done with step 1: DD computation of prior and post SNR
   1784 
   1785   // STEP 2: compute speech/noise likelihood
   1786 
   1787   //compute difference of input spectrum with learned/estimated noise spectrum
   1788   WebRtcNsx_ComputeSpectralDifference(inst, magnU16);
   1789   //compute histograms for determination of parameters (thresholds and weights for features)
   1790   //parameters are extracted once every window time (=inst->modelUpdate)
   1791   //counter update
   1792   inst->cntThresUpdate++;
   1793   flag = (int)(inst->cntThresUpdate == inst->modelUpdate);
   1794   //update histogram
   1795   WebRtcNsx_FeatureParameterExtraction(inst, flag);
   1796   //compute model parameters
   1797   if (flag) {
   1798     inst->cntThresUpdate = 0; // Reset counter
   1799     //update every window:
   1800     // get normalization for spectral difference for next window estimate
   1801 
   1802     // Shift to Q(-2*stages)
   1803     inst->curAvgMagnEnergy >>= STAT_UPDATES;
   1804 
   1805     tmpU32no1 = (inst->curAvgMagnEnergy + inst->timeAvgMagnEnergy + 1) >> 1; //Q(-2*stages)
   1806     // Update featureSpecDiff
   1807     if ((tmpU32no1 != inst->timeAvgMagnEnergy) && (inst->featureSpecDiff) &&
   1808         (inst->timeAvgMagnEnergy > 0)) {
   1809       norm32no1 = 0;
   1810       tmpU32no3 = tmpU32no1;
   1811       while (0xFFFF0000 & tmpU32no3) {
   1812         tmpU32no3 >>= 1;
   1813         norm32no1++;
   1814       }
   1815       tmpU32no2 = inst->featureSpecDiff;
   1816       while (0xFFFF0000 & tmpU32no2) {
   1817         tmpU32no2 >>= 1;
   1818         norm32no1++;
   1819       }
   1820       tmpU32no3 = WEBRTC_SPL_UMUL(tmpU32no3, tmpU32no2);
   1821       tmpU32no3 /= inst->timeAvgMagnEnergy;
   1822       if (WebRtcSpl_NormU32(tmpU32no3) < norm32no1) {
   1823         inst->featureSpecDiff = 0x007FFFFF;
   1824       } else {
   1825         inst->featureSpecDiff = WEBRTC_SPL_MIN(0x007FFFFF,
   1826                                                tmpU32no3 << norm32no1);
   1827       }
   1828     }
   1829 
   1830     inst->timeAvgMagnEnergy = tmpU32no1; // Q(-2*stages)
   1831     inst->curAvgMagnEnergy = 0;
   1832   }
   1833 
   1834   //compute speech/noise probability
   1835   WebRtcNsx_SpeechNoiseProb(inst, nonSpeechProbFinal, priorLocSnr, postLocSnr);
   1836 
   1837   //time-avg parameter for noise update
   1838   gammaNoise = NOISE_UPDATE_Q8; // Q8
   1839 
   1840   maxNoiseU32 = 0;
   1841   postShifts = inst->prevQNoise - qMagn;
   1842   nShifts = inst->prevQMagn - qMagn;
   1843   for (i = 0; i < inst->magnLen; i++) {
   1844     // temporary noise update: use it for speech frames if update value is less than previous
   1845     // the formula has been rewritten into:
   1846     // noiseUpdate = noisePrev[i] + (1 - gammaNoise) * nonSpeechProb * (magn[i] - noisePrev[i])
   1847 
   1848     if (postShifts < 0) {
   1849       tmpU32no2 = magnU16[i] >> -postShifts;  // Q(prevQNoise)
   1850     } else {
   1851       tmpU32no2 = (uint32_t)magnU16[i] << postShifts;  // Q(prevQNoise)
   1852     }
   1853     if (prevNoiseU16[i] > tmpU32no2) {
   1854       sign = -1;
   1855       tmpU32no1 = prevNoiseU16[i] - tmpU32no2;
   1856     } else {
   1857       sign = 1;
   1858       tmpU32no1 = tmpU32no2 - prevNoiseU16[i];
   1859     }
   1860     noiseUpdateU32 = inst->prevNoiseU32[i]; // Q(prevQNoise+11)
   1861     tmpU32no3 = 0;
   1862     if ((tmpU32no1) && (nonSpeechProbFinal[i])) {
   1863       // This value will be used later, if gammaNoise changes
   1864       tmpU32no3 = WEBRTC_SPL_UMUL_32_16(tmpU32no1, nonSpeechProbFinal[i]); // Q(prevQNoise+8)
   1865       if (0x7c000000 & tmpU32no3) {
   1866         // Shifting required before multiplication
   1867         tmpU32no2 = (tmpU32no3 >> 5) * gammaNoise;  // Q(prevQNoise+11)
   1868       } else {
   1869         // We can do shifting after multiplication
   1870         tmpU32no2 = (tmpU32no3 * gammaNoise) >> 5;  // Q(prevQNoise+11)
   1871       }
   1872       if (sign > 0) {
   1873         noiseUpdateU32 += tmpU32no2; // Q(prevQNoise+11)
   1874       } else {
   1875         // This operation is safe. We can never get wrap around, since worst
   1876         // case scenario means magnU16 = 0
   1877         noiseUpdateU32 -= tmpU32no2; // Q(prevQNoise+11)
   1878       }
   1879     }
   1880 
   1881     //increase gamma (i.e., less noise update) for frame likely to be speech
   1882     prevGammaNoise = gammaNoise;
   1883     gammaNoise = NOISE_UPDATE_Q8;
   1884     //time-constant based on speech/noise state
   1885     //increase gamma (i.e., less noise update) for frames likely to be speech
   1886     if (nonSpeechProbFinal[i] < ONE_MINUS_PROB_RANGE_Q8) {
   1887       gammaNoise = GAMMA_NOISE_TRANS_AND_SPEECH_Q8;
   1888     }
   1889 
   1890     if (prevGammaNoise != gammaNoise) {
   1891       // new noise update
   1892       // this line is the same as above, only that the result is stored in a different variable and the gammaNoise
   1893       // has changed
   1894       //
   1895       // noiseUpdate = noisePrev[i] + (1 - gammaNoise) * nonSpeechProb * (magn[i] - noisePrev[i])
   1896 
   1897       if (0x7c000000 & tmpU32no3) {
   1898         // Shifting required before multiplication
   1899         tmpU32no2 = (tmpU32no3 >> 5) * gammaNoise;  // Q(prevQNoise+11)
   1900       } else {
   1901         // We can do shifting after multiplication
   1902         tmpU32no2 = (tmpU32no3 * gammaNoise) >> 5;  // Q(prevQNoise+11)
   1903       }
   1904       if (sign > 0) {
   1905         tmpU32no1 = inst->prevNoiseU32[i] + tmpU32no2; // Q(prevQNoise+11)
   1906       } else {
   1907         tmpU32no1 = inst->prevNoiseU32[i] - tmpU32no2; // Q(prevQNoise+11)
   1908       }
   1909       if (noiseUpdateU32 > tmpU32no1) {
   1910         noiseUpdateU32 = tmpU32no1; // Q(prevQNoise+11)
   1911       }
   1912     }
   1913     noiseU32[i] = noiseUpdateU32; // Q(prevQNoise+11)
   1914     if (noiseUpdateU32 > maxNoiseU32) {
   1915       maxNoiseU32 = noiseUpdateU32;
   1916     }
   1917 
   1918     // conservative noise update
   1919     // // original FLOAT code
   1920     // if (prob_speech < PROB_RANGE) {
   1921     // inst->avgMagnPause[i] = inst->avgMagnPause[i] + (1.0 - gamma_pause)*(magn[i] - inst->avgMagnPause[i]);
   1922     // }
   1923 
   1924     tmp32no2 = WEBRTC_SPL_SHIFT_W32(inst->avgMagnPause[i], -nShifts);
   1925     if (nonSpeechProbFinal[i] > ONE_MINUS_PROB_RANGE_Q8) {
   1926       if (nShifts < 0) {
   1927         tmp32no1 = (int32_t)magnU16[i] - tmp32no2; // Q(qMagn)
   1928         tmp32no1 *= ONE_MINUS_GAMMA_PAUSE_Q8;  // Q(8+prevQMagn+nShifts)
   1929         tmp32no1 = (tmp32no1 + 128) >> 8;  // Q(qMagn).
   1930       } else {
   1931         // In Q(qMagn+nShifts)
   1932         tmp32no1 = ((int32_t)magnU16[i] << nShifts) - inst->avgMagnPause[i];
   1933         tmp32no1 *= ONE_MINUS_GAMMA_PAUSE_Q8;  // Q(8+prevQMagn+nShifts)
   1934         tmp32no1 = (tmp32no1 + (128 << nShifts)) >> (8 + nShifts);  // Q(qMagn).
   1935       }
   1936       tmp32no2 += tmp32no1; // Q(qMagn)
   1937     }
   1938     inst->avgMagnPause[i] = tmp32no2;
   1939   }  // end of frequency loop
   1940 
   1941   norm32no1 = WebRtcSpl_NormU32(maxNoiseU32);
   1942   qNoise = inst->prevQNoise + norm32no1 - 5;
   1943   // done with step 2: noise update
   1944 
   1945   // STEP 3: compute dd update of prior snr and post snr based on new noise estimate
   1946   nShifts = inst->prevQNoise + 11 - qMagn;
   1947   for (i = 0; i < inst->magnLen; i++) {
   1948     // FLOAT code
   1949     // // post and prior SNR
   1950     // curNearSnr = 0.0;
   1951     // if (magn[i] > noise[i])
   1952     // {
   1953     // curNearSnr = magn[i] / (noise[i] + 0.0001) - 1.0;
   1954     // }
   1955     // // DD estimate is sum of two terms: current estimate and previous estimate
   1956     // // directed decision update of snrPrior
   1957     // snrPrior = DD_PR_SNR * prevNearSnr[i] + (1.0 - DD_PR_SNR) * curNearSnr;
   1958     // // gain filter
   1959     // tmpFloat1 = inst->overdrive + snrPrior;
   1960     // tmpFloat2 = snrPrior / tmpFloat1;
   1961     // theFilter[i] = tmpFloat2;
   1962 
   1963     // calculate curNearSnr again, this is necessary because a new noise estimate has been made since then. for the original
   1964     curNearSnr = 0; // Q11
   1965     if (nShifts < 0) {
   1966       // This case is equivalent with magn < noise which implies curNearSnr = 0;
   1967       tmpMagnU32 = (uint32_t)magnU16[i]; // Q(qMagn)
   1968       tmpNoiseU32 = noiseU32[i] << -nShifts;  // Q(qMagn)
   1969     } else if (nShifts > 17) {
   1970       tmpMagnU32 = (uint32_t)magnU16[i] << 17;  // Q(qMagn+17)
   1971       tmpNoiseU32 = noiseU32[i] >> (nShifts - 17);  // Q(qMagn+17)
   1972     } else {
   1973       tmpMagnU32 = (uint32_t)magnU16[i] << nShifts;  // Q(qNoise_prev+11)
   1974       tmpNoiseU32 = noiseU32[i]; // Q(qNoise_prev+11)
   1975     }
   1976     if (tmpMagnU32 > tmpNoiseU32) {
   1977       tmpU32no1 = tmpMagnU32 - tmpNoiseU32; // Q(qCur)
   1978       norm32no2 = WEBRTC_SPL_MIN(11, WebRtcSpl_NormU32(tmpU32no1));
   1979       tmpU32no1 <<= norm32no2;  // Q(qCur+norm32no2)
   1980       tmpU32no2 = tmpNoiseU32 >> (11 - norm32no2);  // Q(qCur+norm32no2-11)
   1981       if (tmpU32no2 > 0) {
   1982         tmpU32no1 /= tmpU32no2;  // Q11
   1983       }
   1984       curNearSnr = WEBRTC_SPL_MIN(satMax, tmpU32no1); // Q11
   1985     }
   1986 
   1987     //directed decision update of priorSnr
   1988     // FLOAT
   1989     // priorSnr = DD_PR_SNR * prevNearSnr + (1.0-DD_PR_SNR) * curNearSnr;
   1990 
   1991     tmpU32no1 = WEBRTC_SPL_UMUL_32_16(prevNearSnr[i], DD_PR_SNR_Q11); // Q22
   1992     tmpU32no2 = WEBRTC_SPL_UMUL_32_16(curNearSnr, ONE_MINUS_DD_PR_SNR_Q11); // Q22
   1993     priorSnr = tmpU32no1 + tmpU32no2; // Q22
   1994 
   1995     //gain filter
   1996     tmpU32no1 = inst->overdrive + ((priorSnr + 8192) >> 14);  // Q8
   1997     assert(inst->overdrive > 0);
   1998     tmpU16no1 = (priorSnr + tmpU32no1 / 2) / tmpU32no1;  // Q14
   1999     inst->noiseSupFilter[i] = WEBRTC_SPL_SAT(16384, tmpU16no1, inst->denoiseBound); // 16384 = Q14(1.0) // Q14
   2000 
   2001     // Weight in the parametric Wiener filter during startup
   2002     if (inst->blockIndex < END_STARTUP_SHORT) {
   2003       // Weight the two suppression filters
   2004       tmpU32no1 = inst->noiseSupFilter[i] * inst->blockIndex;
   2005       tmpU32no2 = noiseSupFilterTmp[i] *
   2006           (END_STARTUP_SHORT - inst->blockIndex);
   2007       tmpU32no1 += tmpU32no2;
   2008       inst->noiseSupFilter[i] = (uint16_t)WebRtcSpl_DivU32U16(tmpU32no1,
   2009                                                                     END_STARTUP_SHORT);
   2010     }
   2011   }  // end of loop over frequencies
   2012   //done with step3
   2013 
   2014   // save noise and magnitude spectrum for next frame
   2015   inst->prevQNoise = qNoise;
   2016   inst->prevQMagn = qMagn;
   2017   if (norm32no1 > 5) {
   2018     for (i = 0; i < inst->magnLen; i++) {
   2019       inst->prevNoiseU32[i] = noiseU32[i] << (norm32no1 - 5);  // Q(qNoise+11)
   2020       inst->prevMagnU16[i] = magnU16[i]; // Q(qMagn)
   2021     }
   2022   } else {
   2023     for (i = 0; i < inst->magnLen; i++) {
   2024       inst->prevNoiseU32[i] = noiseU32[i] >> (5 - norm32no1);  // Q(qNoise+11)
   2025       inst->prevMagnU16[i] = magnU16[i]; // Q(qMagn)
   2026     }
   2027   }
   2028 
   2029   WebRtcNsx_DataSynthesis(inst, outFrame[0]);
   2030 #ifdef NS_FILEDEBUG
   2031   if (fwrite(outframe, sizeof(short),
   2032              inst->blockLen10ms, inst->outfile) != inst->blockLen10ms) {
   2033     assert(false);
   2034   }
   2035 #endif
   2036 
   2037   //for H band:
   2038   // only update data buffer, then apply time-domain gain is applied derived from L band
   2039   if (num_bands > 1) {
   2040     // update analysis buffer for H band
   2041     // append new data to buffer FX
   2042     for (i = 0; i < num_high_bands; ++i) {
   2043       memcpy(inst->dataBufHBFX[i], inst->dataBufHBFX[i] + inst->blockLen10ms,
   2044           (inst->anaLen - inst->blockLen10ms) * sizeof(*inst->dataBufHBFX[i]));
   2045       memcpy(inst->dataBufHBFX[i] + inst->anaLen - inst->blockLen10ms,
   2046           speechFrameHB[i], inst->blockLen10ms * sizeof(*inst->dataBufHBFX[i]));
   2047     }
   2048     // range for averaging low band quantities for H band gain
   2049 
   2050     gainTimeDomainHB = 16384; // 16384 = Q14(1.0)
   2051     //average speech prob from low band
   2052     //average filter gain from low band
   2053     //avg over second half (i.e., 4->8kHz) of freq. spectrum
   2054     tmpU32no1 = 0; // Q12
   2055     tmpU16no1 = 0; // Q8
   2056     for (i = inst->anaLen2 - (inst->anaLen2 >> 2); i < inst->anaLen2; i++) {
   2057       tmpU16no1 += nonSpeechProbFinal[i]; // Q8
   2058       tmpU32no1 += (uint32_t)(inst->noiseSupFilter[i]); // Q14
   2059     }
   2060     assert(inst->stages >= 7);
   2061     avgProbSpeechHB = (4096 - (tmpU16no1 >> (inst->stages - 7)));  // Q12
   2062     avgFilterGainHB = (int16_t)(tmpU32no1 >> (inst->stages - 3));  // Q14
   2063 
   2064     // // original FLOAT code
   2065     // // gain based on speech probability:
   2066     // avg_prob_speech_tt=(float)2.0*avg_prob_speech-(float)1.0;
   2067     // gain_mod=(float)0.5*((float)1.0+(float)tanh(avg_prob_speech_tt)); // between 0 and 1
   2068 
   2069     // gain based on speech probability:
   2070     // original expression: "0.5 * (1 + tanh(2x-1))"
   2071     // avgProbSpeechHB has been anyway saturated to a value between 0 and 1 so the other cases don't have to be dealt with
   2072     // avgProbSpeechHB and gainModHB are in Q12, 3607 = Q12(0.880615234375) which is a zero point of
   2073     // |0.5 * (1 + tanh(2x-1)) - x| - |0.5 * (1 + tanh(2x-1)) - 0.880615234375| meaning that from that point the error of approximating
   2074     // the expression with f(x) = x would be greater than the error of approximating the expression with f(x) = 0.880615234375
   2075     // error: "|0.5 * (1 + tanh(2x-1)) - x| from x=0 to 0.880615234375" -> http://www.wolframalpha.com/input/?i=|0.5+*+(1+%2B+tanh(2x-1))+-+x|+from+x%3D0+to+0.880615234375
   2076     // and:  "|0.5 * (1 + tanh(2x-1)) - 0.880615234375| from x=0.880615234375 to 1" -> http://www.wolframalpha.com/input/?i=+|0.5+*+(1+%2B+tanh(2x-1))+-+0.880615234375|+from+x%3D0.880615234375+to+1
   2077     gainModHB = WEBRTC_SPL_MIN(avgProbSpeechHB, 3607);
   2078 
   2079     // // original FLOAT code
   2080     // //combine gain with low band gain
   2081     // if (avg_prob_speech < (float)0.5) {
   2082     // gain_time_domain_HB=(float)0.5*gain_mod+(float)0.5*avg_filter_gain;
   2083     // }
   2084     // else {
   2085     // gain_time_domain_HB=(float)0.25*gain_mod+(float)0.75*avg_filter_gain;
   2086     // }
   2087 
   2088 
   2089     //combine gain with low band gain
   2090     if (avgProbSpeechHB < 2048) {
   2091       // 2048 = Q12(0.5)
   2092       // the next two lines in float are  "gain_time_domain = 0.5 * gain_mod + 0.5 * avg_filter_gain"; Q2(0.5) = 2 equals one left shift
   2093       gainTimeDomainHB = (gainModHB << 1) + (avgFilterGainHB >> 1); // Q14
   2094     } else {
   2095       // "gain_time_domain = 0.25 * gain_mod + 0.75 * agv_filter_gain;"
   2096       gainTimeDomainHB = (int16_t)((3 * avgFilterGainHB) >> 2);  // 3 = Q2(0.75)
   2097       gainTimeDomainHB += gainModHB; // Q14
   2098     }
   2099     //make sure gain is within flooring range
   2100     gainTimeDomainHB
   2101       = WEBRTC_SPL_SAT(16384, gainTimeDomainHB, (int16_t)(inst->denoiseBound)); // 16384 = Q14(1.0)
   2102 
   2103 
   2104     //apply gain
   2105     for (i = 0; i < num_high_bands; ++i) {
   2106       for (j = 0; j < inst->blockLen10ms; j++) {
   2107         outFrameHB[i][j] = (int16_t)((gainTimeDomainHB *
   2108             inst->dataBufHBFX[i][j]) >> 14);  // Q0
   2109       }
   2110     }
   2111   }  // end of H band gain computation
   2112 }
   2113