Home | History | Annotate | Download | only in source
      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 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   alpha = 0; //warning -DH
     84   a[0] = 1.0;
     85   if (r[0] < LEVINSON_EPS) { /* if r[0] <= 0, set LPC coeff. to zero */
     86     for (i = 0; i < order; i++) {
     87       k[i] = 0;
     88       a[i+1] = 0;
     89     }
     90   } else {
     91     a[1] = k[0] = -r[1]/r[0];
     92     alpha = r[0] + r[1] * k[0];
     93     for (m = 1; m < order; m++){
     94       sum = r[m + 1];
     95       for (i = 0; i < m; i++){
     96         sum += a[i+1] * r[m - i];
     97       }
     98       k[m] = -sum / alpha;
     99       alpha += k[m] * sum;
    100       m_h = (m + 1) >> 1;
    101       for (i = 0; i < m_h; i++){
    102         sum = a[i+1] + k[m] * a[m - i];
    103         a[m - i] += k[m] * a[i+1];
    104         a[i+1] = sum;
    105       }
    106       a[m+1] = k[m];
    107     }
    108   }
    109   return alpha;
    110 }
    111 
    112 
    113 //was static before, but didn't work with MEX file
    114 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   double pitchGains[4]={0,0,0,0};;
    121 
    122   /* Calculate energies of first and second frame halfs */
    123   nrg[0] = 0.0001;
    124   for (k = QLOOKAHEAD/2; k < (FRAMESAMPLES_QUARTER + QLOOKAHEAD) / 2; k++) {
    125     nrg[0] += input[k]*input[k];
    126   }
    127   nrg[1] = 0.0001;
    128   for ( ; k < (FRAMESAMPLES_HALF + QLOOKAHEAD) / 2; k++) {
    129     nrg[1] += input[k]*input[k];
    130   }
    131   nrg[2] = 0.0001;
    132   for ( ; k < (FRAMESAMPLES*3/4 + QLOOKAHEAD) / 2; k++) {
    133     nrg[2] += input[k]*input[k];
    134   }
    135   nrg[3] = 0.0001;
    136   for ( ; k < (FRAMESAMPLES + QLOOKAHEAD) / 2; k++) {
    137     nrg[3] += input[k]*input[k];
    138   }
    139 
    140   /* Calculate average level change */
    141   chng = 0.25 * (fabs(10.0 * log10(nrg[3] / nrg[2])) +
    142                  fabs(10.0 * log10(nrg[2] / nrg[1])) +
    143                  fabs(10.0 * log10(nrg[1] / nrg[0])) +
    144                  fabs(10.0 * log10(nrg[0] / *oldEnergy)));
    145 
    146 
    147   /* Find average pitch gain */
    148   pg = 0.0;
    149   for (k=0; k<4; k++)
    150   {
    151     pitchGains[k] = ((float)pitchGains_Q12[k])/4096;
    152     pg += pitchGains[k];
    153   }
    154   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   *varscale = 0.0 + 1.0 * exp( -1.4 * exp(-200.0 * pg*pg*pg) / (1.0 + 0.4 * chng) );
    161 
    162   *oldEnergy = nrg[3];
    163 }
    164 
    165 void
    166 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   nrg[0] = 0.0001;
    176   for (k = 0; k < (FRAMESAMPLES_QUARTER) / 2; k++) {
    177     nrg[0] += input[k]*input[k];
    178   }
    179   nrg[1] = 0.0001;
    180   for ( ; k < (FRAMESAMPLES_HALF) / 2; k++) {
    181     nrg[1] += input[k]*input[k];
    182   }
    183   nrg[2] = 0.0001;
    184   for ( ; k < (FRAMESAMPLES*3/4) / 2; k++) {
    185     nrg[2] += input[k]*input[k];
    186   }
    187   nrg[3] = 0.0001;
    188   for ( ; k < (FRAMESAMPLES) / 2; k++) {
    189     nrg[3] += input[k]*input[k];
    190   }
    191 
    192   /* Calculate average level change */
    193   chng = 0.25 * (fabs(10.0 * log10(nrg[3] / nrg[2])) +
    194                  fabs(10.0 * log10(nrg[2] / nrg[1])) +
    195                  fabs(10.0 * log10(nrg[1] / nrg[0])) +
    196                  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   *varscale = exp( -1.4 / (1.0 + 0.4 * chng) );
    204 
    205   *oldEnergy = nrg[3];
    206 }
    207 
    208 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   const double HearThresOffset = -28.0;
    227 
    228   /* bandwdith expansion factors for low- and high band */
    229   const double gammaLo = 0.9;
    230   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   const double H_T_H = pow(10.0, 0.05 * HearThresOffset);
    238   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   WebRtcIsac_GetVars(inLo, pitchGains_Q12, &(maskdata->OldEnergy), &varscale);
    242 
    243   /* less-noise-at-low-frequencies factor */
    244   aa = 0.35 * (0.5 + 0.5 * varscale);
    245 
    246   /* replace data in buffer by new look-ahead data */
    247   for (pos1 = 0; pos1 < QLOOKAHEAD; pos1++)
    248     maskdata->DataBufferLo[pos1 + WINLEN - QLOOKAHEAD] = inLo[pos1];
    249 
    250   for (k = 0; k < SUBFRAMES; k++) {
    251 
    252     /* Update input buffer and multiply signal with window */
    253     for (pos1 = 0; pos1 < WINLEN - UPDATE/2; pos1++) {
    254       maskdata->DataBufferLo[pos1] = maskdata->DataBufferLo[pos1 + UPDATE/2];
    255       maskdata->DataBufferHi[pos1] = maskdata->DataBufferHi[pos1 + UPDATE/2];
    256       DataLo[pos1] = maskdata->DataBufferLo[pos1] * kLpcCorrWindow[pos1];
    257       DataHi[pos1] = maskdata->DataBufferHi[pos1] * kLpcCorrWindow[pos1];
    258     }
    259     pos2 = k * UPDATE/2;
    260     for (n = 0; n < UPDATE/2; n++, pos1++) {
    261       maskdata->DataBufferLo[pos1] = inLo[QLOOKAHEAD + pos2];
    262       maskdata->DataBufferHi[pos1] = inHi[pos2++];
    263       DataLo[pos1] = maskdata->DataBufferLo[pos1] * kLpcCorrWindow[pos1];
    264       DataHi[pos1] = maskdata->DataBufferHi[pos1] * kLpcCorrWindow[pos1];
    265     }
    266 
    267     /* Get correlation coefficients */
    268     WebRtcIsac_AutoCorr(corrlo, DataLo, WINLEN, ORDERLO+1); /* computing autocorrelation */
    269     WebRtcIsac_AutoCorr(corrhi, DataHi, WINLEN, ORDERHI);
    270 
    271 
    272     /* less noise for lower frequencies, by filtering/scaling autocorrelation sequences */
    273     corrlo2[0] = (1.0+aa*aa) * corrlo[0] - 2.0*aa * corrlo[1];
    274     tmp = (1.0 + aa*aa);
    275     for (n = 1; n <= ORDERLO; n++) {
    276       corrlo2[n] = tmp * corrlo[n] - aa * (corrlo[n-1] + corrlo[n+1]);
    277     }
    278     tmp = (1.0+aa) * (1.0+aa);
    279     for (n = 0; n <= ORDERHI; n++) {
    280       corrhi[n] = tmp * corrhi[n];
    281     }
    282 
    283     /* add white noise floor */
    284     corrlo2[0] += 1e-6;
    285     corrhi[0] += 1e-6;
    286 
    287 
    288     FwdA = 0.01;
    289     FwdB = 0.01;
    290 
    291     /* recursive filtering of correlation over subframes */
    292     for (n = 0; n <= ORDERLO; n++) {
    293       maskdata->CorrBufLo[n] = FwdA * maskdata->CorrBufLo[n] + corrlo2[n];
    294       corrlo2[n] = ((1.0-FwdA)*FwdB) * maskdata->CorrBufLo[n] + (1.0-FwdB) * corrlo2[n];
    295     }
    296     for (n = 0; n <= ORDERHI; n++) {
    297       maskdata->CorrBufHi[n] = FwdA * maskdata->CorrBufHi[n] + corrhi[n];
    298       corrhi[n] = ((1.0-FwdA)*FwdB) * maskdata->CorrBufHi[n] + (1.0-FwdB) * corrhi[n];
    299     }
    300 
    301     /* compute prediction coefficients */
    302     WebRtcIsac_LevDurb(a_LO, k_veclo, corrlo2, ORDERLO);
    303     WebRtcIsac_LevDurb(a_HI, k_vechi, corrhi, ORDERHI);
    304 
    305     /* bandwidth expansion */
    306     tmp = gammaLo;
    307     for (n = 1; n <= ORDERLO; n++) {
    308       a_LO[n] *= tmp;
    309       tmp *= gammaLo;
    310     }
    311 
    312     /* residual energy */
    313     res_nrg = 0.0;
    314     for (j = 0; j <= ORDERLO; j++) {
    315       for (n = 0; n <= j; n++) {
    316         res_nrg += a_LO[j] * corrlo2[j-n] * a_LO[n];
    317       }
    318       for (n = j+1; n <= ORDERLO; n++) {
    319         res_nrg += a_LO[j] * corrlo2[n-j] * a_LO[n];
    320       }
    321     }
    322 
    323     /* add hearing threshold and compute the gain */
    324     *lo_coeff++ = S_N_R / (sqrt(res_nrg) / varscale + H_T_H);
    325 
    326     /* copy coefficients to output array */
    327     for (n = 1; n <= ORDERLO; n++) {
    328       *lo_coeff++ = a_LO[n];
    329     }
    330 
    331 
    332     /* bandwidth expansion */
    333     tmp = gammaHi;
    334     for (n = 1; n <= ORDERHI; n++) {
    335       a_HI[n] *= tmp;
    336       tmp *= gammaHi;
    337     }
    338 
    339     /* residual energy */
    340     res_nrg = 0.0;
    341     for (j = 0; j <= ORDERHI; j++) {
    342       for (n = 0; n <= j; n++) {
    343         res_nrg += a_HI[j] * corrhi[j-n] * a_HI[n];
    344       }
    345       for (n = j+1; n <= ORDERHI; n++) {
    346         res_nrg += a_HI[j] * corrhi[n-j] * a_HI[n];
    347       }
    348     }
    349 
    350     /* add hearing threshold and compute of the gain */
    351     *hi_coeff++ = S_N_R / (sqrt(res_nrg) / varscale + H_T_H);
    352 
    353     /* copy coefficients to output array */
    354     for (n = 1; n <= ORDERHI; n++) {
    355       *hi_coeff++ = a_HI[n];
    356     }
    357   }
    358 }
    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 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   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   const double gamma = 0.9;
    406 
    407   /* change quallevel depending on pitch gains and level fluctuations */
    408   WebRtcIsac_GetVarsUB(inSignal, &(maskdata->OldEnergy), varscale);
    409 
    410   /* replace data in buffer by new look-ahead data */
    411   for(frameCntr = 0, activeFrameCntr = 0; frameCntr < numSubFrames;
    412       frameCntr++)
    413   {
    414     if(frameCntr == SUBFRAMES)
    415     {
    416       // we are in 16 kHz
    417       varscale++;
    418       WebRtcIsac_GetVarsUB(&inSignal[FRAMESAMPLES_HALF],
    419                           &(maskdata->OldEnergy), varscale);
    420     }
    421     /* Update input buffer and multiply signal with window */
    422     for(pos1 = 0; pos1 < WINLEN - UPDATE/2; pos1++)
    423     {
    424       maskdata->DataBufferLo[pos1] = maskdata->DataBufferLo[pos1 +
    425                                                             UPDATE/2];
    426       data[pos1] = maskdata->DataBufferLo[pos1] * kLpcCorrWindow[pos1];
    427     }
    428     pos2 = frameCntr * UPDATE/2;
    429     for(n = 0; n < UPDATE/2; n++, pos1++, pos2++)
    430     {
    431       maskdata->DataBufferLo[pos1] = inSignal[pos2];
    432       data[pos1] = maskdata->DataBufferLo[pos1] * kLpcCorrWindow[pos1];
    433     }
    434 
    435     /* Get correlation coefficients */
    436     /* computing autocorrelation    */
    437     WebRtcIsac_AutoCorr(corrSubFrame, data, WINLEN, UB_LPC_ORDER+1);
    438     memcpy(corrMat[frameCntr], corrSubFrame,
    439            (UB_LPC_ORDER+1)*sizeof(double));
    440 
    441     criterion1 = ((frameCntr == 0) || (frameCntr == (SUBFRAMES - 1))) &&
    442         (bandwidth == isac12kHz);
    443     criterion2 = (((frameCntr+1) % 4) == 0) &&
    444         (bandwidth == isac16kHz);
    445     if(criterion1 || criterion2)
    446     {
    447       /* add noise */
    448       corrSubFrame[0] += 1e-6;
    449       /* compute prediction coefficients */
    450       WebRtcIsac_LevDurb(aPolynom, reflecCoeff, corrSubFrame,
    451                         UB_LPC_ORDER);
    452 
    453       /* bandwidth expansion */
    454       tmp = gamma;
    455       for (n = 1; n <= UB_LPC_ORDER; n++)
    456       {
    457         *lpCoeff++ = aPolynom[n] * tmp;
    458         tmp *= gamma;
    459       }
    460       activeFrameCntr++;
    461     }
    462   }
    463 }
    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 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   const double HearThresOffset = -28.0;
    501   const double H_T_H = pow(10.0, 0.05 * HearThresOffset);
    502   /* divide by sqrt(12) = 3.46 */
    503   const double S_N_R = pow(10.0, 0.05 * signal_noise_ratio) / 3.46;
    504 
    505   aPolynom[0] = 1;
    506   for(subFrameCntr = 0; subFrameCntr < numVecs; subFrameCntr++)
    507   {
    508     if(subFrameCntr == SUBFRAMES)
    509     {
    510       // we are in second half of a SWB frame. use new varscale
    511       varscale++;
    512     }
    513     memcpy(&aPolynom[1], &filtCoeffVecs[(subFrameCntr * (UB_LPC_ORDER + 1)) +
    514                                         1], sizeof(double) * UB_LPC_ORDER);
    515 
    516     /* residual energy */
    517     res_nrg = 0.0;
    518     for(j = 0; j <= UB_LPC_ORDER; j++)
    519     {
    520       for(n = 0; n <= j; n++)
    521       {
    522         res_nrg += aPolynom[j] * corrMat[subFrameCntr][j-n] *
    523             aPolynom[n];
    524       }
    525       for(n = j+1; n <= UB_LPC_ORDER; n++)
    526       {
    527         res_nrg += aPolynom[j] * corrMat[subFrameCntr][n-j] *
    528             aPolynom[n];
    529       }
    530     }
    531 
    532     /* add hearing threshold and compute the gain */
    533     gain[subFrameCntr] = S_N_R / (sqrt(res_nrg) / *varscale + H_T_H);
    534   }
    535 }
    536