/*
 *  Copyright (c) 2012 The WebRTC project authors. All Rights Reserved.
 *
 *  Use of this source code is governed by a BSD-style license
 *  that can be found in the LICENSE file in the root of the source
 *  tree. An additional intellectual property rights grant can be found
 *  in the file PATENTS.  All contributing project authors may
 *  be found in the AUTHORS file in the root of the source tree.
 */

#include "webrtc/modules/audio_processing/ns/noise_suppression_x.h"

#include <assert.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

#include "webrtc/common_audio/signal_processing/include/real_fft.h"
#include "webrtc/modules/audio_processing/ns/nsx_core.h"
#include "webrtc/system_wrappers/include/cpu_features_wrapper.h"

#if defined(WEBRTC_HAS_NEON)
/* Tables are defined in ARM assembly files. */
extern const int16_t WebRtcNsx_kLogTable[9];
extern const int16_t WebRtcNsx_kCounterDiv[201];
extern const int16_t WebRtcNsx_kLogTableFrac[256];
#else
static const int16_t WebRtcNsx_kLogTable[9] = {
  0, 177, 355, 532, 710, 887, 1065, 1242, 1420
};

static const int16_t WebRtcNsx_kCounterDiv[201] = {
  32767, 16384, 10923, 8192, 6554, 5461, 4681, 4096, 3641, 3277, 2979, 2731,
  2521, 2341, 2185, 2048, 1928, 1820, 1725, 1638, 1560, 1489, 1425, 1365, 1311,
  1260, 1214, 1170, 1130, 1092, 1057, 1024, 993, 964, 936, 910, 886, 862, 840,
  819, 799, 780, 762, 745, 728, 712, 697, 683, 669, 655, 643, 630, 618, 607,
  596, 585, 575, 565, 555, 546, 537, 529, 520, 512, 504, 496, 489, 482, 475,
  468, 462, 455, 449, 443, 437, 431, 426, 420, 415, 410, 405, 400, 395, 390,
  386, 381, 377, 372, 368, 364, 360, 356, 352, 349, 345, 341, 338, 334, 331,
  328, 324, 321, 318, 315, 312, 309, 306, 303, 301, 298, 295, 293, 290, 287,
  285, 282, 280, 278, 275, 273, 271, 269, 266, 264, 262, 260, 258, 256, 254,
  252, 250, 248, 246, 245, 243, 241, 239, 237, 236, 234, 232, 231, 229, 228,
  226, 224, 223, 221, 220, 218, 217, 216, 214, 213, 211, 210, 209, 207, 206,
  205, 204, 202, 201, 200, 199, 197, 196, 195, 194, 193, 192, 191, 189, 188,
  187, 186, 185, 184, 183, 182, 181, 180, 179, 178, 177, 176, 175, 174, 173,
  172, 172, 171, 170, 169, 168, 167, 166, 165, 165, 164, 163
};

static const int16_t WebRtcNsx_kLogTableFrac[256] = {
  0,   1,   3,   4,   6,   7,   9,  10,  11,  13,  14,  16,  17,  18,  20,  21,
  22,  24,  25,  26,  28,  29,  30,  32,  33,  34,  36,  37,  38,  40,  41,  42,
  44,  45,  46,  47,  49,  50,  51,  52,  54,  55,  56,  57,  59,  60,  61,  62,
  63,  65,  66,  67,  68,  69,  71,  72,  73,  74,  75,  77,  78,  79,  80,  81,
  82,  84,  85,  86,  87,  88,  89,  90,  92,  93,  94,  95,  96,  97,  98,  99,
  100, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 116,
  117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131,
  132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146,
  147, 148, 149, 150, 151, 152, 153, 154, 155, 155, 156, 157, 158, 159, 160,
  161, 162, 163, 164, 165, 166, 167, 168, 169, 169, 170, 171, 172, 173, 174,
  175, 176, 177, 178, 178, 179, 180, 181, 182, 183, 184, 185, 185, 186, 187,
  188, 189, 190, 191, 192, 192, 193, 194, 195, 196, 197, 198, 198, 199, 200,
  201, 202, 203, 203, 204, 205, 206, 207, 208, 208, 209, 210, 211, 212, 212,
  213, 214, 215, 216, 216, 217, 218, 219, 220, 220, 221, 222, 223, 224, 224,
  225, 226, 227, 228, 228, 229, 230, 231, 231, 232, 233, 234, 234, 235, 236,
  237, 238, 238, 239, 240, 241, 241, 242, 243, 244, 244, 245, 246, 247, 247,
  248, 249, 249, 250, 251, 252, 252, 253, 254, 255, 255
};
#endif  // WEBRTC_HAS_NEON

// Skip first frequency bins during estimation. (0 <= value < 64)
static const size_t kStartBand = 5;

// hybrib Hanning & flat window
static const int16_t kBlocks80w128x[128] = {
  0,    536,   1072,   1606,   2139,   2669,   3196,   3720,   4240,   4756,   5266,
  5771,   6270,   6762,   7246,   7723,   8192,   8652,   9102,   9543,   9974,  10394,
  10803,  11200,  11585,  11958,  12318,  12665,  12998,  13318,  13623,  13913,  14189,
  14449,  14694,  14924,  15137,  15334,  15515,  15679,  15826,  15956,  16069,  16165,
  16244,  16305,  16349,  16375,  16384,  16384,  16384,  16384,  16384,  16384,  16384,
  16384,  16384,  16384,  16384,  16384,  16384,  16384,  16384,  16384,  16384,  16384,
  16384,  16384,  16384,  16384,  16384,  16384,  16384,  16384,  16384,  16384,  16384,
  16384,  16384,  16384,  16384,  16375,  16349,  16305,  16244,  16165,  16069,  15956,
  15826,  15679,  15515,  15334,  15137,  14924,  14694,  14449,  14189,  13913,  13623,
  13318,  12998,  12665,  12318,  11958,  11585,  11200,  10803,  10394,   9974,   9543,
  9102,   8652,   8192,   7723,   7246,   6762,   6270,   5771,   5266,   4756,   4240,
  3720,   3196,   2669,   2139,   1606,   1072,    536
};

// hybrib Hanning & flat window
static const int16_t kBlocks160w256x[256] = {
  0,   268,   536,   804,  1072,  1339,  1606,  1872,
  2139,  2404,  2669,  2933,  3196,  3459,  3720,  3981,
  4240,  4499,  4756,  5012,  5266,  5520,  5771,  6021,
  6270,  6517,  6762,  7005,  7246,  7486,  7723,  7959,
  8192,  8423,  8652,  8878,  9102,  9324,  9543,  9760,
  9974, 10185, 10394, 10600, 10803, 11003, 11200, 11394,
  11585, 11773, 11958, 12140, 12318, 12493, 12665, 12833,
  12998, 13160, 13318, 13472, 13623, 13770, 13913, 14053,
  14189, 14321, 14449, 14574, 14694, 14811, 14924, 15032,
  15137, 15237, 15334, 15426, 15515, 15599, 15679, 15754,
  15826, 15893, 15956, 16015, 16069, 16119, 16165, 16207,
  16244, 16277, 16305, 16329, 16349, 16364, 16375, 16382,
  16384, 16384, 16384, 16384, 16384, 16384, 16384, 16384,
  16384, 16384, 16384, 16384, 16384, 16384, 16384, 16384,
  16384, 16384, 16384, 16384, 16384, 16384, 16384, 16384,
  16384, 16384, 16384, 16384, 16384, 16384, 16384, 16384,
  16384, 16384, 16384, 16384, 16384, 16384, 16384, 16384,
  16384, 16384, 16384, 16384, 16384, 16384, 16384, 16384,
  16384, 16384, 16384, 16384, 16384, 16384, 16384, 16384,
  16384, 16384, 16384, 16384, 16384, 16384, 16384, 16384,
  16384, 16382, 16375, 16364, 16349, 16329, 16305, 16277,
  16244, 16207, 16165, 16119, 16069, 16015, 15956, 15893,
  15826, 15754, 15679, 15599, 15515, 15426, 15334, 15237,
  15137, 15032, 14924, 14811, 14694, 14574, 14449, 14321,
  14189, 14053, 13913, 13770, 13623, 13472, 13318, 13160,
  12998, 12833, 12665, 12493, 12318, 12140, 11958, 11773,
  11585, 11394, 11200, 11003, 10803, 10600, 10394, 10185,
  9974,  9760,  9543,  9324,  9102,  8878,  8652,  8423,
  8192,  7959,  7723,  7486,  7246,  7005,  6762,  6517,
  6270,  6021,  5771,  5520,  5266,  5012,  4756,  4499,
  4240,  3981,  3720,  3459,  3196,  2933,  2669,  2404,
  2139,  1872,  1606,  1339,  1072,   804,   536,   268
};

// Gain factor1 table: Input value in Q8 and output value in Q13
// original floating point code
//  if (gain > blim) {
//    factor1 = 1.0 + 1.3 * (gain - blim);
//    if (gain * factor1 > 1.0) {
//      factor1 = 1.0 / gain;
//    }
//  } else {
//    factor1 = 1.0;
//  }
static const int16_t kFactor1Table[257] = {
  8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8233, 8274, 8315, 8355, 8396, 8436, 8475, 8515, 8554, 8592, 8631, 8669,
  8707, 8745, 8783, 8820, 8857, 8894, 8931, 8967, 9003, 9039, 9075, 9111, 9146, 9181,
  9216, 9251, 9286, 9320, 9354, 9388, 9422, 9456, 9489, 9523, 9556, 9589, 9622, 9655,
  9687, 9719, 9752, 9784, 9816, 9848, 9879, 9911, 9942, 9973, 10004, 10035, 10066,
  10097, 10128, 10158, 10188, 10218, 10249, 10279, 10308, 10338, 10368, 10397, 10426,
  10456, 10485, 10514, 10543, 10572, 10600, 10629, 10657, 10686, 10714, 10742, 10770,
  10798, 10826, 10854, 10882, 10847, 10810, 10774, 10737, 10701, 10666, 10631, 10596,
  10562, 10527, 10494, 10460, 10427, 10394, 10362, 10329, 10297, 10266, 10235, 10203,
  10173, 10142, 10112, 10082, 10052, 10023, 9994, 9965, 9936, 9908, 9879, 9851, 9824,
  9796, 9769, 9742, 9715, 9689, 9662, 9636, 9610, 9584, 9559, 9534, 9508, 9484, 9459,
  9434, 9410, 9386, 9362, 9338, 9314, 9291, 9268, 9245, 9222, 9199, 9176, 9154, 9132,
  9110, 9088, 9066, 9044, 9023, 9002, 8980, 8959, 8939, 8918, 8897, 8877, 8857, 8836,
  8816, 8796, 8777, 8757, 8738, 8718, 8699, 8680, 8661, 8642, 8623, 8605, 8586, 8568,
  8550, 8532, 8514, 8496, 8478, 8460, 8443, 8425, 8408, 8391, 8373, 8356, 8339, 8323,
  8306, 8289, 8273, 8256, 8240, 8224, 8208, 8192
};

// For Factor2 tables
// original floating point code
// if (gain > blim) {
//   factor2 = 1.0;
// } else {
//   factor2 = 1.0 - 0.3 * (blim - gain);
//   if (gain <= inst->denoiseBound) {
//     factor2 = 1.0 - 0.3 * (blim - inst->denoiseBound);
//   }
// }
//
// Gain factor table: Input value in Q8 and output value in Q13
static const int16_t kFactor2Aggressiveness1[257] = {
  7577, 7577, 7577, 7577, 7577, 7577,
  7577, 7577, 7577, 7577, 7577, 7577, 7577, 7577, 7577, 7577, 7577, 7596, 7614, 7632,
  7650, 7667, 7683, 7699, 7715, 7731, 7746, 7761, 7775, 7790, 7804, 7818, 7832, 7845,
  7858, 7871, 7884, 7897, 7910, 7922, 7934, 7946, 7958, 7970, 7982, 7993, 8004, 8016,
  8027, 8038, 8049, 8060, 8070, 8081, 8091, 8102, 8112, 8122, 8132, 8143, 8152, 8162,
  8172, 8182, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192
};

// Gain factor table: Input value in Q8 and output value in Q13
static const int16_t kFactor2Aggressiveness2[257] = {
  7270, 7270, 7270, 7270, 7270, 7306,
  7339, 7369, 7397, 7424, 7448, 7472, 7495, 7517, 7537, 7558, 7577, 7596, 7614, 7632,
  7650, 7667, 7683, 7699, 7715, 7731, 7746, 7761, 7775, 7790, 7804, 7818, 7832, 7845,
  7858, 7871, 7884, 7897, 7910, 7922, 7934, 7946, 7958, 7970, 7982, 7993, 8004, 8016,
  8027, 8038, 8049, 8060, 8070, 8081, 8091, 8102, 8112, 8122, 8132, 8143, 8152, 8162,
  8172, 8182, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192
};

