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. */
24extern const int16_t WebRtcNsx_kLogTable[9];
25extern const int16_t WebRtcNsx_kCounterDiv[201];
26extern const int16_t WebRtcNsx_kLogTableFrac[256];
27#else
28static const int16_t WebRtcNsx_kLogTable[9] = {
29  0, 177, 355, 532, 710, 887, 1065, 1242, 1420
30};
31
32static 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
49static 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)
71static const size_t kStartBand = 5;
72
73// hybrib Hanning & flat window
74static 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
90static 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//  }
135static 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
170static 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
193static 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
216static 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
240static 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
254static 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
268static 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
290static 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.
303static 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
334static 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.
456static 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.
477static 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.
491static 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.
523static 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.
542static 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.
553NoiseEstimation WebRtcNsx_NoiseEstimation;
554PrepareSpectrum WebRtcNsx_PrepareSpectrum;
555SynthesisUpdate WebRtcNsx_SynthesisUpdate;
556AnalysisUpdate WebRtcNsx_AnalysisUpdate;
557Denormalize WebRtcNsx_Denormalize;
558NormalizeRealBuffer WebRtcNsx_NormalizeRealBuffer;
559
560#if (defined WEBRTC_DETECT_NEON || defined WEBRTC_HAS_NEON)
561// Initialize function pointers for ARM Neon platform.
562static 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.
572static 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
583void 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
628int32_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
783int 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
818void 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
1019void 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
1088void 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)
1182void 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
1420void 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
1500void 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