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 #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, rshift)      \
 135 {                                                       \
 136    pre = (MUL16(are, bre) - MUL16(aim, bim)) >> rshift; \
 137    pim = (MUL16(are, bim) + MUL16(bre, aim)) >> rshift; \
 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, 15);
 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], 15);
 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], 0);
 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, unsigned int len)
 256 {
 257     dsp->apply_window_int16(output, input, window, len);
 258 }
 259
 260
 261 /**
 262  * Calculate the log2() of the maximum absolute value in an array.
 263  * @param tab input array
 264  * @param n   number of values in the array
 265  * @return    log2(max(abs(tab[])))
 266  */
 267 static int log2_tab(AC3EncodeContext *s, int16_t *src, int len)
 268 {
 269     int v = s->ac3dsp.ac3_max_msb_abs_int16(src, len);
 270     return av_log2(v);
 271 }
 272
 273
 274 /**
 275  * Normalize the input samples to use the maximum available precision.
 276  * This assumes signed 16-bit input samples.
 277  *
 278  * @return exponent shift
 279  */
 280 static int normalize_samples(AC3EncodeContext *s)
 281 {
 282     int v = 14 - log2_tab(s, s->windowed_samples, AC3_WINDOW_SIZE);
 283     if (v > 0)
 284         s->ac3dsp.ac3_lshift_int16(s->windowed_samples, AC3_WINDOW_SIZE, v);
 285     /* +6 to right-shift from 31-bit to 25-bit */
 286     return v + 6;
 287 }
 288
 289
 290 /**
 291  * Scale MDCT coefficients to 25-bit signed fixed-point.
 292  */
 293 static void scale_coefficients(AC3EncodeContext *s)
 294 {
 295     int blk, ch;
 296
 297     for (blk = 0; blk < AC3_MAX_BLOCKS; blk++) {
 298         AC3Block *block = &s->blocks[blk];
 299         for (ch = 0; ch < s->channels; ch++) {
 300             s->ac3dsp.ac3_rshift_int32(block->mdct_coef[ch], AC3_MAX_COEFS,
 301                                        block->coeff_shift[ch]);
 302         }
 303     }
 304 }
 305
 306
 307 #ifdef TEST
 308 /*************************************************************************/
 309 /* TEST */
 310
 311 #include "libavutil/lfg.h"
 312
 313 #define MDCT_NBITS 9
 314 #define MDCT_SAMPLES (1 << MDCT_NBITS)
 315 #define FN (MDCT_SAMPLES/4)
 316
 317
 318 static void fft_test(AC3MDCTContext *mdct, AVLFG *lfg)
 319 {
 320     IComplex in[FN], in1[FN];
 321     int k, n, i;
 322     float sum_re, sum_im, a;
 323
 324     for (i = 0; i < FN; i++) {
 325         in[i].re = av_lfg_get(lfg) % 65535 - 32767;
 326         in[i].im = av_lfg_get(lfg) % 65535 - 32767;
 327         in1[i]   = in[i];
 328     }
 329     fft(mdct, in, 7);
 330
 331     /* do it by hand */
 332     for (k = 0; k < FN; k++) {
 333         sum_re = 0;
 334         sum_im = 0;
 335         for (n = 0; n < FN; n++) {
 336             a = -2 * M_PI * (n * k) / FN;
 337             sum_re += in1[n].re * cos(a) - in1[n].im * sin(a);
 338             sum_im += in1[n].re * sin(a) + in1[n].im * cos(a);
 339         }
 340         av_log(NULL, AV_LOG_DEBUG, "%3d: %6d,%6d %6.0f,%6.0f\n",
 341                k, in[k].re, in[k].im, sum_re / FN, sum_im / FN);
 342     }
 343 }
 344
 345
 346 static void mdct_test(AC3MDCTContext *mdct, AVLFG *lfg)
 347 {
 348     int16_t input[MDCT_SAMPLES];
 349     int32_t output[AC3_MAX_COEFS];
 350     float input1[MDCT_SAMPLES];
 351     float output1[AC3_MAX_COEFS];
 352     float s, a, err, e, emax;
 353     int i, k, n;
 354
 355     for (i = 0; i < MDCT_SAMPLES; i++) {
 356         input[i]  = (av_lfg_get(lfg) % 65535 - 32767) * 9 / 10;
 357         input1[i] = input[i];
 358     }
 359
 360     mdct512(mdct, output, input);
 361
 362     /* do it by hand */
 363     for (k = 0; k < AC3_MAX_COEFS; k++) {
 364         s = 0;
 365         for (n = 0; n < MDCT_SAMPLES; n++) {
 366             a = (2*M_PI*(2*n+1+MDCT_SAMPLES/2)*(2*k+1) / (4 * MDCT_SAMPLES));
 367             s += input1[n] * cos(a);
 368         }
 369         output1[k] = -2 * s / MDCT_SAMPLES;
 370     }
 371
 372     err  = 0;
 373     emax = 0;
 374     for (i = 0; i < AC3_MAX_COEFS; i++) {
 375         av_log(NULL, AV_LOG_DEBUG, "%3d: %7d %7.0f\n", i, output[i], output1[i]);
 376         e = output[i] - output1[i];
 377         if (e > emax)
 378             emax = e;
 379         err += e * e;
 380     }
 381     av_log(NULL, AV_LOG_DEBUG, "err2=%f emax=%f\n", err / AC3_MAX_COEFS, emax);
 382 }
 383
 384
 385 int main(void)
 386 {
 387     AVLFG lfg;
 388     AC3MDCTContext mdct;
 389
 390     mdct.avctx = NULL;
 391     av_log_set_level(AV_LOG_DEBUG);
 392     mdct_init(&mdct, 9);
 393
 394     fft_test(&mdct, &lfg);
 395     mdct_test(&mdct, &lfg);
 396
 397     return 0;
 398 }
 399 #endif /* TEST */
 400
 401
 402 AVCodec ff_ac3_fixed_encoder = {
 403     "ac3_fixed",
 404     AVMEDIA_TYPE_AUDIO,
 405     CODEC_ID_AC3,
 406     sizeof(AC3EncodeContext),
 407     ac3_encode_init,
 408     ac3_encode_frame,
 409     ac3_encode_close,
 410     NULL,
 411     .sample_fmts = (const enum AVSampleFormat[]){AV_SAMPLE_FMT_S16,AV_SAMPLE_FMT_NONE},
 412     .long_name = NULL_IF_CONFIG_SMALL("ATSC A/52A (AC-3)"),
 413     .channel_layouts = ac3_channel_layouts,
 414 };