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