1 /*
2 * Copyright (c) 2011 The WebRTC project authors. All Rights Reserved.
3 *
4 * Use of this source code is governed by a BSD-style license
5 * that can be found in the LICENSE file in the root of the source
6 * tree. An additional intellectual property rights grant can be found
7 * in the file PATENTS. All contributing project authors may
8 * be found in the AUTHORS file in the root of the source tree.
9 */
10
11 /* digital_agc.c
12 *
13 */
14
15 #include "webrtc/modules/audio_processing/agc/legacy/digital_agc.h"
16
17 #include <assert.h>
18 #include <string.h>
19 #ifdef WEBRTC_AGC_DEBUG_DUMP
20 #include <stdio.h>
21 #endif
22
23 #include "webrtc/modules/audio_processing/agc/legacy/gain_control.h"
24
25 // To generate the gaintable, copy&paste the following lines to a Matlab window:
26 // MaxGain = 6; MinGain = 0; CompRatio = 3; Knee = 1;
27 // zeros = 0:31; lvl = 2.^(1-zeros);
28 // A = -10*log10(lvl) * (CompRatio - 1) / CompRatio;
29 // B = MaxGain - MinGain;
30 // gains = round(2^16*10.^(0.05 * (MinGain + B * ( log(exp(-Knee*A)+exp(-Knee*B)) - log(1+exp(-Knee*B)) ) / log(1/(1+exp(Knee*B))))));
31 // fprintf(1, '\t%i, %i, %i, %i,\n', gains);
32 // % Matlab code for plotting the gain and input/output level characteristic (copy/paste the following 3 lines):
33 // in = 10*log10(lvl); out = 20*log10(gains/65536);
34 // subplot(121); plot(in, out); axis([-30, 0, -5, 20]); grid on; xlabel('Input (dB)'); ylabel('Gain (dB)');
35 // subplot(122); plot(in, in+out); axis([-30, 0, -30, 5]); grid on; xlabel('Input (dB)'); ylabel('Output (dB)');
36 // zoom on;
37
38 // Generator table for y=log2(1+e^x) in Q8.
39 enum { kGenFuncTableSize = 128 };
40 static const uint16_t kGenFuncTable[kGenFuncTableSize] = {
41 256, 485, 786, 1126, 1484, 1849, 2217, 2586,
42 2955, 3324, 3693, 4063, 4432, 4801, 5171, 5540,
43 5909, 6279, 6648, 7017, 7387, 7756, 8125, 8495,
44 8864, 9233, 9603, 9972, 10341, 10711, 11080, 11449,
45 11819, 12188, 12557, 12927, 13296, 13665, 14035, 14404,
46 14773, 15143, 15512, 15881, 16251, 16620, 16989, 17359,
47 17728, 18097, 18466, 18836, 19205, 19574, 19944, 20313,
48 20682, 21052, 21421, 21790, 22160, 22529, 22898, 23268,
49 23637, 24006, 24376, 24745, 25114, 25484, 25853, 26222,
50 26592, 26961, 27330, 27700, 28069, 28438, 28808, 29177,
51 29546, 29916, 30285, 30654, 31024, 31393, 31762, 32132,
52 32501, 32870, 33240, 33609, 33978, 34348, 34717, 35086,
53 35456, 35825, 36194, 36564, 36933, 37302, 37672, 38041,
54 38410, 38780, 39149, 39518, 39888, 40257, 40626, 40996,
55 41365, 41734, 42104, 42473, 42842, 43212, 43581, 43950,
56 44320, 44689, 45058, 45428, 45797, 46166, 46536, 46905
57 };
58
59 static const int16_t kAvgDecayTime = 250; // frames; < 3000
60
WebRtcAgc_CalculateGainTable(int32_t * gainTable,int16_t digCompGaindB,int16_t targetLevelDbfs,uint8_t limiterEnable,int16_t analogTarget)61 int32_t WebRtcAgc_CalculateGainTable(int32_t *gainTable, // Q16
62 int16_t digCompGaindB, // Q0
63 int16_t targetLevelDbfs,// Q0
64 uint8_t limiterEnable,
65 int16_t analogTarget) // Q0
66 {
67 // This function generates the compressor gain table used in the fixed digital part.
68 uint32_t tmpU32no1, tmpU32no2, absInLevel, logApprox;
69 int32_t inLevel, limiterLvl;
70 int32_t tmp32, tmp32no1, tmp32no2, numFIX, den, y32;
71 const uint16_t kLog10 = 54426; // log2(10) in Q14
72 const uint16_t kLog10_2 = 49321; // 10*log10(2) in Q14
73 const uint16_t kLogE_1 = 23637; // log2(e) in Q14
74 uint16_t constMaxGain;
75 uint16_t tmpU16, intPart, fracPart;
76 const int16_t kCompRatio = 3;
77 const int16_t kSoftLimiterLeft = 1;
78 int16_t limiterOffset = 0; // Limiter offset
79 int16_t limiterIdx, limiterLvlX;
80 int16_t constLinApprox, zeroGainLvl, maxGain, diffGain;
81 int16_t i, tmp16, tmp16no1;
82 int zeros, zerosScale;
83
84 // Constants
85 // kLogE_1 = 23637; // log2(e) in Q14
86 // kLog10 = 54426; // log2(10) in Q14
87 // kLog10_2 = 49321; // 10*log10(2) in Q14
88
89 // Calculate maximum digital gain and zero gain level
90 tmp32no1 = (digCompGaindB - analogTarget) * (kCompRatio - 1);
91 tmp16no1 = analogTarget - targetLevelDbfs;
92 tmp16no1 += WebRtcSpl_DivW32W16ResW16(tmp32no1 + (kCompRatio >> 1), kCompRatio);
93 maxGain = WEBRTC_SPL_MAX(tmp16no1, (analogTarget - targetLevelDbfs));
94 tmp32no1 = maxGain * kCompRatio;
95 zeroGainLvl = digCompGaindB;
96 zeroGainLvl -= WebRtcSpl_DivW32W16ResW16(tmp32no1 + ((kCompRatio - 1) >> 1),
97 kCompRatio - 1);
98 if ((digCompGaindB <= analogTarget) && (limiterEnable))
99 {
100 zeroGainLvl += (analogTarget - digCompGaindB + kSoftLimiterLeft);
101 limiterOffset = 0;
102 }
103
104 // Calculate the difference between maximum gain and gain at 0dB0v:
105 // diffGain = maxGain + (compRatio-1)*zeroGainLvl/compRatio
106 // = (compRatio-1)*digCompGaindB/compRatio
107 tmp32no1 = digCompGaindB * (kCompRatio - 1);
108 diffGain = WebRtcSpl_DivW32W16ResW16(tmp32no1 + (kCompRatio >> 1), kCompRatio);
109 if (diffGain < 0 || diffGain >= kGenFuncTableSize)
110 {
111 assert(0);
112 return -1;
113 }
114
115 // Calculate the limiter level and index:
116 // limiterLvlX = analogTarget - limiterOffset
117 // limiterLvl = targetLevelDbfs + limiterOffset/compRatio
118 limiterLvlX = analogTarget - limiterOffset;
119 limiterIdx =
120 2 + WebRtcSpl_DivW32W16ResW16((int32_t)limiterLvlX << 13, kLog10_2 / 2);
121 tmp16no1 = WebRtcSpl_DivW32W16ResW16(limiterOffset + (kCompRatio >> 1), kCompRatio);
122 limiterLvl = targetLevelDbfs + tmp16no1;
123
124 // Calculate (through table lookup):
125 // constMaxGain = log2(1+2^(log2(e)*diffGain)); (in Q8)
126 constMaxGain = kGenFuncTable[diffGain]; // in Q8
127
128 // Calculate a parameter used to approximate the fractional part of 2^x with a
129 // piecewise linear function in Q14:
130 // constLinApprox = round(3/2*(4*(3-2*sqrt(2))/(log(2)^2)-0.5)*2^14);
131 constLinApprox = 22817; // in Q14
132
133 // Calculate a denominator used in the exponential part to convert from dB to linear scale:
134 // den = 20*constMaxGain (in Q8)
135 den = WEBRTC_SPL_MUL_16_U16(20, constMaxGain); // in Q8
136
137 for (i = 0; i < 32; i++)
138 {
139 // Calculate scaled input level (compressor):
140 // inLevel = fix((-constLog10_2*(compRatio-1)*(1-i)+fix(compRatio/2))/compRatio)
141 tmp16 = (int16_t)((kCompRatio - 1) * (i - 1)); // Q0
142 tmp32 = WEBRTC_SPL_MUL_16_U16(tmp16, kLog10_2) + 1; // Q14
143 inLevel = WebRtcSpl_DivW32W16(tmp32, kCompRatio); // Q14
144
145 // Calculate diffGain-inLevel, to map using the genFuncTable
146 inLevel = ((int32_t)diffGain << 14) - inLevel; // Q14
147
148 // Make calculations on abs(inLevel) and compensate for the sign afterwards.
149 absInLevel = (uint32_t)WEBRTC_SPL_ABS_W32(inLevel); // Q14
150
151 // LUT with interpolation
152 intPart = (uint16_t)(absInLevel >> 14);
153 fracPart = (uint16_t)(absInLevel & 0x00003FFF); // extract the fractional part
154 tmpU16 = kGenFuncTable[intPart + 1] - kGenFuncTable[intPart]; // Q8
155 tmpU32no1 = tmpU16 * fracPart; // Q22
156 tmpU32no1 += (uint32_t)kGenFuncTable[intPart] << 14; // Q22
157 logApprox = tmpU32no1 >> 8; // Q14
158 // Compensate for negative exponent using the relation:
159 // log2(1 + 2^-x) = log2(1 + 2^x) - x
160 if (inLevel < 0)
161 {
162 zeros = WebRtcSpl_NormU32(absInLevel);
163 zerosScale = 0;
164 if (zeros < 15)
165 {
166 // Not enough space for multiplication
167 tmpU32no2 = absInLevel >> (15 - zeros); // Q(zeros-1)
168 tmpU32no2 = WEBRTC_SPL_UMUL_32_16(tmpU32no2, kLogE_1); // Q(zeros+13)
169 if (zeros < 9)
170 {
171 zerosScale = 9 - zeros;
172 tmpU32no1 >>= zerosScale; // Q(zeros+13)
173 } else
174 {
175 tmpU32no2 >>= zeros - 9; // Q22
176 }
177 } else
178 {
179 tmpU32no2 = WEBRTC_SPL_UMUL_32_16(absInLevel, kLogE_1); // Q28
180 tmpU32no2 >>= 6; // Q22
181 }
182 logApprox = 0;
183 if (tmpU32no2 < tmpU32no1)
184 {
185 logApprox = (tmpU32no1 - tmpU32no2) >> (8 - zerosScale); //Q14
186 }
187 }
188 numFIX = (maxGain * constMaxGain) << 6; // Q14
189 numFIX -= (int32_t)logApprox * diffGain; // Q14
190
191 // Calculate ratio
192 // Shift |numFIX| as much as possible.
193 // Ensure we avoid wrap-around in |den| as well.
194 if (numFIX > (den >> 8)) // |den| is Q8.
195 {
196 zeros = WebRtcSpl_NormW32(numFIX);
197 } else
198 {
199 zeros = WebRtcSpl_NormW32(den) + 8;
200 }
201 numFIX <<= zeros; // Q(14+zeros)
202
203 // Shift den so we end up in Qy1
204 tmp32no1 = WEBRTC_SPL_SHIFT_W32(den, zeros - 8); // Q(zeros)
205 if (numFIX < 0)
206 {
207 numFIX -= tmp32no1 / 2;
208 } else
209 {
210 numFIX += tmp32no1 / 2;
211 }
212 y32 = numFIX / tmp32no1; // in Q14
213 if (limiterEnable && (i < limiterIdx))
214 {
215 tmp32 = WEBRTC_SPL_MUL_16_U16(i - 1, kLog10_2); // Q14
216 tmp32 -= limiterLvl << 14; // Q14
217 y32 = WebRtcSpl_DivW32W16(tmp32 + 10, 20);
218 }
219 if (y32 > 39000)
220 {
221 tmp32 = (y32 >> 1) * kLog10 + 4096; // in Q27
222 tmp32 >>= 13; // In Q14.
223 } else
224 {
225 tmp32 = y32 * kLog10 + 8192; // in Q28
226 tmp32 >>= 14; // In Q14.
227 }
228 tmp32 += 16 << 14; // in Q14 (Make sure final output is in Q16)
229
230 // Calculate power
231 if (tmp32 > 0)
232 {
233 intPart = (int16_t)(tmp32 >> 14);
234 fracPart = (uint16_t)(tmp32 & 0x00003FFF); // in Q14
235 if ((fracPart >> 13) != 0)
236 {
237 tmp16 = (2 << 14) - constLinApprox;
238 tmp32no2 = (1 << 14) - fracPart;
239 tmp32no2 *= tmp16;
240 tmp32no2 >>= 13;
241 tmp32no2 = (1 << 14) - tmp32no2;
242 } else
243 {
244 tmp16 = constLinApprox - (1 << 14);
245 tmp32no2 = (fracPart * tmp16) >> 13;
246 }
247 fracPart = (uint16_t)tmp32no2;
248 gainTable[i] =
249 (1 << intPart) + WEBRTC_SPL_SHIFT_W32(fracPart, intPart - 14);
250 } else
251 {
252 gainTable[i] = 0;
253 }
254 }
255
256 return 0;
257 }
258
WebRtcAgc_InitDigital(DigitalAgc * stt,int16_t agcMode)259 int32_t WebRtcAgc_InitDigital(DigitalAgc* stt, int16_t agcMode) {
260 if (agcMode == kAgcModeFixedDigital)
261 {
262 // start at minimum to find correct gain faster
263 stt->capacitorSlow = 0;
264 } else
265 {
266 // start out with 0 dB gain
267 stt->capacitorSlow = 134217728; // (int32_t)(0.125f * 32768.0f * 32768.0f);
268 }
269 stt->capacitorFast = 0;
270 stt->gain = 65536;
271 stt->gatePrevious = 0;
272 stt->agcMode = agcMode;
273 #ifdef WEBRTC_AGC_DEBUG_DUMP
274 stt->frameCounter = 0;
275 #endif
276
277 // initialize VADs
278 WebRtcAgc_InitVad(&stt->vadNearend);
279 WebRtcAgc_InitVad(&stt->vadFarend);
280
281 return 0;
282 }
283
WebRtcAgc_AddFarendToDigital(DigitalAgc * stt,const int16_t * in_far,size_t nrSamples)284 int32_t WebRtcAgc_AddFarendToDigital(DigitalAgc* stt,
285 const int16_t* in_far,
286 size_t nrSamples) {
287 assert(stt != NULL);
288 // VAD for far end
289 WebRtcAgc_ProcessVad(&stt->vadFarend, in_far, nrSamples);
290
291 return 0;
292 }
293
WebRtcAgc_ProcessDigital(DigitalAgc * stt,const int16_t * const * in_near,size_t num_bands,int16_t * const * out,uint32_t FS,int16_t lowlevelSignal)294 int32_t WebRtcAgc_ProcessDigital(DigitalAgc* stt,
295 const int16_t* const* in_near,
296 size_t num_bands,
297 int16_t* const* out,
298 uint32_t FS,
299 int16_t lowlevelSignal) {
300 // array for gains (one value per ms, incl start & end)
301 int32_t gains[11];
302
303 int32_t out_tmp, tmp32;
304 int32_t env[10];
305 int32_t max_nrg;
306 int32_t cur_level;
307 int32_t gain32, delta;
308 int16_t logratio;
309 int16_t lower_thr, upper_thr;
310 int16_t zeros = 0, zeros_fast, frac = 0;
311 int16_t decay;
312 int16_t gate, gain_adj;
313 int16_t k;
314 size_t n, i, L;
315 int16_t L2; // samples/subframe
316
317 // determine number of samples per ms
318 if (FS == 8000)
319 {
320 L = 8;
321 L2 = 3;
322 } else if (FS == 16000 || FS == 32000 || FS == 48000)
323 {
324 L = 16;
325 L2 = 4;
326 } else
327 {
328 return -1;
329 }
330
331 for (i = 0; i < num_bands; ++i)
332 {
333 if (in_near[i] != out[i])
334 {
335 // Only needed if they don't already point to the same place.
336 memcpy(out[i], in_near[i], 10 * L * sizeof(in_near[i][0]));
337 }
338 }
339 // VAD for near end
340 logratio = WebRtcAgc_ProcessVad(&stt->vadNearend, out[0], L * 10);
341
342 // Account for far end VAD
343 if (stt->vadFarend.counter > 10)
344 {
345 tmp32 = 3 * logratio;
346 logratio = (int16_t)((tmp32 - stt->vadFarend.logRatio) >> 2);
347 }
348
349 // Determine decay factor depending on VAD
350 // upper_thr = 1.0f;
351 // lower_thr = 0.25f;
352 upper_thr = 1024; // Q10
353 lower_thr = 0; // Q10
354 if (logratio > upper_thr)
355 {
356 // decay = -2^17 / DecayTime; -> -65
357 decay = -65;
358 } else if (logratio < lower_thr)
359 {
360 decay = 0;
361 } else
362 {
363 // decay = (int16_t)(((lower_thr - logratio)
364 // * (2^27/(DecayTime*(upper_thr-lower_thr)))) >> 10);
365 // SUBSTITUTED: 2^27/(DecayTime*(upper_thr-lower_thr)) -> 65
366 tmp32 = (lower_thr - logratio) * 65;
367 decay = (int16_t)(tmp32 >> 10);
368 }
369
370 // adjust decay factor for long silence (detected as low standard deviation)
371 // This is only done in the adaptive modes
372 if (stt->agcMode != kAgcModeFixedDigital)
373 {
374 if (stt->vadNearend.stdLongTerm < 4000)
375 {
376 decay = 0;
377 } else if (stt->vadNearend.stdLongTerm < 8096)
378 {
379 // decay = (int16_t)(((stt->vadNearend.stdLongTerm - 4000) * decay) >> 12);
380 tmp32 = (stt->vadNearend.stdLongTerm - 4000) * decay;
381 decay = (int16_t)(tmp32 >> 12);
382 }
383
384 if (lowlevelSignal != 0)
385 {
386 decay = 0;
387 }
388 }
389 #ifdef WEBRTC_AGC_DEBUG_DUMP
390 stt->frameCounter++;
391 fprintf(stt->logFile,
392 "%5.2f\t%d\t%d\t%d\t",
393 (float)(stt->frameCounter) / 100,
394 logratio,
395 decay,
396 stt->vadNearend.stdLongTerm);
397 #endif
398 // Find max amplitude per sub frame
399 // iterate over sub frames
400 for (k = 0; k < 10; k++)
401 {
402 // iterate over samples
403 max_nrg = 0;
404 for (n = 0; n < L; n++)
405 {
406 int32_t nrg = out[0][k * L + n] * out[0][k * L + n];
407 if (nrg > max_nrg)
408 {
409 max_nrg = nrg;
410 }
411 }
412 env[k] = max_nrg;
413 }
414
415 // Calculate gain per sub frame
416 gains[0] = stt->gain;
417 for (k = 0; k < 10; k++)
418 {
419 // Fast envelope follower
420 // decay time = -131000 / -1000 = 131 (ms)
421 stt->capacitorFast = AGC_SCALEDIFF32(-1000, stt->capacitorFast, stt->capacitorFast);
422 if (env[k] > stt->capacitorFast)
423 {
424 stt->capacitorFast = env[k];
425 }
426 // Slow envelope follower
427 if (env[k] > stt->capacitorSlow)
428 {
429 // increase capacitorSlow
430 stt->capacitorSlow
431 = AGC_SCALEDIFF32(500, (env[k] - stt->capacitorSlow), stt->capacitorSlow);
432 } else
433 {
434 // decrease capacitorSlow
435 stt->capacitorSlow
436 = AGC_SCALEDIFF32(decay, stt->capacitorSlow, stt->capacitorSlow);
437 }
438
439 // use maximum of both capacitors as current level
440 if (stt->capacitorFast > stt->capacitorSlow)
441 {
442 cur_level = stt->capacitorFast;
443 } else
444 {
445 cur_level = stt->capacitorSlow;
446 }
447 // Translate signal level into gain, using a piecewise linear approximation
448 // find number of leading zeros
449 zeros = WebRtcSpl_NormU32((uint32_t)cur_level);
450 if (cur_level == 0)
451 {
452 zeros = 31;
453 }
454 tmp32 = (cur_level << zeros) & 0x7FFFFFFF;
455 frac = (int16_t)(tmp32 >> 19); // Q12.
456 tmp32 = (stt->gainTable[zeros-1] - stt->gainTable[zeros]) * frac;
457 gains[k + 1] = stt->gainTable[zeros] + (tmp32 >> 12);
458 #ifdef WEBRTC_AGC_DEBUG_DUMP
459 if (k == 0) {
460 fprintf(stt->logFile,
461 "%d\t%d\t%d\t%d\t%d\n",
462 env[0],
463 cur_level,
464 stt->capacitorFast,
465 stt->capacitorSlow,
466 zeros);
467 }
468 #endif
469 }
470
471 // Gate processing (lower gain during absence of speech)
472 zeros = (zeros << 9) - (frac >> 3);
473 // find number of leading zeros
474 zeros_fast = WebRtcSpl_NormU32((uint32_t)stt->capacitorFast);
475 if (stt->capacitorFast == 0)
476 {
477 zeros_fast = 31;
478 }
479 tmp32 = (stt->capacitorFast << zeros_fast) & 0x7FFFFFFF;
480 zeros_fast <<= 9;
481 zeros_fast -= (int16_t)(tmp32 >> 22);
482
483 gate = 1000 + zeros_fast - zeros - stt->vadNearend.stdShortTerm;
484
485 if (gate < 0)
486 {
487 stt->gatePrevious = 0;
488 } else
489 {
490 tmp32 = stt->gatePrevious * 7;
491 gate = (int16_t)((gate + tmp32) >> 3);
492 stt->gatePrevious = gate;
493 }
494 // gate < 0 -> no gate
495 // gate > 2500 -> max gate
496 if (gate > 0)
497 {
498 if (gate < 2500)
499 {
500 gain_adj = (2500 - gate) >> 5;
501 } else
502 {
503 gain_adj = 0;
504 }
505 for (k = 0; k < 10; k++)
506 {
507 if ((gains[k + 1] - stt->gainTable[0]) > 8388608)
508 {
509 // To prevent wraparound
510 tmp32 = (gains[k + 1] - stt->gainTable[0]) >> 8;
511 tmp32 *= 178 + gain_adj;
512 } else
513 {
514 tmp32 = (gains[k+1] - stt->gainTable[0]) * (178 + gain_adj);
515 tmp32 >>= 8;
516 }
517 gains[k + 1] = stt->gainTable[0] + tmp32;
518 }
519 }
520
521 // Limit gain to avoid overload distortion
522 for (k = 0; k < 10; k++)
523 {
524 // To prevent wrap around
525 zeros = 10;
526 if (gains[k + 1] > 47453132)
527 {
528 zeros = 16 - WebRtcSpl_NormW32(gains[k + 1]);
529 }
530 gain32 = (gains[k + 1] >> zeros) + 1;
531 gain32 *= gain32;
532 // check for overflow
533 while (AGC_MUL32((env[k] >> 12) + 1, gain32)
534 > WEBRTC_SPL_SHIFT_W32((int32_t)32767, 2 * (1 - zeros + 10)))
535 {
536 // multiply by 253/256 ==> -0.1 dB
537 if (gains[k + 1] > 8388607)
538 {
539 // Prevent wrap around
540 gains[k + 1] = (gains[k+1] / 256) * 253;
541 } else
542 {
543 gains[k + 1] = (gains[k+1] * 253) / 256;
544 }
545 gain32 = (gains[k + 1] >> zeros) + 1;
546 gain32 *= gain32;
547 }
548 }
549 // gain reductions should be done 1 ms earlier than gain increases
550 for (k = 1; k < 10; k++)
551 {
552 if (gains[k] > gains[k + 1])
553 {
554 gains[k] = gains[k + 1];
555 }
556 }
557 // save start gain for next frame
558 stt->gain = gains[10];
559
560 // Apply gain
561 // handle first sub frame separately
562 delta = (gains[1] - gains[0]) << (4 - L2);
563 gain32 = gains[0] << 4;
564 // iterate over samples
565 for (n = 0; n < L; n++)
566 {
567 for (i = 0; i < num_bands; ++i)
568 {
569 tmp32 = out[i][n] * ((gain32 + 127) >> 7);
570 out_tmp = tmp32 >> 16;
571 if (out_tmp > 4095)
572 {
573 out[i][n] = (int16_t)32767;
574 } else if (out_tmp < -4096)
575 {
576 out[i][n] = (int16_t)-32768;
577 } else
578 {
579 tmp32 = out[i][n] * (gain32 >> 4);
580 out[i][n] = (int16_t)(tmp32 >> 16);
581 }
582 }
583 //
584
585 gain32 += delta;
586 }
587 // iterate over subframes
588 for (k = 1; k < 10; k++)
589 {
590 delta = (gains[k+1] - gains[k]) << (4 - L2);
591 gain32 = gains[k] << 4;
592 // iterate over samples
593 for (n = 0; n < L; n++)
594 {
595 for (i = 0; i < num_bands; ++i)
596 {
597 tmp32 = out[i][k * L + n] * (gain32 >> 4);
598 out[i][k * L + n] = (int16_t)(tmp32 >> 16);
599 }
600 gain32 += delta;
601 }
602 }
603
604 return 0;
605 }
606
WebRtcAgc_InitVad(AgcVad * state)607 void WebRtcAgc_InitVad(AgcVad* state) {
608 int16_t k;
609
610 state->HPstate = 0; // state of high pass filter
611 state->logRatio = 0; // log( P(active) / P(inactive) )
612 // average input level (Q10)
613 state->meanLongTerm = 15 << 10;
614
615 // variance of input level (Q8)
616 state->varianceLongTerm = 500 << 8;
617
618 state->stdLongTerm = 0; // standard deviation of input level in dB
619 // short-term average input level (Q10)
620 state->meanShortTerm = 15 << 10;
621
622 // short-term variance of input level (Q8)
623 state->varianceShortTerm = 500 << 8;
624
625 state->stdShortTerm = 0; // short-term standard deviation of input level in dB
626 state->counter = 3; // counts updates
627 for (k = 0; k < 8; k++)
628 {
629 // downsampling filter
630 state->downState[k] = 0;
631 }
632 }
633
WebRtcAgc_ProcessVad(AgcVad * state,const int16_t * in,size_t nrSamples)634 int16_t WebRtcAgc_ProcessVad(AgcVad* state, // (i) VAD state
635 const int16_t* in, // (i) Speech signal
636 size_t nrSamples) // (i) number of samples
637 {
638 int32_t out, nrg, tmp32, tmp32b;
639 uint16_t tmpU16;
640 int16_t k, subfr, tmp16;
641 int16_t buf1[8];
642 int16_t buf2[4];
643 int16_t HPstate;
644 int16_t zeros, dB;
645
646 // process in 10 sub frames of 1 ms (to save on memory)
647 nrg = 0;
648 HPstate = state->HPstate;
649 for (subfr = 0; subfr < 10; subfr++)
650 {
651 // downsample to 4 kHz
652 if (nrSamples == 160)
653 {
654 for (k = 0; k < 8; k++)
655 {
656 tmp32 = (int32_t)in[2 * k] + (int32_t)in[2 * k + 1];
657 tmp32 >>= 1;
658 buf1[k] = (int16_t)tmp32;
659 }
660 in += 16;
661
662 WebRtcSpl_DownsampleBy2(buf1, 8, buf2, state->downState);
663 } else
664 {
665 WebRtcSpl_DownsampleBy2(in, 8, buf2, state->downState);
666 in += 8;
667 }
668
669 // high pass filter and compute energy
670 for (k = 0; k < 4; k++)
671 {
672 out = buf2[k] + HPstate;
673 tmp32 = 600 * out;
674 HPstate = (int16_t)((tmp32 >> 10) - buf2[k]);
675 nrg += (out * out) >> 6;
676 }
677 }
678 state->HPstate = HPstate;
679
680 // find number of leading zeros
681 if (!(0xFFFF0000 & nrg))
682 {
683 zeros = 16;
684 } else
685 {
686 zeros = 0;
687 }
688 if (!(0xFF000000 & (nrg << zeros)))
689 {
690 zeros += 8;
691 }
692 if (!(0xF0000000 & (nrg << zeros)))
693 {
694 zeros += 4;
695 }
696 if (!(0xC0000000 & (nrg << zeros)))
697 {
698 zeros += 2;
699 }
700 if (!(0x80000000 & (nrg << zeros)))
701 {
702 zeros += 1;
703 }
704
705 // energy level (range {-32..30}) (Q10)
706 dB = (15 - zeros) << 11;
707
708 // Update statistics
709
710 if (state->counter < kAvgDecayTime)
711 {
712 // decay time = AvgDecTime * 10 ms
713 state->counter++;
714 }
715
716 // update short-term estimate of mean energy level (Q10)
717 tmp32 = state->meanShortTerm * 15 + dB;
718 state->meanShortTerm = (int16_t)(tmp32 >> 4);
719
720 // update short-term estimate of variance in energy level (Q8)
721 tmp32 = (dB * dB) >> 12;
722 tmp32 += state->varianceShortTerm * 15;
723 state->varianceShortTerm = tmp32 / 16;
724
725 // update short-term estimate of standard deviation in energy level (Q10)
726 tmp32 = state->meanShortTerm * state->meanShortTerm;
727 tmp32 = (state->varianceShortTerm << 12) - tmp32;
728 state->stdShortTerm = (int16_t)WebRtcSpl_Sqrt(tmp32);
729
730 // update long-term estimate of mean energy level (Q10)
731 tmp32 = state->meanLongTerm * state->counter + dB;
732 state->meanLongTerm = WebRtcSpl_DivW32W16ResW16(
733 tmp32, WebRtcSpl_AddSatW16(state->counter, 1));
734
735 // update long-term estimate of variance in energy level (Q8)
736 tmp32 = (dB * dB) >> 12;
737 tmp32 += state->varianceLongTerm * state->counter;
738 state->varianceLongTerm = WebRtcSpl_DivW32W16(
739 tmp32, WebRtcSpl_AddSatW16(state->counter, 1));
740
741 // update long-term estimate of standard deviation in energy level (Q10)
742 tmp32 = state->meanLongTerm * state->meanLongTerm;
743 tmp32 = (state->varianceLongTerm << 12) - tmp32;
744 state->stdLongTerm = (int16_t)WebRtcSpl_Sqrt(tmp32);
745
746 // update voice activity measure (Q10)
747 tmp16 = 3 << 12;
748 // TODO(bjornv): (dB - state->meanLongTerm) can overflow, e.g., in
749 // ApmTest.Process unit test. Previously the macro WEBRTC_SPL_MUL_16_16()
750 // was used, which did an intermediate cast to (int16_t), hence losing
751 // significant bits. This cause logRatio to max out positive, rather than
752 // negative. This is a bug, but has very little significance.
753 tmp32 = tmp16 * (int16_t)(dB - state->meanLongTerm);
754 tmp32 = WebRtcSpl_DivW32W16(tmp32, state->stdLongTerm);
755 tmpU16 = (13 << 12);
756 tmp32b = WEBRTC_SPL_MUL_16_U16(state->logRatio, tmpU16);
757 tmp32 += tmp32b >> 10;
758
759 state->logRatio = (int16_t)(tmp32 >> 6);
760
761 // limit
762 if (state->logRatio > 2048)
763 {
764 state->logRatio = 2048;
765 }
766 if (state->logRatio < -2048)
767 {
768 state->logRatio = -2048;
769 }
770
771 return state->logRatio; // Q10
772 }
773