// Gain factor table: Input value in Q8 and output value in Q13
static const int16_t kFactor2Aggressiveness3[257] = {
  7184, 7184, 7184, 7229, 7270, 7306,
  7339, 7369, 7397, 7424, 7448, 7472, 7495, 7517, 7537, 7558, 7577, 7596, 7614, 7632,
  7650, 7667, 7683, 7699, 7715, 7731, 7746, 7761, 7775, 7790, 7804, 7818, 7832, 7845,
  7858, 7871, 7884, 7897, 7910, 7922, 7934, 7946, 7958, 7970, 7982, 7993, 8004, 8016,
  8027, 8038, 8049, 8060, 8070, 8081, 8091, 8102, 8112, 8122, 8132, 8143, 8152, 8162,
  8172, 8182, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
  8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192
};

// sum of log2(i) from table index to inst->anaLen2 in Q5
// Note that the first table value is invalid, since log2(0) = -infinity
static const int16_t kSumLogIndex[66] = {
  0,  22917,  22917,  22885,  22834,  22770,  22696,  22613,
  22524,  22428,  22326,  22220,  22109,  21994,  21876,  21754,
  21629,  21501,  21370,  21237,  21101,  20963,  20822,  20679,
  20535,  20388,  20239,  20089,  19937,  19783,  19628,  19470,
  19312,  19152,  18991,  18828,  18664,  18498,  18331,  18164,
  17994,  17824,  17653,  17480,  17306,  17132,  16956,  16779,
  16602,  16423,  16243,  16063,  15881,  15699,  15515,  15331,
  15146,  14960,  14774,  14586,  14398,  14209,  14019,  13829,
  13637,  13445
};

// sum of log2(i)^2 from table index to inst->anaLen2 in Q2
// Note that the first table value is invalid, since log2(0) = -infinity
static const int16_t kSumSquareLogIndex[66] = {
  0,  16959,  16959,  16955,  16945,  16929,  16908,  16881,
  16850,  16814,  16773,  16729,  16681,  16630,  16575,  16517,
  16456,  16392,  16325,  16256,  16184,  16109,  16032,  15952,
  15870,  15786,  15700,  15612,  15521,  15429,  15334,  15238,
  15140,  15040,  14938,  14834,  14729,  14622,  14514,  14404,
  14292,  14179,  14064,  13947,  13830,  13710,  13590,  13468,
  13344,  13220,  13094,  12966,  12837,  12707,  12576,  12444,
  12310,  12175,  12039,  11902,  11763,  11624,  11483,  11341,
  11198,  11054
};

// log2(table index) in Q12
// Note that the first table value is invalid, since log2(0) = -infinity
static const int16_t kLogIndex[129] = {
  0,      0,   4096,   6492,   8192,   9511,  10588,  11499,
  12288,  12984,  13607,  14170,  14684,  15157,  15595,  16003,
  16384,  16742,  17080,  17400,  17703,  17991,  18266,  18529,
  18780,  19021,  19253,  19476,  19691,  19898,  20099,  20292,
  20480,  20662,  20838,  21010,  21176,  21338,  21496,  21649,
  21799,  21945,  22087,  22226,  22362,  22495,  22625,  22752,
  22876,  22998,  23117,  23234,  23349,  23462,  23572,  23680,
  23787,  23892,  23994,  24095,  24195,  24292,  24388,  24483,
  24576,  24668,  24758,  24847,  24934,  25021,  25106,  25189,
  25272,  25354,  25434,  25513,  25592,  25669,  25745,  25820,
  25895,  25968,  26041,  26112,  26183,  26253,  26322,  26390,
  26458,  26525,  26591,  26656,  26721,  26784,  26848,  26910,
  26972,  27033,  27094,  27154,  27213,  27272,  27330,  27388,
  27445,  27502,  27558,  27613,  27668,  27722,  27776,  27830,
  27883,  27935,  27988,  28039,  28090,  28141,  28191,  28241,
  28291,  28340,  28388,  28437,  28484,  28532,  28579,  28626,
  28672
};

// determinant of estimation matrix in Q0 corresponding to the log2 tables above
// Note that the first table value is invalid, since log2(0) = -infinity
static const int16_t kDeterminantEstMatrix[66] = {
  0,  29814,  25574,  22640,  20351,  18469,  16873,  15491,
  14277,  13199,  12233,  11362,  10571,   9851,   9192,   8587,
  8030,   7515,   7038,   6596,   6186,   5804,   5448,   5115,
  4805,   4514,   4242,   3988,   3749,   3524,   3314,   3116,
  2930,   2755,   2590,   2435,   2289,   2152,   2022,   1900,
  1785,   1677,   1575,   1478,   1388,   1302,   1221,   1145,
  1073,   1005,    942,    881,    825,    771,    721,    674,
  629,    587,    547,    510,    475,    442,    411,    382,
  355,    330
};

// Update the noise estimation information.
static void UpdateNoiseEstimate(NoiseSuppressionFixedC* inst, int offset) {
  int32_t tmp32no1 = 0;
  int32_t tmp32no2 = 0;
  int16_t tmp16 = 0;
  const int16_t kExp2Const = 11819; // Q13

  size_t i = 0;

  tmp16 = WebRtcSpl_MaxValueW16(inst->noiseEstLogQuantile + offset,
                                   inst->magnLen);
  // Guarantee a Q-domain as high as possible and still fit in int16
  inst->qNoise = 14 - (int) WEBRTC_SPL_MUL_16_16_RSFT_WITH_ROUND(
                   kExp2Const, tmp16, 21);
  for (i = 0; i < inst->magnLen; i++) {
    // inst->quantile[i]=exp(inst->lquantile[offset+i]);
    // in Q21
    tmp32no2 = kExp2Const * inst->noiseEstLogQuantile[offset + i];
    tmp32no1 = (0x00200000 | (tmp32no2 & 0x001FFFFF)); // 2^21 + frac
    tmp16 = (int16_t)(tmp32no2 >> 21);
    tmp16 -= 21;// shift 21 to get result in Q0
    tmp16 += (int16_t) inst->qNoise; //shift to get result in Q(qNoise)
    if (tmp16 < 0) {
      tmp32no1 >>= -tmp16;
    } else {
      tmp32no1 <<= tmp16;
    }
    inst->noiseEstQuantile[i] = WebRtcSpl_SatW32ToW16(tmp32no1);
  }
}

// Noise Estimation
static void NoiseEstimationC(NoiseSuppressionFixedC* inst,
                             uint16_t* magn,
                             uint32_t* noise,
                             int16_t* q_noise) {
  int16_t lmagn[HALF_ANAL_BLOCKL], counter, countDiv;
  int16_t countProd, delta, zeros, frac;
  int16_t log2, tabind, logval, tmp16, tmp16no1, tmp16no2;
  const int16_t log2_const = 22713; // Q15
  const int16_t width_factor = 21845;

  size_t i, s, offset;

  tabind = inst->stages - inst->normData;
  assert(tabind < 9);
  assert(tabind > -9);
  if (tabind < 0) {
    logval = -WebRtcNsx_kLogTable[-tabind];
  } else {
    logval = WebRtcNsx_kLogTable[tabind];
  }

  // lmagn(i)=log(magn(i))=log(2)*log2(magn(i))
  // magn is in Q(-stages), and the real lmagn values are:
  // real_lmagn(i)=log(magn(i)*2^stages)=log(magn(i))+log(2^stages)
  // lmagn in Q8
  for (i = 0; i < inst->magnLen; i++) {
    if (magn[i]) {
      zeros = WebRtcSpl_NormU32((uint32_t)magn[i]);
      frac = (int16_t)((((uint32_t)magn[i] << zeros)
                              & 0x7FFFFFFF) >> 23);
      // log2(magn(i))
      assert(frac < 256);
      log2 = (int16_t)(((31 - zeros) << 8)
                             + WebRtcNsx_kLogTableFrac[frac]);
      // log2(magn(i))*log(2)
      lmagn[i] = (int16_t)((log2 * log2_const) >> 15);
      // + log(2^stages)
      lmagn[i] += logval;
    } else {
      lmagn[i] = logval;//0;
    }
  }

  // loop over simultaneous estimates
  for (s = 0; s < SIMULT; s++) {
    offset = s * inst->magnLen;

    // Get counter values from state
    counter = inst->noiseEstCounter[s];
    assert(counter < 201);
    countDiv = WebRtcNsx_kCounterDiv[counter];
    countProd = (int16_t)(counter * countDiv);

    // quant_est(...)
    for (i = 0; i < inst->magnLen; i++) {
      // compute delta
      if (inst->noiseEstDensity[offset + i] > 512) {
        // Get the value for delta by shifting intead of dividing.
        int factor = WebRtcSpl_NormW16(inst->noiseEstDensity[offset + i]);
        delta = (int16_t)(FACTOR_Q16 >> (14 - factor));
      } else {
        delta = FACTOR_Q7;
        if (inst->blockIndex < END_STARTUP_LONG) {
          // Smaller step size during startup. This prevents from using
          // unrealistic values causing overflow.
          delta = FACTOR_Q7_STARTUP;
        }
      }

      // update log quantile estimate
      tmp16 = (int16_t)((delta * countDiv) >> 14);
      if (lmagn[i] > inst->noiseEstLogQuantile[offset + i]) {
        // +=QUANTILE*delta/(inst->counter[s]+1) QUANTILE=0.25, =1 in Q2
        // CounterDiv=1/(inst->counter[s]+1) in Q15
        tmp16 += 2;
        inst->noiseEstLogQuantile[offset + i] += tmp16 / 4;
      } else {
        tmp16 += 1;
        // *(1-QUANTILE), in Q2 QUANTILE=0.25, 1-0.25=0.75=3 in Q2
        // TODO(bjornv): investigate why we need to truncate twice.
        tmp16no2 = (int16_t)((tmp16 / 2) * 3 / 2);
        inst->noiseEstLogQuantile[offset + i] -= tmp16no2;
        if (inst->noiseEstLogQuantile[offset + i] < logval) {
          // This is the smallest fixed point representation we can
          // have, hence we limit the output.
          inst->noiseEstLogQuantile[offset + i] = logval;
        }
      }

      // update density estimate
      if (WEBRTC_SPL_ABS_W16(lmagn[i] - inst->noiseEstLogQuantile[offset + i])
          < WIDTH_Q8) {
        tmp16no1 = (int16_t)WEBRTC_SPL_MUL_16_16_RSFT_WITH_ROUND(
                     inst->noiseEstDensity[offset + i], countProd, 15);
        tmp16no2 = (int16_t)WEBRTC_SPL_MUL_16_16_RSFT_WITH_ROUND(
                     width_factor, countDiv, 15);
        inst->noiseEstDensity[offset + i] = tmp16no1 + tmp16no2;
      }
    }  // end loop over magnitude spectrum

    if (counter >= END_STARTUP_LONG) {
      inst->noiseEstCounter[s] = 0;
      if (inst->blockIndex >= END_STARTUP_LONG) {
        UpdateNoiseEstimate(inst, offset);
      }
    }
    inst->noiseEstCounter[s]++;

  }  // end loop over simultaneous estimates

  // Sequentially update the noise during startup
  if (inst->blockIndex < END_STARTUP_LONG) {
    UpdateNoiseEstimate(inst, offset);
  }

  for (i = 0; i < inst->magnLen; i++) {
    noise[i] = (uint32_t)(inst->noiseEstQuantile[i]); // Q(qNoise)
  }
  (*q_noise) = (int16_t)inst->qNoise;
}

// Filter the data in the frequency domain, and create spectrum.
static void PrepareSpectrumC(NoiseSuppressionFixedC* inst, int16_t* freq_buf) {
  size_t i = 0, j = 0;

  for (i = 0; i < inst->magnLen; i++) {
    inst->real[i] = (int16_t)((inst->real[i] *
        (int16_t)(inst->noiseSupFilter[i])) >> 14);  // Q(normData-stages)
    inst->imag[i] = (int16_t)((inst->imag[i] *
        (int16_t)(inst->noiseSupFilter[i])) >> 14);  // Q(normData-stages)
  }

  freq_buf[0] = inst->real[0];
  freq_buf[1] = -inst->imag[0];
  for (i = 1, j = 2; i < inst->anaLen2; i += 1, j += 2) {
    freq_buf[j] = inst->real[i];
    freq_buf[j + 1] = -inst->imag[i];
  }
  freq_buf[inst->anaLen] = inst->real[inst->anaLen2];
  freq_buf[inst->anaLen + 1] = -inst->imag[inst->anaLen2];
}

// Denormalize the real-valued signal |in|, the output from inverse FFT.
static void DenormalizeC(NoiseSuppressionFixedC* inst,
                         int16_t* in,
                         int factor) {
  size_t i = 0;
  int32_t tmp32 = 0;
  for (i = 0; i < inst->anaLen; i += 1) {
    tmp32 = WEBRTC_SPL_SHIFT_W32((int32_t)in[i],
                                 factor - inst->normData);
    inst->real[i] = WebRtcSpl_SatW32ToW16(tmp32); // Q0
  }
}

