git.sesse.net Git - ffmpeg/blob - 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 #include "ac3enc.c"
  30
  31
  32 /** Scale a float value by 2^15, convert to an integer, and clip to range -32767..32767. */
  33 #define FIX15(a) av_clip(SCALE_FLOAT(a, 15), -32767, 32767)
  34
  35
  36 /**
  37  * Finalize MDCT and free allocated memory.
  38  */
  39 static av_cold void mdct_end(AC3MDCTContext *mdct)
  40 {
  41     mdct->nbits = 0;
  42     av_freep(&mdct->costab);
  43     av_freep(&mdct->sintab);
  44     av_freep(&mdct->xcos1);
  45     av_freep(&mdct->xsin1);
  46     av_freep(&mdct->rot_tmp);
  47     av_freep(&mdct->cplx_tmp);
  48 }
  49
  50
  51 /**
  52  * Initialize FFT tables.
  53  * @param ln log2(FFT size)
  54  */
  55 static av_cold int fft_init(AVCodecContext *avctx, AC3MDCTContext *mdct, int ln)
  56 {
  57     int i, n, n2;
  58     float alpha;
  59
  60     n  = 1 << ln;
  61     n2 = n >> 1;
  62
  63     FF_ALLOC_OR_GOTO(avctx, mdct->costab, n2 * sizeof(*mdct->costab), fft_alloc_fail);
  64     FF_ALLOC_OR_GOTO(avctx, mdct->sintab, n2 * sizeof(*mdct->sintab), fft_alloc_fail);
  65
  66     for (i = 0; i < n2; i++) {
  67         alpha     = 2.0 * M_PI * i / n;
  68         mdct->costab[i] = FIX15(cos(alpha));
  69         mdct->sintab[i] = FIX15(sin(alpha));
  70     }
  71
  72     return 0;
  73 fft_alloc_fail:
  74     mdct_end(mdct);
  75     return AVERROR(ENOMEM);
  76 }
  77
  78
  79 /**
  80  * Initialize MDCT tables.
  81  * @param nbits log2(MDCT size)
  82  */
  83 static av_cold int mdct_init(AVCodecContext *avctx, AC3MDCTContext *mdct,
  84                              int nbits)
  85 {
  86     int i, n, n4, ret;
  87
  88     n  = 1 << nbits;
  89     n4 = n >> 2;
  90
  91     mdct->nbits = nbits;
  92
  93     ret = fft_init(avctx, mdct, nbits - 2);
  94     if (ret)
  95         return ret;
  96
  97     mdct->window = ff_ac3_window;
  98
  99     FF_ALLOC_OR_GOTO(avctx, mdct->xcos1,    n4 * sizeof(*mdct->xcos1),    mdct_alloc_fail);
 100     FF_ALLOC_OR_GOTO(avctx, mdct->xsin1,    n4 * sizeof(*mdct->xsin1),    mdct_alloc_fail);
 101     FF_ALLOC_OR_GOTO(avctx, mdct->rot_tmp,  n  * sizeof(*mdct->rot_tmp),  mdct_alloc_fail);
 102     FF_ALLOC_OR_GOTO(avctx, mdct->cplx_tmp, n4 * sizeof(*mdct->cplx_tmp), mdct_alloc_fail);
 103
 104     for (i = 0; i < n4; i++) {
 105         float alpha = 2.0 * M_PI * (i + 1.0 / 8.0) / n;
 106         mdct->xcos1[i] = FIX15(-cos(alpha));
 107         mdct->xsin1[i] = FIX15(-sin(alpha));
 108     }
 109
 110     return 0;
 111 mdct_alloc_fail:
 112     mdct_end(mdct);
 113     return AVERROR(ENOMEM);
 114 }
 115
 116
 117 /** Butterfly op */
 118 #define BF(pre, pim, qre, qim, pre1, pim1, qre1, qim1)  \
 119 {                                                       \
 120   int ax, ay, bx, by;                                   \
 121   bx  = pre1;                                           \
 122   by  = pim1;                                           \
 123   ax  = qre1;                                           \
 124   ay  = qim1;                                           \
 125   pre = (bx + ax) >> 1;                                 \
 126   pim = (by + ay) >> 1;                                 \
 127   qre = (bx - ax) >> 1;                                 \
 128   qim = (by - ay) >> 1;                                 \
 129 }
 130
 131
 132 /** Complex multiply */
 133 #define CMUL(pre, pim, are, aim, bre, bim)              \
 134 {                                                       \
 135    pre = (MUL16(are, bre) - MUL16(aim, bim)) >> 15;     \
 136    pim = (MUL16(are, bim) + MUL16(bre, aim)) >> 15;     \
 137 }
 138
 139
 140 /**
 141  * Calculate a 2^n point complex FFT on 2^ln points.
 142  * @param z  complex input/output samples
 143  * @param ln log2(FFT size)
 144  */
 145 static void fft(AC3MDCTContext *mdct, IComplex *z, int ln)
 146 {
 147     int j, l, np, np2;
 148     int nblocks, nloops;
 149     register IComplex *p,*q;
 150     int tmp_re, tmp_im;
 151
 152     np = 1 << ln;
 153
 154     /* reverse */
 155     for (j = 0; j < np; j++) {
 156         int k = av_reverse[j] >> (8 - ln);
 157         if (k < j)
 158             FFSWAP(IComplex, z[k], z[j]);
 159     }
 160
 161     /* pass 0 */
 162
 163     p = &z[0];
 164     j = np >> 1;
 165     do {
 166         BF(p[0].re, p[0].im, p[1].re, p[1].im,
 167            p[0].re, p[0].im, p[1].re, p[1].im);
 168         p += 2;
 169     } while (--j);
 170
 171     /* pass 1 */
 172
 173     p = &z[0];
 174     j = np >> 2;
 175     do {
 176         BF(p[0].re, p[0].im, p[2].re,  p[2].im,
 177            p[0].re, p[0].im, p[2].re,  p[2].im);
 178         BF(p[1].re, p[1].im, p[3].re,  p[3].im,
 179            p[1].re, p[1].im, p[3].im, -p[3].re);
 180         p+=4;
 181     } while (--j);
 182
 183     /* pass 2 .. ln-1 */
 184
 185     nblocks = np >> 3;
 186     nloops  =  1 << 2;
 187     np2     = np >> 1;
 188     do {
 189         p = z;
 190         q = z + nloops;
 191         for (j = 0; j < nblocks; j++) {
 192             BF(p->re, p->im, q->re, q->im,
 193                p->re, p->im, q->re, q->im);
 194             p++;
 195             q++;
 196             for(l = nblocks; l < np2; l += nblocks) {
 197                 CMUL(tmp_re, tmp_im, mdct->costab[l], -mdct->sintab[l], q->re, q->im);
 198                 BF(p->re, p->im, q->re,  q->im,
 199                    p->re, p->im, tmp_re, tmp_im);
 200                 p++;
 201                 q++;
 202             }
 203             p += nloops;
 204             q += nloops;
 205         }
 206         nblocks = nblocks >> 1;
 207         nloops  = nloops  << 1;
 208     } while (nblocks);
 209 }
 210
 211
 212 /**
 213  * Calculate a 512-point MDCT
 214  * @param out 256 output frequency coefficients
 215  * @param in  512 windowed input audio samples
 216  */
 217 static void mdct512(AC3MDCTContext *mdct, int32_t *out, int16_t *in)
 218 {
 219     int i, re, im, n, n2, n4;
 220     int16_t *rot = mdct->rot_tmp;
 221     IComplex *x  = mdct->cplx_tmp;
 222
 223     n  = 1 << mdct->nbits;
 224     n2 = n >> 1;
 225     n4 = n >> 2;
 226
 227     /* shift to simplify computations */
 228     for (i = 0; i <n4; i++)
 229         rot[i] = -in[i + 3*n4];
 230     memcpy(&rot[n4], &in[0], 3*n4*sizeof(*in));
 231
 232     /* pre rotation */
 233     for (i = 0; i < n4; i++) {
 234         re =  ((int)rot[   2*i] - (int)rot[ n-1-2*i]) >> 1;
 235         im = -((int)rot[n2+2*i] - (int)rot[n2-1-2*i]) >> 1;
 236         CMUL(x[i].re, x[i].im, re, im, -mdct->xcos1[i], mdct->xsin1[i]);
 237     }
 238
 239     fft(mdct, x, mdct->nbits - 2);
 240
 241     /* post rotation */
 242     for (i = 0; i < n4; i++) {
 243         re = x[i].re;
 244         im = x[i].im;
 245         CMUL(out[n2-1-2*i], out[2*i], re, im, mdct->xsin1[i], mdct->xcos1[i]);
 246     }
 247 }
 248
 249
 250 /**
 251  * Apply KBD window to input samples prior to MDCT.
 252  */
 253 static void apply_window(int16_t *output, const int16_t *input,
 254                          const int16_t *window, int n)
 255 {
 256     int i;
 257     int n2 = n >> 1;
 258
 259     for (i = 0; i < n2; i++) {
 260         output[i]     = MUL16(input[i],     window[i]) >> 15;
 261         output[n-i-1] = MUL16(input[n-i-1], window[i]) >> 15;
 262     }
 263 }
 264
 265
 266 /**
 267  * Calculate the log2() of the maximum absolute value in an array.
 268  * @param tab input array
 269  * @param n   number of values in the array
 270  * @return    log2(max(abs(tab[])))
 271  */
 272 static int log2_tab(int16_t *tab, int n)
 273 {
 274     int i, v;
 275
 276     v = 0;
 277     for (i = 0; i < n; i++)
 278         v |= abs(tab[i]);
 279
 280     return av_log2(v);
 281 }
 282
 283
 284 /**
 285  * Left-shift each value in an array by a specified amount.
 286  * @param tab    input array
 287  * @param n      number of values in the array
 288  * @param lshift left shift amount. a negative value means right shift.
 289  */
 290 static void lshift_tab(int16_t *tab, int n, int lshift)
 291 {
 292     int i;
 293
 294     if (lshift > 0) {
 295         for (i = 0; i < n; i++)
 296             tab[i] <<= lshift;
 297     } else if (lshift < 0) {
 298         lshift = -lshift;
 299         for (i = 0; i < n; i++)
 300             tab[i] >>= lshift;
 301     }
 302 }
 303
 304
 305 /**
 306  * Normalize the input samples to use the maximum available precision.
 307  * This assumes signed 16-bit input samples. Exponents are reduced by 9 to
 308  * match the 24-bit internal precision for MDCT coefficients.
 309  *
 310  * @return exponent shift
 311  */
 312 static int normalize_samples(AC3EncodeContext *s)
 313 {
 314     int v = 14 - log2_tab(s->windowed_samples, AC3_WINDOW_SIZE);
 315     v = FFMAX(0, v);
 316     lshift_tab(s->windowed_samples, AC3_WINDOW_SIZE, v);
 317     return v - 9;
 318 }
 319
 320
 321 #ifdef TEST
 322 /*************************************************************************/
 323 /* TEST */
 324
 325 #include "libavutil/lfg.h"
 326
 327 #define MDCT_NBITS 9
 328 #define MDCT_SAMPLES (1 << MDCT_NBITS)
 329 #define FN (MDCT_SAMPLES/4)
 330
 331
 332 static void fft_test(AC3MDCTContext *mdct, AVLFG *lfg)
 333 {
 334     IComplex in[FN], in1[FN];
 335     int k, n, i;
 336     float sum_re, sum_im, a;
 337
 338     for (i = 0; i < FN; i++) {
 339         in[i].re = av_lfg_get(lfg) % 65535 - 32767;
 340         in[i].im = av_lfg_get(lfg) % 65535 - 32767;
 341         in1[i]   = in[i];
 342     }
 343     fft(mdct, in, 7);
 344
 345     /* do it by hand */
 346     for (k = 0; k < FN; k++) {
 347         sum_re = 0;
 348         sum_im = 0;
 349         for (n = 0; n < FN; n++) {
 350             a = -2 * M_PI * (n * k) / FN;
 351             sum_re += in1[n].re * cos(a) - in1[n].im * sin(a);
 352             sum_im += in1[n].re * sin(a) + in1[n].im * cos(a);
 353         }
 354         av_log(NULL, AV_LOG_DEBUG, "%3d: %6d,%6d %6.0f,%6.0f\n",
 355                k, in[k].re, in[k].im, sum_re / FN, sum_im / FN);
 356     }
 357 }
 358
 359
 360 static void mdct_test(AC3MDCTContext *mdct, AVLFG *lfg)
 361 {
 362     int16_t input[MDCT_SAMPLES];
 363     int32_t output[AC3_MAX_COEFS];
 364     float input1[MDCT_SAMPLES];
 365     float output1[AC3_MAX_COEFS];
 366     float s, a, err, e, emax;
 367     int i, k, n;
 368
 369     for (i = 0; i < MDCT_SAMPLES; i++) {
 370         input[i]  = (av_lfg_get(lfg) % 65535 - 32767) * 9 / 10;
 371         input1[i] = input[i];
 372     }
 373
 374     mdct512(mdct, output, input);
 375
 376     /* do it by hand */
 377     for (k = 0; k < AC3_MAX_COEFS; k++) {
 378         s = 0;
 379         for (n = 0; n < MDCT_SAMPLES; n++) {
 380             a = (2*M_PI*(2*n+1+MDCT_SAMPLES/2)*(2*k+1) / (4 * MDCT_SAMPLES));
 381             s += input1[n] * cos(a);
 382         }
 383         output1[k] = -2 * s / MDCT_SAMPLES;
 384     }
 385
 386     err  = 0;
 387     emax = 0;
 388     for (i = 0; i < AC3_MAX_COEFS; i++) {
 389         av_log(NULL, AV_LOG_DEBUG, "%3d: %7d %7.0f\n", i, output[i], output1[i]);
 390         e = output[i] - output1[i];
 391         if (e > emax)
 392             emax = e;
 393         err += e * e;
 394     }
 395     av_log(NULL, AV_LOG_DEBUG, "err2=%f emax=%f\n", err / AC3_MAX_COEFS, emax);
 396 }
 397
 398
 399 int main(void)
 400 {
 401     AVLFG lfg;
 402     AC3MDCTContext mdct;
 403
 404     mdct.avctx = NULL;
 405     av_log_set_level(AV_LOG_DEBUG);
 406     mdct_init(&mdct, 9);
 407
 408     fft_test(&mdct, &lfg);
 409     mdct_test(&mdct, &lfg);
 410
 411     return 0;
 412 }
 413 #endif /* TEST */
 414
 415
 416 AVCodec ac3_encoder = {
 417     "ac3",
 418     AVMEDIA_TYPE_AUDIO,
 419     CODEC_ID_AC3,
 420     sizeof(AC3EncodeContext),
 421     ac3_encode_init,
 422     ac3_encode_frame,
 423     ac3_encode_close,
 424     NULL,
 425     .sample_fmts = (const enum AVSampleFormat[]){AV_SAMPLE_FMT_S16,AV_SAMPLE_FMT_NONE},
 426     .long_name = NULL_IF_CONFIG_SMALL("ATSC A/52A (AC-3)"),
 427     .channel_layouts = ac3_channel_layouts,
 428 };