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