// For the noise supression process, synthesis, read out fully processed
// segment, and update synthesis buffer.
static void SynthesisUpdateC(NoiseSuppressionFixedC* inst,
                             int16_t* out_frame,
                             int16_t gain_factor) {
  size_t i = 0;
  int16_t tmp16a = 0;
  int16_t tmp16b = 0;
  int32_t tmp32 = 0;

  // synthesis
  for (i = 0; i < inst->anaLen; i++) {
    tmp16a = (int16_t)WEBRTC_SPL_MUL_16_16_RSFT_WITH_ROUND(
                 inst->window[i], inst->real[i], 14); // Q0, window in Q14
    tmp32 = WEBRTC_SPL_MUL_16_16_RSFT_WITH_ROUND(tmp16a, gain_factor, 13); // Q0
    // Down shift with rounding
    tmp16b = WebRtcSpl_SatW32ToW16(tmp32); // Q0
    inst->synthesisBuffer[i] = WebRtcSpl_AddSatW16(inst->synthesisBuffer[i],
                                                   tmp16b); // Q0
  }

  // read out fully processed segment
  for (i = 0; i < inst->blockLen10ms; i++) {
    out_frame[i] = inst->synthesisBuffer[i]; // Q0
  }

  // update synthesis buffer
  memcpy(inst->synthesisBuffer, inst->synthesisBuffer + inst->blockLen10ms,
      (inst->anaLen - inst->blockLen10ms) * sizeof(*inst->synthesisBuffer));
  WebRtcSpl_ZerosArrayW16(inst->synthesisBuffer
      + inst->anaLen - inst->blockLen10ms, inst->blockLen10ms);
}

// Update analysis buffer for lower band, and window data before FFT.
static void AnalysisUpdateC(NoiseSuppressionFixedC* inst,
                            int16_t* out,
                            int16_t* new_speech) {
  size_t i = 0;

  // For lower band update analysis buffer.
  memcpy(inst->analysisBuffer, inst->analysisBuffer + inst->blockLen10ms,
      (inst->anaLen - inst->blockLen10ms) * sizeof(*inst->analysisBuffer));
  memcpy(inst->analysisBuffer + inst->anaLen - inst->blockLen10ms, new_speech,
      inst->blockLen10ms * sizeof(*inst->analysisBuffer));

  // Window data before FFT.
  for (i = 0; i < inst->anaLen; i++) {
    out[i] = (int16_t)WEBRTC_SPL_MUL_16_16_RSFT_WITH_ROUND(
               inst->window[i], inst->analysisBuffer[i], 14); // Q0
  }
}

// Normalize the real-valued signal |in|, the input to forward FFT.
static void NormalizeRealBufferC(NoiseSuppressionFixedC* inst,
                                 const int16_t* in,
                                 int16_t* out) {
  size_t i = 0;
  assert(inst->normData >= 0);
  for (i = 0; i < inst->anaLen; ++i) {
    out[i] = in[i] << inst->normData;  // Q(normData)
  }
}

// Declare function pointers.
NoiseEstimation WebRtcNsx_NoiseEstimation;
PrepareSpectrum WebRtcNsx_PrepareSpectrum;
SynthesisUpdate WebRtcNsx_SynthesisUpdate;
AnalysisUpdate WebRtcNsx_AnalysisUpdate;
Denormalize WebRtcNsx_Denormalize;
NormalizeRealBuffer WebRtcNsx_NormalizeRealBuffer;

#if defined(WEBRTC_HAS_NEON)
// Initialize function pointers for ARM Neon platform.
static void WebRtcNsx_InitNeon(void) {
  WebRtcNsx_NoiseEstimation = WebRtcNsx_NoiseEstimationNeon;
  WebRtcNsx_PrepareSpectrum = WebRtcNsx_PrepareSpectrumNeon;
  WebRtcNsx_SynthesisUpdate = WebRtcNsx_SynthesisUpdateNeon;
  WebRtcNsx_AnalysisUpdate = WebRtcNsx_AnalysisUpdateNeon;
}
#endif

#if defined(MIPS32_LE)
// Initialize function pointers for MIPS platform.
static void WebRtcNsx_InitMips(void) {
  WebRtcNsx_PrepareSpectrum = WebRtcNsx_PrepareSpectrum_mips;
  WebRtcNsx_SynthesisUpdate = WebRtcNsx_SynthesisUpdate_mips;
  WebRtcNsx_AnalysisUpdate = WebRtcNsx_AnalysisUpdate_mips;
  WebRtcNsx_NormalizeRealBuffer = WebRtcNsx_NormalizeRealBuffer_mips;
#if defined(MIPS_DSP_R1_LE)
  WebRtcNsx_Denormalize = WebRtcNsx_Denormalize_mips;
#endif
}
#endif

void WebRtcNsx_CalcParametricNoiseEstimate(NoiseSuppressionFixedC* inst,
                                           int16_t pink_noise_exp_avg,
                                           int32_t pink_noise_num_avg,
                                           int freq_index,
                                           uint32_t* noise_estimate,
                                           uint32_t* noise_estimate_avg) {
  int32_t tmp32no1 = 0;
  int32_t tmp32no2 = 0;

  int16_t int_part = 0;
  int16_t frac_part = 0;

  // Use pink noise estimate
  // noise_estimate = 2^(pinkNoiseNumerator + pinkNoiseExp * log2(j))
  assert(freq_index >= 0);
  assert(freq_index < 129);
  tmp32no2 = (pink_noise_exp_avg * kLogIndex[freq_index]) >> 15;  // Q11
  tmp32no1 = pink_noise_num_avg - tmp32no2; // Q11

  // Calculate output: 2^tmp32no1
  // Output in Q(minNorm-stages)
  tmp32no1 += (inst->minNorm - inst->stages) << 11;
  if (tmp32no1 > 0) {
    int_part = (int16_t)(tmp32no1 >> 11);
    frac_part = (int16_t)(tmp32no1 & 0x000007ff); // Q11
    // Piecewise linear approximation of 'b' in
    // 2^(int_part+frac_part) = 2^int_part * (1 + b)
    // 'b' is given in Q11 and below stored in frac_part.
    if (frac_part >> 10) {
      // Upper fractional part
      tmp32no2 = (2048 - frac_part) * 1244;  // Q21
      tmp32no2 = 2048 - (tmp32no2 >> 10);
    } else {
      // Lower fractional part
      tmp32no2 = (frac_part * 804) >> 10;
    }
    // Shift fractional part to Q(minNorm-stages)
    tmp32no2 = WEBRTC_SPL_SHIFT_W32(tmp32no2, int_part - 11);
    *noise_estimate_avg = (1 << int_part) + (uint32_t)tmp32no2;
    // Scale up to initMagnEst, which is not block averaged
    *noise_estimate = (*noise_estimate_avg) * (uint32_t)(inst->blockIndex + 1);
  }
}

// Initialize state
int32_t WebRtcNsx_InitCore(NoiseSuppressionFixedC* inst, uint32_t fs) {
  int i;

  //check for valid pointer
  if (inst == NULL) {
    return -1;
  }
  //

  // Initialization of struct
  if (fs == 8000 || fs == 16000 || fs == 32000 || fs == 48000) {
    inst->fs = fs;
  } else {
    return -1;
  }

  if (fs == 8000) {
    inst->blockLen10ms = 80;
    inst->anaLen = 128;
    inst->stages = 7;
    inst->window = kBlocks80w128x;
    inst->thresholdLogLrt = 131072; //default threshold for LRT feature
    inst->maxLrt = 0x0040000;
    inst->minLrt = 52429;
  } else {
    inst->blockLen10ms = 160;
    inst->anaLen = 256;
    inst->stages = 8;
    inst->window = kBlocks160w256x;
    inst->thresholdLogLrt = 212644; //default threshold for LRT feature
    inst->maxLrt = 0x0080000;
    inst->minLrt = 104858;
  }
  inst->anaLen2 = inst->anaLen / 2;
  inst->magnLen = inst->anaLen2 + 1;

  if (inst->real_fft != NULL) {
    WebRtcSpl_FreeRealFFT(inst->real_fft);
  }
  inst->real_fft = WebRtcSpl_CreateRealFFT(inst->stages);
  if (inst->real_fft == NULL) {
    return -1;
  }

  WebRtcSpl_ZerosArrayW16(inst->analysisBuffer, ANAL_BLOCKL_MAX);
  WebRtcSpl_ZerosArrayW16(inst->synthesisBuffer, ANAL_BLOCKL_MAX);

  // for HB processing
  WebRtcSpl_ZerosArrayW16(inst->dataBufHBFX[0],
                          NUM_HIGH_BANDS_MAX * ANAL_BLOCKL_MAX);
  // for quantile noise estimation
  WebRtcSpl_ZerosArrayW16(inst->noiseEstQuantile, HALF_ANAL_BLOCKL);
  for (i = 0; i < SIMULT * HALF_ANAL_BLOCKL; i++) {
    inst->noiseEstLogQuantile[i] = 2048; // Q8
    inst->noiseEstDensity[i] = 153; // Q9
  }
  for (i = 0; i < SIMULT; i++) {
    inst->noiseEstCounter[i] = (int16_t)(END_STARTUP_LONG * (i + 1)) / SIMULT;
  }

  // Initialize suppression filter with ones
  WebRtcSpl_MemSetW16((int16_t*)inst->noiseSupFilter, 16384, HALF_ANAL_BLOCKL);

  // Set the aggressiveness: default
  inst->aggrMode = 0;

  //initialize variables for new method
  inst->priorNonSpeechProb = 8192; // Q14(0.5) prior probability for speech/noise
  for (i = 0; i < HALF_ANAL_BLOCKL; i++) {
    inst->prevMagnU16[i] = 0;
    inst->prevNoiseU32[i] = 0; //previous noise-spectrum
    inst->logLrtTimeAvgW32[i] = 0; //smooth LR ratio
    inst->avgMagnPause[i] = 0; //conservative noise spectrum estimate
    inst->initMagnEst[i] = 0; //initial average magnitude spectrum
  }

  //feature quantities
  inst->thresholdSpecDiff = 50; //threshold for difference feature: determined on-line
  inst->thresholdSpecFlat = 20480; //threshold for flatness: determined on-line
  inst->featureLogLrt = inst->thresholdLogLrt; //average LRT factor (= threshold)
  inst->featureSpecFlat = inst->thresholdSpecFlat; //spectral flatness (= threshold)
  inst->featureSpecDiff = inst->thresholdSpecDiff; //spectral difference (= threshold)
  inst->weightLogLrt = 6; //default weighting par for LRT feature
  inst->weightSpecFlat = 0; //default weighting par for spectral flatness feature
  inst->weightSpecDiff = 0; //default weighting par for spectral difference feature

  inst->curAvgMagnEnergy = 0; //window time-average of input magnitude spectrum
  inst->timeAvgMagnEnergy = 0; //normalization for spectral difference
  inst->timeAvgMagnEnergyTmp = 0; //normalization for spectral difference

  //histogram quantities: used to estimate/update thresholds for features
  WebRtcSpl_ZerosArrayW16(inst->histLrt, HIST_PAR_EST);
  WebRtcSpl_ZerosArrayW16(inst->histSpecDiff, HIST_PAR_EST);
  WebRtcSpl_ZerosArrayW16(inst->histSpecFlat, HIST_PAR_EST);

  inst->blockIndex = -1; //frame counter

  //inst->modelUpdate    = 500;   //window for update
  inst->modelUpdate = (1 << STAT_UPDATES); //window for update
  inst->cntThresUpdate = 0; //counter feature thresholds updates

  inst->sumMagn = 0;
  inst->magnEnergy = 0;
  inst->prevQMagn = 0;
  inst->qNoise = 0;
  inst->prevQNoise = 0;

  inst->energyIn = 0;
  inst->scaleEnergyIn = 0;

  inst->whiteNoiseLevel = 0;
  inst->pinkNoiseNumerator = 0;
  inst->pinkNoiseExp = 0;
  inst->minNorm = 15; // Start with full scale
  inst->zeroInputSignal = 0;

  //default mode
  WebRtcNsx_set_policy_core(inst, 0);

#ifdef NS_FILEDEBUG
  inst->infile = fopen("indebug.pcm", "wb");
  inst->outfile = fopen("outdebug.pcm", "wb");
  inst->file1 = fopen("file1.pcm", "wb");
  inst->file2 = fopen("file2.pcm", "wb");
  inst->file3 = fopen("file3.pcm", "wb");
  inst->file4 = fopen("file4.pcm", "wb");
  inst->file5 = fopen("file5.pcm", "wb");
#endif

  // Initialize function pointers.
  WebRtcNsx_NoiseEstimation = NoiseEstimationC;
  WebRtcNsx_PrepareSpectrum = PrepareSpectrumC;
  WebRtcNsx_SynthesisUpdate = SynthesisUpdateC;
  WebRtcNsx_AnalysisUpdate = AnalysisUpdateC;
  WebRtcNsx_Denormalize = DenormalizeC;
  WebRtcNsx_NormalizeRealBuffer = NormalizeRealBufferC;

#if defined(WEBRTC_HAS_NEON)
  WebRtcNsx_InitNeon();
#endif

#if defined(MIPS32_LE)
  WebRtcNsx_InitMips();
#endif

  inst->initFlag = 1;

  return 0;
}

