1 /*
2 * Copyright (c) 2013 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 <assert.h>
12
13 #include "webrtc/modules/audio_processing/ns/include/noise_suppression_x.h"
14 #include "webrtc/modules/audio_processing/ns/nsx_core.h"
15 #include "webrtc/modules/audio_processing/ns/nsx_defines.h"
16
17 static const int16_t kIndicatorTable[17] = {
18 0, 2017, 3809, 5227, 6258, 6963, 7424, 7718,
19 7901, 8014, 8084, 8126, 8152, 8168, 8177, 8183, 8187
20 };
21
22 // Compute speech/noise probability
23 // speech/noise probability is returned in: probSpeechFinal
24 //snrLocPrior is the prior SNR for each frequency (in Q11)
25 //snrLocPost is the post SNR for each frequency (in Q11)
WebRtcNsx_SpeechNoiseProb(NoiseSuppressionFixedC * inst,uint16_t * nonSpeechProbFinal,uint32_t * priorLocSnr,uint32_t * postLocSnr)26 void WebRtcNsx_SpeechNoiseProb(NoiseSuppressionFixedC* inst,
27 uint16_t* nonSpeechProbFinal,
28 uint32_t* priorLocSnr,
29 uint32_t* postLocSnr) {
30 uint32_t zeros, num, den, tmpU32no1, tmpU32no2, tmpU32no3;
31 int32_t invLrtFX, indPriorFX, tmp32, tmp32no1, tmp32no2, besselTmpFX32;
32 int32_t frac32, logTmp;
33 int32_t logLrtTimeAvgKsumFX;
34 int16_t indPriorFX16;
35 int16_t tmp16, tmp16no1, tmp16no2, tmpIndFX, tableIndex, frac, intPart;
36 size_t i;
37 int normTmp, normTmp2, nShifts;
38
39 // compute feature based on average LR factor
40 // this is the average over all frequencies of the smooth log LRT
41 logLrtTimeAvgKsumFX = 0;
42 for (i = 0; i < inst->magnLen; i++) {
43 besselTmpFX32 = (int32_t)postLocSnr[i]; // Q11
44 normTmp = WebRtcSpl_NormU32(postLocSnr[i]);
45 num = postLocSnr[i] << normTmp; // Q(11+normTmp)
46 if (normTmp > 10) {
47 den = priorLocSnr[i] << (normTmp - 11); // Q(normTmp)
48 } else {
49 den = priorLocSnr[i] >> (11 - normTmp); // Q(normTmp)
50 }
51 if (den > 0) {
52 besselTmpFX32 -= num / den; // Q11
53 } else {
54 besselTmpFX32 = 0;
55 }
56
57 // inst->logLrtTimeAvg[i] += LRT_TAVG * (besselTmp - log(snrLocPrior)
58 // - inst->logLrtTimeAvg[i]);
59 // Here, LRT_TAVG = 0.5
60 zeros = WebRtcSpl_NormU32(priorLocSnr[i]);
61 frac32 = (int32_t)(((priorLocSnr[i] << zeros) & 0x7FFFFFFF) >> 19);
62 tmp32 = (frac32 * frac32 * -43) >> 19;
63 tmp32 += ((int16_t)frac32 * 5412) >> 12;
64 frac32 = tmp32 + 37;
65 // tmp32 = log2(priorLocSnr[i])
66 tmp32 = (int32_t)(((31 - zeros) << 12) + frac32) - (11 << 12); // Q12
67 logTmp = (tmp32 * 178) >> 8; // log2(priorLocSnr[i])*log(2)
68 // tmp32no1 = LRT_TAVG * (log(snrLocPrior) + inst->logLrtTimeAvg[i]) in Q12.
69 tmp32no1 = (logTmp + inst->logLrtTimeAvgW32[i]) / 2;
70 inst->logLrtTimeAvgW32[i] += (besselTmpFX32 - tmp32no1); // Q12
71
72 logLrtTimeAvgKsumFX += inst->logLrtTimeAvgW32[i]; // Q12
73 }
74 inst->featureLogLrt = (logLrtTimeAvgKsumFX * BIN_SIZE_LRT) >>
75 (inst->stages + 11);
76
77 // done with computation of LR factor
78
79 //
80 //compute the indicator functions
81 //
82
83 // average LRT feature
84 // FLOAT code
85 // indicator0 = 0.5 * (tanh(widthPrior *
86 // (logLrtTimeAvgKsum - threshPrior0)) + 1.0);
87 tmpIndFX = 16384; // Q14(1.0)
88 tmp32no1 = logLrtTimeAvgKsumFX - inst->thresholdLogLrt; // Q12
89 nShifts = 7 - inst->stages; // WIDTH_PR_MAP_SHIFT - inst->stages + 5;
90 //use larger width in tanh map for pause regions
91 if (tmp32no1 < 0) {
92 tmpIndFX = 0;
93 tmp32no1 = -tmp32no1;
94 //widthPrior = widthPrior * 2.0;
95 nShifts++;
96 }
97 tmp32no1 = WEBRTC_SPL_SHIFT_W32(tmp32no1, nShifts); // Q14
98 // compute indicator function: sigmoid map
99 tableIndex = (int16_t)(tmp32no1 >> 14);
100 if ((tableIndex < 16) && (tableIndex >= 0)) {
101 tmp16no2 = kIndicatorTable[tableIndex];
102 tmp16no1 = kIndicatorTable[tableIndex + 1] - kIndicatorTable[tableIndex];
103 frac = (int16_t)(tmp32no1 & 0x00003fff); // Q14
104 tmp16no2 += (int16_t)((tmp16no1 * frac) >> 14);
105 if (tmpIndFX == 0) {
106 tmpIndFX = 8192 - tmp16no2; // Q14
107 } else {
108 tmpIndFX = 8192 + tmp16no2; // Q14
109 }
110 }
111 indPriorFX = inst->weightLogLrt * tmpIndFX; // 6*Q14
112
113 //spectral flatness feature
114 if (inst->weightSpecFlat) {
115 tmpU32no1 = WEBRTC_SPL_UMUL(inst->featureSpecFlat, 400); // Q10
116 tmpIndFX = 16384; // Q14(1.0)
117 //use larger width in tanh map for pause regions
118 tmpU32no2 = inst->thresholdSpecFlat - tmpU32no1; //Q10
119 nShifts = 4;
120 if (inst->thresholdSpecFlat < tmpU32no1) {
121 tmpIndFX = 0;
122 tmpU32no2 = tmpU32no1 - inst->thresholdSpecFlat;
123 //widthPrior = widthPrior * 2.0;
124 nShifts++;
125 }
126 tmpU32no1 = WebRtcSpl_DivU32U16(tmpU32no2 << nShifts, 25); // Q14
127 // compute indicator function: sigmoid map
128 // FLOAT code
129 // indicator1 = 0.5 * (tanh(sgnMap * widthPrior *
130 // (threshPrior1 - tmpFloat1)) + 1.0);
131 tableIndex = (int16_t)(tmpU32no1 >> 14);
132 if (tableIndex < 16) {
133 tmp16no2 = kIndicatorTable[tableIndex];
134 tmp16no1 = kIndicatorTable[tableIndex + 1] - kIndicatorTable[tableIndex];
135 frac = (int16_t)(tmpU32no1 & 0x00003fff); // Q14
136 tmp16no2 += (int16_t)((tmp16no1 * frac) >> 14);
137 if (tmpIndFX) {
138 tmpIndFX = 8192 + tmp16no2; // Q14
139 } else {
140 tmpIndFX = 8192 - tmp16no2; // Q14
141 }
142 }
143 indPriorFX += inst->weightSpecFlat * tmpIndFX; // 6*Q14
144 }
145
146 //for template spectral-difference
147 if (inst->weightSpecDiff) {
148 tmpU32no1 = 0;
149 if (inst->featureSpecDiff) {
150 normTmp = WEBRTC_SPL_MIN(20 - inst->stages,
151 WebRtcSpl_NormU32(inst->featureSpecDiff));
152 assert(normTmp >= 0);
153 tmpU32no1 = inst->featureSpecDiff << normTmp; // Q(normTmp-2*stages)
154 tmpU32no2 = inst->timeAvgMagnEnergy >> (20 - inst->stages - normTmp);
155 if (tmpU32no2 > 0) {
156 // Q(20 - inst->stages)
157 tmpU32no1 /= tmpU32no2;
158 } else {
159 tmpU32no1 = (uint32_t)(0x7fffffff);
160 }
161 }
162 tmpU32no3 = (inst->thresholdSpecDiff << 17) / 25;
163 tmpU32no2 = tmpU32no1 - tmpU32no3;
164 nShifts = 1;
165 tmpIndFX = 16384; // Q14(1.0)
166 //use larger width in tanh map for pause regions
167 if (tmpU32no2 & 0x80000000) {
168 tmpIndFX = 0;
169 tmpU32no2 = tmpU32no3 - tmpU32no1;
170 //widthPrior = widthPrior * 2.0;
171 nShifts--;
172 }
173 tmpU32no1 = tmpU32no2 >> nShifts;
174 // compute indicator function: sigmoid map
175 /* FLOAT code
176 indicator2 = 0.5 * (tanh(widthPrior * (tmpFloat1 - threshPrior2)) + 1.0);
177 */
178 tableIndex = (int16_t)(tmpU32no1 >> 14);
179 if (tableIndex < 16) {
180 tmp16no2 = kIndicatorTable[tableIndex];
181 tmp16no1 = kIndicatorTable[tableIndex + 1] - kIndicatorTable[tableIndex];
182 frac = (int16_t)(tmpU32no1 & 0x00003fff); // Q14
183 tmp16no2 += (int16_t)WEBRTC_SPL_MUL_16_16_RSFT_WITH_ROUND(
184 tmp16no1, frac, 14);
185 if (tmpIndFX) {
186 tmpIndFX = 8192 + tmp16no2;
187 } else {
188 tmpIndFX = 8192 - tmp16no2;
189 }
190 }
191 indPriorFX += inst->weightSpecDiff * tmpIndFX; // 6*Q14
192 }
193
194 //combine the indicator function with the feature weights
195 // FLOAT code
196 // indPrior = 1 - (weightIndPrior0 * indicator0 + weightIndPrior1 *
197 // indicator1 + weightIndPrior2 * indicator2);
198 indPriorFX16 = WebRtcSpl_DivW32W16ResW16(98307 - indPriorFX, 6); // Q14
199 // done with computing indicator function
200
201 //compute the prior probability
202 // FLOAT code
203 // inst->priorNonSpeechProb += PRIOR_UPDATE *
204 // (indPriorNonSpeech - inst->priorNonSpeechProb);
205 tmp16 = indPriorFX16 - inst->priorNonSpeechProb; // Q14
206 inst->priorNonSpeechProb += (int16_t)((PRIOR_UPDATE_Q14 * tmp16) >> 14);
207
208 //final speech probability: combine prior model with LR factor:
209
210 memset(nonSpeechProbFinal, 0, sizeof(uint16_t) * inst->magnLen);
211
212 if (inst->priorNonSpeechProb > 0) {
213 for (i = 0; i < inst->magnLen; i++) {
214 // FLOAT code
215 // invLrt = exp(inst->logLrtTimeAvg[i]);
216 // invLrt = inst->priorSpeechProb * invLrt;
217 // nonSpeechProbFinal[i] = (1.0 - inst->priorSpeechProb) /
218 // (1.0 - inst->priorSpeechProb + invLrt);
219 // invLrt = (1.0 - inst->priorNonSpeechProb) * invLrt;
220 // nonSpeechProbFinal[i] = inst->priorNonSpeechProb /
221 // (inst->priorNonSpeechProb + invLrt);
222 if (inst->logLrtTimeAvgW32[i] < 65300) {
223 tmp32no1 = (inst->logLrtTimeAvgW32[i] * 23637) >> 14; // Q12
224 intPart = (int16_t)(tmp32no1 >> 12);
225 if (intPart < -8) {
226 intPart = -8;
227 }
228 frac = (int16_t)(tmp32no1 & 0x00000fff); // Q12
229
230 // Quadratic approximation of 2^frac
231 tmp32no2 = (frac * frac * 44) >> 19; // Q12.
232 tmp32no2 += (frac * 84) >> 7; // Q12
233 invLrtFX = (1 << (8 + intPart)) +
234 WEBRTC_SPL_SHIFT_W32(tmp32no2, intPart - 4); // Q8
235
236 normTmp = WebRtcSpl_NormW32(invLrtFX);
237 normTmp2 = WebRtcSpl_NormW16((16384 - inst->priorNonSpeechProb));
238 if (normTmp + normTmp2 >= 7) {
239 if (normTmp + normTmp2 < 15) {
240 invLrtFX >>= 15 - normTmp2 - normTmp;
241 // Q(normTmp+normTmp2-7)
242 tmp32no1 = invLrtFX * (16384 - inst->priorNonSpeechProb);
243 // Q(normTmp+normTmp2+7)
244 invLrtFX = WEBRTC_SPL_SHIFT_W32(tmp32no1, 7 - normTmp - normTmp2);
245 // Q14
246 } else {
247 tmp32no1 = invLrtFX * (16384 - inst->priorNonSpeechProb);
248 // Q22
249 invLrtFX = tmp32no1 >> 8; // Q14.
250 }
251
252 tmp32no1 = (int32_t)inst->priorNonSpeechProb << 8; // Q22
253
254 nonSpeechProbFinal[i] = tmp32no1 /
255 (inst->priorNonSpeechProb + invLrtFX); // Q8
256 }
257 }
258 }
259 }
260 }
261
262