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