int WebRtcNsx_set_policy_core(NoiseSuppressionFixedC* inst, int mode) {
  // allow for modes:0,1,2,3
  if (mode < 0 || mode > 3) {
    return -1;
  }

  inst->aggrMode = mode;
  if (mode == 0) {
    inst->overdrive = 256; // Q8(1.0)
    inst->denoiseBound = 8192; // Q14(0.5)
    inst->gainMap = 0; // No gain compensation
  } else if (mode == 1) {
    inst->overdrive = 256; // Q8(1.0)
    inst->denoiseBound = 4096; // Q14(0.25)
    inst->factor2Table = kFactor2Aggressiveness1;
    inst->gainMap = 1;
  } else if (mode == 2) {
    inst->overdrive = 282; // ~= Q8(1.1)
    inst->denoiseBound = 2048; // Q14(0.125)
    inst->factor2Table = kFactor2Aggressiveness2;
    inst->gainMap = 1;
  } else if (mode == 3) {
    inst->overdrive = 320; // Q8(1.25)
    inst->denoiseBound = 1475; // ~= Q14(0.09)
    inst->factor2Table = kFactor2Aggressiveness3;
    inst->gainMap = 1;
  }
  return 0;
}

// Extract thresholds for feature parameters
// histograms are computed over some window_size (given by window_pars)
// thresholds and weights are extracted every window
// flag 0 means update histogram only, flag 1 means compute the thresholds/weights
// threshold and weights are returned in: inst->priorModelPars
void WebRtcNsx_FeatureParameterExtraction(NoiseSuppressionFixedC* inst,
                                          int flag) {
  uint32_t tmpU32;
  uint32_t histIndex;
  uint32_t posPeak1SpecFlatFX, posPeak2SpecFlatFX;
  uint32_t posPeak1SpecDiffFX, posPeak2SpecDiffFX;

  int32_t tmp32;
  int32_t fluctLrtFX, thresFluctLrtFX;
  int32_t avgHistLrtFX, avgSquareHistLrtFX, avgHistLrtComplFX;

  int16_t j;
  int16_t numHistLrt;

  int i;
  int useFeatureSpecFlat, useFeatureSpecDiff, featureSum;
  int maxPeak1, maxPeak2;
  int weightPeak1SpecFlat, weightPeak2SpecFlat;
  int weightPeak1SpecDiff, weightPeak2SpecDiff;

  //update histograms
  if (!flag) {
    // LRT
    // Type casting to UWord32 is safe since negative values will not be wrapped to larger
    // values than HIST_PAR_EST
    histIndex = (uint32_t)(inst->featureLogLrt);
    if (histIndex < HIST_PAR_EST) {
      inst->histLrt[histIndex]++;
    }
    // Spectral flatness
    // (inst->featureSpecFlat*20)>>10 = (inst->featureSpecFlat*5)>>8
    histIndex = (inst->featureSpecFlat * 5) >> 8;
    if (histIndex < HIST_PAR_EST) {
      inst->histSpecFlat[histIndex]++;
    }
    // Spectral difference
    histIndex = HIST_PAR_EST;
    if (inst->timeAvgMagnEnergy > 0) {
      // Guard against division by zero
      // If timeAvgMagnEnergy == 0 we have no normalizing statistics and
      // therefore can't update the histogram
      histIndex = ((inst->featureSpecDiff * 5) >> inst->stages) /
          inst->timeAvgMagnEnergy;
    }
    if (histIndex < HIST_PAR_EST) {
      inst->histSpecDiff[histIndex]++;
    }
  }

  // extract parameters for speech/noise probability
  if (flag) {
    useFeatureSpecDiff = 1;
    //for LRT feature:
    // compute the average over inst->featureExtractionParams.rangeAvgHistLrt
    avgHistLrtFX = 0;
    avgSquareHistLrtFX = 0;
    numHistLrt = 0;
    for (i = 0; i < BIN_SIZE_LRT; i++) {
      j = (2 * i + 1);
      tmp32 = inst->histLrt[i] * j;
      avgHistLrtFX += tmp32;
      numHistLrt += inst->histLrt[i];
      avgSquareHistLrtFX += tmp32 * j;
    }
    avgHistLrtComplFX = avgHistLrtFX;
    for (; i < HIST_PAR_EST; i++) {
      j = (2 * i + 1);
      tmp32 = inst->histLrt[i] * j;
      avgHistLrtComplFX += tmp32;
      avgSquareHistLrtFX += tmp32 * j;
    }
    fluctLrtFX = avgSquareHistLrtFX * numHistLrt -
        avgHistLrtFX * avgHistLrtComplFX;
    thresFluctLrtFX = THRES_FLUCT_LRT * numHistLrt;
    // get threshold for LRT feature:
    tmpU32 = (FACTOR_1_LRT_DIFF * (uint32_t)avgHistLrtFX);
    if ((fluctLrtFX < thresFluctLrtFX) || (numHistLrt == 0) ||
        (tmpU32 > (uint32_t)(100 * numHistLrt))) {
      //very low fluctuation, so likely noise
      inst->thresholdLogLrt = inst->maxLrt;
    } else {
      tmp32 = (int32_t)((tmpU32 << (9 + inst->stages)) / numHistLrt /
                              25);
      // check if value is within min/max range
      inst->thresholdLogLrt = WEBRTC_SPL_SAT(inst->maxLrt,
                                             tmp32,
                                             inst->minLrt);
    }
    if (fluctLrtFX < thresFluctLrtFX) {
      // Do not use difference feature if fluctuation of LRT feature is very low:
      // most likely just noise state
      useFeatureSpecDiff = 0;
    }

    // for spectral flatness and spectral difference: compute the main peaks of histogram
    maxPeak1 = 0;
    maxPeak2 = 0;
    posPeak1SpecFlatFX = 0;
    posPeak2SpecFlatFX = 0;
    weightPeak1SpecFlat = 0;
    weightPeak2SpecFlat = 0;

    // peaks for flatness
    for (i = 0; i < HIST_PAR_EST; i++) {
      if (inst->histSpecFlat[i] > maxPeak1) {
        // Found new "first" peak
        maxPeak2 = maxPeak1;
        weightPeak2SpecFlat = weightPeak1SpecFlat;
        posPeak2SpecFlatFX = posPeak1SpecFlatFX;

        maxPeak1 = inst->histSpecFlat[i];
        weightPeak1SpecFlat = inst->histSpecFlat[i];
        posPeak1SpecFlatFX = (uint32_t)(2 * i + 1);
      } else if (inst->histSpecFlat[i] > maxPeak2) {
        // Found new "second" peak
        maxPeak2 = inst->histSpecFlat[i];
        weightPeak2SpecFlat = inst->histSpecFlat[i];
        posPeak2SpecFlatFX = (uint32_t)(2 * i + 1);
      }
    }

    // for spectral flatness feature
    useFeatureSpecFlat = 1;
    // merge the two peaks if they are close
    if ((posPeak1SpecFlatFX - posPeak2SpecFlatFX < LIM_PEAK_SPACE_FLAT_DIFF)
        && (weightPeak2SpecFlat * LIM_PEAK_WEIGHT_FLAT_DIFF > weightPeak1SpecFlat)) {
      weightPeak1SpecFlat += weightPeak2SpecFlat;
      posPeak1SpecFlatFX = (posPeak1SpecFlatFX + posPeak2SpecFlatFX) >> 1;
    }
    //reject if weight of peaks is not large enough, or peak value too small
    if (weightPeak1SpecFlat < THRES_WEIGHT_FLAT_DIFF || posPeak1SpecFlatFX
        < THRES_PEAK_FLAT) {
      useFeatureSpecFlat = 0;
    } else { // if selected, get the threshold
      // compute the threshold and check if value is within min/max range
      inst->thresholdSpecFlat = WEBRTC_SPL_SAT(MAX_FLAT_Q10, FACTOR_2_FLAT_Q10
                                               * posPeak1SpecFlatFX, MIN_FLAT_Q10); //Q10
    }
    // done with flatness feature

    if (useFeatureSpecDiff) {
      //compute two peaks for spectral difference
      maxPeak1 = 0;
      maxPeak2 = 0;
      posPeak1SpecDiffFX = 0;
      posPeak2SpecDiffFX = 0;
      weightPeak1SpecDiff = 0;
      weightPeak2SpecDiff = 0;
      // peaks for spectral difference
      for (i = 0; i < HIST_PAR_EST; i++) {
        if (inst->histSpecDiff[i] > maxPeak1) {
          // Found new "first" peak
          maxPeak2 = maxPeak1;
          weightPeak2SpecDiff = weightPeak1SpecDiff;
          posPeak2SpecDiffFX = posPeak1SpecDiffFX;

          maxPeak1 = inst->histSpecDiff[i];
          weightPeak1SpecDiff = inst->histSpecDiff[i];
          posPeak1SpecDiffFX = (uint32_t)(2 * i + 1);
        } else if (inst->histSpecDiff[i] > maxPeak2) {
          // Found new "second" peak
          maxPeak2 = inst->histSpecDiff[i];
          weightPeak2SpecDiff = inst->histSpecDiff[i];
          posPeak2SpecDiffFX = (uint32_t)(2 * i + 1);
        }
      }

      // merge the two peaks if they are close
      if ((posPeak1SpecDiffFX - posPeak2SpecDiffFX < LIM_PEAK_SPACE_FLAT_DIFF)
          && (weightPeak2SpecDiff * LIM_PEAK_WEIGHT_FLAT_DIFF > weightPeak1SpecDiff)) {
        weightPeak1SpecDiff += weightPeak2SpecDiff;
        posPeak1SpecDiffFX = (posPeak1SpecDiffFX + posPeak2SpecDiffFX) >> 1;
      }
      // get the threshold value and check if value is within min/max range
      inst->thresholdSpecDiff = WEBRTC_SPL_SAT(MAX_DIFF, FACTOR_1_LRT_DIFF
                                               * posPeak1SpecDiffFX, MIN_DIFF); //5x bigger
      //reject if weight of peaks is not large enough
      if (weightPeak1SpecDiff < THRES_WEIGHT_FLAT_DIFF) {
        useFeatureSpecDiff = 0;
      }
      // done with spectral difference feature
    }

    // select the weights between the features
    // inst->priorModelPars[4] is weight for LRT: always selected
    featureSum = 6 / (1 + useFeatureSpecFlat + useFeatureSpecDiff);
    inst->weightLogLrt = featureSum;
    inst->weightSpecFlat = useFeatureSpecFlat * featureSum;
    inst->weightSpecDiff = useFeatureSpecDiff * featureSum;

    // set histograms to zero for next update
    WebRtcSpl_ZerosArrayW16(inst->histLrt, HIST_PAR_EST);
    WebRtcSpl_ZerosArrayW16(inst->histSpecDiff, HIST_PAR_EST);
    WebRtcSpl_ZerosArrayW16(inst->histSpecFlat, HIST_PAR_EST);
  }  // end of flag == 1
}


// Compute spectral flatness on input spectrum
// magn is the magnitude spectrum
// spectral flatness is returned in inst->featureSpecFlat
void WebRtcNsx_ComputeSpectralFlatness(NoiseSuppressionFixedC* inst,
                                       uint16_t* magn) {
  uint32_t tmpU32;
  uint32_t avgSpectralFlatnessNum, avgSpectralFlatnessDen;

  int32_t tmp32;
  int32_t currentSpectralFlatness, logCurSpectralFlatness;

  int16_t zeros, frac, intPart;

  size_t i;

  // for flatness
  avgSpectralFlatnessNum = 0;
  avgSpectralFlatnessDen = inst->sumMagn - (uint32_t)magn[0]; // Q(normData-stages)

  // compute log of ratio of the geometric to arithmetic mean: check for log(0) case
  // flatness = exp( sum(log(magn[i]))/N - log(sum(magn[i])/N) )
  //          = exp( sum(log(magn[i]))/N ) * N / sum(magn[i])
  //          = 2^( sum(log2(magn[i]))/N - (log2(sum(magn[i])) - log2(N)) ) [This is used]
  for (i = 1; i < inst->magnLen; i++) {
    // First bin is excluded from spectrum measures. Number of bins is now a power of 2
    if (magn[i]) {
      zeros = WebRtcSpl_NormU32((uint32_t)magn[i]);
      frac = (int16_t)(((uint32_t)((uint32_t)(magn[i]) << zeros)
                              & 0x7FFFFFFF) >> 23);
      // log2(magn(i))
      assert(frac < 256);
      tmpU32 = (uint32_t)(((31 - zeros) << 8)
                                + WebRtcNsx_kLogTableFrac[frac]); // Q8
      avgSpectralFlatnessNum += tmpU32; // Q8
    } else {
      //if at least one frequency component is zero, treat separately
      tmpU32 = WEBRTC_SPL_UMUL_32_16(inst->featureSpecFlat, SPECT_FLAT_TAVG_Q14); // Q24
      inst->featureSpecFlat -= tmpU32 >> 14;  // Q10
      return;
    }
  }
  //ratio and inverse log: check for case of log(0)
  zeros = WebRtcSpl_NormU32(avgSpectralFlatnessDen);
  frac = (int16_t)(((avgSpectralFlatnessDen << zeros) & 0x7FFFFFFF) >> 23);
  // log2(avgSpectralFlatnessDen)
  assert(frac < 256);
  tmp32 = (int32_t)(((31 - zeros) << 8) + WebRtcNsx_kLogTableFrac[frac]); // Q8
  logCurSpectralFlatness = (int32_t)avgSpectralFlatnessNum;
  logCurSpectralFlatness += ((int32_t)(inst->stages - 1) << (inst->stages + 7)); // Q(8+stages-1)
  logCurSpectralFlatness -= (tmp32 << (inst->stages - 1));
  logCurSpectralFlatness <<= (10 - inst->stages);  // Q17
  tmp32 = (int32_t)(0x00020000 | (WEBRTC_SPL_ABS_W32(logCurSpectralFlatness)
                                        & 0x0001FFFF)); //Q17
  intPart = 7 - (logCurSpectralFlatness >> 17);  // Add 7 for output in Q10.
  if (intPart > 0) {
    currentSpectralFlatness = tmp32 >> intPart;
  } else {
    currentSpectralFlatness = tmp32 << -intPart;
  }

  //time average update of spectral flatness feature
  tmp32 = currentSpectralFlatness - (int32_t)inst->featureSpecFlat; // Q10
  tmp32 *= SPECT_FLAT_TAVG_Q14;  // Q24
  inst->featureSpecFlat += tmp32 >> 14;  // Q10
  // done with flatness feature
}


