1 /*
2 * Copyright (c) 2011 The WebRTC project authors. All Rights Reserved.
3 *
4 * Use of this source code is governed by a BSD-style license
5 * that can be found in the LICENSE file in the root of the source
6 * tree. An additional intellectual property rights grant can be found
7 * in the file PATENTS. All contributing project authors may
8 * be found in the AUTHORS file in the root of the source tree.
9 */
10
11 #include <math.h>
12 #include <string.h>
13
14 #include "modules/audio_coding/codecs/isac/main/source/lpc_analysis.h"
15 #include "modules/audio_coding/codecs/isac/main/source/settings.h"
16 #include "modules/audio_coding/codecs/isac/main/source/codec.h"
17 #include "modules/audio_coding/codecs/isac/main/source/entropy_coding.h"
18 #include "modules/audio_coding/codecs/isac/main/source/filter_functions.h"
19 #include "modules/audio_coding/codecs/isac/main/source/isac_vad.h"
20
21 /* window */
22 /* Matlab generation code:
23 * t = (1:256)/257; r = 1-(1-t).^.45; w = sin(r*pi).^3; w = w/sum(w); plot((1:256)/8, w); grid;
24 * for k=1:16, fprintf(1, '%.8f, ', w(k*16 + (-15:0))); fprintf(1, '\n'); end
25 */
26 static const double kLpcCorrWindow[WINLEN] = {
27 0.00000000, 0.00000001, 0.00000004, 0.00000010, 0.00000020,
28 0.00000035, 0.00000055, 0.00000083, 0.00000118, 0.00000163,
29 0.00000218, 0.00000283, 0.00000361, 0.00000453, 0.00000558, 0.00000679,
30 0.00000817, 0.00000973, 0.00001147, 0.00001342, 0.00001558,
31 0.00001796, 0.00002058, 0.00002344, 0.00002657, 0.00002997,
32 0.00003365, 0.00003762, 0.00004190, 0.00004651, 0.00005144, 0.00005673,
33 0.00006236, 0.00006837, 0.00007476, 0.00008155, 0.00008875,
34 0.00009636, 0.00010441, 0.00011290, 0.00012186, 0.00013128,
35 0.00014119, 0.00015160, 0.00016252, 0.00017396, 0.00018594, 0.00019846,
36 0.00021155, 0.00022521, 0.00023946, 0.00025432, 0.00026978,
37 0.00028587, 0.00030260, 0.00031998, 0.00033802, 0.00035674,
38 0.00037615, 0.00039626, 0.00041708, 0.00043863, 0.00046092, 0.00048396,
39 0.00050775, 0.00053233, 0.00055768, 0.00058384, 0.00061080,
40 0.00063858, 0.00066720, 0.00069665, 0.00072696, 0.00075813,
41 0.00079017, 0.00082310, 0.00085692, 0.00089164, 0.00092728, 0.00096384,
42 0.00100133, 0.00103976, 0.00107914, 0.00111947, 0.00116077,
43 0.00120304, 0.00124630, 0.00129053, 0.00133577, 0.00138200,
44 0.00142924, 0.00147749, 0.00152676, 0.00157705, 0.00162836, 0.00168070,
45 0.00173408, 0.00178850, 0.00184395, 0.00190045, 0.00195799,
46 0.00201658, 0.00207621, 0.00213688, 0.00219860, 0.00226137,
47 0.00232518, 0.00239003, 0.00245591, 0.00252284, 0.00259079, 0.00265977,
48 0.00272977, 0.00280078, 0.00287280, 0.00294582, 0.00301984,
49 0.00309484, 0.00317081, 0.00324774, 0.00332563, 0.00340446,
50 0.00348421, 0.00356488, 0.00364644, 0.00372889, 0.00381220, 0.00389636,
51 0.00398135, 0.00406715, 0.00415374, 0.00424109, 0.00432920,
52 0.00441802, 0.00450754, 0.00459773, 0.00468857, 0.00478001,
53 0.00487205, 0.00496464, 0.00505775, 0.00515136, 0.00524542, 0.00533990,
54 0.00543476, 0.00552997, 0.00562548, 0.00572125, 0.00581725,
55 0.00591342, 0.00600973, 0.00610612, 0.00620254, 0.00629895,
56 0.00639530, 0.00649153, 0.00658758, 0.00668341, 0.00677894, 0.00687413,
57 0.00696891, 0.00706322, 0.00715699, 0.00725016, 0.00734266,
58 0.00743441, 0.00752535, 0.00761540, 0.00770449, 0.00779254,
59 0.00787947, 0.00796519, 0.00804963, 0.00813270, 0.00821431, 0.00829437,
60 0.00837280, 0.00844949, 0.00852436, 0.00859730, 0.00866822,
61 0.00873701, 0.00880358, 0.00886781, 0.00892960, 0.00898884,
62 0.00904542, 0.00909923, 0.00915014, 0.00919805, 0.00924283, 0.00928436,
63 0.00932252, 0.00935718, 0.00938821, 0.00941550, 0.00943890,
64 0.00945828, 0.00947351, 0.00948446, 0.00949098, 0.00949294,
65 0.00949020, 0.00948262, 0.00947005, 0.00945235, 0.00942938, 0.00940099,
66 0.00936704, 0.00932738, 0.00928186, 0.00923034, 0.00917268,
67 0.00910872, 0.00903832, 0.00896134, 0.00887763, 0.00878706,
68 0.00868949, 0.00858478, 0.00847280, 0.00835343, 0.00822653, 0.00809199,
69 0.00794970, 0.00779956, 0.00764145, 0.00747530, 0.00730103,
70 0.00711857, 0.00692787, 0.00672888, 0.00652158, 0.00630597,
71 0.00608208, 0.00584994, 0.00560962, 0.00536124, 0.00510493, 0.00484089,
72 0.00456935, 0.00429062, 0.00400505, 0.00371310, 0.00341532,
73 0.00311238, 0.00280511, 0.00249452, 0.00218184, 0.00186864,
74 0.00155690, 0.00124918, 0.00094895, 0.00066112, 0.00039320, 0.00015881
75 };
76
WebRtcIsac_GetVars(const double * input,const int16_t * pitchGains_Q12,double * oldEnergy,double * varscale)77 static void WebRtcIsac_GetVars(const double* input,
78 const int16_t* pitchGains_Q12,
79 double* oldEnergy,
80 double* varscale) {
81 double nrg[4], chng, pg;
82 int k;
83
84 double pitchGains[4]={0,0,0,0};;
85
86 /* Calculate energies of first and second frame halfs */
87 nrg[0] = 0.0001;
88 for (k = QLOOKAHEAD/2; k < (FRAMESAMPLES_QUARTER + QLOOKAHEAD) / 2; k++) {
89 nrg[0] += input[k]*input[k];
90 }
91 nrg[1] = 0.0001;
92 for ( ; k < (FRAMESAMPLES_HALF + QLOOKAHEAD) / 2; k++) {
93 nrg[1] += input[k]*input[k];
94 }
95 nrg[2] = 0.0001;
96 for ( ; k < (FRAMESAMPLES*3/4 + QLOOKAHEAD) / 2; k++) {
97 nrg[2] += input[k]*input[k];
98 }
99 nrg[3] = 0.0001;
100 for ( ; k < (FRAMESAMPLES + QLOOKAHEAD) / 2; k++) {
101 nrg[3] += input[k]*input[k];
102 }
103
104 /* Calculate average level change */
105 chng = 0.25 * (fabs(10.0 * log10(nrg[3] / nrg[2])) +
106 fabs(10.0 * log10(nrg[2] / nrg[1])) +
107 fabs(10.0 * log10(nrg[1] / nrg[0])) +
108 fabs(10.0 * log10(nrg[0] / *oldEnergy)));
109
110
111 /* Find average pitch gain */
112 pg = 0.0;
113 for (k=0; k<4; k++)
114 {
115 pitchGains[k] = ((float)pitchGains_Q12[k])/4096;
116 pg += pitchGains[k];
117 }
118 pg *= 0.25;
119
120 /* If pitch gain is low and energy constant - increase noise level*/
121 /* Matlab code:
122 pg = 0:.01:.45; plot(pg, 0.0 + 1.0 * exp( -1.0 * exp(-200.0 * pg.*pg.*pg) / (1.0 + 0.4 * 0) ))
123 */
124 *varscale = 0.0 + 1.0 * exp( -1.4 * exp(-200.0 * pg*pg*pg) / (1.0 + 0.4 * chng) );
125
126 *oldEnergy = nrg[3];
127 }
128
WebRtcIsac_GetVarsUB(const double * input,double * oldEnergy,double * varscale)129 static void WebRtcIsac_GetVarsUB(const double* input,
130 double* oldEnergy,
131 double* varscale) {
132 double nrg[4], chng;
133 int k;
134
135 /* Calculate energies of first and second frame halfs */
136 nrg[0] = 0.0001;
137 for (k = 0; k < (FRAMESAMPLES_QUARTER) / 2; k++) {
138 nrg[0] += input[k]*input[k];
139 }
140 nrg[1] = 0.0001;
141 for ( ; k < (FRAMESAMPLES_HALF) / 2; k++) {
142 nrg[1] += input[k]*input[k];
143 }
144 nrg[2] = 0.0001;
145 for ( ; k < (FRAMESAMPLES*3/4) / 2; k++) {
146 nrg[2] += input[k]*input[k];
147 }
148 nrg[3] = 0.0001;
149 for ( ; k < (FRAMESAMPLES) / 2; k++) {
150 nrg[3] += input[k]*input[k];
151 }
152
153 /* Calculate average level change */
154 chng = 0.25 * (fabs(10.0 * log10(nrg[3] / nrg[2])) +
155 fabs(10.0 * log10(nrg[2] / nrg[1])) +
156 fabs(10.0 * log10(nrg[1] / nrg[0])) +
157 fabs(10.0 * log10(nrg[0] / *oldEnergy)));
158
159
160 /* If pitch gain is low and energy constant - increase noise level*/
161 /* Matlab code:
162 pg = 0:.01:.45; plot(pg, 0.0 + 1.0 * exp( -1.0 * exp(-200.0 * pg.*pg.*pg) / (1.0 + 0.4 * 0) ))
163 */
164 *varscale = exp( -1.4 / (1.0 + 0.4 * chng) );
165
166 *oldEnergy = nrg[3];
167 }
168
WebRtcIsac_GetLpcCoefLb(double * inLo,double * inHi,MaskFiltstr * maskdata,double signal_noise_ratio,const int16_t * pitchGains_Q12,double * lo_coeff,double * hi_coeff)169 void WebRtcIsac_GetLpcCoefLb(double *inLo, double *inHi, MaskFiltstr *maskdata,
170 double signal_noise_ratio, const int16_t *pitchGains_Q12,
171 double *lo_coeff, double *hi_coeff)
172 {
173 int k, n, j, pos1, pos2;
174 double varscale;
175
176 double DataLo[WINLEN], DataHi[WINLEN];
177 double corrlo[ORDERLO+2], corrlo2[ORDERLO+1];
178 double corrhi[ORDERHI+1];
179 double k_veclo[ORDERLO], k_vechi[ORDERHI];
180
181 double a_LO[ORDERLO+1], a_HI[ORDERHI+1];
182 double tmp, res_nrg;
183
184 double FwdA, FwdB;
185
186 /* hearing threshold level in dB; higher value gives more noise */
187 const double HearThresOffset = -28.0;
188
189 /* bandwdith expansion factors for low- and high band */
190 const double gammaLo = 0.9;
191 const double gammaHi = 0.8;
192
193 /* less-noise-at-low-frequencies factor */
194 double aa;
195
196
197 /* convert from dB to signal level */
198 const double H_T_H = pow(10.0, 0.05 * HearThresOffset);
199 double S_N_R = pow(10.0, 0.05 * signal_noise_ratio) / 3.46; /* divide by sqrt(12) */
200
201 /* change quallevel depending on pitch gains and level fluctuations */
202 WebRtcIsac_GetVars(inLo, pitchGains_Q12, &(maskdata->OldEnergy), &varscale);
203
204 /* less-noise-at-low-frequencies factor */
205 aa = 0.35 * (0.5 + 0.5 * varscale);
206
207 /* replace data in buffer by new look-ahead data */
208 for (pos1 = 0; pos1 < QLOOKAHEAD; pos1++)
209 maskdata->DataBufferLo[pos1 + WINLEN - QLOOKAHEAD] = inLo[pos1];
210
211 for (k = 0; k < SUBFRAMES; k++) {
212
213 /* Update input buffer and multiply signal with window */
214 for (pos1 = 0; pos1 < WINLEN - UPDATE/2; pos1++) {
215 maskdata->DataBufferLo[pos1] = maskdata->DataBufferLo[pos1 + UPDATE/2];
216 maskdata->DataBufferHi[pos1] = maskdata->DataBufferHi[pos1 + UPDATE/2];
217 DataLo[pos1] = maskdata->DataBufferLo[pos1] * kLpcCorrWindow[pos1];
218 DataHi[pos1] = maskdata->DataBufferHi[pos1] * kLpcCorrWindow[pos1];
219 }
220 pos2 = k * UPDATE/2;
221 for (n = 0; n < UPDATE/2; n++, pos1++) {
222 maskdata->DataBufferLo[pos1] = inLo[QLOOKAHEAD + pos2];
223 maskdata->DataBufferHi[pos1] = inHi[pos2++];
224 DataLo[pos1] = maskdata->DataBufferLo[pos1] * kLpcCorrWindow[pos1];
225 DataHi[pos1] = maskdata->DataBufferHi[pos1] * kLpcCorrWindow[pos1];
226 }
227
228 /* Get correlation coefficients */
229 WebRtcIsac_AutoCorr(corrlo, DataLo, WINLEN, ORDERLO+1); /* computing autocorrelation */
230 WebRtcIsac_AutoCorr(corrhi, DataHi, WINLEN, ORDERHI);
231
232
233 /* less noise for lower frequencies, by filtering/scaling autocorrelation sequences */
234 corrlo2[0] = (1.0+aa*aa) * corrlo[0] - 2.0*aa * corrlo[1];
235 tmp = (1.0 + aa*aa);
236 for (n = 1; n <= ORDERLO; n++) {
237 corrlo2[n] = tmp * corrlo[n] - aa * (corrlo[n-1] + corrlo[n+1]);
238 }
239 tmp = (1.0+aa) * (1.0+aa);
240 for (n = 0; n <= ORDERHI; n++) {
241 corrhi[n] = tmp * corrhi[n];
242 }
243
244 /* add white noise floor */
245 corrlo2[0] += 1e-6;
246 corrhi[0] += 1e-6;
247
248
249 FwdA = 0.01;
250 FwdB = 0.01;
251
252 /* recursive filtering of correlation over subframes */
253 for (n = 0; n <= ORDERLO; n++) {
254 maskdata->CorrBufLo[n] = FwdA * maskdata->CorrBufLo[n] + corrlo2[n];
255 corrlo2[n] = ((1.0-FwdA)*FwdB) * maskdata->CorrBufLo[n] + (1.0-FwdB) * corrlo2[n];
256 }
257 for (n = 0; n <= ORDERHI; n++) {
258 maskdata->CorrBufHi[n] = FwdA * maskdata->CorrBufHi[n] + corrhi[n];
259 corrhi[n] = ((1.0-FwdA)*FwdB) * maskdata->CorrBufHi[n] + (1.0-FwdB) * corrhi[n];
260 }
261
262 /* compute prediction coefficients */
263 WebRtcIsac_LevDurb(a_LO, k_veclo, corrlo2, ORDERLO);
264 WebRtcIsac_LevDurb(a_HI, k_vechi, corrhi, ORDERHI);
265
266 /* bandwidth expansion */
267 tmp = gammaLo;
268 for (n = 1; n <= ORDERLO; n++) {
269 a_LO[n] *= tmp;
270 tmp *= gammaLo;
271 }
272
273 /* residual energy */
274 res_nrg = 0.0;
275 for (j = 0; j <= ORDERLO; j++) {
276 for (n = 0; n <= j; n++) {
277 res_nrg += a_LO[j] * corrlo2[j-n] * a_LO[n];
278 }
279 for (n = j+1; n <= ORDERLO; n++) {
280 res_nrg += a_LO[j] * corrlo2[n-j] * a_LO[n];
281 }
282 }
283
284 /* add hearing threshold and compute the gain */
285 *lo_coeff++ = S_N_R / (sqrt(res_nrg) / varscale + H_T_H);
286
287 /* copy coefficients to output array */
288 for (n = 1; n <= ORDERLO; n++) {
289 *lo_coeff++ = a_LO[n];
290 }
291
292
293 /* bandwidth expansion */
294 tmp = gammaHi;
295 for (n = 1; n <= ORDERHI; n++) {
296 a_HI[n] *= tmp;
297 tmp *= gammaHi;
298 }
299
300 /* residual energy */
301 res_nrg = 0.0;
302 for (j = 0; j <= ORDERHI; j++) {
303 for (n = 0; n <= j; n++) {
304 res_nrg += a_HI[j] * corrhi[j-n] * a_HI[n];
305 }
306 for (n = j+1; n <= ORDERHI; n++) {
307 res_nrg += a_HI[j] * corrhi[n-j] * a_HI[n];
308 }
309 }
310
311 /* add hearing threshold and compute of the gain */
312 *hi_coeff++ = S_N_R / (sqrt(res_nrg) / varscale + H_T_H);
313
314 /* copy coefficients to output array */
315 for (n = 1; n <= ORDERHI; n++) {
316 *hi_coeff++ = a_HI[n];
317 }
318 }
319 }
320
321
322
323 /******************************************************************************
324 * WebRtcIsac_GetLpcCoefUb()
325 *
326 * Compute LP coefficients and correlation coefficients. At 12 kHz LP
327 * coefficients of the first and the last sub-frame is computed. At 16 kHz
328 * LP coefficients of 4th, 8th and 12th sub-frames are computed. We always
329 * compute correlation coefficients of all sub-frames.
330 *
331 * Inputs:
332 * -inSignal : Input signal
333 * -maskdata : a structure keeping signal from previous frame.
334 * -bandwidth : specifies if the codec is in 0-16 kHz mode or
335 * 0-12 kHz mode.
336 *
337 * Outputs:
338 * -lpCoeff : pointer to a buffer where A-polynomials are
339 * written to (first coeff is 1 and it is not
340 * written)
341 * -corrMat : a matrix where correlation coefficients of each
342 * sub-frame are written to one row.
343 * -varscale : a scale used to compute LPC gains.
344 */
345 void
WebRtcIsac_GetLpcCoefUb(double * inSignal,MaskFiltstr * maskdata,double * lpCoeff,double corrMat[][UB_LPC_ORDER+1],double * varscale,int16_t bandwidth)346 WebRtcIsac_GetLpcCoefUb(
347 double* inSignal,
348 MaskFiltstr* maskdata,
349 double* lpCoeff,
350 double corrMat[][UB_LPC_ORDER + 1],
351 double* varscale,
352 int16_t bandwidth)
353 {
354 int frameCntr, activeFrameCntr, n, pos1, pos2;
355 int16_t criterion1;
356 int16_t criterion2;
357 int16_t numSubFrames = SUBFRAMES * (1 + (bandwidth == isac16kHz));
358 double data[WINLEN];
359 double corrSubFrame[UB_LPC_ORDER+2];
360 double reflecCoeff[UB_LPC_ORDER];
361
362 double aPolynom[UB_LPC_ORDER+1];
363 double tmp;
364
365 /* bandwdith expansion factors */
366 const double gamma = 0.9;
367
368 /* change quallevel depending on pitch gains and level fluctuations */
369 WebRtcIsac_GetVarsUB(inSignal, &(maskdata->OldEnergy), varscale);
370
371 /* replace data in buffer by new look-ahead data */
372 for(frameCntr = 0, activeFrameCntr = 0; frameCntr < numSubFrames;
373 frameCntr++)
374 {
375 if(frameCntr == SUBFRAMES)
376 {
377 // we are in 16 kHz
378 varscale++;
379 WebRtcIsac_GetVarsUB(&inSignal[FRAMESAMPLES_HALF],
380 &(maskdata->OldEnergy), varscale);
381 }
382 /* Update input buffer and multiply signal with window */
383 for(pos1 = 0; pos1 < WINLEN - UPDATE/2; pos1++)
384 {
385 maskdata->DataBufferLo[pos1] = maskdata->DataBufferLo[pos1 +
386 UPDATE/2];
387 data[pos1] = maskdata->DataBufferLo[pos1] * kLpcCorrWindow[pos1];
388 }
389 pos2 = frameCntr * UPDATE/2;
390 for(n = 0; n < UPDATE/2; n++, pos1++, pos2++)
391 {
392 maskdata->DataBufferLo[pos1] = inSignal[pos2];
393 data[pos1] = maskdata->DataBufferLo[pos1] * kLpcCorrWindow[pos1];
394 }
395
396 /* Get correlation coefficients */
397 /* computing autocorrelation */
398 WebRtcIsac_AutoCorr(corrSubFrame, data, WINLEN, UB_LPC_ORDER+1);
399 memcpy(corrMat[frameCntr], corrSubFrame,
400 (UB_LPC_ORDER+1)*sizeof(double));
401
402 criterion1 = ((frameCntr == 0) || (frameCntr == (SUBFRAMES - 1))) &&
403 (bandwidth == isac12kHz);
404 criterion2 = (((frameCntr+1) % 4) == 0) &&
405 (bandwidth == isac16kHz);
406 if(criterion1 || criterion2)
407 {
408 /* add noise */
409 corrSubFrame[0] += 1e-6;
410 /* compute prediction coefficients */
411 WebRtcIsac_LevDurb(aPolynom, reflecCoeff, corrSubFrame,
412 UB_LPC_ORDER);
413
414 /* bandwidth expansion */
415 tmp = gamma;
416 for (n = 1; n <= UB_LPC_ORDER; n++)
417 {
418 *lpCoeff++ = aPolynom[n] * tmp;
419 tmp *= gamma;
420 }
421 activeFrameCntr++;
422 }
423 }
424 }
425
426
427
428 /******************************************************************************
429 * WebRtcIsac_GetLpcGain()
430 *
431 * Compute the LPC gains for each sub-frame, given the LPC of each sub-frame
432 * and the corresponding correlation coefficients.
433 *
434 * Inputs:
435 * -signal_noise_ratio : the desired SNR in dB.
436 * -numVecs : number of sub-frames
437 * -corrMat : a matrix of correlation coefficients where
438 * each row is a set of correlation coefficients of
439 * one sub-frame.
440 * -varscale : a scale computed when WebRtcIsac_GetLpcCoefUb()
441 * is called.
442 *
443 * Outputs:
444 * -gain : pointer to a buffer where LP gains are written.
445 *
446 */
447 void
WebRtcIsac_GetLpcGain(double signal_noise_ratio,const double * filtCoeffVecs,int numVecs,double * gain,double corrMat[][UB_LPC_ORDER+1],const double * varscale)448 WebRtcIsac_GetLpcGain(
449 double signal_noise_ratio,
450 const double* filtCoeffVecs,
451 int numVecs,
452 double* gain,
453 double corrMat[][UB_LPC_ORDER + 1],
454 const double* varscale)
455 {
456 int16_t j, n;
457 int16_t subFrameCntr;
458 double aPolynom[ORDERLO + 1];
459 double res_nrg;
460
461 const double HearThresOffset = -28.0;
462 const double H_T_H = pow(10.0, 0.05 * HearThresOffset);
463 /* divide by sqrt(12) = 3.46 */
464 const double S_N_R = pow(10.0, 0.05 * signal_noise_ratio) / 3.46;
465
466 aPolynom[0] = 1;
467 for(subFrameCntr = 0; subFrameCntr < numVecs; subFrameCntr++)
468 {
469 if(subFrameCntr == SUBFRAMES)
470 {
471 // we are in second half of a SWB frame. use new varscale
472 varscale++;
473 }
474 memcpy(&aPolynom[1], &filtCoeffVecs[(subFrameCntr * (UB_LPC_ORDER + 1)) +
475 1], sizeof(double) * UB_LPC_ORDER);
476
477 /* residual energy */
478 res_nrg = 0.0;
479 for(j = 0; j <= UB_LPC_ORDER; j++)
480 {
481 for(n = 0; n <= j; n++)
482 {
483 res_nrg += aPolynom[j] * corrMat[subFrameCntr][j-n] *
484 aPolynom[n];
485 }
486 for(n = j+1; n <= UB_LPC_ORDER; n++)
487 {
488 res_nrg += aPolynom[j] * corrMat[subFrameCntr][n-j] *
489 aPolynom[n];
490 }
491 }
492
493 /* add hearing threshold and compute the gain */
494 gain[subFrameCntr] = S_N_R / (sqrt(res_nrg) / *varscale + H_T_H);
495 }
496 }
497