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