// Compute the difference measure between input spectrum and a template/learned noise spectrum
// magn_tmp is the input spectrum
// the reference/template spectrum is  inst->magn_avg_pause[i]
// returns (normalized) spectral difference in inst->featureSpecDiff
void WebRtcNsx_ComputeSpectralDifference(NoiseSuppressionFixedC* inst,
                                         uint16_t* magnIn) {
  // This is to be calculated:
  // avgDiffNormMagn = var(magnIn) - cov(magnIn, magnAvgPause)^2 / var(magnAvgPause)

  uint32_t tmpU32no1, tmpU32no2;
  uint32_t varMagnUFX, varPauseUFX, avgDiffNormMagnUFX;

  int32_t tmp32no1, tmp32no2;
  int32_t avgPauseFX, avgMagnFX, covMagnPauseFX;
  int32_t maxPause, minPause;

  int16_t tmp16no1;

  size_t i;
  int norm32, nShifts;

  avgPauseFX = 0;
  maxPause = 0;
  minPause = inst->avgMagnPause[0]; // Q(prevQMagn)
  // compute average quantities
  for (i = 0; i < inst->magnLen; i++) {
    // Compute mean of magn_pause
    avgPauseFX += inst->avgMagnPause[i]; // in Q(prevQMagn)
    maxPause = WEBRTC_SPL_MAX(maxPause, inst->avgMagnPause[i]);
    minPause = WEBRTC_SPL_MIN(minPause, inst->avgMagnPause[i]);
  }
  // normalize by replacing div of "inst->magnLen" with "inst->stages-1" shifts
  avgPauseFX >>= inst->stages - 1;
  avgMagnFX = inst->sumMagn >> (inst->stages - 1);
  // Largest possible deviation in magnPause for (co)var calculations
  tmp32no1 = WEBRTC_SPL_MAX(maxPause - avgPauseFX, avgPauseFX - minPause);
  // Get number of shifts to make sure we don't get wrap around in varPause
  nShifts = WEBRTC_SPL_MAX(0, 10 + inst->stages - WebRtcSpl_NormW32(tmp32no1));

  varMagnUFX = 0;
  varPauseUFX = 0;
  covMagnPauseFX = 0;
  for (i = 0; i < inst->magnLen; i++) {
    // Compute var and cov of magn and magn_pause
    tmp16no1 = (int16_t)((int32_t)magnIn[i] - avgMagnFX);
    tmp32no2 = inst->avgMagnPause[i] - avgPauseFX;
    varMagnUFX += (uint32_t)(tmp16no1 * tmp16no1);  // Q(2*qMagn)
    tmp32no1 = tmp32no2 * tmp16no1;  // Q(prevQMagn+qMagn)
    covMagnPauseFX += tmp32no1; // Q(prevQMagn+qMagn)
    tmp32no1 = tmp32no2 >> nShifts;  // Q(prevQMagn-minPause).
    varPauseUFX += tmp32no1 * tmp32no1;  // Q(2*(prevQMagn-minPause))
  }
  //update of average magnitude spectrum: Q(-2*stages) and averaging replaced by shifts
  inst->curAvgMagnEnergy +=
      inst->magnEnergy >> (2 * inst->normData + inst->stages - 1);

  avgDiffNormMagnUFX = varMagnUFX; // Q(2*qMagn)
  if ((varPauseUFX) && (covMagnPauseFX)) {
    tmpU32no1 = (uint32_t)WEBRTC_SPL_ABS_W32(covMagnPauseFX); // Q(prevQMagn+qMagn)
    norm32 = WebRtcSpl_NormU32(tmpU32no1) - 16;
    if (norm32 > 0) {
      tmpU32no1 <<= norm32;  // Q(prevQMagn+qMagn+norm32)
    } else {
      tmpU32no1 >>= -norm32;  // Q(prevQMagn+qMagn+norm32)
    }
    tmpU32no2 = WEBRTC_SPL_UMUL(tmpU32no1, tmpU32no1); // Q(2*(prevQMagn+qMagn-norm32))

    nShifts += norm32;
    nShifts <<= 1;
    if (nShifts < 0) {
      varPauseUFX >>= (-nShifts); // Q(2*(qMagn+norm32+minPause))
      nShifts = 0;
    }
    if (varPauseUFX > 0) {
      // Q(2*(qMagn+norm32-16+minPause))
      tmpU32no1 = tmpU32no2 / varPauseUFX;
      tmpU32no1 >>= nShifts;

      // Q(2*qMagn)
      avgDiffNormMagnUFX -= WEBRTC_SPL_MIN(avgDiffNormMagnUFX, tmpU32no1);
    } else {
      avgDiffNormMagnUFX = 0;
    }
  }
  //normalize and compute time average update of difference feature
  tmpU32no1 = avgDiffNormMagnUFX >> (2 * inst->normData);
  if (inst->featureSpecDiff > tmpU32no1) {
    tmpU32no2 = WEBRTC_SPL_UMUL_32_16(inst->featureSpecDiff - tmpU32no1,
                                      SPECT_DIFF_TAVG_Q8); // Q(8-2*stages)
    inst->featureSpecDiff -= tmpU32no2 >> 8;  // Q(-2*stages)
  } else {
    tmpU32no2 = WEBRTC_SPL_UMUL_32_16(tmpU32no1 - inst->featureSpecDiff,
                                      SPECT_DIFF_TAVG_Q8); // Q(8-2*stages)
    inst->featureSpecDiff += tmpU32no2 >> 8;  // Q(-2*stages)
  }
}

// Transform input (speechFrame) to frequency domain magnitude (magnU16)
void WebRtcNsx_DataAnalysis(NoiseSuppressionFixedC* inst,
                            short* speechFrame,
                            uint16_t* magnU16) {
  uint32_t tmpU32no1;

  int32_t   tmp_1_w32 = 0;
  int32_t   tmp_2_w32 = 0;
  int32_t   sum_log_magn = 0;
  int32_t   sum_log_i_log_magn = 0;

  uint16_t  sum_log_magn_u16 = 0;
  uint16_t  tmp_u16 = 0;

  int16_t   sum_log_i = 0;
  int16_t   sum_log_i_square = 0;
  int16_t   frac = 0;
  int16_t   log2 = 0;
  int16_t   matrix_determinant = 0;
  int16_t   maxWinData;

  size_t i, j;
  int zeros;
  int net_norm = 0;
  int right_shifts_in_magnU16 = 0;
  int right_shifts_in_initMagnEst = 0;

  int16_t winData_buff[ANAL_BLOCKL_MAX * 2 + 16];
  int16_t realImag_buff[ANAL_BLOCKL_MAX * 2 + 16];

  // Align the structures to 32-byte boundary for the FFT function.
  int16_t* winData = (int16_t*) (((uintptr_t)winData_buff + 31) & ~31);
  int16_t* realImag = (int16_t*) (((uintptr_t) realImag_buff + 31) & ~31);

  // Update analysis buffer for lower band, and window data before FFT.
  WebRtcNsx_AnalysisUpdate(inst, winData, speechFrame);

  // Get input energy
  inst->energyIn =
      WebRtcSpl_Energy(winData, inst->anaLen, &inst->scaleEnergyIn);

  // Reset zero input flag
  inst->zeroInputSignal = 0;
  // Acquire norm for winData
  maxWinData = WebRtcSpl_MaxAbsValueW16(winData, inst->anaLen);
  inst->normData = WebRtcSpl_NormW16(maxWinData);
  if (maxWinData == 0) {
    // Treat zero input separately.
    inst->zeroInputSignal = 1;
    return;
  }

  // Determine the net normalization in the frequency domain
  net_norm = inst->stages - inst->normData;
  // Track lowest normalization factor and use it to prevent wrap around in shifting
  right_shifts_in_magnU16 = inst->normData - inst->minNorm;
  right_shifts_in_initMagnEst = WEBRTC_SPL_MAX(-right_shifts_in_magnU16, 0);
  inst->minNorm -= right_shifts_in_initMagnEst;
  right_shifts_in_magnU16 = WEBRTC_SPL_MAX(right_shifts_in_magnU16, 0);

  // create realImag as winData interleaved with zeros (= imag. part), normalize it
  WebRtcNsx_NormalizeRealBuffer(inst, winData, realImag);

  // FFT output will be in winData[].
  WebRtcSpl_RealForwardFFT(inst->real_fft, realImag, winData);

  inst->imag[0] = 0; // Q(normData-stages)
  inst->imag[inst->anaLen2] = 0;
  inst->real[0] = winData[0]; // Q(normData-stages)
  inst->real[inst->anaLen2] = winData[inst->anaLen];
  // Q(2*(normData-stages))
  inst->magnEnergy = (uint32_t)(inst->real[0] * inst->real[0]);
  inst->magnEnergy += (uint32_t)(inst->real[inst->anaLen2] *
                                 inst->real[inst->anaLen2]);
  magnU16[0] = (uint16_t)WEBRTC_SPL_ABS_W16(inst->real[0]); // Q(normData-stages)
  magnU16[inst->anaLen2] = (uint16_t)WEBRTC_SPL_ABS_W16(inst->real[inst->anaLen2]);
  inst->sumMagn = (uint32_t)magnU16[0]; // Q(normData-stages)
  inst->sumMagn += (uint32_t)magnU16[inst->anaLen2];

  if (inst->blockIndex >= END_STARTUP_SHORT) {
    for (i = 1, j = 2; i < inst->anaLen2; i += 1, j += 2) {
      inst->real[i] = winData[j];
      inst->imag[i] = -winData[j + 1];
      // magnitude spectrum
      // energy in Q(2*(normData-stages))
      tmpU32no1 = (uint32_t)(winData[j] * winData[j]);
      tmpU32no1 += (uint32_t)(winData[j + 1] * winData[j + 1]);
      inst->magnEnergy += tmpU32no1; // Q(2*(normData-stages))

      magnU16[i] = (uint16_t)WebRtcSpl_SqrtFloor(tmpU32no1); // Q(normData-stages)
      inst->sumMagn += (uint32_t)magnU16[i]; // Q(normData-stages)
    }
  } else {
    //
    // Gather information during startup for noise parameter estimation
    //

    // Switch initMagnEst to Q(minNorm-stages)
    inst->initMagnEst[0] >>= right_shifts_in_initMagnEst;
    inst->initMagnEst[inst->anaLen2] >>= right_shifts_in_initMagnEst;

    // Update initMagnEst with magnU16 in Q(minNorm-stages).
    inst->initMagnEst[0] += magnU16[0] >> right_shifts_in_magnU16;
    inst->initMagnEst[inst->anaLen2] +=
        magnU16[inst->anaLen2] >> right_shifts_in_magnU16;

    log2 = 0;
    if (magnU16[inst->anaLen2]) {
      // Calculate log2(magnU16[inst->anaLen2])
      zeros = WebRtcSpl_NormU32((uint32_t)magnU16[inst->anaLen2]);
      frac = (int16_t)((((uint32_t)magnU16[inst->anaLen2] << zeros) &
                              0x7FFFFFFF) >> 23); // Q8
      // log2(magnU16(i)) in Q8
      assert(frac < 256);
      log2 = (int16_t)(((31 - zeros) << 8) + WebRtcNsx_kLogTableFrac[frac]);
    }

    sum_log_magn = (int32_t)log2; // Q8
    // sum_log_i_log_magn in Q17
    sum_log_i_log_magn = (kLogIndex[inst->anaLen2] * log2) >> 3;

    for (i = 1, j = 2; i < inst->anaLen2; i += 1, j += 2) {
      inst->real[i] = winData[j];
      inst->imag[i] = -winData[j + 1];
      // magnitude spectrum
      // energy in Q(2*(normData-stages))
      tmpU32no1 = (uint32_t)(winData[j] * winData[j]);
      tmpU32no1 += (uint32_t)(winData[j + 1] * winData[j + 1]);
      inst->magnEnergy += tmpU32no1; // Q(2*(normData-stages))

      magnU16[i] = (uint16_t)WebRtcSpl_SqrtFloor(tmpU32no1); // Q(normData-stages)
      inst->sumMagn += (uint32_t)magnU16[i]; // Q(normData-stages)

      // Switch initMagnEst to Q(minNorm-stages)
      inst->initMagnEst[i] >>= right_shifts_in_initMagnEst;

      // Update initMagnEst with magnU16 in Q(minNorm-stages).
      inst->initMagnEst[i] += magnU16[i] >> right_shifts_in_magnU16;

      if (i >= kStartBand) {
        // For pink noise estimation. Collect data neglecting lower frequency band
        log2 = 0;
        if (magnU16[i]) {
          zeros = WebRtcSpl_NormU32((uint32_t)magnU16[i]);
          frac = (int16_t)((((uint32_t)magnU16[i] << zeros) &
                                  0x7FFFFFFF) >> 23);
          // log2(magnU16(i)) in Q8
          assert(frac < 256);
          log2 = (int16_t)(((31 - zeros) << 8)
                                 + WebRtcNsx_kLogTableFrac[frac]);
        }
        sum_log_magn += (int32_t)log2; // Q8
        // sum_log_i_log_magn in Q17
        sum_log_i_log_magn += (kLogIndex[i] * log2) >> 3;
      }
    }

    //
    //compute simplified noise model during startup
    //

    // Estimate White noise

    // Switch whiteNoiseLevel to Q(minNorm-stages)
    inst->whiteNoiseLevel >>= right_shifts_in_initMagnEst;

    // Update the average magnitude spectrum, used as noise estimate.
    tmpU32no1 = WEBRTC_SPL_UMUL_32_16(inst->sumMagn, inst->overdrive);
    tmpU32no1 >>= inst->stages + 8;

    // Replacing division above with 'stages' shifts
    // Shift to same Q-domain as whiteNoiseLevel
    tmpU32no1 >>= right_shifts_in_magnU16;
    // This operation is safe from wrap around as long as END_STARTUP_SHORT < 128
    assert(END_STARTUP_SHORT < 128);
    inst->whiteNoiseLevel += tmpU32no1; // Q(minNorm-stages)

    // Estimate Pink noise parameters
    // Denominator used in both parameter estimates.
    // The value is only dependent on the size of the frequency band (kStartBand)
    // and to reduce computational complexity stored in a table (kDeterminantEstMatrix[])
    assert(kStartBand < 66);
    matrix_determinant = kDeterminantEstMatrix[kStartBand]; // Q0
    sum_log_i = kSumLogIndex[kStartBand]; // Q5
    sum_log_i_square = kSumSquareLogIndex[kStartBand]; // Q2
    if (inst->fs == 8000) {
      // Adjust values to shorter blocks in narrow band.
      tmp_1_w32 = (int32_t)matrix_determinant;
      tmp_1_w32 += (kSumLogIndex[65] * sum_log_i) >> 9;
      tmp_1_w32 -= (kSumLogIndex[65] * kSumLogIndex[65]) >> 10;
      tmp_1_w32 -= (int32_t)sum_log_i_square << 4;
      tmp_1_w32 -= ((inst->magnLen - kStartBand) * kSumSquareLogIndex[65]) >> 2;
      matrix_determinant = (int16_t)tmp_1_w32;
      sum_log_i -= kSumLogIndex[65]; // Q5
      sum_log_i_square -= kSumSquareLogIndex[65]; // Q2
    }

    // Necessary number of shifts to fit sum_log_magn in a word16
    zeros = 16 - WebRtcSpl_NormW32(sum_log_magn);
    if (zeros < 0) {
      zeros = 0;
    }
    tmp_1_w32 = sum_log_magn << 1;  // Q9
    sum_log_magn_u16 = (uint16_t)(tmp_1_w32 >> zeros);  // Q(9-zeros).

    // Calculate and update pinkNoiseNumerator. Result in Q11.
    tmp_2_w32 = WEBRTC_SPL_MUL_16_U16(sum_log_i_square, sum_log_magn_u16); // Q(11-zeros)
    tmpU32no1 = sum_log_i_log_magn >> 12;  // Q5

    // Shift the largest value of sum_log_i and tmp32no3 before multiplication
    tmp_u16 = ((uint16_t)sum_log_i << 1);  // Q6
    if ((uint32_t)sum_log_i > tmpU32no1) {
      tmp_u16 >>= zeros;
    } else {
      tmpU32no1 >>= zeros;
    }
    tmp_2_w32 -= (int32_t)WEBRTC_SPL_UMUL_32_16(tmpU32no1, tmp_u16); // Q(11-zeros)
    matrix_determinant >>= zeros;  // Q(-zeros)
    tmp_2_w32 = WebRtcSpl_DivW32W16(tmp_2_w32, matrix_determinant); // Q11
    tmp_2_w32 += (int32_t)net_norm << 11;  // Q11
    if (tmp_2_w32 < 0) {
      tmp_2_w32 = 0;
    }
    inst->pinkNoiseNumerator += tmp_2_w32; // Q11

    // Calculate and update pinkNoiseExp. Result in Q14.
    tmp_2_w32 = WEBRTC_SPL_MUL_16_U16(sum_log_i, sum_log_magn_u16); // Q(14-zeros)
    tmp_1_w32 = sum_log_i_log_magn >> (3 + zeros);
    tmp_1_w32 *= inst->magnLen - kStartBand;
    tmp_2_w32 -= tmp_1_w32; // Q(14-zeros)
    if (tmp_2_w32 > 0) {
      // If the exponential parameter is negative force it to zero, which means a
      // flat spectrum.
      tmp_1_w32 = WebRtcSpl_DivW32W16(tmp_2_w32, matrix_determinant); // Q14
      inst->pinkNoiseExp += WEBRTC_SPL_SAT(16384, tmp_1_w32, 0); // Q14
    }
  }
}

