2 * G.729, G729 Annex D postfilter
3 * Copyright (c) 2008 Vladimir Voroshilov
5 * This file is part of FFmpeg.
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.
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.
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
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"
37 * short interpolation filter (of length 33, according to spec)
38 * for computing signal with non-integer delay
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,
46 * long interpolation filter (of length 129, according to spec)
47 * for computing signal with non-integer delay
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,
61 * formant_pp_factor_num_pow[i] = FORMANT_PP_FACTOR_NUM^(i+1)
63 static const int16_t formant_pp_factor_num_pow
[10]= {
65 18022, 9912, 5451, 2998, 1649, 907, 499, 274, 151, 83
69 * formant_pp_factor_den_pow[i] = FORMANT_PP_FACTOR_DEN^(i+1)
71 static const int16_t formant_pp_factor_den_pow
[10] = {
73 22938, 16057, 11240, 7868, 5508, 3856, 2699, 1889, 1322, 925
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
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.
86 static void residual_filter(int16_t* out
, const int16_t* filter_coeffs
, const int16_t* in
,
91 for (n
= subframe_size
- 1; n
>= 0; n
--) {
93 for (i
= 0; i
< 10; i
++)
94 sum
+= filter_coeffs
[i
] * in
[n
- i
- 1];
96 out
[n
] = in
[n
] + (sum
>> 12);
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
108 * \return 0 if long-term prediction gain is less than 3dB, 1 - otherwise
110 static int16_t long_term_filter(AudioDSPContext
*adsp
, int pitch_delay_int
,
111 const int16_t* residual
, int16_t *residual_filt
,
121 int corr_int_num
, corr_int_den
;
126 int16_t gain_num
,gain_den
; //selected signal's gain numerator and denominator
127 int16_t sh_gain_num
, sh_gain_den
;
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
;
133 int16_t best_delay_int
, best_delay_frac
;
135 int16_t delayed_signal_offset
;
136 int lt_filt_factor_a
, lt_filt_factor_b
;
138 int16_t * selected_signal
;
139 const int16_t * selected_signal_const
; //Necessary to avoid compiler warning
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];
146 for(i
=0; i
<subframe_size
+ RES_PREV_DATA_SIZE
; i
++)
147 tmp
|= FFABS(residual
[i
]);
152 shift
= av_log2(tmp
) - 11;
155 for (i
= 0; i
< subframe_size
+ RES_PREV_DATA_SIZE
; i
++)
156 sig_scaled
[i
] = residual
[i
] >> shift
;
158 for (i
= 0; i
< subframe_size
+ RES_PREV_DATA_SIZE
; i
++)
159 sig_scaled
[i
] = residual
[i
] << -shift
;
161 /* Start of best delay searching code */
164 ener
= adsp
->scalarproduct_int16(sig_scaled
+ RES_PREV_DATA_SIZE
,
165 sig_scaled
+ RES_PREV_DATA_SIZE
,
168 sh_ener
= FFMAX(av_log2(ener
) - 14, 0);
170 /* Search for best pitch delay.
172 sum{ r(n) * r(k,n) ] }^2
173 R'(k)^2 := -------------------------
174 sum{ r(k,n) * r(k,n) }
177 R(T) := sum{ r(n) * r(n-T) ] }
181 r(n-T) is integer delayed signal with delay T
182 r(k,n) is non-integer delayed signal with integer delay best_delay
183 and fractional delay k */
185 /* Find integer delay best_delay which maximizes correlation R(T).
187 This is also equals to numerator of R'(0),
188 since the fine search (second step) is done with 1/8
189 precision around best_delay. */
191 best_delay_int
= pitch_delay_int
- 1;
192 for (i
= pitch_delay_int
- 1; i
<= pitch_delay_int
+ 1; i
++) {
193 sum
= adsp
->scalarproduct_int16(sig_scaled
+ RES_PREV_DATA_SIZE
,
194 sig_scaled
+ RES_PREV_DATA_SIZE
- i
,
196 if (sum
> corr_int_num
) {
202 /* Compute denominator of pseudo-normalized correlation R'(0). */
203 corr_int_den
= adsp
->scalarproduct_int16(sig_scaled
- best_delay_int
+ RES_PREV_DATA_SIZE
,
204 sig_scaled
- best_delay_int
+ RES_PREV_DATA_SIZE
,
207 /* Compute signals with non-integer delay k (with 1/8 precision),
208 where k is in [0;6] range.
209 Entire delay is qual to best_delay+(k+1)/8
210 This is archieved by applying an interpolation filter of
211 legth 33 to source signal. */
212 for (k
= 0; k
< ANALYZED_FRAC_DELAYS
; k
++) {
213 ff_acelp_interpolate(&delayed_signal
[k
][0],
214 &sig_scaled
[RES_PREV_DATA_SIZE
- best_delay_int
],
215 ff_g729_interp_filt_short
,
216 ANALYZED_FRAC_DELAYS
+1,
222 /* Compute denominator of pseudo-normalized correlation R'(k).
224 corr_den[k][0] is square root of R'(k) denominator, for int(T) == int(T0)
225 corr_den[k][1] is square root of R'(k) denominator, for int(T) == int(T0)+1
227 Also compute maximum value of above denominators over all k. */
229 for (k
= 0; k
< ANALYZED_FRAC_DELAYS
; k
++) {
230 sum
= adsp
->scalarproduct_int16(&delayed_signal
[k
][1],
231 &delayed_signal
[k
][1],
233 corr_den
[k
][0] = sum
+ delayed_signal
[k
][0 ] * delayed_signal
[k
][0 ];
234 corr_den
[k
][1] = sum
+ delayed_signal
[k
][subframe_size
] * delayed_signal
[k
][subframe_size
];
236 tmp
= FFMAX3(tmp
, corr_den
[k
][0], corr_den
[k
][1]);
239 sh_gain_den
= av_log2(tmp
) - 14;
240 if (sh_gain_den
>= 0) {
242 sh_gain_num
= FFMAX(sh_gain_den
, sh_ener
);
243 /* Loop through all k and find delay that maximizes
245 Search is done in [int(T0)-1; intT(0)+1] range
246 with 1/8 precision. */
247 delayed_signal_offset
= 1;
249 gain_den
= corr_int_den
>> sh_gain_den
;
250 gain_num
= corr_int_num
>> sh_gain_num
;
251 gain_num_square
= gain_num
* gain_num
;
252 for (k
= 0; k
< ANALYZED_FRAC_DELAYS
; k
++) {
253 for (i
= 0; i
< 2; i
++) {
254 int16_t gain_num_short
, gain_den_short
;
255 int gain_num_short_square
;
256 /* Compute numerator of pseudo-normalized
257 correlation R'(k). */
258 sum
= adsp
->scalarproduct_int16(&delayed_signal
[k
][i
],
259 sig_scaled
+ RES_PREV_DATA_SIZE
,
261 gain_num_short
= FFMAX(sum
>> sh_gain_num
, 0);
264 gain_num_short_square gain_num_square
265 R'(T)^2 = -----------------------, max R'(T)^2= --------------
268 gain_num_short_square
= gain_num_short
* gain_num_short
;
269 gain_den_short
= corr_den
[k
][i
] >> sh_gain_den
;
271 tmp
= MULL(gain_num_short_square
, gain_den
, FRAC_BITS
);
272 tmp2
= MULL(gain_num_square
, gain_den_short
, FRAC_BITS
);
274 // R'(T)^2 > max R'(T)^2
276 gain_num
= gain_num_short
;
277 gain_den
= gain_den_short
;
278 gain_num_square
= gain_num_short_square
;
279 delayed_signal_offset
= i
;
280 best_delay_frac
= k
+ 1;
290 L64_temp0
= (int64_t)gain_num_square
<< ((sh_gain_num
<< 1) + 1);
291 L64_temp1
= ((int64_t)gain_den
* ener
) << (sh_gain_den
+ sh_ener
);
292 if (L64_temp0
< L64_temp1
)
294 } // if(sh_gain_den >= 0)
295 } // if(corr_int_num)
297 /* End of best delay searching code */
300 memcpy(residual_filt
, residual
+ RES_PREV_DATA_SIZE
, subframe_size
* sizeof(int16_t));
302 /* Long-term prediction gain is less than 3dB. Long-term postfilter is disabled. */
305 if (best_delay_frac
) {
306 /* Recompute delayed signal with an interpolation filter of length 129. */
307 ff_acelp_interpolate(residual_filt
,
308 &sig_scaled
[RES_PREV_DATA_SIZE
- best_delay_int
+ delayed_signal_offset
],
309 ff_g729_interp_filt_long
,
310 ANALYZED_FRAC_DELAYS
+ 1,
314 /* Compute R'(k) correlation's numerator. */
315 sum
= adsp
->scalarproduct_int16(residual_filt
,
316 sig_scaled
+ RES_PREV_DATA_SIZE
,
321 sh_gain_long_num
= 0;
323 tmp
= FFMAX(av_log2(sum
) - 14, 0);
326 sh_gain_long_num
= tmp
;
329 /* Compute R'(k) correlation's denominator. */
330 sum
= adsp
->scalarproduct_int16(residual_filt
, residual_filt
, subframe_size
);
332 tmp
= FFMAX(av_log2(sum
) - 14, 0);
335 sh_gain_long_den
= tmp
;
337 /* Select between original and delayed signal.
338 Delayed signal will be selected if it increases R'(k)
340 L_temp0
= gain_num
* gain_num
;
341 L_temp0
= MULL(L_temp0
, gain_long_den
, FRAC_BITS
);
343 L_temp1
= gain_long_num
* gain_long_num
;
344 L_temp1
= MULL(L_temp1
, gain_den
, FRAC_BITS
);
346 tmp
= ((sh_gain_long_num
- sh_gain_num
) << 1) - (sh_gain_long_den
- sh_gain_den
);
352 /* Check if longer filter increases the values of R'(k). */
353 if (L_temp1
> L_temp0
) {
354 /* Select long filter. */
355 selected_signal
= residual_filt
;
356 gain_num
= gain_long_num
;
357 gain_den
= gain_long_den
;
358 sh_gain_num
= sh_gain_long_num
;
359 sh_gain_den
= sh_gain_long_den
;
361 /* Select short filter. */
362 selected_signal
= &delayed_signal
[best_delay_frac
-1][delayed_signal_offset
];
364 /* Rescale selected signal to original value. */
366 for (i
= 0; i
< subframe_size
; i
++)
367 selected_signal
[i
] <<= shift
;
369 for (i
= 0; i
< subframe_size
; i
++)
370 selected_signal
[i
] >>= -shift
;
372 /* necessary to avoid compiler warning */
373 selected_signal_const
= selected_signal
;
374 } // if(best_delay_frac)
376 selected_signal_const
= residual
+ RES_PREV_DATA_SIZE
- (best_delay_int
+ 1 - delayed_signal_offset
);
378 tmp
= sh_gain_num
- sh_gain_den
;
384 if (gain_num
> gain_den
)
385 lt_filt_factor_a
= MIN_LT_FILT_FACTOR_A
;
389 lt_filt_factor_a
= (gain_den
<< 15) / (gain_den
+ gain_num
);
392 L64_temp0
= (((int64_t)gain_num
) << sh_gain_num
) >> 1;
393 L64_temp1
= ((int64_t)gain_den
) << sh_gain_den
;
394 lt_filt_factor_a
= FFMAX((L64_temp1
<< 15) / (L64_temp1
+ L64_temp0
), MIN_LT_FILT_FACTOR_A
);
397 /* Filter through selected filter. */
398 lt_filt_factor_b
= 32767 - lt_filt_factor_a
+ 1;
400 ff_acelp_weighted_vector_sum(residual_filt
, residual
+ RES_PREV_DATA_SIZE
,
401 selected_signal_const
,
402 lt_filt_factor_a
, lt_filt_factor_b
,
403 1<<14, 15, subframe_size
);
405 // Long-term prediction gain is larger than 3dB.
410 * \brief Calculate reflection coefficient for tilt compensation filter (4.2.3).
411 * \param dsp initialized DSP context
412 * \param lp_gn (3.12) coefficients of A(z/FORMANT_PP_FACTOR_NUM) filter
413 * \param lp_gd (3.12) coefficients of A(z/FORMANT_PP_FACTOR_DEN) filter
414 * \param speech speech to update
415 * \param subframe_size size of subframe
417 * \return (3.12) reflection coefficient
419 * \remark The routine also calculates the gain term for the short-term
420 * filter (gf) and multiplies the speech data by 1/gf.
422 * \note All members of lp_gn, except 10-19 must be equal to zero.
424 static int16_t get_tilt_comp(AudioDSPContext
*adsp
, int16_t *lp_gn
,
425 const int16_t *lp_gd
, int16_t* speech
,
428 int rh1
,rh0
; // (3.12)
433 lp_gn
[10] = 4096; //1.0 in (3.12)
435 /* Apply 1/A(z/FORMANT_PP_FACTOR_DEN) filter to hf. */
436 ff_celp_lp_synthesis_filter(lp_gn
+ 11, lp_gd
+ 1, lp_gn
+ 11, 22, 10, 0, 0, 0x800);
437 /* Now lp_gn (starting with 10) contains impulse response
438 of A(z/FORMANT_PP_FACTOR_NUM)/A(z/FORMANT_PP_FACTOR_DEN) filter. */
440 rh0
= adsp
->scalarproduct_int16(lp_gn
+ 10, lp_gn
+ 10, 20);
441 rh1
= adsp
->scalarproduct_int16(lp_gn
+ 10, lp_gn
+ 11, 20);
443 /* downscale to avoid overflow */
444 temp
= av_log2(rh0
) - 14;
450 if (FFABS(rh1
) > rh0
|| !rh0
)
454 for (i
= 0; i
< 20; i
++)
455 gain_term
+= FFABS(lp_gn
[i
+ 10]);
456 gain_term
>>= 2; // (3.12) -> (5.10)
458 if (gain_term
> 0x400) { // 1.0 in (5.10)
459 temp
= 0x2000000 / gain_term
; // 1.0/gain_term in (0.15)
460 for (i
= 0; i
< subframe_size
; i
++)
461 speech
[i
] = (speech
[i
] * temp
+ 0x4000) >> 15;
464 return -(rh1
<< 15) / rh0
;
468 * \brief Apply tilt compensation filter (4.2.3).
469 * \param res_pst [in/out] residual signal (partially filtered)
470 * \param k1 (3.12) reflection coefficient
471 * \param subframe_size size of subframe
472 * \param ht_prev_data previous data for 4.2.3, equation 86
474 * \return new value for ht_prev_data
476 static int16_t apply_tilt_comp(int16_t* out
, int16_t* res_pst
, int refl_coeff
,
477 int subframe_size
, int16_t ht_prev_data
)
484 if (refl_coeff
> 0) {
485 gt
= (refl_coeff
* G729_TILT_FACTOR_PLUS
+ 0x4000) >> 15;
486 fact
= 0x4000; // 0.5 in (0.15)
489 gt
= (refl_coeff
* G729_TILT_FACTOR_MINUS
+ 0x4000) >> 15;
490 fact
= 0x800; // 0.5 in (3.12)
493 ga
= (fact
<< 15) / av_clip_int16(32768 - FFABS(gt
));
496 /* Apply tilt compensation filter to signal. */
497 tmp
= res_pst
[subframe_size
- 1];
499 for (i
= subframe_size
- 1; i
>= 1; i
--) {
500 tmp2
= (res_pst
[i
] << 15) + ((gt
* res_pst
[i
-1]) << 1);
501 tmp2
= (tmp2
+ 0x4000) >> 15;
503 tmp2
= (tmp2
* ga
* 2 + fact
) >> sh_fact
;
506 tmp2
= (res_pst
[0] << 15) + ((gt
* ht_prev_data
) << 1);
507 tmp2
= (tmp2
+ 0x4000) >> 15;
508 tmp2
= (tmp2
* ga
* 2 + fact
) >> sh_fact
;
514 void ff_g729_postfilter(AudioDSPContext
*adsp
, int16_t* ht_prev_data
, int* voicing
,
515 const int16_t *lp_filter_coeffs
, int pitch_delay_int
,
516 int16_t* residual
, int16_t* res_filter_data
,
517 int16_t* pos_filter_data
, int16_t *speech
, int subframe_size
)
519 int16_t residual_filt_buf
[SUBFRAME_SIZE
+11];
520 int16_t lp_gn
[33]; // (3.12)
521 int16_t lp_gd
[11]; // (3.12)
525 /* Zero-filling is necessary for tilt-compensation filter. */
526 memset(lp_gn
, 0, 33 * sizeof(int16_t));
528 /* Calculate A(z/FORMANT_PP_FACTOR_NUM) filter coefficients. */
529 for (i
= 0; i
< 10; i
++)
530 lp_gn
[i
+ 11] = (lp_filter_coeffs
[i
+ 1] * formant_pp_factor_num_pow
[i
] + 0x4000) >> 15;
532 /* Calculate A(z/FORMANT_PP_FACTOR_DEN) filter coefficients. */
533 for (i
= 0; i
< 10; i
++)
534 lp_gd
[i
+ 1] = (lp_filter_coeffs
[i
+ 1] * formant_pp_factor_den_pow
[i
] + 0x4000) >> 15;
536 /* residual signal calculation (one-half of short-term postfilter) */
537 memcpy(speech
- 10, res_filter_data
, 10 * sizeof(int16_t));
538 residual_filter(residual
+ RES_PREV_DATA_SIZE
, lp_gn
+ 11, speech
, subframe_size
);
539 /* Save data to use it in the next subframe. */
540 memcpy(res_filter_data
, speech
+ subframe_size
- 10, 10 * sizeof(int16_t));
542 /* long-term filter. If long-term prediction gain is larger than 3dB (returned value is
543 nonzero) then declare current subframe as periodic. */
544 *voicing
= FFMAX(*voicing
, long_term_filter(adsp
, pitch_delay_int
,
545 residual
, residual_filt_buf
+ 10,
548 /* shift residual for using in next subframe */
549 memmove(residual
, residual
+ subframe_size
, RES_PREV_DATA_SIZE
* sizeof(int16_t));
551 /* short-term filter tilt compensation */
552 tilt_comp_coeff
= get_tilt_comp(adsp
, lp_gn
, lp_gd
, residual_filt_buf
+ 10, subframe_size
);
554 /* Apply second half of short-term postfilter: 1/A(z/FORMANT_PP_FACTOR_DEN) */
555 ff_celp_lp_synthesis_filter(pos_filter_data
+ 10, lp_gd
+ 1,
556 residual_filt_buf
+ 10,
557 subframe_size
, 10, 0, 0, 0x800);
558 memcpy(pos_filter_data
, pos_filter_data
+ subframe_size
, 10 * sizeof(int16_t));
560 *ht_prev_data
= apply_tilt_comp(speech
, pos_filter_data
+ 10, tilt_comp_coeff
,
561 subframe_size
, *ht_prev_data
);
565 * \brief Adaptive gain control (4.2.4)
566 * \param gain_before gain of speech before applying postfilters
567 * \param gain_after gain of speech after applying postfilters
568 * \param speech [in/out] signal buffer
569 * \param subframe_size length of subframe
570 * \param gain_prev (3.12) previous value of gain coefficient
572 * \return (3.12) last value of gain coefficient
574 int16_t ff_g729_adaptive_gain_control(int gain_before
, int gain_after
, int16_t *speech
,
575 int subframe_size
, int16_t gain_prev
)
579 int exp_before
, exp_after
;
581 if(!gain_after
&& gain_before
)
586 exp_before
= 14 - av_log2(gain_before
);
587 gain_before
= bidir_sal(gain_before
, exp_before
);
589 exp_after
= 14 - av_log2(gain_after
);
590 gain_after
= bidir_sal(gain_after
, exp_after
);
592 if (gain_before
< gain_after
) {
593 gain
= (gain_before
<< 15) / gain_after
;
594 gain
= bidir_sal(gain
, exp_after
- exp_before
- 1);
596 gain
= ((gain_before
- gain_after
) << 14) / gain_after
+ 0x4000;
597 gain
= bidir_sal(gain
, exp_after
- exp_before
);
599 gain
= (gain
* G729_AGC_FAC1
+ 0x4000) >> 15; // gain * (1-0.9875)
603 for (n
= 0; n
< subframe_size
; n
++) {
604 // gain_prev = gain + 0.9875 * gain_prev
605 gain_prev
= (G729_AGC_FACTOR
* gain_prev
+ 0x4000) >> 15;
606 gain_prev
= av_clip_int16(gain
+ gain_prev
);
607 speech
[n
] = av_clip_int16((speech
[n
] * gain_prev
+ 0x2000) >> 14);