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) 23static const int kStartBand = 5; 24 25// Constants to compensate for shifting signal log(2^shifts). 26const WebRtc_Word16 WebRtcNsx_kLogTable[9] = { 27 0, 177, 355, 532, 710, 887, 1065, 1242, 1420 28}; 29 30const 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 46const 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 65static 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 196static 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 202static 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 218static 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// } 263static 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 298static 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 321static 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 344static 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 368static 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 382static 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 396static 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 418static 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. 431NoiseEstimation WebRtcNsx_NoiseEstimation; 432PrepareSpectrum WebRtcNsx_PrepareSpectrum; 433SynthesisUpdate WebRtcNsx_SynthesisUpdate; 434AnalysisUpdate WebRtcNsx_AnalysisUpdate; 435Denormalize WebRtcNsx_Denormalize; 436CreateComplexBuffer WebRtcNsx_CreateComplexBuffer; 437 438// Update the noise estimation information. 439static 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 471static 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. 594static 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. 619static __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. 631static 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. 664static 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. 685static __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 695void 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 741WebRtc_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 892int 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 927void 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 1127void 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 1197void 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) 1292void 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) 1522void 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 1768void 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 1844int 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