void WebRtcNsx_DataSynthesis(NoiseSuppressionFixedC* inst, short* outFrame) {
  int32_t energyOut;

  int16_t realImag_buff[ANAL_BLOCKL_MAX * 2 + 16];
  int16_t rfft_out_buff[ANAL_BLOCKL_MAX * 2 + 16];

  // Align the structures to 32-byte boundary for the FFT function.
  int16_t* realImag = (int16_t*) (((uintptr_t)realImag_buff + 31) & ~31);
  int16_t* rfft_out = (int16_t*) (((uintptr_t) rfft_out_buff + 31) & ~31);

  int16_t tmp16no1, tmp16no2;
  int16_t energyRatio;
  int16_t gainFactor, gainFactor1, gainFactor2;

  size_t i;
  int outCIFFT;
  int scaleEnergyOut = 0;

  if (inst->zeroInputSignal) {
    // synthesize the special case of zero input
    // read out fully processed segment
    for (i = 0; i < inst->blockLen10ms; i++) {
      outFrame[i] = inst->synthesisBuffer[i]; // Q0
    }
    // update synthesis buffer
    memcpy(inst->synthesisBuffer, inst->synthesisBuffer + inst->blockLen10ms,
        (inst->anaLen - inst->blockLen10ms) * sizeof(*inst->synthesisBuffer));
    WebRtcSpl_ZerosArrayW16(inst->synthesisBuffer + inst->anaLen - inst->blockLen10ms,
                            inst->blockLen10ms);
    return;
  }

  // Filter the data in the frequency domain, and create spectrum.
  WebRtcNsx_PrepareSpectrum(inst, realImag);

  // Inverse FFT output will be in rfft_out[].
  outCIFFT = WebRtcSpl_RealInverseFFT(inst->real_fft, realImag, rfft_out);

  WebRtcNsx_Denormalize(inst, rfft_out, outCIFFT);

  //scale factor: only do it after END_STARTUP_LONG time
  gainFactor = 8192; // 8192 = Q13(1.0)
  if (inst->gainMap == 1 &&
      inst->blockIndex > END_STARTUP_LONG &&
      inst->energyIn > 0) {
    // Q(-scaleEnergyOut)
    energyOut = WebRtcSpl_Energy(inst->real, inst->anaLen, &scaleEnergyOut);
    if (scaleEnergyOut == 0 && !(energyOut & 0x7f800000)) {
      energyOut = WEBRTC_SPL_SHIFT_W32(energyOut, 8 + scaleEnergyOut
                                       - inst->scaleEnergyIn);
    } else {
      // |energyIn| is currently in Q(|scaleEnergyIn|), but to later on end up
      // with an |energyRatio| in Q8 we need to change the Q-domain to
      // Q(-8-scaleEnergyOut).
      inst->energyIn >>= 8 + scaleEnergyOut - inst->scaleEnergyIn;
    }

    assert(inst->energyIn > 0);
    energyRatio = (energyOut + inst->energyIn / 2) / inst->energyIn;  // Q8
    // Limit the ratio to [0, 1] in Q8, i.e., [0, 256]
    energyRatio = WEBRTC_SPL_SAT(256, energyRatio, 0);

    // all done in lookup tables now
    assert(energyRatio < 257);
    gainFactor1 = kFactor1Table[energyRatio]; // Q8
    gainFactor2 = inst->factor2Table[energyRatio]; // Q8

    //combine both scales with speech/noise prob: note prior (priorSpeechProb) is not frequency dependent

    // factor = inst->priorSpeechProb*factor1 + (1.0-inst->priorSpeechProb)*factor2; // original code
    tmp16no1 = (int16_t)(((16384 - inst->priorNonSpeechProb) * gainFactor1) >>
        14);  // in Q13, where 16384 = Q14(1.0)
    tmp16no2 = (int16_t)((inst->priorNonSpeechProb * gainFactor2) >> 14);
    gainFactor = tmp16no1 + tmp16no2; // Q13
  }  // out of flag_gain_map==1

  // Synthesis, read out fully processed segment, and update synthesis buffer.
  WebRtcNsx_SynthesisUpdate(inst, outFrame, gainFactor);
}

