git.sesse.net Git - ffmpeg/blob - libavcodec/g729postfilter.c

   1 /*
   2  * G.729, G729 Annex D postfilter
   3  * Copyright (c) 2008 Vladimir Voroshilov
   4  *
   5  * This file is part of FFmpeg.
   6  *
   7  * FFmpeg is free software; you can redistribute it and/or
   8  * modify it under the terms of the GNU Lesser General Public
   9  * License as published by the Free Software Foundation; either
  10  * version 2.1 of the License, or (at your option) any later version.
  11  *
  12  * FFmpeg is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15  * Lesser General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU Lesser General Public
  18  * License along with FFmpeg; if not, write to the Free Software
  19  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  20  */
  21 #include <inttypes.h>
  22 #include <limits.h>
  23
  24 #include "avcodec.h"
  25 #include "g729.h"
  26 #include "acelp_pitch_delay.h"
  27 #include "g729postfilter.h"
  28 #include "celp_math.h"
  29 #include "acelp_filters.h"
  30 #include "acelp_vectors.h"
  31 #include "celp_filters.h"
  32
  33 #define FRAC_BITS 15
  34 #include "mathops.h"
  35
  36 /**
  37  * short interpolation filter (of length 33, according to spec)
  38  * for computing signal with non-integer delay
  39  */
  40 static const int16_t ff_g729_interp_filt_short[(ANALYZED_FRAC_DELAYS+1)*SHORT_INT_FILT_LEN] = {
  41       0, 31650, 28469, 23705, 18050, 12266,  7041,  2873,
  42       0, -1597, -2147, -1992, -1492,  -933,  -484,  -188,
  43 };
  44
  45 /**
  46  * long interpolation filter (of length 129, according to spec)
  47  * for computing signal with non-integer delay
  48  */
  49 static const int16_t ff_g729_interp_filt_long[(ANALYZED_FRAC_DELAYS+1)*LONG_INT_FILT_LEN] = {
  50    0, 31915, 29436, 25569, 20676, 15206,  9639,  4439,
  51    0, -3390, -5579, -6549, -6414, -5392, -3773, -1874,
  52    0,  1595,  2727,  3303,  3319,  2850,  2030,  1023,
  53    0,  -887, -1527, -1860, -1876, -1614, -1150,  -579,
  54    0,   501,   859,  1041,  1044,   892,   631,   315,
  55    0,  -266,  -453,  -543,  -538,  -455,  -317,  -156,
  56    0,   130,   218,   258,   253,   212,   147,    72,
  57    0,   -59,  -101,  -122,  -123,  -106,   -77,   -40,
  58 };
  59
  60 /**
  61  * formant_pp_factor_num_pow[i] = FORMANT_PP_FACTOR_NUM^(i+1)
  62  */
  63 static const int16_t formant_pp_factor_num_pow[10]= {
  64   /* (0.15) */
  65   18022, 9912, 5451, 2998, 1649, 907, 499, 274, 151, 83
  66 };
  67
  68 /**
  69  * formant_pp_factor_den_pow[i] = FORMANT_PP_FACTOR_DEN^(i+1)
  70  */
  71 static const int16_t formant_pp_factor_den_pow[10] = {
  72   /* (0.15) */
  73   22938, 16057, 11240, 7868, 5508, 3856, 2699, 1889, 1322, 925
  74 };
  75
  76 /**
  77  * \brief Residual signal calculation (4.2.1 if G.729)
  78  * \param out [out] output data filtered through A(z/FORMANT_PP_FACTOR_NUM)
  79  * \param filter_coeffs (3.12) A(z/FORMANT_PP_FACTOR_NUM) filter coefficients
  80  * \param in input speech data to process
  81  * \param subframe_size size of one subframe
  82  *
  83  * \note in buffer must contain 10 items of previous speech data before top of the buffer
  84  * \remark It is safe to pass the same buffer for input and output.
  85  */
  86 static void residual_filter(int16_t* out, const int16_t* filter_coeffs, const int16_t* in,
  87                             int subframe_size)
  88 {
  89     int i, n;
  90
  91     for (n = subframe_size - 1; n >= 0; n--) {
  92         int sum = 0x800;
  93         for (i = 0; i < 10; i++)
  94             sum += filter_coeffs[i] * in[n - i - 1];
  95
  96         out[n] = in[n] + (sum >> 12);
  97     }
  98 }
  99
 100 /**
 101  * \brief long-term postfilter (4.2.1)
 102  * \param dsp initialized DSP context
 103  * \param pitch_delay_int integer part of the pitch delay in the first subframe
 104  * \param residual filtering input data
 105  * \param residual_filt [out] speech signal with applied A(z/FORMANT_PP_FACTOR_NUM) filter
 106  * \param subframe_size size of subframe
 107  *
 108  * \return 0 if long-term prediction gain is less than 3dB, 1 -  otherwise
 109  */
 110 static int16_t long_term_filter(AudioDSPContext *adsp, int pitch_delay_int,
 111                                 const int16_t* residual, int16_t *residual_filt,
 112                                 int subframe_size)
 113 {
 114     int i, k, tmp, tmp2;
 115     int sum;
 116     int L_temp0;
 117     int L_temp1;
 118     int64_t L64_temp0;
 119     int64_t L64_temp1;
 120     int16_t shift;
 121     int corr_int_num, corr_int_den;
 122
 123     int ener;
 124     int16_t sh_ener;
 125
 126     int16_t gain_num,gain_den; //selected signal's gain numerator and denominator
 127     int16_t sh_gain_num, sh_gain_den;
 128     int gain_num_square;
 129
 130     int16_t gain_long_num,gain_long_den; //filtered through long interpolation filter signal's gain numerator and denominator
 131     int16_t sh_gain_long_num, sh_gain_long_den;
 132
 133     int16_t best_delay_int, best_delay_frac;
 134
 135     int16_t delayed_signal_offset;
 136     int lt_filt_factor_a, lt_filt_factor_b;
 137
 138     int16_t * selected_signal;
 139     const int16_t * selected_signal_const; //Necessary to avoid compiler warning
 140
 141     int16_t sig_scaled[SUBFRAME_SIZE + RES_PREV_DATA_SIZE];
 142     int16_t delayed_signal[ANALYZED_FRAC_DELAYS][SUBFRAME_SIZE+1];
 143     int corr_den[ANALYZED_FRAC_DELAYS][2];
 144
 145     tmp = 0;
 146     for(i=0; i<subframe_size + RES_PREV_DATA_SIZE; i++)
 147         tmp |= FFABS(residual[i]);
 148
 149     if(!tmp)
 150         shift = 3;
 151     else
 152         shift = av_log2(tmp) - 11;
 153
 154     if (shift > 0)
 155         for (i = 0; i < subframe_size + RES_PREV_DATA_SIZE; i++)
 156             sig_scaled[i] = residual[i] >> shift;
 157     else
 158         for (i = 0; i < subframe_size + RES_PREV_DATA_SIZE; i++)
 159             sig_scaled[i] = (unsigned)residual[i] << -shift;
 160
 161     /* Start of best delay searching code */
 162     gain_num = 0;
 163
 164     ener = adsp->scalarproduct_int16(sig_scaled + RES_PREV_DATA_SIZE,
 165                                     sig_scaled + RES_PREV_DATA_SIZE,
 166                                     subframe_size);
 167     if (ener) {
 168         sh_ener = av_log2(ener) - 14;
 169         sh_ener = FFMAX(sh_ener, 0);
 170         ener >>= sh_ener;
 171         /* Search for best pitch delay.
 172
 173                        sum{ r(n) * r(k,n) ] }^2
 174            R'(k)^2 := -------------------------
 175                        sum{ r(k,n) * r(k,n) }
 176
 177
 178            R(T)    :=  sum{ r(n) * r(n-T) ] }
 179
 180
 181            where
 182            r(n-T) is integer delayed signal with delay T
 183            r(k,n) is non-integer delayed signal with integer delay best_delay
 184            and fractional delay k */
 185
 186         /* Find integer delay best_delay which maximizes correlation R(T).
 187
 188            This is also equals to numerator of R'(0),
 189            since the fine search (second step) is done with 1/8
 190            precision around best_delay. */
 191         corr_int_num = 0;
 192         best_delay_int = pitch_delay_int - 1;
 193         for (i = pitch_delay_int - 1; i <= pitch_delay_int + 1; i++) {
 194             sum = adsp->scalarproduct_int16(sig_scaled + RES_PREV_DATA_SIZE,
 195                                            sig_scaled + RES_PREV_DATA_SIZE - i,
 196                                            subframe_size);
 197             if (sum > corr_int_num) {
 198                 corr_int_num = sum;
 199                 best_delay_int = i;
 200             }
 201         }
 202         if (corr_int_num) {
 203             /* Compute denominator of pseudo-normalized correlation R'(0). */
 204             corr_int_den = adsp->scalarproduct_int16(sig_scaled + RES_PREV_DATA_SIZE - best_delay_int,
 205                                                      sig_scaled + RES_PREV_DATA_SIZE - best_delay_int,
 206                                                     subframe_size);
 207
 208             /* Compute signals with non-integer delay k (with 1/8 precision),
 209                where k is in [0;6] range.
 210                Entire delay is qual to best_delay+(k+1)/8
 211                This is archieved by applying an interpolation filter of
 212                legth 33 to source signal. */
 213             for (k = 0; k < ANALYZED_FRAC_DELAYS; k++) {
 214                 ff_acelp_interpolate(&delayed_signal[k][0],
 215                                      &sig_scaled[RES_PREV_DATA_SIZE - best_delay_int],
 216                                      ff_g729_interp_filt_short,
 217                                      ANALYZED_FRAC_DELAYS+1,
 218                                      8 - k - 1,
 219                                      SHORT_INT_FILT_LEN,
 220                                      subframe_size + 1);
 221             }
 222
 223             /* Compute denominator of pseudo-normalized correlation R'(k).
 224
 225                  corr_den[k][0] is square root of R'(k) denominator, for int(T) == int(T0)
 226                  corr_den[k][1] is square root of R'(k) denominator, for int(T) == int(T0)+1
 227
 228               Also compute maximum value of above denominators over all k. */
 229             tmp = corr_int_den;
 230             for (k = 0; k < ANALYZED_FRAC_DELAYS; k++) {
 231                 sum = adsp->scalarproduct_int16(&delayed_signal[k][1],
 232                                                &delayed_signal[k][1],
 233                                                subframe_size - 1);
 234                 corr_den[k][0] = sum + delayed_signal[k][0            ] * delayed_signal[k][0            ];
 235                 corr_den[k][1] = sum + delayed_signal[k][subframe_size] * delayed_signal[k][subframe_size];
 236
 237                 tmp = FFMAX3(tmp, corr_den[k][0], corr_den[k][1]);
 238             }
 239
 240             sh_gain_den = av_log2(tmp) - 14;
 241             if (sh_gain_den >= 0) {
 242
 243                 sh_gain_num =  FFMAX(sh_gain_den, sh_ener);
 244                 /* Loop through all k and find delay that maximizes
 245                    R'(k) correlation.
 246                    Search is done in [int(T0)-1; intT(0)+1] range
 247                    with 1/8 precision. */
 248                 delayed_signal_offset = 1;
 249                 best_delay_frac = 0;
 250                 gain_den = corr_int_den >> sh_gain_den;
 251                 gain_num = corr_int_num >> sh_gain_num;
 252                 gain_num_square = gain_num * gain_num;
 253                 for (k = 0; k < ANALYZED_FRAC_DELAYS; k++) {
 254                     for (i = 0; i < 2; i++) {
 255                         int16_t gain_num_short, gain_den_short;
 256                         int gain_num_short_square;
 257                         /* Compute numerator of pseudo-normalized
 258                            correlation R'(k). */
 259                         sum = adsp->scalarproduct_int16(&delayed_signal[k][i],
 260                                                        sig_scaled + RES_PREV_DATA_SIZE,
 261                                                        subframe_size);
 262                         gain_num_short = FFMAX(sum >> sh_gain_num, 0);
 263
 264                         /*
 265                                       gain_num_short_square                gain_num_square
 266                            R'(T)^2 = -----------------------, max R'(T)^2= --------------
 267                                            den                                 gain_den
 268                         */
 269                         gain_num_short_square = gain_num_short * gain_num_short;
 270                         gain_den_short = corr_den[k][i] >> sh_gain_den;
 271
 272                         tmp = MULL(gain_num_short_square, gain_den, FRAC_BITS);
 273                         tmp2 = MULL(gain_num_square, gain_den_short, FRAC_BITS);
 274
 275                         // R'(T)^2 > max R'(T)^2
 276                         if (tmp > tmp2) {
 277                             gain_num = gain_num_short;
 278                             gain_den = gain_den_short;
 279                             gain_num_square = gain_num_short_square;
 280                             delayed_signal_offset = i;
 281                             best_delay_frac = k + 1;
 282                         }
 283                     }
 284                 }
 285
 286                 /*
 287                        R'(T)^2
 288                   2 * --------- < 1
 289                         R(0)
 290                 */
 291                 L64_temp0 =  (int64_t)gain_num_square  << ((sh_gain_num << 1) + 1);
 292                 L64_temp1 = ((int64_t)gain_den * ener) << (sh_gain_den + sh_ener);
 293                 if (L64_temp0 < L64_temp1)
 294                     gain_num = 0;
 295             } // if(sh_gain_den >= 0)
 296         } // if(corr_int_num)
 297     } // if(ener)
 298     /* End of best delay searching code  */
 299
 300     if (!gain_num) {
 301         memcpy(residual_filt, residual + RES_PREV_DATA_SIZE, subframe_size * sizeof(int16_t));
 302
 303         /* Long-term prediction gain is less than 3dB. Long-term postfilter is disabled. */
 304         return 0;
 305     }
 306     if (best_delay_frac) {
 307         /* Recompute delayed signal with an interpolation filter of length 129. */
 308         ff_acelp_interpolate(residual_filt,
 309                              &sig_scaled[RES_PREV_DATA_SIZE - best_delay_int + delayed_signal_offset],
 310                              ff_g729_interp_filt_long,
 311                              ANALYZED_FRAC_DELAYS + 1,
 312                              8 - best_delay_frac,
 313                              LONG_INT_FILT_LEN,
 314                              subframe_size + 1);
 315         /* Compute R'(k) correlation's numerator. */
 316         sum = adsp->scalarproduct_int16(residual_filt,
 317                                        sig_scaled + RES_PREV_DATA_SIZE,
 318                                        subframe_size);
 319
 320         if (sum < 0) {
 321             gain_long_num = 0;
 322             sh_gain_long_num = 0;
 323         } else {
 324             tmp = av_log2(sum) - 14;
 325             tmp = FFMAX(tmp, 0);
 326             sum >>= tmp;
 327             gain_long_num = sum;
 328             sh_gain_long_num = tmp;
 329         }
 330
 331         /* Compute R'(k) correlation's denominator. */
 332         sum = adsp->scalarproduct_int16(residual_filt, residual_filt, subframe_size);
 333
 334         tmp = av_log2(sum) - 14;
 335         tmp = FFMAX(tmp, 0);
 336         sum >>= tmp;
 337         gain_long_den = sum;
 338         sh_gain_long_den = tmp;
 339
 340         /* Select between original and delayed signal.
 341            Delayed signal will be selected if it increases R'(k)
 342            correlation. */
 343         L_temp0 = gain_num * gain_num;
 344         L_temp0 = MULL(L_temp0, gain_long_den, FRAC_BITS);
 345
 346         L_temp1 = gain_long_num * gain_long_num;
 347         L_temp1 = MULL(L_temp1, gain_den, FRAC_BITS);
 348
 349         tmp = ((sh_gain_long_num - sh_gain_num) * 2) - (sh_gain_long_den - sh_gain_den);
 350         if (tmp > 0)
 351             L_temp0 >>= tmp;
 352         else
 353             L_temp1 >>= -tmp;
 354
 355         /* Check if longer filter increases the values of R'(k). */
 356         if (L_temp1 > L_temp0) {
 357             /* Select long filter. */
 358             selected_signal = residual_filt;
 359             gain_num = gain_long_num;
 360             gain_den = gain_long_den;
 361             sh_gain_num = sh_gain_long_num;
 362             sh_gain_den = sh_gain_long_den;
 363         } else
 364             /* Select short filter. */
 365             selected_signal = &delayed_signal[best_delay_frac-1][delayed_signal_offset];
 366
 367         /* Rescale selected signal to original value. */
 368         if (shift > 0)
 369             for (i = 0; i < subframe_size; i++)
 370                 selected_signal[i] *= 1 << shift;
 371         else
 372             for (i = 0; i < subframe_size; i++)
 373                 selected_signal[i] >>= -shift;
 374
 375         /* necessary to avoid compiler warning */
 376         selected_signal_const = selected_signal;
 377     } // if(best_delay_frac)
 378     else
 379         selected_signal_const = residual + RES_PREV_DATA_SIZE - (best_delay_int + 1 - delayed_signal_offset);
 380 #ifdef G729_BITEXACT
 381     tmp = sh_gain_num - sh_gain_den;
 382     if (tmp > 0)
 383         gain_den >>= tmp;
 384     else
 385         gain_num >>= -tmp;
 386
 387     if (gain_num > gain_den)
 388         lt_filt_factor_a = MIN_LT_FILT_FACTOR_A;
 389     else {
 390         gain_num >>= 2;
 391         gain_den >>= 1;
 392         lt_filt_factor_a = (gain_den << 15) / (gain_den + gain_num);
 393     }
 394 #else
 395     L64_temp0 = (((int64_t)gain_num) << sh_gain_num) >> 1;
 396     L64_temp1 = ((int64_t)gain_den) << sh_gain_den;
 397     lt_filt_factor_a = FFMAX((L64_temp1 << 15) / (L64_temp1 + L64_temp0), MIN_LT_FILT_FACTOR_A);
 398 #endif
 399
 400     /* Filter through selected filter. */
 401     lt_filt_factor_b = 32767 - lt_filt_factor_a + 1;
 402
 403     ff_acelp_weighted_vector_sum(residual_filt, residual + RES_PREV_DATA_SIZE,
 404                                  selected_signal_const,
 405                                  lt_filt_factor_a, lt_filt_factor_b,
 406                                  1<<14, 15, subframe_size);
 407
 408     // Long-term prediction gain is larger than 3dB.
 409     return 1;
 410 }
 411
 412 /**
 413  * \brief Calculate reflection coefficient for tilt compensation filter (4.2.3).
 414  * \param dsp initialized DSP context
 415  * \param lp_gn (3.12) coefficients of A(z/FORMANT_PP_FACTOR_NUM) filter
 416  * \param lp_gd (3.12) coefficients of A(z/FORMANT_PP_FACTOR_DEN) filter
 417  * \param speech speech to update
 418  * \param subframe_size size of subframe
 419  *
 420  * \return (3.12) reflection coefficient
 421  *
 422  * \remark The routine also calculates the gain term for the short-term
 423  *         filter (gf) and multiplies the speech data by 1/gf.
 424  *
 425  * \note All members of lp_gn, except 10-19 must be equal to zero.
 426  */
 427 static int16_t get_tilt_comp(AudioDSPContext *adsp, int16_t *lp_gn,
 428                              const int16_t *lp_gd, int16_t* speech,
 429                              int subframe_size)
 430 {
 431     int rh1,rh0; // (3.12)
 432     int temp;
 433     int i;
 434     int gain_term;
 435
 436     lp_gn[10] = 4096; //1.0 in (3.12)
 437
 438     /* Apply 1/A(z/FORMANT_PP_FACTOR_DEN) filter to hf. */
 439     ff_celp_lp_synthesis_filter(lp_gn + 11, lp_gd + 1, lp_gn + 11, 22, 10, 0, 0, 0x800);
 440     /* Now lp_gn (starting with 10) contains impulse response
 441        of A(z/FORMANT_PP_FACTOR_NUM)/A(z/FORMANT_PP_FACTOR_DEN) filter. */
 442
 443     rh0 = adsp->scalarproduct_int16(lp_gn + 10, lp_gn + 10, 20);
 444     rh1 = adsp->scalarproduct_int16(lp_gn + 10, lp_gn + 11, 20);
 445
 446     /* downscale to avoid overflow */
 447     temp = av_log2(rh0) - 14;
 448     if (temp > 0) {
 449         rh0 >>= temp;
 450         rh1 >>= temp;
 451     }
 452
 453     if (FFABS(rh1) > rh0 || !rh0)
 454         return 0;
 455
 456     gain_term = 0;
 457     for (i = 0; i < 20; i++)
 458         gain_term += FFABS(lp_gn[i + 10]);
 459     gain_term >>= 2; // (3.12) -> (5.10)
 460
 461     if (gain_term > 0x400) { // 1.0 in (5.10)
 462         temp = 0x2000000 / gain_term; // 1.0/gain_term in (0.15)
 463         for (i = 0; i < subframe_size; i++)
 464             speech[i] = (speech[i] * temp + 0x4000) >> 15;
 465     }
 466
 467     return -(rh1 * (1 << 15)) / rh0;
 468 }
 469
 470 /**
 471  * \brief Apply tilt compensation filter (4.2.3).
 472  * \param res_pst [in/out] residual signal (partially filtered)
 473  * \param k1 (3.12) reflection coefficient
 474  * \param subframe_size size of subframe
 475  * \param ht_prev_data previous data for 4.2.3, equation 86
 476  *
 477  * \return new value for ht_prev_data
 478 */
 479 static int16_t apply_tilt_comp(int16_t* out, int16_t* res_pst, int refl_coeff,
 480                                int subframe_size, int16_t ht_prev_data)
 481 {
 482     int tmp, tmp2;
 483     int i;
 484     int gt, ga;
 485     int fact, sh_fact;
 486
 487     if (refl_coeff > 0) {
 488         gt = (refl_coeff * G729_TILT_FACTOR_PLUS + 0x4000) >> 15;
 489         fact = 0x2000; // 0.5 in (0.15)
 490         sh_fact = 14;
 491     } else {
 492         gt = (refl_coeff * G729_TILT_FACTOR_MINUS + 0x4000) >> 15;
 493         fact = 0x400; // 0.5 in (3.12)
 494         sh_fact = 11;
 495     }
 496     ga = (fact << 16) / av_clip_int16(32768 - FFABS(gt));
 497     gt >>= 1;
 498
 499     /* Apply tilt compensation filter to signal. */
 500     tmp = res_pst[subframe_size - 1];
 501
 502     for (i = subframe_size - 1; i >= 1; i--) {
 503         tmp2 = (gt * res_pst[i-1]) * 2 + 0x4000;
 504         tmp2 = res_pst[i] + (tmp2 >> 15);
 505
 506         tmp2 = (tmp2 * ga + fact) >> sh_fact;
 507         out[i] = tmp2;
 508     }
 509     tmp2 = (gt * ht_prev_data) * 2 + 0x4000;
 510     tmp2 = res_pst[0] + (tmp2 >> 15);
 511     tmp2 = (tmp2 * ga + fact) >> sh_fact;
 512     out[0] = tmp2;
 513
 514     return tmp;
 515 }
 516
 517 void ff_g729_postfilter(AudioDSPContext *adsp, int16_t* ht_prev_data, int* voicing,
 518                      const int16_t *lp_filter_coeffs, int pitch_delay_int,
 519                      int16_t* residual, int16_t* res_filter_data,
 520                      int16_t* pos_filter_data, int16_t *speech, int subframe_size)
 521 {
 522     int16_t residual_filt_buf[SUBFRAME_SIZE+11];
 523     int16_t lp_gn[33]; // (3.12)
 524     int16_t lp_gd[11]; // (3.12)
 525     int tilt_comp_coeff;
 526     int i;
 527
 528     /* Zero-filling is necessary for tilt-compensation filter. */
 529     memset(lp_gn, 0, 33 * sizeof(int16_t));
 530
 531     /* Calculate A(z/FORMANT_PP_FACTOR_NUM) filter coefficients. */
 532     for (i = 0; i < 10; i++)
 533         lp_gn[i + 11] = (lp_filter_coeffs[i + 1] * formant_pp_factor_num_pow[i] + 0x4000) >> 15;
 534
 535     /* Calculate A(z/FORMANT_PP_FACTOR_DEN) filter coefficients. */
 536     for (i = 0; i < 10; i++)
 537         lp_gd[i + 1] = (lp_filter_coeffs[i + 1] * formant_pp_factor_den_pow[i] + 0x4000) >> 15;
 538
 539     /* residual signal calculation (one-half of short-term postfilter) */
 540     memcpy(speech - 10, res_filter_data, 10 * sizeof(int16_t));
 541     residual_filter(residual + RES_PREV_DATA_SIZE, lp_gn + 11, speech, subframe_size);
 542     /* Save data to use it in the next subframe. */
 543     memcpy(res_filter_data, speech + subframe_size - 10, 10 * sizeof(int16_t));
 544
 545     /* long-term filter. If long-term prediction gain is larger than 3dB (returned value is
 546        nonzero) then declare current subframe as periodic. */
 547     i = long_term_filter(adsp, pitch_delay_int,
 548                                                 residual, residual_filt_buf + 10,
 549                                                 subframe_size);
 550     *voicing = FFMAX(*voicing, i);
 551
 552     /* shift residual for using in next subframe */
 553     memmove(residual, residual + subframe_size, RES_PREV_DATA_SIZE * sizeof(int16_t));
 554
 555     /* short-term filter tilt compensation */
 556     tilt_comp_coeff = get_tilt_comp(adsp, lp_gn, lp_gd, residual_filt_buf + 10, subframe_size);
 557
 558     /* Apply second half of short-term postfilter: 1/A(z/FORMANT_PP_FACTOR_DEN) */
 559     ff_celp_lp_synthesis_filter(pos_filter_data + 10, lp_gd + 1,
 560                                 residual_filt_buf + 10,
 561                                 subframe_size, 10, 0, 0, 0x800);
 562     memcpy(pos_filter_data, pos_filter_data + subframe_size, 10 * sizeof(int16_t));
 563
 564     *ht_prev_data = apply_tilt_comp(speech, pos_filter_data + 10, tilt_comp_coeff,
 565                                     subframe_size, *ht_prev_data);
 566 }
 567
 568 /**
 569  * \brief Adaptive gain control (4.2.4)
 570  * \param gain_before gain of speech before applying postfilters
 571  * \param gain_after  gain of speech after applying postfilters
 572  * \param speech [in/out] signal buffer
 573  * \param subframe_size length of subframe
 574  * \param gain_prev (3.12) previous value of gain coefficient
 575  *
 576  * \return (3.12) last value of gain coefficient
 577  */
 578 int16_t ff_g729_adaptive_gain_control(int gain_before, int gain_after, int16_t *speech,
 579                                    int subframe_size, int16_t gain_prev)
 580 {
 581     int gain; // (3.12)
 582     int n;
 583     int exp_before, exp_after;
 584
 585     if(!gain_after && gain_before)
 586         return 0;
 587
 588     if (gain_before) {
 589
 590         exp_before  = 14 - av_log2(gain_before);
 591         gain_before = bidir_sal(gain_before, exp_before);
 592
 593         exp_after  = 14 - av_log2(gain_after);
 594         gain_after = bidir_sal(gain_after, exp_after);
 595
 596         if (gain_before < gain_after) {
 597             gain = (gain_before << 15) / gain_after;
 598             gain = bidir_sal(gain, exp_after - exp_before - 1);
 599         } else {
 600             gain = ((gain_before - gain_after) << 14) / gain_after + 0x4000;
 601             gain = bidir_sal(gain, exp_after - exp_before);
 602         }
 603         gain = av_clip_int16(gain);
 604         gain = (gain * G729_AGC_FAC1 + 0x4000) >> 15; // gain * (1-0.9875)
 605     } else
 606         gain = 0;
 607
 608     for (n = 0; n < subframe_size; n++) {
 609         // gain_prev = gain + 0.9875 * gain_prev
 610         gain_prev = (G729_AGC_FACTOR * gain_prev + 0x4000) >> 15;
 611         gain_prev = av_clip_int16(gain + gain_prev);
 612         speech[n] = av_clip_int16((speech[n] * gain_prev + 0x2000) >> 14);
 613     }
 614     return gain_prev;
 615 }