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