void WebRtcNsx_ProcessCore(NoiseSuppressionFixedC* inst,
                           const short* const* speechFrame,
                           int num_bands,
                           short* const* outFrame) {
  // main routine for noise suppression

  uint32_t tmpU32no1, tmpU32no2, tmpU32no3;
  uint32_t satMax, maxNoiseU32;
  uint32_t tmpMagnU32, tmpNoiseU32;
  uint32_t nearMagnEst;
  uint32_t noiseUpdateU32;
  uint32_t noiseU32[HALF_ANAL_BLOCKL];
  uint32_t postLocSnr[HALF_ANAL_BLOCKL];
  uint32_t priorLocSnr[HALF_ANAL_BLOCKL];
  uint32_t prevNearSnr[HALF_ANAL_BLOCKL];
  uint32_t curNearSnr;
  uint32_t priorSnr;
  uint32_t noise_estimate = 0;
  uint32_t noise_estimate_avg = 0;
  uint32_t numerator = 0;

  int32_t tmp32no1, tmp32no2;
  int32_t pink_noise_num_avg = 0;

  uint16_t tmpU16no1;
  uint16_t magnU16[HALF_ANAL_BLOCKL];
  uint16_t prevNoiseU16[HALF_ANAL_BLOCKL];
  uint16_t nonSpeechProbFinal[HALF_ANAL_BLOCKL];
  uint16_t gammaNoise, prevGammaNoise;
  uint16_t noiseSupFilterTmp[HALF_ANAL_BLOCKL];

  int16_t qMagn, qNoise;
  int16_t avgProbSpeechHB, gainModHB, avgFilterGainHB, gainTimeDomainHB;
  int16_t pink_noise_exp_avg = 0;

  size_t i, j;
  int nShifts, postShifts;
  int norm32no1, norm32no2;
  int flag, sign;
  int q_domain_to_use = 0;

  // Code for ARMv7-Neon platform assumes the following:
  assert(inst->anaLen > 0);
  assert(inst->anaLen2 > 0);
  assert(inst->anaLen % 16 == 0);
  assert(inst->anaLen2 % 8 == 0);
  assert(inst->blockLen10ms > 0);
  assert(inst->blockLen10ms % 16 == 0);
  assert(inst->magnLen == inst->anaLen2 + 1);

#ifdef NS_FILEDEBUG
  if (fwrite(spframe, sizeof(short),
             inst->blockLen10ms, inst->infile) != inst->blockLen10ms) {
    assert(false);
  }
#endif

  // Check that initialization has been done
  assert(inst->initFlag == 1);
  assert((num_bands - 1) <= NUM_HIGH_BANDS_MAX);

  const short* const* speechFrameHB = NULL;
  short* const* outFrameHB = NULL;
  size_t num_high_bands = 0;
  if (num_bands > 1) {
    speechFrameHB = &speechFrame[1];
    outFrameHB = &outFrame[1];
    num_high_bands = (size_t)(num_bands - 1);
  }

  // Store speechFrame and transform to frequency domain
  WebRtcNsx_DataAnalysis(inst, (short*)speechFrame[0], magnU16);

  if (inst->zeroInputSignal) {
    WebRtcNsx_DataSynthesis(inst, outFrame[0]);

    if (num_bands > 1) {
      // update analysis buffer for H band
      // append new data to buffer FX
      for (i = 0; i < num_high_bands; ++i) {
        int block_shift = inst->anaLen - inst->blockLen10ms;
        memcpy(inst->dataBufHBFX[i], inst->dataBufHBFX[i] + inst->blockLen10ms,
            block_shift * sizeof(*inst->dataBufHBFX[i]));
        memcpy(inst->dataBufHBFX[i] + block_shift, speechFrameHB[i],
            inst->blockLen10ms * sizeof(*inst->dataBufHBFX[i]));
        for (j = 0; j < inst->blockLen10ms; j++) {
          outFrameHB[i][j] = inst->dataBufHBFX[i][j]; // Q0
        }
      }
    }  // end of H band gain computation
    return;
  }

  // Update block index when we have something to process
  inst->blockIndex++;
  //

  // Norm of magn
  qMagn = inst->normData - inst->stages;

  // Compute spectral flatness on input spectrum
  WebRtcNsx_ComputeSpectralFlatness(inst, magnU16);

  // quantile noise estimate
  WebRtcNsx_NoiseEstimation(inst, magnU16, noiseU32, &qNoise);

  //noise estimate from previous frame
  for (i = 0; i < inst->magnLen; i++) {
    prevNoiseU16[i] = (uint16_t)(inst->prevNoiseU32[i] >> 11);  // Q(prevQNoise)
  }

  if (inst->blockIndex < END_STARTUP_SHORT) {
    // Noise Q-domain to be used later; see description at end of section.
    q_domain_to_use = WEBRTC_SPL_MIN((int)qNoise, inst->minNorm - inst->stages);

    // Calculate frequency independent parts in parametric noise estimate and calculate
    // the estimate for the lower frequency band (same values for all frequency bins)
    if (inst->pinkNoiseExp) {
      pink_noise_exp_avg = (int16_t)WebRtcSpl_DivW32W16(inst->pinkNoiseExp,
                                                              (int16_t)(inst->blockIndex + 1)); // Q14
      pink_noise_num_avg = WebRtcSpl_DivW32W16(inst->pinkNoiseNumerator,
                                               (int16_t)(inst->blockIndex + 1)); // Q11
      WebRtcNsx_CalcParametricNoiseEstimate(inst,
                                            pink_noise_exp_avg,
                                            pink_noise_num_avg,
                                            kStartBand,
                                            &noise_estimate,
                                            &noise_estimate_avg);
    } else {
      // Use white noise estimate if we have poor pink noise parameter estimates
      noise_estimate = inst->whiteNoiseLevel; // Q(minNorm-stages)
      noise_estimate_avg = noise_estimate / (inst->blockIndex + 1); // Q(minNorm-stages)
    }
    for (i = 0; i < inst->magnLen; i++) {
      // Estimate the background noise using the pink noise parameters if permitted
      if ((inst->pinkNoiseExp) && (i >= kStartBand)) {
        // Reset noise_estimate
        noise_estimate = 0;
        noise_estimate_avg = 0;
        // Calculate the parametric noise estimate for current frequency bin
        WebRtcNsx_CalcParametricNoiseEstimate(inst,
                                              pink_noise_exp_avg,
                                              pink_noise_num_avg,
                                              i,
                                              &noise_estimate,
                                              &noise_estimate_avg);
      }
      // Calculate parametric Wiener filter
      noiseSupFilterTmp[i] = inst->denoiseBound;
      if (inst->initMagnEst[i]) {
        // numerator = (initMagnEst - noise_estimate * overdrive)
        // Result in Q(8+minNorm-stages)
        tmpU32no1 = WEBRTC_SPL_UMUL_32_16(noise_estimate, inst->overdrive);
        numerator = inst->initMagnEst[i] << 8;
        if (numerator > tmpU32no1) {
          // Suppression filter coefficient larger than zero, so calculate.
          numerator -= tmpU32no1;

          // Determine number of left shifts in numerator for best accuracy after
          // division
          nShifts = WebRtcSpl_NormU32(numerator);
          nShifts = WEBRTC_SPL_SAT(6, nShifts, 0);

          // Shift numerator to Q(nShifts+8+minNorm-stages)
          numerator <<= nShifts;

          // Shift denominator to Q(nShifts-6+minNorm-stages)
          tmpU32no1 = inst->initMagnEst[i] >> (6 - nShifts);
          if (tmpU32no1 == 0) {
            // This is only possible if numerator = 0, in which case
            // we don't need any division.
            tmpU32no1 = 1;
          }
          tmpU32no2 = numerator / tmpU32no1;  // Q14
          noiseSupFilterTmp[i] = (uint16_t)WEBRTC_SPL_SAT(16384, tmpU32no2,
              (uint32_t)(inst->denoiseBound)); // Q14
        }
      }
      // Weight quantile noise 'noiseU32' with modeled noise 'noise_estimate_avg'
      // 'noiseU32 is in Q(qNoise) and 'noise_estimate' in Q(minNorm-stages)
      // To guarantee that we do not get wrap around when shifting to the same domain
      // we use the lowest one. Furthermore, we need to save 6 bits for the weighting.
      // 'noise_estimate_avg' can handle this operation by construction, but 'noiseU32'
      // may not.

      // Shift 'noiseU32' to 'q_domain_to_use'
      tmpU32no1 = noiseU32[i] >> (qNoise - q_domain_to_use);
      // Shift 'noise_estimate_avg' to 'q_domain_to_use'
      tmpU32no2 = noise_estimate_avg >>
          (inst->minNorm - inst->stages - q_domain_to_use);
      // Make a simple check to see if we have enough room for weighting 'tmpU32no1'
      // without wrap around
      nShifts = 0;
      if (tmpU32no1 & 0xfc000000) {
        tmpU32no1 >>= 6;
        tmpU32no2 >>= 6;
        nShifts = 6;
      }
      tmpU32no1 *= inst->blockIndex;
      tmpU32no2 *= (END_STARTUP_SHORT - inst->blockIndex);
      // Add them together and divide by startup length
      noiseU32[i] = WebRtcSpl_DivU32U16(tmpU32no1 + tmpU32no2, END_STARTUP_SHORT);
      // Shift back if necessary
      noiseU32[i] <<= nShifts;
    }
    // Update new Q-domain for 'noiseU32'
    qNoise = q_domain_to_use;
  }
  // compute average signal during END_STARTUP_LONG time:
  // used to normalize spectral difference measure
  if (inst->blockIndex < END_STARTUP_LONG) {
    // substituting division with shift ending up in Q(-2*stages)
    inst->timeAvgMagnEnergyTmp +=
        inst->magnEnergy >> (2 * inst->normData + inst->stages - 1);
    inst->timeAvgMagnEnergy = WebRtcSpl_DivU32U16(inst->timeAvgMagnEnergyTmp,
                                                  inst->blockIndex + 1);
  }

  //start processing at frames == converged+1
  // STEP 1: compute prior and post SNR based on quantile noise estimates

  // compute direct decision (DD) estimate of prior SNR: needed for new method
  satMax = (uint32_t)1048575;// Largest possible value without getting overflow despite shifting 12 steps
  postShifts = 6 + qMagn - qNoise;
  nShifts = 5 - inst->prevQMagn + inst->prevQNoise;
  for (i = 0; i < inst->magnLen; i++) {
    // FLOAT:
    // post SNR
    // postLocSnr[i] = 0.0;
    // if (magn[i] > noise[i])
    // {
    //   postLocSnr[i] = magn[i] / (noise[i] + 0.0001);
    // }
    // // previous post SNR
    // // previous estimate: based on previous frame with gain filter (smooth is previous filter)
    //
    // prevNearSnr[i] = inst->prevMagnU16[i] / (inst->noisePrev[i] + 0.0001) * (inst->smooth[i]);
    //
    // // DD estimate is sum of two terms: current estimate and previous estimate
    // // directed decision update of priorSnr (or we actually store [2*priorSnr+1])
    //
    // priorLocSnr[i] = DD_PR_SNR * prevNearSnr[i] + (1.0 - DD_PR_SNR) * (postLocSnr[i] - 1.0);

    // calculate post SNR: output in Q11
    postLocSnr[i] = 2048; // 1.0 in Q11
    tmpU32no1 = (uint32_t)magnU16[i] << 6;  // Q(6+qMagn)
    if (postShifts < 0) {
      tmpU32no2 = noiseU32[i] >> -postShifts;  // Q(6+qMagn)
    } else {
      tmpU32no2 = noiseU32[i] << postShifts;  // Q(6+qMagn)
    }
    if (tmpU32no1 > tmpU32no2) {
      // Current magnitude larger than noise
      tmpU32no1 <<= 11;  // Q(17+qMagn)
      if (tmpU32no2 > 0) {
        tmpU32no1 /= tmpU32no2;  // Q11
        postLocSnr[i] = WEBRTC_SPL_MIN(satMax, tmpU32no1); // Q11
      } else {
        postLocSnr[i] = satMax;
      }
    }

    // calculate prevNearSnr[i] and save for later instead of recalculating it later
    // |nearMagnEst| in Q(prevQMagn + 14)
    nearMagnEst = inst->prevMagnU16[i] * inst->noiseSupFilter[i];
    tmpU32no1 = nearMagnEst << 3;  // Q(prevQMagn+17)
    tmpU32no2 = inst->prevNoiseU32[i] >> nShifts;  // Q(prevQMagn+6)

    if (tmpU32no2 > 0) {
      tmpU32no1 /= tmpU32no2;  // Q11
      tmpU32no1 = WEBRTC_SPL_MIN(satMax, tmpU32no1); // Q11
    } else {
      tmpU32no1 = satMax; // Q11
    }
    prevNearSnr[i] = tmpU32no1; // Q11

    //directed decision update of priorSnr
    tmpU32no1 = WEBRTC_SPL_UMUL_32_16(prevNearSnr[i], DD_PR_SNR_Q11); // Q22
    tmpU32no2 = WEBRTC_SPL_UMUL_32_16(postLocSnr[i] - 2048, ONE_MINUS_DD_PR_SNR_Q11); // Q22
    priorSnr = tmpU32no1 + tmpU32no2 + 512; // Q22 (added 512 for rounding)
    // priorLocSnr = 1 + 2*priorSnr
    priorLocSnr[i] = 2048 + (priorSnr >> 10);  // Q11
  }  // end of loop over frequencies
  // done with step 1: DD computation of prior and post SNR

  // STEP 2: compute speech/noise likelihood

  //compute difference of input spectrum with learned/estimated noise spectrum
  WebRtcNsx_ComputeSpectralDifference(inst, magnU16);
  //compute histograms for determination of parameters (thresholds and weights for features)
  //parameters are extracted once every window time (=inst->modelUpdate)
  //counter update
  inst->cntThresUpdate++;
  flag = (int)(inst->cntThresUpdate == inst->modelUpdate);
  //update histogram
  WebRtcNsx_FeatureParameterExtraction(inst, flag);
  //compute model parameters
  if (flag) {
    inst->cntThresUpdate = 0; // Reset counter
    //update every window:
    // get normalization for spectral difference for next window estimate

    // Shift to Q(-2*stages)
    inst->curAvgMagnEnergy >>= STAT_UPDATES;

    tmpU32no1 = (inst->curAvgMagnEnergy + inst->timeAvgMagnEnergy + 1) >> 1; //Q(-2*stages)
    // Update featureSpecDiff
    if ((tmpU32no1 != inst->timeAvgMagnEnergy) && (inst->featureSpecDiff) &&
        (inst->timeAvgMagnEnergy > 0)) {
      norm32no1 = 0;
      tmpU32no3 = tmpU32no1;
      while (0xFFFF0000 & tmpU32no3) {
        tmpU32no3 >>= 1;
        norm32no1++;
      }
      tmpU32no2 = inst->featureSpecDiff;
      while (0xFFFF0000 & tmpU32no2) {
        tmpU32no2 >>= 1;
        norm32no1++;
      }
      tmpU32no3 = WEBRTC_SPL_UMUL(tmpU32no3, tmpU32no2);
      tmpU32no3 /= inst->timeAvgMagnEnergy;
      if (WebRtcSpl_NormU32(tmpU32no3) < norm32no1) {
        inst->featureSpecDiff = 0x007FFFFF;
      } else {
        inst->featureSpecDiff = WEBRTC_SPL_MIN(0x007FFFFF,
                                               tmpU32no3 << norm32no1);
      }
    }

    inst->timeAvgMagnEnergy = tmpU32no1; // Q(-2*stages)
    inst->curAvgMagnEnergy = 0;
  }

  //compute speech/noise probability
  WebRtcNsx_SpeechNoiseProb(inst, nonSpeechProbFinal, priorLocSnr, postLocSnr);

  //time-avg parameter for noise update
  gammaNoise = NOISE_UPDATE_Q8; // Q8

  maxNoiseU32 = 0;
  postShifts = inst->prevQNoise - qMagn;
  nShifts = inst->prevQMagn - qMagn;
  for (i = 0; i < inst->magnLen; i++) {
    // temporary noise update: use it for speech frames if update value is less than previous
    // the formula has been rewritten into:
    // noiseUpdate = noisePrev[i] + (1 - gammaNoise) * nonSpeechProb * (magn[i] - noisePrev[i])

    if (postShifts < 0) {
      tmpU32no2 = magnU16[i] >> -postShifts;  // Q(prevQNoise)
    } else {
      tmpU32no2 = (uint32_t)magnU16[i] << postShifts;  // Q(prevQNoise)
    }
    if (prevNoiseU16[i] > tmpU32no2) {
      sign = -1;
      tmpU32no1 = prevNoiseU16[i] - tmpU32no2;
    } else {
      sign = 1;
      tmpU32no1 = tmpU32no2 - prevNoiseU16[i];
    }
    noiseUpdateU32 = inst->prevNoiseU32[i]; // Q(prevQNoise+11)
    tmpU32no3 = 0;
    if ((tmpU32no1) && (nonSpeechProbFinal[i])) {
      // This value will be used later, if gammaNoise changes
      tmpU32no3 = WEBRTC_SPL_UMUL_32_16(tmpU32no1, nonSpeechProbFinal[i]); // Q(prevQNoise+8)
      if (0x7c000000 & tmpU32no3) {
        // Shifting required before multiplication
        tmpU32no2 = (tmpU32no3 >> 5) * gammaNoise;  // Q(prevQNoise+11)
      } else {
        // We can do shifting after multiplication
        tmpU32no2 = (tmpU32no3 * gammaNoise) >> 5;  // Q(prevQNoise+11)
      }
      if (sign > 0) {
        noiseUpdateU32 += tmpU32no2; // Q(prevQNoise+11)
      } else {
        // This operation is safe. We can never get wrap around, since worst
        // case scenario means magnU16 = 0
        noiseUpdateU32 -= tmpU32no2; // Q(prevQNoise+11)
      }
    }

    //increase gamma (i.e., less noise update) for frame likely to be speech
    prevGammaNoise = gammaNoise;
    gammaNoise = NOISE_UPDATE_Q8;
    //time-constant based on speech/noise state
    //increase gamma (i.e., less noise update) for frames likely to be speech
    if (nonSpeechProbFinal[i] < ONE_MINUS_PROB_RANGE_Q8) {
      gammaNoise = GAMMA_NOISE_TRANS_AND_SPEECH_Q8;
    }

    if (prevGammaNoise != gammaNoise) {
      // new noise update
      // this line is the same as above, only that the result is stored in a different variable and the gammaNoise
      // has changed
      //
      // noiseUpdate = noisePrev[i] + (1 - gammaNoise) * nonSpeechProb * (magn[i] - noisePrev[i])

      if (0x7c000000 & tmpU32no3) {
        // Shifting required before multiplication
        tmpU32no2 = (tmpU32no3 >> 5) * gammaNoise;  // Q(prevQNoise+11)
      } else {
        // We can do shifting after multiplication
        tmpU32no2 = (tmpU32no3 * gammaNoise) >> 5;  // Q(prevQNoise+11)
      }
      if (sign > 0) {
        tmpU32no1 = inst->prevNoiseU32[i] + tmpU32no2; // Q(prevQNoise+11)
      } else {
        tmpU32no1 = inst->prevNoiseU32[i] - tmpU32no2; // Q(prevQNoise+11)
      }
      if (noiseUpdateU32 > tmpU32no1) {
        noiseUpdateU32 = tmpU32no1; // Q(prevQNoise+11)
      }
    }
    noiseU32[i] = noiseUpdateU32; // Q(prevQNoise+11)
    if (noiseUpdateU32 > maxNoiseU32) {
      maxNoiseU32 = noiseUpdateU32;
    }

    // conservative noise update
    // // original FLOAT code
    // if (prob_speech < PROB_RANGE) {
    // inst->avgMagnPause[i] = inst->avgMagnPause[i] + (1.0 - gamma_pause)*(magn[i] - inst->avgMagnPause[i]);
    // }

    tmp32no2 = WEBRTC_SPL_SHIFT_W32(inst->avgMagnPause[i], -nShifts);
    if (nonSpeechProbFinal[i] > ONE_MINUS_PROB_RANGE_Q8) {
      if (nShifts < 0) {
        tmp32no1 = (int32_t)magnU16[i] - tmp32no2; // Q(qMagn)
        tmp32no1 *= ONE_MINUS_GAMMA_PAUSE_Q8;  // Q(8+prevQMagn+nShifts)
        tmp32no1 = (tmp32no1 + 128) >> 8;  // Q(qMagn).
      } else {
        // In Q(qMagn+nShifts)
        tmp32no1 = ((int32_t)magnU16[i] << nShifts) - inst->avgMagnPause[i];
        tmp32no1 *= ONE_MINUS_GAMMA_PAUSE_Q8;  // Q(8+prevQMagn+nShifts)
        tmp32no1 = (tmp32no1 + (128 << nShifts)) >> (8 + nShifts);  // Q(qMagn).
      }
      tmp32no2 += tmp32no1; // Q(qMagn)
    }
    inst->avgMagnPause[i] = tmp32no2;
  }  // end of frequency loop

  norm32no1 = WebRtcSpl_NormU32(maxNoiseU32);
  qNoise = inst->prevQNoise + norm32no1 - 5;
  // done with step 2: noise update

  // STEP 3: compute dd update of prior snr and post snr based on new noise estimate
  nShifts = inst->prevQNoise + 11 - qMagn;
  for (i = 0; i < inst->magnLen; i++) {
    // FLOAT code
    // // post and prior SNR
    // curNearSnr = 0.0;
    // if (magn[i] > noise[i])
    // {
    // curNearSnr = magn[i] / (noise[i] + 0.0001) - 1.0;
    // }
    // // DD estimate is sum of two terms: current estimate and previous estimate
    // // directed decision update of snrPrior
    // snrPrior = DD_PR_SNR * prevNearSnr[i] + (1.0 - DD_PR_SNR) * curNearSnr;
    // // gain filter
    // tmpFloat1 = inst->overdrive + snrPrior;
    // tmpFloat2 = snrPrior / tmpFloat1;
    // theFilter[i] = tmpFloat2;

    // calculate curNearSnr again, this is necessary because a new noise estimate has been made since then. for the original
    curNearSnr = 0; // Q11
    if (nShifts < 0) {
      // This case is equivalent with magn < noise which implies curNearSnr = 0;
      tmpMagnU32 = (uint32_t)magnU16[i]; // Q(qMagn)
      tmpNoiseU32 = noiseU32[i] << -nShifts;  // Q(qMagn)
    } else if (nShifts > 17) {
      tmpMagnU32 = (uint32_t)magnU16[i] << 17;  // Q(qMagn+17)
      tmpNoiseU32 = noiseU32[i] >> (nShifts - 17);  // Q(qMagn+17)
    } else {
      tmpMagnU32 = (uint32_t)magnU16[i] << nShifts;  // Q(qNoise_prev+11)
      tmpNoiseU32 = noiseU32[i]; // Q(qNoise_prev+11)
    }
    if (tmpMagnU32 > tmpNoiseU32) {
      tmpU32no1 = tmpMagnU32 - tmpNoiseU32; // Q(qCur)
      norm32no2 = WEBRTC_SPL_MIN(11, WebRtcSpl_NormU32(tmpU32no1));
      tmpU32no1 <<= norm32no2;  // Q(qCur+norm32no2)
      tmpU32no2 = tmpNoiseU32 >> (11 - norm32no2);  // Q(qCur+norm32no2-11)
      if (tmpU32no2 > 0) {
        tmpU32no1 /= tmpU32no2;  // Q11
      }
      curNearSnr = WEBRTC_SPL_MIN(satMax, tmpU32no1); // Q11
    }

    //directed decision update of priorSnr
    // FLOAT
    // priorSnr = DD_PR_SNR * prevNearSnr + (1.0-DD_PR_SNR) * curNearSnr;

    tmpU32no1 = WEBRTC_SPL_UMUL_32_16(prevNearSnr[i], DD_PR_SNR_Q11); // Q22
    tmpU32no2 = WEBRTC_SPL_UMUL_32_16(curNearSnr, ONE_MINUS_DD_PR_SNR_Q11); // Q22
    priorSnr = tmpU32no1 + tmpU32no2; // Q22

    //gain filter
    tmpU32no1 = inst->overdrive + ((priorSnr + 8192) >> 14);  // Q8
    assert(inst->overdrive > 0);
    tmpU16no1 = (priorSnr + tmpU32no1 / 2) / tmpU32no1;  // Q14
    inst->noiseSupFilter[i] = WEBRTC_SPL_SAT(16384, tmpU16no1, inst->denoiseBound); // 16384 = Q14(1.0) // Q14

    // Weight in the parametric Wiener filter during startup
    if (inst->blockIndex < END_STARTUP_SHORT) {
      // Weight the two suppression filters
      tmpU32no1 = inst->noiseSupFilter[i] * inst->blockIndex;
      tmpU32no2 = noiseSupFilterTmp[i] *
          (END_STARTUP_SHORT - inst->blockIndex);
      tmpU32no1 += tmpU32no2;
      inst->noiseSupFilter[i] = (uint16_t)WebRtcSpl_DivU32U16(tmpU32no1,
                                                                    END_STARTUP_SHORT);
    }
  }  // end of loop over frequencies
  //done with step3

  // save noise and magnitude spectrum for next frame
  inst->prevQNoise = qNoise;
  inst->prevQMagn = qMagn;
  if (norm32no1 > 5) {
    for (i = 0; i < inst->magnLen; i++) {
      inst->prevNoiseU32[i] = noiseU32[i] << (norm32no1 - 5);  // Q(qNoise+11)
      inst->prevMagnU16[i] = magnU16[i]; // Q(qMagn)
    }
  } else {
    for (i = 0; i < inst->magnLen; i++) {
      inst->prevNoiseU32[i] = noiseU32[i] >> (5 - norm32no1);  // Q(qNoise+11)
      inst->prevMagnU16[i] = magnU16[i]; // Q(qMagn)
    }
  }

  WebRtcNsx_DataSynthesis(inst, outFrame[0]);
#ifdef NS_FILEDEBUG
  if (fwrite(outframe, sizeof(short),
             inst->blockLen10ms, inst->outfile) != inst->blockLen10ms) {
    assert(false);
  }
#endif

  //for H band:
  // only update data buffer, then apply time-domain gain is applied derived from L band
  if (num_bands > 1) {
    // update analysis buffer for H band
    // append new data to buffer FX
    for (i = 0; i < num_high_bands; ++i) {
      memcpy(inst->dataBufHBFX[i], inst->dataBufHBFX[i] + inst->blockLen10ms,
          (inst->anaLen - inst->blockLen10ms) * sizeof(*inst->dataBufHBFX[i]));
      memcpy(inst->dataBufHBFX[i] + inst->anaLen - inst->blockLen10ms,
          speechFrameHB[i], inst->blockLen10ms * sizeof(*inst->dataBufHBFX[i]));
    }
    // range for averaging low band quantities for H band gain

    gainTimeDomainHB = 16384; // 16384 = Q14(1.0)
    //average speech prob from low band
    //average filter gain from low band
    //avg over second half (i.e., 4->8kHz) of freq. spectrum
    tmpU32no1 = 0; // Q12
    tmpU16no1 = 0; // Q8
    for (i = inst->anaLen2 - (inst->anaLen2 >> 2); i < inst->anaLen2; i++) {
      tmpU16no1 += nonSpeechProbFinal[i]; // Q8
      tmpU32no1 += (uint32_t)(inst->noiseSupFilter[i]); // Q14
    }
    assert(inst->stages >= 7);
    avgProbSpeechHB = (4096 - (tmpU16no1 >> (inst->stages - 7)));  // Q12
    avgFilterGainHB = (int16_t)(tmpU32no1 >> (inst->stages - 3));  // Q14

    // // original FLOAT code
    // // gain based on speech probability:
    // avg_prob_speech_tt=(float)2.0*avg_prob_speech-(float)1.0;
    // gain_mod=(float)0.5*((float)1.0+(float)tanh(avg_prob_speech_tt)); // between 0 and 1

    // gain based on speech probability:
    // original expression: "0.5 * (1 + tanh(2x-1))"
    // avgProbSpeechHB has been anyway saturated to a value between 0 and 1 so the other cases don't have to be dealt with
    // avgProbSpeechHB and gainModHB are in Q12, 3607 = Q12(0.880615234375) which is a zero point of
    // |0.5 * (1 + tanh(2x-1)) - x| - |0.5 * (1 + tanh(2x-1)) - 0.880615234375| meaning that from that point the error of approximating
    // the expression with f(x) = x would be greater than the error of approximating the expression with f(x) = 0.880615234375
    // 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
    // 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
    gainModHB = WEBRTC_SPL_MIN(avgProbSpeechHB, 3607);

    // // original FLOAT code
    // //combine gain with low band gain
    // if (avg_prob_speech < (float)0.5) {
    // gain_time_domain_HB=(float)0.5*gain_mod+(float)0.5*avg_filter_gain;
    // }
    // else {
    // gain_time_domain_HB=(float)0.25*gain_mod+(float)0.75*avg_filter_gain;
    // }


    //combine gain with low band gain
    if (avgProbSpeechHB < 2048) {
      // 2048 = Q12(0.5)
      // 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
      gainTimeDomainHB = (gainModHB << 1) + (avgFilterGainHB >> 1); // Q14
    } else {
      // "gain_time_domain = 0.25 * gain_mod + 0.75 * agv_filter_gain;"
      gainTimeDomainHB = (int16_t)((3 * avgFilterGainHB) >> 2);  // 3 = Q2(0.75)
      gainTimeDomainHB += gainModHB; // Q14
    }
    //make sure gain is within flooring range
    gainTimeDomainHB
      = WEBRTC_SPL_SAT(16384, gainTimeDomainHB, (int16_t)(inst->denoiseBound)); // 16384 = Q14(1.0)


    //apply gain
    for (i = 0; i < num_high_bands; ++i) {
      for (j = 0; j < inst->blockLen10ms; j++) {
        outFrameHB[i][j] = (int16_t)((gainTimeDomainHB *
            inst->dataBufHBFX[i][j]) >> 14);  // Q0
      }
    }
  }  // end of H band gain computation
}