ac3enc: shift coefficients to 24-bit following MDCT rather than using an exponent...
[ffmpeg.git] / libavcodec / ac3enc_fixed.c
1 /*
2  * The simplest AC-3 encoder
3  * Copyright (c) 2000 Fabrice Bellard
4  * Copyright (c) 2006-2010 Justin Ruggles <justin.ruggles@gmail.com>
5  * Copyright (c) 2006-2010 Prakash Punnoor <prakash@punnoor.de>
6  *
7  * This file is part of FFmpeg.
8  *
9  * FFmpeg is free software; you can redistribute it and/or
10  * modify it under the terms of the GNU Lesser General Public
11  * License as published by the Free Software Foundation; either
12  * version 2.1 of the License, or (at your option) any later version.
13  *
14  * FFmpeg is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17  * Lesser General Public License for more details.
18  *
19  * You should have received a copy of the GNU Lesser General Public
20  * License along with FFmpeg; if not, write to the Free Software
21  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
22  */
23
24 /**
25  * @file
26  * fixed-point AC-3 encoder.
27  */
28
29 #undef CONFIG_AC3ENC_FLOAT
30 #include "ac3enc.c"
31
32
33 /** Scale a float value by 2^15, convert to an integer, and clip to range -32767..32767. */
34 #define FIX15(a) av_clip(SCALE_FLOAT(a, 15), -32767, 32767)
35
36
37 /**
38  * Finalize MDCT and free allocated memory.
39  */
40 static av_cold void mdct_end(AC3MDCTContext *mdct)
41 {
42     mdct->nbits = 0;
43     av_freep(&mdct->costab);
44     av_freep(&mdct->sintab);
45     av_freep(&mdct->xcos1);
46     av_freep(&mdct->xsin1);
47     av_freep(&mdct->rot_tmp);
48     av_freep(&mdct->cplx_tmp);
49 }
50
51
52 /**
53  * Initialize FFT tables.
54  * @param ln log2(FFT size)
55  */
56 static av_cold int fft_init(AVCodecContext *avctx, AC3MDCTContext *mdct, int ln)
57 {
58     int i, n, n2;
59     float alpha;
60
61     n  = 1 << ln;
62     n2 = n >> 1;
63
64     FF_ALLOC_OR_GOTO(avctx, mdct->costab, n2 * sizeof(*mdct->costab), fft_alloc_fail);
65     FF_ALLOC_OR_GOTO(avctx, mdct->sintab, n2 * sizeof(*mdct->sintab), fft_alloc_fail);
66
67     for (i = 0; i < n2; i++) {
68         alpha     = 2.0 * M_PI * i / n;
69         mdct->costab[i] = FIX15(cos(alpha));
70         mdct->sintab[i] = FIX15(sin(alpha));
71     }
72
73     return 0;
74 fft_alloc_fail:
75     mdct_end(mdct);
76     return AVERROR(ENOMEM);
77 }
78
79
80 /**
81  * Initialize MDCT tables.
82  * @param nbits log2(MDCT size)
83  */
84 static av_cold int mdct_init(AVCodecContext *avctx, AC3MDCTContext *mdct,
85                              int nbits)
86 {
87     int i, n, n4, ret;
88
89     n  = 1 << nbits;
90     n4 = n >> 2;
91
92     mdct->nbits = nbits;
93
94     ret = fft_init(avctx, mdct, nbits - 2);
95     if (ret)
96         return ret;
97
98     mdct->window = ff_ac3_window;
99
100     FF_ALLOC_OR_GOTO(avctx, mdct->xcos1,    n4 * sizeof(*mdct->xcos1),    mdct_alloc_fail);
101     FF_ALLOC_OR_GOTO(avctx, mdct->xsin1,    n4 * sizeof(*mdct->xsin1),    mdct_alloc_fail);
102     FF_ALLOC_OR_GOTO(avctx, mdct->rot_tmp,  n  * sizeof(*mdct->rot_tmp),  mdct_alloc_fail);
103     FF_ALLOC_OR_GOTO(avctx, mdct->cplx_tmp, n4 * sizeof(*mdct->cplx_tmp), mdct_alloc_fail);
104
105     for (i = 0; i < n4; i++) {
106         float alpha = 2.0 * M_PI * (i + 1.0 / 8.0) / n;
107         mdct->xcos1[i] = FIX15(-cos(alpha));
108         mdct->xsin1[i] = FIX15(-sin(alpha));
109     }
110
111     return 0;
112 mdct_alloc_fail:
113     mdct_end(mdct);
114     return AVERROR(ENOMEM);
115 }
116
117
118 /** Butterfly op */
119 #define BF(pre, pim, qre, qim, pre1, pim1, qre1, qim1)  \
120 {                                                       \
121   int ax, ay, bx, by;                                   \
122   bx  = pre1;                                           \
123   by  = pim1;                                           \
124   ax  = qre1;                                           \
125   ay  = qim1;                                           \
126   pre = (bx + ax) >> 1;                                 \
127   pim = (by + ay) >> 1;                                 \
128   qre = (bx - ax) >> 1;                                 \
129   qim = (by - ay) >> 1;                                 \
130 }
131
132
133 /** Complex multiply */
134 #define CMUL(pre, pim, are, aim, bre, bim)              \
135 {                                                       \
136    pre = (MUL16(are, bre) - MUL16(aim, bim)) >> 15;     \
137    pim = (MUL16(are, bim) + MUL16(bre, aim)) >> 15;     \
138 }
139
140
141 /**
142  * Calculate a 2^n point complex FFT on 2^ln points.
143  * @param z  complex input/output samples
144  * @param ln log2(FFT size)
145  */
146 static void fft(AC3MDCTContext *mdct, IComplex *z, int ln)
147 {
148     int j, l, np, np2;
149     int nblocks, nloops;
150     register IComplex *p,*q;
151     int tmp_re, tmp_im;
152
153     np = 1 << ln;
154
155     /* reverse */
156     for (j = 0; j < np; j++) {
157         int k = av_reverse[j] >> (8 - ln);
158         if (k < j)
159             FFSWAP(IComplex, z[k], z[j]);
160     }
161
162     /* pass 0 */
163
164     p = &z[0];
165     j = np >> 1;
166     do {
167         BF(p[0].re, p[0].im, p[1].re, p[1].im,
168            p[0].re, p[0].im, p[1].re, p[1].im);
169         p += 2;
170     } while (--j);
171
172     /* pass 1 */
173
174     p = &z[0];
175     j = np >> 2;
176     do {
177         BF(p[0].re, p[0].im, p[2].re,  p[2].im,
178            p[0].re, p[0].im, p[2].re,  p[2].im);
179         BF(p[1].re, p[1].im, p[3].re,  p[3].im,
180            p[1].re, p[1].im, p[3].im, -p[3].re);
181         p+=4;
182     } while (--j);
183
184     /* pass 2 .. ln-1 */
185
186     nblocks = np >> 3;
187     nloops  =  1 << 2;
188     np2     = np >> 1;
189     do {
190         p = z;
191         q = z + nloops;
192         for (j = 0; j < nblocks; j++) {
193             BF(p->re, p->im, q->re, q->im,
194                p->re, p->im, q->re, q->im);
195             p++;
196             q++;
197             for(l = nblocks; l < np2; l += nblocks) {
198                 CMUL(tmp_re, tmp_im, mdct->costab[l], -mdct->sintab[l], q->re, q->im);
199                 BF(p->re, p->im, q->re,  q->im,
200                    p->re, p->im, tmp_re, tmp_im);
201                 p++;
202                 q++;
203             }
204             p += nloops;
205             q += nloops;
206         }
207         nblocks = nblocks >> 1;
208         nloops  = nloops  << 1;
209     } while (nblocks);
210 }
211
212
213 /**
214  * Calculate a 512-point MDCT
215  * @param out 256 output frequency coefficients
216  * @param in  512 windowed input audio samples
217  */
218 static void mdct512(AC3MDCTContext *mdct, int32_t *out, int16_t *in)
219 {
220     int i, re, im, n, n2, n4;
221     int16_t *rot = mdct->rot_tmp;
222     IComplex *x  = mdct->cplx_tmp;
223
224     n  = 1 << mdct->nbits;
225     n2 = n >> 1;
226     n4 = n >> 2;
227
228     /* shift to simplify computations */
229     for (i = 0; i <n4; i++)
230         rot[i] = -in[i + 3*n4];
231     memcpy(&rot[n4], &in[0], 3*n4*sizeof(*in));
232
233     /* pre rotation */
234     for (i = 0; i < n4; i++) {
235         re =  ((int)rot[   2*i] - (int)rot[ n-1-2*i]) >> 1;
236         im = -((int)rot[n2+2*i] - (int)rot[n2-1-2*i]) >> 1;
237         CMUL(x[i].re, x[i].im, re, im, -mdct->xcos1[i], mdct->xsin1[i]);
238     }
239
240     fft(mdct, x, mdct->nbits - 2);
241
242     /* post rotation */
243     for (i = 0; i < n4; i++) {
244         re = x[i].re;
245         im = x[i].im;
246         CMUL(out[n2-1-2*i], out[2*i], re, im, mdct->xsin1[i], mdct->xcos1[i]);
247     }
248 }
249
250
251 /**
252  * Apply KBD window to input samples prior to MDCT.
253  */
254 static void apply_window(DSPContext *dsp, int16_t *output, const int16_t *input,
255                          const int16_t *window, int n)
256 {
257     int i;
258     int n2 = n >> 1;
259
260     for (i = 0; i < n2; i++) {
261         output[i]     = MUL16(input[i],     window[i]) >> 15;
262         output[n-i-1] = MUL16(input[n-i-1], window[i]) >> 15;
263     }
264 }
265
266
267 /**
268  * Calculate the log2() of the maximum absolute value in an array.
269  * @param tab input array
270  * @param n   number of values in the array
271  * @return    log2(max(abs(tab[])))
272  */
273 static int log2_tab(AC3EncodeContext *s, int16_t *src, int len)
274 {
275     int v = s->ac3dsp.ac3_max_msb_abs_int16(src, len);
276     return av_log2(v);
277 }
278
279
280 /**
281  * Left-shift each value in an array by a specified amount.
282  * @param tab    input array
283  * @param n      number of values in the array
284  * @param lshift left shift amount
285  */
286 static void lshift_tab(int16_t *tab, int n, unsigned int lshift)
287 {
288     int i;
289
290     if (lshift > 0) {
291         for (i = 0; i < n; i++)
292             tab[i] <<= lshift;
293     }
294 }
295
296
297 /**
298  * Shift each value in an array by a specified amount.
299  * @param src    input array
300  * @param n      number of values in the array
301  * @param shift  shift amount (negative=right, positive=left)
302  */
303 static void shift_int32(int32_t *src, int n, int shift)
304 {
305     int i;
306
307     if (shift > 0) {
308         for (i = 0; i < n; i++)
309             src[i] <<= shift;
310     } else if (shift < 0) {
311         shift = -shift;
312         for (i = 0; i < n; i++)
313             src[i] >>= shift;
314     }
315 }
316
317
318 /**
319  * Normalize the input samples to use the maximum available precision.
320  * This assumes signed 16-bit input samples.
321  *
322  * @return coefficient shift
323  */
324 static int normalize_samples(AC3EncodeContext *s)
325 {
326     int v = 14 - log2_tab(s, s->windowed_samples, AC3_WINDOW_SIZE);
327     lshift_tab(s->windowed_samples, AC3_WINDOW_SIZE, v);
328     return 9 - v;
329 }
330
331
332 /**
333  * Scale MDCT coefficients to 24-bit fixed-point.
334  */
335 static void scale_coefficients(AC3EncodeContext *s)
336 {
337     int blk, ch;
338
339     for (blk = 0; blk < AC3_MAX_BLOCKS; blk++) {
340         AC3Block *block = &s->blocks[blk];
341         for (ch = 0; ch < s->channels; ch++) {
342             shift_int32(block->mdct_coef[ch], AC3_MAX_COEFS,
343                         block->coeff_shift[ch]);
344         }
345     }
346 }
347
348
349 #ifdef TEST
350 /*************************************************************************/
351 /* TEST */
352
353 #include "libavutil/lfg.h"
354
355 #define MDCT_NBITS 9
356 #define MDCT_SAMPLES (1 << MDCT_NBITS)
357 #define FN (MDCT_SAMPLES/4)
358
359
360 static void fft_test(AC3MDCTContext *mdct, AVLFG *lfg)
361 {
362     IComplex in[FN], in1[FN];
363     int k, n, i;
364     float sum_re, sum_im, a;
365
366     for (i = 0; i < FN; i++) {
367         in[i].re = av_lfg_get(lfg) % 65535 - 32767;
368         in[i].im = av_lfg_get(lfg) % 65535 - 32767;
369         in1[i]   = in[i];
370     }
371     fft(mdct, in, 7);
372
373     /* do it by hand */
374     for (k = 0; k < FN; k++) {
375         sum_re = 0;
376         sum_im = 0;
377         for (n = 0; n < FN; n++) {
378             a = -2 * M_PI * (n * k) / FN;
379             sum_re += in1[n].re * cos(a) - in1[n].im * sin(a);
380             sum_im += in1[n].re * sin(a) + in1[n].im * cos(a);
381         }
382         av_log(NULL, AV_LOG_DEBUG, "%3d: %6d,%6d %6.0f,%6.0f\n",
383                k, in[k].re, in[k].im, sum_re / FN, sum_im / FN);
384     }
385 }
386
387
388 static void mdct_test(AC3MDCTContext *mdct, AVLFG *lfg)
389 {
390     int16_t input[MDCT_SAMPLES];
391     int32_t output[AC3_MAX_COEFS];
392     float input1[MDCT_SAMPLES];
393     float output1[AC3_MAX_COEFS];
394     float s, a, err, e, emax;
395     int i, k, n;
396
397     for (i = 0; i < MDCT_SAMPLES; i++) {
398         input[i]  = (av_lfg_get(lfg) % 65535 - 32767) * 9 / 10;
399         input1[i] = input[i];
400     }
401
402     mdct512(mdct, output, input);
403
404     /* do it by hand */
405     for (k = 0; k < AC3_MAX_COEFS; k++) {
406         s = 0;
407         for (n = 0; n < MDCT_SAMPLES; n++) {
408             a = (2*M_PI*(2*n+1+MDCT_SAMPLES/2)*(2*k+1) / (4 * MDCT_SAMPLES));
409             s += input1[n] * cos(a);
410         }
411         output1[k] = -2 * s / MDCT_SAMPLES;
412     }
413
414     err  = 0;
415     emax = 0;
416     for (i = 0; i < AC3_MAX_COEFS; i++) {
417         av_log(NULL, AV_LOG_DEBUG, "%3d: %7d %7.0f\n", i, output[i], output1[i]);
418         e = output[i] - output1[i];
419         if (e > emax)
420             emax = e;
421         err += e * e;
422     }
423     av_log(NULL, AV_LOG_DEBUG, "err2=%f emax=%f\n", err / AC3_MAX_COEFS, emax);
424 }
425
426
427 int main(void)
428 {
429     AVLFG lfg;
430     AC3MDCTContext mdct;
431
432     mdct.avctx = NULL;
433     av_log_set_level(AV_LOG_DEBUG);
434     mdct_init(&mdct, 9);
435
436     fft_test(&mdct, &lfg);
437     mdct_test(&mdct, &lfg);
438
439     return 0;
440 }
441 #endif /* TEST */
442
443
444 AVCodec ff_ac3_fixed_encoder = {
445     "ac3_fixed",
446     AVMEDIA_TYPE_AUDIO,
447     CODEC_ID_AC3,
448     sizeof(AC3EncodeContext),
449     ac3_encode_init,
450     ac3_encode_frame,
451     ac3_encode_close,
452     NULL,
453     .sample_fmts = (const enum AVSampleFormat[]){AV_SAMPLE_FMT_S16,AV_SAMPLE_FMT_NONE},
454     .long_name = NULL_IF_CONFIG_SMALL("ATSC A/52A (AC-3)"),
455     .channel_layouts = ac3_channel_layouts,
456 };