Browse code

G.729 postfilter

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