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