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