git.sesse.net Git - ffmpeg/blob - libavcodec/ac3enc.c

   1 /*
   2  * The simplest AC3 encoder
   3  * Copyright (c) 2000 Fabrice Bellard.
   4  *
   5  * This library is free software; you can redistribute it and/or
   6  * modify it under the terms of the GNU Lesser General Public
   7  * License as published by the Free Software Foundation; either
   8  * version 2 of the License, or (at your option) any later version.
   9  *
  10  * This library is distributed in the hope that it will be useful,
  11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * Lesser General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU Lesser General Public
  16  * License along with this library; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  18  */
  19
  20 /**
  21  * @file ac3enc.c
  22  * The simplest AC3 encoder.
  23  */
  24 //#define DEBUG
  25 //#define DEBUG_BITALLOC
  26 #include "avcodec.h"
  27 #include "bitstream.h"
  28 #include "crc.h"
  29 #include "ac3.h"
  30
  31 typedef struct AC3EncodeContext {
  32     PutBitContext pb;
  33     int nb_channels;
  34     int nb_all_channels;
  35     int lfe_channel;
  36     int bit_rate;
  37     unsigned int sample_rate;
  38     unsigned int bsid;
  39     unsigned int frame_size_min; /* minimum frame size in case rounding is necessary */
  40     unsigned int frame_size; /* current frame size in words */
  41     unsigned int bits_written;
  42     unsigned int samples_written;
  43     int halfratecod;
  44     unsigned int frmsizecod;
  45     unsigned int fscod; /* frequency */
  46     unsigned int acmod;
  47     int lfe;
  48     unsigned int bsmod;
  49     short last_samples[AC3_MAX_CHANNELS][256];
  50     unsigned int chbwcod[AC3_MAX_CHANNELS];
  51     int nb_coefs[AC3_MAX_CHANNELS];
  52
  53     /* bitrate allocation control */
  54     int sgaincod, sdecaycod, fdecaycod, dbkneecod, floorcod;
  55     AC3BitAllocParameters bit_alloc;
  56     int csnroffst;
  57     int fgaincod[AC3_MAX_CHANNELS];
  58     int fsnroffst[AC3_MAX_CHANNELS];
  59     /* mantissa encoding */
  60     int mant1_cnt, mant2_cnt, mant4_cnt;
  61 } AC3EncodeContext;
  62
  63 #include "ac3tab.h"
  64
  65 #define MDCT_NBITS 9
  66 #define N         (1 << MDCT_NBITS)
  67
  68 /* new exponents are sent if their Norm 1 exceed this number */
  69 #define EXP_DIFF_THRESHOLD 1000
  70
  71 static void fft_init(int ln);
  72
  73 static inline int16_t fix15(float a)
  74 {
  75     int v;
  76     v = (int)(a * (float)(1 << 15));
  77     if (v < -32767)
  78         v = -32767;
  79     else if (v > 32767)
  80         v = 32767;
  81     return v;
  82 }
  83
  84
  85 /**
  86  * Generate a Kaiser Window.
  87  */
  88 static void k_window_init(int alpha, double *window, int n, int iter)
  89 {
  90     int j, k;
  91     double a, x;
  92     a = alpha * M_PI / n;
  93     a = a*a;
  94     for(k=0; k<n; k++) {
  95         x = k * (n - k) * a;
  96         window[k] = 1.0;
  97         for(j=iter; j>0; j--) {
  98             window[k] = (window[k] * x / (j*j)) + 1.0;
  99         }
 100     }
 101 }
 102
 103 /**
 104  * Generate a Kaiser-Bessel Derived Window.
 105  * @param alpha  determines window shape
 106  * @param window array to fill with window values
 107  * @param iter   number of iterations to use in BesselI0
 108  */
 109
 110 static void kbd_window_init(int alpha, int16_t *out_window, int iter)
 111 {
 112     int k, n2;
 113     double kwindow[256];
 114     double window[256];
 115
 116     n2 = 256;
 117     k_window_init(alpha, kwindow, n2, iter);
 118     window[0] = kwindow[0];
 119     for(k=1; k<n2; k++) {
 120         window[k] = window[k-1] + kwindow[k];
 121     }
 122     for(k=0; k<n2; k++) {
 123         window[k] = sqrt(window[k] / (window[n2-1]+1));
 124         //out_window[k] = round(window[k]*((1<<15)-1)); //enable this for a rounded window
 125         out_window[k] = fix15(window[k]);
 126     }
 127 }
 128
 129 static inline int calc_lowcomp1(int a, int b0, int b1)
 130 {
 131     if ((b0 + 256) == b1) {
 132         a = 384 ;
 133     } else if (b0 > b1) {
 134         a = a - 64;
 135         if (a < 0) a=0;
 136     }
 137     return a;
 138 }
 139
 140 static inline int calc_lowcomp(int a, int b0, int b1, int bin)
 141 {
 142     if (bin < 7) {
 143         if ((b0 + 256) == b1) {
 144             a = 384 ;
 145         } else if (b0 > b1) {
 146             a = a - 64;
 147             if (a < 0) a=0;
 148         }
 149     } else if (bin < 20) {
 150         if ((b0 + 256) == b1) {
 151             a = 320 ;
 152         } else if (b0 > b1) {
 153             a= a - 64;
 154             if (a < 0) a=0;
 155         }
 156     } else {
 157         a = a - 128;
 158         if (a < 0) a=0;
 159     }
 160     return a;
 161 }
 162
 163 /* AC3 bit allocation. The algorithm is the one described in the AC3
 164    spec. */
 165 void ac3_parametric_bit_allocation(AC3BitAllocParameters *s, uint8_t *bap,
 166                                    int8_t *exp, int start, int end,
 167                                    int snroffset, int fgain, int is_lfe,
 168                                    int deltbae,int deltnseg,
 169                                    uint8_t *deltoffst, uint8_t *deltlen, uint8_t *deltba)
 170 {
 171     int bin,i,j,k,end1,v,v1,bndstrt,bndend,lowcomp,begin;
 172     int fastleak,slowleak,address,tmp;
 173     int16_t psd[256]; /* scaled exponents */
 174     int16_t bndpsd[50]; /* interpolated exponents */
 175     int16_t excite[50]; /* excitation */
 176     int16_t mask[50];   /* masking value */
 177
 178     /* exponent mapping to PSD */
 179     for(bin=start;bin<end;bin++) {
 180         psd[bin]=(3072 - (exp[bin] << 7));
 181     }
 182
 183     /* PSD integration */
 184     j=start;
 185     k=masktab[start];
 186     do {
 187         v=psd[j];
 188         j++;
 189         end1=bndtab[k+1];
 190         if (end1 > end) end1=end;
 191         for(i=j;i<end1;i++) {
 192             int c,adr;
 193             /* logadd */
 194             v1=psd[j];
 195             c=v-v1;
 196             if (c >= 0) {
 197                 adr=c >> 1;
 198                 if (adr > 255) adr=255;
 199                 v=v + latab[adr];
 200             } else {
 201                 adr=(-c) >> 1;
 202                 if (adr > 255) adr=255;
 203                 v=v1 + latab[adr];
 204             }
 205             j++;
 206         }
 207         bndpsd[k]=v;
 208         k++;
 209     } while (end > bndtab[k]);
 210
 211     /* excitation function */
 212     bndstrt = masktab[start];
 213     bndend = masktab[end-1] + 1;
 214
 215     if (bndstrt == 0) {
 216         lowcomp = 0;
 217         lowcomp = calc_lowcomp1(lowcomp, bndpsd[0], bndpsd[1]) ;
 218         excite[0] = bndpsd[0] - fgain - lowcomp ;
 219         lowcomp = calc_lowcomp1(lowcomp, bndpsd[1], bndpsd[2]) ;
 220         excite[1] = bndpsd[1] - fgain - lowcomp ;
 221         begin = 7 ;
 222         for (bin = 2; bin < 7; bin++) {
 223             if (!(is_lfe && bin == 6))
 224                 lowcomp = calc_lowcomp1(lowcomp, bndpsd[bin], bndpsd[bin+1]) ;
 225             fastleak = bndpsd[bin] - fgain ;
 226             slowleak = bndpsd[bin] - s->sgain ;
 227             excite[bin] = fastleak - lowcomp ;
 228             if (!(is_lfe && bin == 6)) {
 229                 if (bndpsd[bin] <= bndpsd[bin+1]) {
 230                     begin = bin + 1 ;
 231                     break ;
 232                 }
 233             }
 234         }
 235
 236         end1=bndend;
 237         if (end1 > 22) end1=22;
 238
 239         for (bin = begin; bin < end1; bin++) {
 240             if (!(is_lfe && bin == 6))
 241                 lowcomp = calc_lowcomp(lowcomp, bndpsd[bin], bndpsd[bin+1], bin) ;
 242
 243             fastleak -= s->fdecay ;
 244             v = bndpsd[bin] - fgain;
 245             if (fastleak < v) fastleak = v;
 246
 247             slowleak -= s->sdecay ;
 248             v = bndpsd[bin] - s->sgain;
 249             if (slowleak < v) slowleak = v;
 250
 251             v=fastleak - lowcomp;
 252             if (slowleak > v) v=slowleak;
 253
 254             excite[bin] = v;
 255         }
 256         begin = 22;
 257     } else {
 258         /* coupling channel */
 259         begin = bndstrt;
 260
 261         fastleak = (s->cplfleak << 8) + 768;
 262         slowleak = (s->cplsleak << 8) + 768;
 263     }
 264
 265     for (bin = begin; bin < bndend; bin++) {
 266         fastleak -= s->fdecay ;
 267         v = bndpsd[bin] - fgain;
 268         if (fastleak < v) fastleak = v;
 269         slowleak -= s->sdecay ;
 270         v = bndpsd[bin] - s->sgain;
 271         if (slowleak < v) slowleak = v;
 272
 273         v=fastleak;
 274         if (slowleak > v) v = slowleak;
 275         excite[bin] = v;
 276     }
 277
 278     /* compute masking curve */
 279
 280     for (bin = bndstrt; bin < bndend; bin++) {
 281         v1 = excite[bin];
 282         tmp = s->dbknee - bndpsd[bin];
 283         if (tmp > 0) {
 284             v1 += tmp >> 2;
 285         }
 286         v=hth[bin >> s->halfratecod][s->fscod];
 287         if (v1 > v) v=v1;
 288         mask[bin] = v;
 289     }
 290
 291     /* delta bit allocation */
 292
 293     if (deltbae == 0 || deltbae == 1) {
 294         int band, seg, delta;
 295         band = 0 ;
 296         for (seg = 0; seg < deltnseg; seg++) {
 297             band += deltoffst[seg] ;
 298             if (deltba[seg] >= 4) {
 299                 delta = (deltba[seg] - 3) << 7;
 300             } else {
 301                 delta = (deltba[seg] - 4) << 7;
 302             }
 303             for (k = 0; k < deltlen[seg]; k++) {
 304                 mask[band] += delta ;
 305                 band++ ;
 306             }
 307         }
 308     }
 309
 310     /* compute bit allocation */
 311
 312     i = start ;
 313     j = masktab[start] ;
 314     do {
 315         v=mask[j];
 316         v -= snroffset ;
 317         v -= s->floor ;
 318         if (v < 0) v = 0;
 319         v &= 0x1fe0 ;
 320         v += s->floor ;
 321
 322         end1=bndtab[j] + bndsz[j];
 323         if (end1 > end) end1=end;
 324
 325         for (k = i; k < end1; k++) {
 326             address = (psd[i] - v) >> 5 ;
 327             if (address < 0) address=0;
 328             else if (address > 63) address=63;
 329             bap[i] = baptab[address];
 330             i++;
 331         }
 332     } while (end > bndtab[j++]) ;
 333 }
 334
 335 typedef struct IComplex {
 336     short re,im;
 337 } IComplex;
 338
 339 static void fft_init(int ln)
 340 {
 341     int i, j, m, n;
 342     float alpha;
 343
 344     n = 1 << ln;
 345
 346     for(i=0;i<(n/2);i++) {
 347         alpha = 2 * M_PI * (float)i / (float)n;
 348         costab[i] = fix15(cos(alpha));
 349         sintab[i] = fix15(sin(alpha));
 350     }
 351
 352     for(i=0;i<n;i++) {
 353         m=0;
 354         for(j=0;j<ln;j++) {
 355             m |= ((i >> j) & 1) << (ln-j-1);
 356         }
 357         fft_rev[i]=m;
 358     }
 359 }
 360
 361 /* butter fly op */
 362 #define BF(pre, pim, qre, qim, pre1, pim1, qre1, qim1) \
 363 {\
 364   int ax, ay, bx, by;\
 365   bx=pre1;\
 366   by=pim1;\
 367   ax=qre1;\
 368   ay=qim1;\
 369   pre = (bx + ax) >> 1;\
 370   pim = (by + ay) >> 1;\
 371   qre = (bx - ax) >> 1;\
 372   qim = (by - ay) >> 1;\
 373 }
 374
 375 #define MUL16(a,b) ((a) * (b))
 376
 377 #define CMUL(pre, pim, are, aim, bre, bim) \
 378 {\
 379    pre = (MUL16(are, bre) - MUL16(aim, bim)) >> 15;\
 380    pim = (MUL16(are, bim) + MUL16(bre, aim)) >> 15;\
 381 }
 382
 383
 384 /* do a 2^n point complex fft on 2^ln points. */
 385 static void fft(IComplex *z, int ln)
 386 {
 387     int        j, l, np, np2;
 388     int        nblocks, nloops;
 389     register IComplex *p,*q;
 390     int tmp_re, tmp_im;
 391
 392     np = 1 << ln;
 393
 394     /* reverse */
 395     for(j=0;j<np;j++) {
 396         int k;
 397         IComplex tmp;
 398         k = fft_rev[j];
 399         if (k < j) {
 400             tmp = z[k];
 401             z[k] = z[j];
 402             z[j] = tmp;
 403         }
 404     }
 405
 406     /* pass 0 */
 407
 408     p=&z[0];
 409     j=(np >> 1);
 410     do {
 411         BF(p[0].re, p[0].im, p[1].re, p[1].im,
 412            p[0].re, p[0].im, p[1].re, p[1].im);
 413         p+=2;
 414     } while (--j != 0);
 415
 416     /* pass 1 */
 417
 418     p=&z[0];
 419     j=np >> 2;
 420     do {
 421         BF(p[0].re, p[0].im, p[2].re, p[2].im,
 422            p[0].re, p[0].im, p[2].re, p[2].im);
 423         BF(p[1].re, p[1].im, p[3].re, p[3].im,
 424            p[1].re, p[1].im, p[3].im, -p[3].re);
 425         p+=4;
 426     } while (--j != 0);
 427
 428     /* pass 2 .. ln-1 */
 429
 430     nblocks = np >> 3;
 431     nloops = 1 << 2;
 432     np2 = np >> 1;
 433     do {
 434         p = z;
 435         q = z + nloops;
 436         for (j = 0; j < nblocks; ++j) {
 437
 438             BF(p->re, p->im, q->re, q->im,
 439                p->re, p->im, q->re, q->im);
 440
 441             p++;
 442             q++;
 443             for(l = nblocks; l < np2; l += nblocks) {
 444                 CMUL(tmp_re, tmp_im, costab[l], -sintab[l], q->re, q->im);
 445                 BF(p->re, p->im, q->re, q->im,
 446                    p->re, p->im, tmp_re, tmp_im);
 447                 p++;
 448                 q++;
 449             }
 450             p += nloops;
 451             q += nloops;
 452         }
 453         nblocks = nblocks >> 1;
 454         nloops = nloops << 1;
 455     } while (nblocks != 0);
 456 }
 457
 458 /* do a 512 point mdct */
 459 static void mdct512(int32_t *out, int16_t *in)
 460 {
 461     int i, re, im, re1, im1;
 462     int16_t rot[N];
 463     IComplex x[N/4];
 464
 465     /* shift to simplify computations */
 466     for(i=0;i<N/4;i++)
 467         rot[i] = -in[i + 3*N/4];
 468     for(i=N/4;i<N;i++)
 469         rot[i] = in[i - N/4];
 470
 471     /* pre rotation */
 472     for(i=0;i<N/4;i++) {
 473         re = ((int)rot[2*i] - (int)rot[N-1-2*i]) >> 1;
 474         im = -((int)rot[N/2+2*i] - (int)rot[N/2-1-2*i]) >> 1;
 475         CMUL(x[i].re, x[i].im, re, im, -xcos1[i], xsin1[i]);
 476     }
 477
 478     fft(x, MDCT_NBITS - 2);
 479
 480     /* post rotation */
 481     for(i=0;i<N/4;i++) {
 482         re = x[i].re;
 483         im = x[i].im;
 484         CMUL(re1, im1, re, im, xsin1[i], xcos1[i]);
 485         out[2*i] = im1;
 486         out[N/2-1-2*i] = re1;
 487     }
 488 }
 489
 490 /* XXX: use another norm ? */
 491 static int calc_exp_diff(uint8_t *exp1, uint8_t *exp2, int n)
 492 {
 493     int sum, i;
 494     sum = 0;
 495     for(i=0;i<n;i++) {
 496         sum += abs(exp1[i] - exp2[i]);
 497     }
 498     return sum;
 499 }
 500
 501 static void compute_exp_strategy(uint8_t exp_strategy[NB_BLOCKS][AC3_MAX_CHANNELS],
 502                                  uint8_t exp[NB_BLOCKS][AC3_MAX_CHANNELS][N/2],
 503                                  int ch, int is_lfe)
 504 {
 505     int i, j;
 506     int exp_diff;
 507
 508     /* estimate if the exponent variation & decide if they should be
 509        reused in the next frame */
 510     exp_strategy[0][ch] = EXP_NEW;
 511     for(i=1;i<NB_BLOCKS;i++) {
 512         exp_diff = calc_exp_diff(exp[i][ch], exp[i-1][ch], N/2);
 513 #ifdef DEBUG
 514         av_log(NULL, AV_LOG_DEBUG, "exp_diff=%d\n", exp_diff);
 515 #endif
 516         if (exp_diff > EXP_DIFF_THRESHOLD)
 517             exp_strategy[i][ch] = EXP_NEW;
 518         else
 519             exp_strategy[i][ch] = EXP_REUSE;
 520     }
 521     if (is_lfe)
 522         return;
 523
 524     /* now select the encoding strategy type : if exponents are often
 525        recoded, we use a coarse encoding */
 526     i = 0;
 527     while (i < NB_BLOCKS) {
 528         j = i + 1;
 529         while (j < NB_BLOCKS && exp_strategy[j][ch] == EXP_REUSE)
 530             j++;
 531         switch(j - i) {
 532         case 1:
 533             exp_strategy[i][ch] = EXP_D45;
 534             break;
 535         case 2:
 536         case 3:
 537             exp_strategy[i][ch] = EXP_D25;
 538             break;
 539         default:
 540             exp_strategy[i][ch] = EXP_D15;
 541             break;
 542         }
 543         i = j;
 544     }
 545 }
 546
 547 /* set exp[i] to min(exp[i], exp1[i]) */
 548 static void exponent_min(uint8_t exp[N/2], uint8_t exp1[N/2], int n)
 549 {
 550     int i;
 551
 552     for(i=0;i<n;i++) {
 553         if (exp1[i] < exp[i])
 554             exp[i] = exp1[i];
 555     }
 556 }
 557
 558 /* update the exponents so that they are the ones the decoder will
 559    decode. Return the number of bits used to code the exponents */
 560 static int encode_exp(uint8_t encoded_exp[N/2],
 561                       uint8_t exp[N/2],
 562                       int nb_exps,
 563                       int exp_strategy)
 564 {
 565     int group_size, nb_groups, i, j, k, exp_min;
 566     uint8_t exp1[N/2];
 567
 568     switch(exp_strategy) {
 569     case EXP_D15:
 570         group_size = 1;
 571         break;
 572     case EXP_D25:
 573         group_size = 2;
 574         break;
 575     default:
 576     case EXP_D45:
 577         group_size = 4;
 578         break;
 579     }
 580     nb_groups = ((nb_exps + (group_size * 3) - 4) / (3 * group_size)) * 3;
 581
 582     /* for each group, compute the minimum exponent */
 583     exp1[0] = exp[0]; /* DC exponent is handled separately */
 584     k = 1;
 585     for(i=1;i<=nb_groups;i++) {
 586         exp_min = exp[k];
 587         assert(exp_min >= 0 && exp_min <= 24);
 588         for(j=1;j<group_size;j++) {
 589             if (exp[k+j] < exp_min)
 590                 exp_min = exp[k+j];
 591         }
 592         exp1[i] = exp_min;
 593         k += group_size;
 594     }
 595
 596     /* constraint for DC exponent */
 597     if (exp1[0] > 15)
 598         exp1[0] = 15;
 599
 600     /* Decrease the delta between each groups to within 2
 601      * so that they can be differentially encoded */
 602     for (i=1;i<=nb_groups;i++)
 603         exp1[i] = FFMIN(exp1[i], exp1[i-1] + 2);
 604     for (i=nb_groups-1;i>=0;i--)
 605         exp1[i] = FFMIN(exp1[i], exp1[i+1] + 2);
 606
 607     /* now we have the exponent values the decoder will see */
 608     encoded_exp[0] = exp1[0];
 609     k = 1;
 610     for(i=1;i<=nb_groups;i++) {
 611         for(j=0;j<group_size;j++) {
 612             encoded_exp[k+j] = exp1[i];
 613         }
 614         k += group_size;
 615     }
 616
 617 #if defined(DEBUG)
 618     av_log(NULL, AV_LOG_DEBUG, "exponents: strategy=%d\n", exp_strategy);
 619     for(i=0;i<=nb_groups * group_size;i++) {
 620         av_log(NULL, AV_LOG_DEBUG, "%d ", encoded_exp[i]);
 621     }
 622     av_log(NULL, AV_LOG_DEBUG, "\n");
 623 #endif
 624
 625     return 4 + (nb_groups / 3) * 7;
 626 }
 627
 628 /* return the size in bits taken by the mantissa */
 629 static int compute_mantissa_size(AC3EncodeContext *s, uint8_t *m, int nb_coefs)
 630 {
 631     int bits, mant, i;
 632
 633     bits = 0;
 634     for(i=0;i<nb_coefs;i++) {
 635         mant = m[i];
 636         switch(mant) {
 637         case 0:
 638             /* nothing */
 639             break;
 640         case 1:
 641             /* 3 mantissa in 5 bits */
 642             if (s->mant1_cnt == 0)
 643                 bits += 5;
 644             if (++s->mant1_cnt == 3)
 645                 s->mant1_cnt = 0;
 646             break;
 647         case 2:
 648             /* 3 mantissa in 7 bits */
 649             if (s->mant2_cnt == 0)
 650                 bits += 7;
 651             if (++s->mant2_cnt == 3)
 652                 s->mant2_cnt = 0;
 653             break;
 654         case 3:
 655             bits += 3;
 656             break;
 657         case 4:
 658             /* 2 mantissa in 7 bits */
 659             if (s->mant4_cnt == 0)
 660                 bits += 7;
 661             if (++s->mant4_cnt == 2)
 662                 s->mant4_cnt = 0;
 663             break;
 664         case 14:
 665             bits += 14;
 666             break;
 667         case 15:
 668             bits += 16;
 669             break;
 670         default:
 671             bits += mant - 1;
 672             break;
 673         }
 674     }
 675     return bits;
 676 }
 677
 678
 679 static int bit_alloc(AC3EncodeContext *s,
 680                      uint8_t bap[NB_BLOCKS][AC3_MAX_CHANNELS][N/2],
 681                      uint8_t encoded_exp[NB_BLOCKS][AC3_MAX_CHANNELS][N/2],
 682                      uint8_t exp_strategy[NB_BLOCKS][AC3_MAX_CHANNELS],
 683                      int frame_bits, int csnroffst, int fsnroffst)
 684 {
 685     int i, ch;
 686
 687     /* compute size */
 688     for(i=0;i<NB_BLOCKS;i++) {
 689         s->mant1_cnt = 0;
 690         s->mant2_cnt = 0;
 691         s->mant4_cnt = 0;
 692         for(ch=0;ch<s->nb_all_channels;ch++) {
 693             ac3_parametric_bit_allocation(&s->bit_alloc,
 694                                           bap[i][ch], (int8_t *)encoded_exp[i][ch],
 695                                           0, s->nb_coefs[ch],
 696                                           (((csnroffst-15) << 4) +
 697                                            fsnroffst) << 2,
 698                                           fgaintab[s->fgaincod[ch]],
 699                                           ch == s->lfe_channel,
 700                                           2, 0, NULL, NULL, NULL);
 701             frame_bits += compute_mantissa_size(s, bap[i][ch],
 702                                                  s->nb_coefs[ch]);
 703         }
 704     }
 705 #if 0
 706     printf("csnr=%d fsnr=%d frame_bits=%d diff=%d\n",
 707            csnroffst, fsnroffst, frame_bits,
 708            16 * s->frame_size - ((frame_bits + 7) & ~7));
 709 #endif
 710     return 16 * s->frame_size - frame_bits;
 711 }
 712
 713 #define SNR_INC1 4
 714
 715 static int compute_bit_allocation(AC3EncodeContext *s,
 716                                   uint8_t bap[NB_BLOCKS][AC3_MAX_CHANNELS][N/2],
 717                                   uint8_t encoded_exp[NB_BLOCKS][AC3_MAX_CHANNELS][N/2],
 718                                   uint8_t exp_strategy[NB_BLOCKS][AC3_MAX_CHANNELS],
 719                                   int frame_bits)
 720 {
 721     int i, ch;
 722     int csnroffst, fsnroffst;
 723     uint8_t bap1[NB_BLOCKS][AC3_MAX_CHANNELS][N/2];
 724     static int frame_bits_inc[8] = { 0, 0, 2, 2, 2, 4, 2, 4 };
 725
 726     /* init default parameters */
 727     s->sdecaycod = 2;
 728     s->fdecaycod = 1;
 729     s->sgaincod = 1;
 730     s->dbkneecod = 2;
 731     s->floorcod = 4;
 732     for(ch=0;ch<s->nb_all_channels;ch++)
 733         s->fgaincod[ch] = 4;
 734
 735     /* compute real values */
 736     s->bit_alloc.fscod = s->fscod;
 737     s->bit_alloc.halfratecod = s->halfratecod;
 738     s->bit_alloc.sdecay = sdecaytab[s->sdecaycod] >> s->halfratecod;
 739     s->bit_alloc.fdecay = fdecaytab[s->fdecaycod] >> s->halfratecod;
 740     s->bit_alloc.sgain = sgaintab[s->sgaincod];
 741     s->bit_alloc.dbknee = dbkneetab[s->dbkneecod];
 742     s->bit_alloc.floor = floortab[s->floorcod];
 743
 744     /* header size */
 745     frame_bits += 65;
 746     // if (s->acmod == 2)
 747     //    frame_bits += 2;
 748     frame_bits += frame_bits_inc[s->acmod];
 749
 750     /* audio blocks */
 751     for(i=0;i<NB_BLOCKS;i++) {
 752         frame_bits += s->nb_channels * 2 + 2; /* blksw * c, dithflag * c, dynrnge, cplstre */
 753         if (s->acmod == 2) {
 754             frame_bits++; /* rematstr */
 755             if(i==0) frame_bits += 4;
 756         }
 757         frame_bits += 2 * s->nb_channels; /* chexpstr[2] * c */
 758         if (s->lfe)
 759             frame_bits++; /* lfeexpstr */
 760         for(ch=0;ch<s->nb_channels;ch++) {
 761             if (exp_strategy[i][ch] != EXP_REUSE)
 762                 frame_bits += 6 + 2; /* chbwcod[6], gainrng[2] */
 763         }
 764         frame_bits++; /* baie */
 765         frame_bits++; /* snr */
 766         frame_bits += 2; /* delta / skip */
 767     }
 768     frame_bits++; /* cplinu for block 0 */
 769     /* bit alloc info */
 770     /* sdcycod[2], fdcycod[2], sgaincod[2], dbpbcod[2], floorcod[3] */
 771     /* csnroffset[6] */
 772     /* (fsnoffset[4] + fgaincod[4]) * c */
 773     frame_bits += 2*4 + 3 + 6 + s->nb_all_channels * (4 + 3);
 774
 775     /* auxdatae, crcrsv */
 776     frame_bits += 2;
 777
 778     /* CRC */
 779     frame_bits += 16;
 780
 781     /* now the big work begins : do the bit allocation. Modify the snr
 782        offset until we can pack everything in the requested frame size */
 783
 784     csnroffst = s->csnroffst;
 785     while (csnroffst >= 0 &&
 786            bit_alloc(s, bap, encoded_exp, exp_strategy, frame_bits, csnroffst, 0) < 0)
 787         csnroffst -= SNR_INC1;
 788     if (csnroffst < 0) {
 789         av_log(NULL, AV_LOG_ERROR, "Bit allocation failed, try increasing the bitrate, -ab 384 for example!\n");
 790         return -1;
 791     }
 792     while ((csnroffst + SNR_INC1) <= 63 &&
 793            bit_alloc(s, bap1, encoded_exp, exp_strategy, frame_bits,
 794                      csnroffst + SNR_INC1, 0) >= 0) {
 795         csnroffst += SNR_INC1;
 796         memcpy(bap, bap1, sizeof(bap1));
 797     }
 798     while ((csnroffst + 1) <= 63 &&
 799            bit_alloc(s, bap1, encoded_exp, exp_strategy, frame_bits, csnroffst + 1, 0) >= 0) {
 800         csnroffst++;
 801         memcpy(bap, bap1, sizeof(bap1));
 802     }
 803
 804     fsnroffst = 0;
 805     while ((fsnroffst + SNR_INC1) <= 15 &&
 806            bit_alloc(s, bap1, encoded_exp, exp_strategy, frame_bits,
 807                      csnroffst, fsnroffst + SNR_INC1) >= 0) {
 808         fsnroffst += SNR_INC1;
 809         memcpy(bap, bap1, sizeof(bap1));
 810     }
 811     while ((fsnroffst + 1) <= 15 &&
 812            bit_alloc(s, bap1, encoded_exp, exp_strategy, frame_bits,
 813                      csnroffst, fsnroffst + 1) >= 0) {
 814         fsnroffst++;
 815         memcpy(bap, bap1, sizeof(bap1));
 816     }
 817
 818     s->csnroffst = csnroffst;
 819     for(ch=0;ch<s->nb_all_channels;ch++)
 820         s->fsnroffst[ch] = fsnroffst;
 821 #if defined(DEBUG_BITALLOC)
 822     {
 823         int j;
 824
 825         for(i=0;i<6;i++) {
 826             for(ch=0;ch<s->nb_all_channels;ch++) {
 827                 printf("Block #%d Ch%d:\n", i, ch);
 828                 printf("bap=");
 829                 for(j=0;j<s->nb_coefs[ch];j++) {
 830                     printf("%d ",bap[i][ch][j]);
 831                 }
 832                 printf("\n");
 833             }
 834         }
 835     }
 836 #endif
 837     return 0;
 838 }
 839
 840 void ac3_common_init(void)
 841 {
 842     int i, j, k, l, v;
 843     /* compute bndtab and masktab from bandsz */
 844     k = 0;
 845     l = 0;
 846     for(i=0;i<50;i++) {
 847         bndtab[i] = l;
 848         v = bndsz[i];
 849         for(j=0;j<v;j++) masktab[k++]=i;
 850         l += v;
 851     }
 852     bndtab[50] = 0;
 853 }
 854
 855
 856 static int AC3_encode_init(AVCodecContext *avctx)
 857 {
 858     int freq = avctx->sample_rate;
 859     int bitrate = avctx->bit_rate;
 860     int channels = avctx->channels;
 861     AC3EncodeContext *s = avctx->priv_data;
 862     int i, j, ch;
 863     float alpha;
 864     static const uint8_t acmod_defs[6] = {
 865         0x01, /* C */
 866         0x02, /* L R */
 867         0x03, /* L C R */
 868         0x06, /* L R SL SR */
 869         0x07, /* L C R SL SR */
 870         0x07, /* L C R SL SR (+LFE) */
 871     };
 872
 873     avctx->frame_size = AC3_FRAME_SIZE;
 874
 875     /* number of channels */
 876     if (channels < 1 || channels > 6)
 877         return -1;
 878     s->acmod = acmod_defs[channels - 1];
 879     s->lfe = (channels == 6) ? 1 : 0;
 880     s->nb_all_channels = channels;
 881     s->nb_channels = channels > 5 ? 5 : channels;
 882     s->lfe_channel = s->lfe ? 5 : -1;
 883
 884     /* frequency */
 885     for(i=0;i<3;i++) {
 886         for(j=0;j<3;j++)
 887             if ((ac3_freqs[j] >> i) == freq)
 888                 goto found;
 889     }
 890     return -1;
 891  found:
 892     s->sample_rate = freq;
 893     s->halfratecod = i;
 894     s->fscod = j;
 895     s->bsid = 8 + s->halfratecod;
 896     s->bsmod = 0; /* complete main audio service */
 897
 898     /* bitrate & frame size */
 899     bitrate /= 1000;
 900     for(i=0;i<19;i++) {
 901         if ((ac3_bitratetab[i] >> s->halfratecod) == bitrate)
 902             break;
 903     }
 904     if (i == 19)
 905         return -1;
 906     s->bit_rate = bitrate;
 907     s->frmsizecod = i << 1;
 908     s->frame_size_min = (bitrate * 1000 * AC3_FRAME_SIZE) / (freq * 16);
 909     s->bits_written = 0;
 910     s->samples_written = 0;
 911     s->frame_size = s->frame_size_min;
 912
 913     /* bit allocation init */
 914     for(ch=0;ch<s->nb_channels;ch++) {
 915         /* bandwidth for each channel */
 916         /* XXX: should compute the bandwidth according to the frame
 917            size, so that we avoid anoying high freq artefacts */
 918         s->chbwcod[ch] = 50; /* sample bandwidth as mpeg audio layer 2 table 0 */
 919         s->nb_coefs[ch] = ((s->chbwcod[ch] + 12) * 3) + 37;
 920     }
 921     if (s->lfe) {
 922         s->nb_coefs[s->lfe_channel] = 7; /* fixed */
 923     }
 924     /* initial snr offset */
 925     s->csnroffst = 40;
 926
 927     ac3_common_init();
 928
 929     /* mdct init */
 930     kbd_window_init(5.0, ac3_window, 50);
 931
 932     fft_init(MDCT_NBITS - 2);
 933     for(i=0;i<N/4;i++) {
 934         alpha = 2 * M_PI * (i + 1.0 / 8.0) / (float)N;
 935         xcos1[i] = fix15(-cos(alpha));
 936         xsin1[i] = fix15(-sin(alpha));
 937     }
 938
 939     avctx->coded_frame= avcodec_alloc_frame();
 940     avctx->coded_frame->key_frame= 1;
 941
 942     return 0;
 943 }
 944
 945 /* output the AC3 frame header */
 946 static void output_frame_header(AC3EncodeContext *s, unsigned char *frame)
 947 {
 948     init_put_bits(&s->pb, frame, AC3_MAX_CODED_FRAME_SIZE);
 949
 950     put_bits(&s->pb, 16, 0x0b77); /* frame header */
 951     put_bits(&s->pb, 16, 0); /* crc1: will be filled later */
 952     put_bits(&s->pb, 2, s->fscod);
 953     put_bits(&s->pb, 6, s->frmsizecod + (s->frame_size - s->frame_size_min));
 954     put_bits(&s->pb, 5, s->bsid);
 955     put_bits(&s->pb, 3, s->bsmod);
 956     put_bits(&s->pb, 3, s->acmod);
 957     if ((s->acmod & 0x01) && s->acmod != 0x01)
 958         put_bits(&s->pb, 2, 1); /* XXX -4.5 dB */
 959     if (s->acmod & 0x04)
 960         put_bits(&s->pb, 2, 1); /* XXX -6 dB */
 961     if (s->acmod == 0x02)
 962         put_bits(&s->pb, 2, 0); /* surround not indicated */
 963     put_bits(&s->pb, 1, s->lfe); /* LFE */
 964     put_bits(&s->pb, 5, 31); /* dialog norm: -31 db */
 965     put_bits(&s->pb, 1, 0); /* no compression control word */
 966     put_bits(&s->pb, 1, 0); /* no lang code */
 967     put_bits(&s->pb, 1, 0); /* no audio production info */
 968     put_bits(&s->pb, 1, 0); /* no copyright */
 969     put_bits(&s->pb, 1, 1); /* original bitstream */
 970     put_bits(&s->pb, 1, 0); /* no time code 1 */
 971     put_bits(&s->pb, 1, 0); /* no time code 2 */
 972     put_bits(&s->pb, 1, 0); /* no addtional bit stream info */
 973 }
 974
 975 /* symetric quantization on 'levels' levels */
 976 static inline int sym_quant(int c, int e, int levels)
 977 {
 978     int v;
 979
 980     if (c >= 0) {
 981         v = (levels * (c << e)) >> 24;
 982         v = (v + 1) >> 1;
 983         v = (levels >> 1) + v;
 984     } else {
 985         v = (levels * ((-c) << e)) >> 24;
 986         v = (v + 1) >> 1;
 987         v = (levels >> 1) - v;
 988     }
 989     assert (v >= 0 && v < levels);
 990     return v;
 991 }
 992
 993 /* asymetric quantization on 2^qbits levels */
 994 static inline int asym_quant(int c, int e, int qbits)
 995 {
 996     int lshift, m, v;
 997
 998     lshift = e + qbits - 24;
 999     if (lshift >= 0)
1000         v = c << lshift;
1001     else
1002         v = c >> (-lshift);
1003     /* rounding */
1004     v = (v + 1) >> 1;
1005     m = (1 << (qbits-1));
1006     if (v >= m)
1007         v = m - 1;
1008     assert(v >= -m);
1009     return v & ((1 << qbits)-1);
1010 }
1011
1012 /* Output one audio block. There are NB_BLOCKS audio blocks in one AC3
1013    frame */
1014 static void output_audio_block(AC3EncodeContext *s,
1015                                uint8_t exp_strategy[AC3_MAX_CHANNELS],
1016                                uint8_t encoded_exp[AC3_MAX_CHANNELS][N/2],
1017                                uint8_t bap[AC3_MAX_CHANNELS][N/2],
1018                                int32_t mdct_coefs[AC3_MAX_CHANNELS][N/2],
1019                                int8_t global_exp[AC3_MAX_CHANNELS],
1020                                int block_num)
1021 {
1022     int ch, nb_groups, group_size, i, baie, rbnd;
1023     uint8_t *p;
1024     uint16_t qmant[AC3_MAX_CHANNELS][N/2];
1025     int exp0, exp1;
1026     int mant1_cnt, mant2_cnt, mant4_cnt;
1027     uint16_t *qmant1_ptr, *qmant2_ptr, *qmant4_ptr;
1028     int delta0, delta1, delta2;
1029
1030     for(ch=0;ch<s->nb_channels;ch++)
1031         put_bits(&s->pb, 1, 0); /* 512 point MDCT */
1032     for(ch=0;ch<s->nb_channels;ch++)
1033         put_bits(&s->pb, 1, 1); /* no dither */
1034     put_bits(&s->pb, 1, 0); /* no dynamic range */
1035     if (block_num == 0) {
1036         /* for block 0, even if no coupling, we must say it. This is a
1037            waste of bit :-) */
1038         put_bits(&s->pb, 1, 1); /* coupling strategy present */
1039         put_bits(&s->pb, 1, 0); /* no coupling strategy */
1040     } else {
1041         put_bits(&s->pb, 1, 0); /* no new coupling strategy */
1042     }
1043
1044     if (s->acmod == 2)
1045       {
1046         if(block_num==0)
1047           {
1048             /* first block must define rematrixing (rematstr)  */
1049             put_bits(&s->pb, 1, 1);
1050
1051             /* dummy rematrixing rematflg(1:4)=0 */
1052             for (rbnd=0;rbnd<4;rbnd++)
1053               put_bits(&s->pb, 1, 0);
1054           }
1055         else
1056           {
1057             /* no matrixing (but should be used in the future) */
1058             put_bits(&s->pb, 1, 0);
1059           }
1060       }
1061
1062 #if defined(DEBUG)
1063     {
1064       static int count = 0;
1065       av_log(NULL, AV_LOG_DEBUG, "Block #%d (%d)\n", block_num, count++);
1066     }
1067 #endif
1068     /* exponent strategy */
1069     for(ch=0;ch<s->nb_channels;ch++) {
1070         put_bits(&s->pb, 2, exp_strategy[ch]);
1071     }
1072
1073     if (s->lfe) {
1074         put_bits(&s->pb, 1, exp_strategy[s->lfe_channel]);
1075     }
1076
1077     for(ch=0;ch<s->nb_channels;ch++) {
1078         if (exp_strategy[ch] != EXP_REUSE)
1079             put_bits(&s->pb, 6, s->chbwcod[ch]);
1080     }
1081
1082     /* exponents */
1083     for (ch = 0; ch < s->nb_all_channels; ch++) {
1084         switch(exp_strategy[ch]) {
1085         case EXP_REUSE:
1086             continue;
1087         case EXP_D15:
1088             group_size = 1;
1089             break;
1090         case EXP_D25:
1091             group_size = 2;
1092             break;
1093         default:
1094         case EXP_D45:
1095             group_size = 4;
1096             break;
1097         }
1098         nb_groups = (s->nb_coefs[ch] + (group_size * 3) - 4) / (3 * group_size);
1099         p = encoded_exp[ch];
1100
1101         /* first exponent */
1102         exp1 = *p++;
1103         put_bits(&s->pb, 4, exp1);
1104
1105         /* next ones are delta encoded */
1106         for(i=0;i<nb_groups;i++) {
1107             /* merge three delta in one code */
1108             exp0 = exp1;
1109             exp1 = p[0];
1110             p += group_size;
1111             delta0 = exp1 - exp0 + 2;
1112
1113             exp0 = exp1;
1114             exp1 = p[0];
1115             p += group_size;
1116             delta1 = exp1 - exp0 + 2;
1117
1118             exp0 = exp1;
1119             exp1 = p[0];
1120             p += group_size;
1121             delta2 = exp1 - exp0 + 2;
1122
1123             put_bits(&s->pb, 7, ((delta0 * 5 + delta1) * 5) + delta2);
1124         }
1125
1126         if (ch != s->lfe_channel)
1127             put_bits(&s->pb, 2, 0); /* no gain range info */
1128     }
1129
1130     /* bit allocation info */
1131     baie = (block_num == 0);
1132     put_bits(&s->pb, 1, baie);
1133     if (baie) {
1134         put_bits(&s->pb, 2, s->sdecaycod);
1135         put_bits(&s->pb, 2, s->fdecaycod);
1136         put_bits(&s->pb, 2, s->sgaincod);
1137         put_bits(&s->pb, 2, s->dbkneecod);
1138         put_bits(&s->pb, 3, s->floorcod);
1139     }
1140
1141     /* snr offset */
1142     put_bits(&s->pb, 1, baie); /* always present with bai */
1143     if (baie) {
1144         put_bits(&s->pb, 6, s->csnroffst);
1145         for(ch=0;ch<s->nb_all_channels;ch++) {
1146             put_bits(&s->pb, 4, s->fsnroffst[ch]);
1147             put_bits(&s->pb, 3, s->fgaincod[ch]);
1148         }
1149     }
1150
1151     put_bits(&s->pb, 1, 0); /* no delta bit allocation */
1152     put_bits(&s->pb, 1, 0); /* no data to skip */
1153
1154     /* mantissa encoding : we use two passes to handle the grouping. A
1155        one pass method may be faster, but it would necessitate to
1156        modify the output stream. */
1157
1158     /* first pass: quantize */
1159     mant1_cnt = mant2_cnt = mant4_cnt = 0;
1160     qmant1_ptr = qmant2_ptr = qmant4_ptr = NULL;
1161
1162     for (ch = 0; ch < s->nb_all_channels; ch++) {
1163         int b, c, e, v;
1164
1165         for(i=0;i<s->nb_coefs[ch];i++) {
1166             c = mdct_coefs[ch][i];
1167             e = encoded_exp[ch][i] - global_exp[ch];
1168             b = bap[ch][i];
1169             switch(b) {
1170             case 0:
1171                 v = 0;
1172                 break;
1173             case 1:
1174                 v = sym_quant(c, e, 3);
1175                 switch(mant1_cnt) {
1176                 case 0:
1177                     qmant1_ptr = &qmant[ch][i];
1178                     v = 9 * v;
1179                     mant1_cnt = 1;
1180                     break;
1181                 case 1:
1182                     *qmant1_ptr += 3 * v;
1183                     mant1_cnt = 2;
1184                     v = 128;
1185                     break;
1186                 default:
1187                     *qmant1_ptr += v;
1188                     mant1_cnt = 0;
1189                     v = 128;
1190                     break;
1191                 }
1192                 break;
1193             case 2:
1194                 v = sym_quant(c, e, 5);
1195                 switch(mant2_cnt) {
1196                 case 0:
1197                     qmant2_ptr = &qmant[ch][i];
1198                     v = 25 * v;
1199                     mant2_cnt = 1;
1200                     break;
1201                 case 1:
1202                     *qmant2_ptr += 5 * v;
1203                     mant2_cnt = 2;
1204                     v = 128;
1205                     break;
1206                 default:
1207                     *qmant2_ptr += v;
1208                     mant2_cnt = 0;
1209                     v = 128;
1210                     break;
1211                 }
1212                 break;
1213             case 3:
1214                 v = sym_quant(c, e, 7);
1215                 break;
1216             case 4:
1217                 v = sym_quant(c, e, 11);
1218                 switch(mant4_cnt) {
1219                 case 0:
1220                     qmant4_ptr = &qmant[ch][i];
1221                     v = 11 * v;
1222                     mant4_cnt = 1;
1223                     break;
1224                 default:
1225                     *qmant4_ptr += v;
1226                     mant4_cnt = 0;
1227                     v = 128;
1228                     break;
1229                 }
1230                 break;
1231             case 5:
1232                 v = sym_quant(c, e, 15);
1233                 break;
1234             case 14:
1235                 v = asym_quant(c, e, 14);
1236                 break;
1237             case 15:
1238                 v = asym_quant(c, e, 16);
1239                 break;
1240             default:
1241                 v = asym_quant(c, e, b - 1);
1242                 break;
1243             }
1244             qmant[ch][i] = v;
1245         }
1246     }
1247
1248     /* second pass : output the values */
1249     for (ch = 0; ch < s->nb_all_channels; ch++) {
1250         int b, q;
1251
1252         for(i=0;i<s->nb_coefs[ch];i++) {
1253             q = qmant[ch][i];
1254             b = bap[ch][i];
1255             switch(b) {
1256             case 0:
1257                 break;
1258             case 1:
1259                 if (q != 128)
1260                     put_bits(&s->pb, 5, q);
1261                 break;
1262             case 2:
1263                 if (q != 128)
1264                     put_bits(&s->pb, 7, q);
1265                 break;
1266             case 3:
1267                 put_bits(&s->pb, 3, q);
1268                 break;
1269             case 4:
1270                 if (q != 128)
1271                     put_bits(&s->pb, 7, q);
1272                 break;
1273             case 14:
1274                 put_bits(&s->pb, 14, q);
1275                 break;
1276             case 15:
1277                 put_bits(&s->pb, 16, q);
1278                 break;
1279             default:
1280                 put_bits(&s->pb, b - 1, q);
1281                 break;
1282             }
1283         }
1284     }
1285 }
1286
1287 #define CRC16_POLY ((1 << 0) | (1 << 2) | (1 << 15) | (1 << 16))
1288
1289 static unsigned int mul_poly(unsigned int a, unsigned int b, unsigned int poly)
1290 {
1291     unsigned int c;
1292
1293     c = 0;
1294     while (a) {
1295         if (a & 1)
1296             c ^= b;
1297         a = a >> 1;
1298         b = b << 1;
1299         if (b & (1 << 16))
1300             b ^= poly;
1301     }
1302     return c;
1303 }
1304
1305 static unsigned int pow_poly(unsigned int a, unsigned int n, unsigned int poly)
1306 {
1307     unsigned int r;
1308     r = 1;
1309     while (n) {
1310         if (n & 1)
1311             r = mul_poly(r, a, poly);
1312         a = mul_poly(a, a, poly);
1313         n >>= 1;
1314     }
1315     return r;
1316 }
1317
1318
1319 /* compute log2(max(abs(tab[]))) */
1320 static int log2_tab(int16_t *tab, int n)
1321 {
1322     int i, v;
1323
1324     v = 0;
1325     for(i=0;i<n;i++) {
1326         v |= abs(tab[i]);
1327     }
1328     return av_log2(v);
1329 }
1330
1331 static void lshift_tab(int16_t *tab, int n, int lshift)
1332 {
1333     int i;
1334
1335     if (lshift > 0) {
1336         for(i=0;i<n;i++) {
1337             tab[i] <<= lshift;
1338         }
1339     } else if (lshift < 0) {
1340         lshift = -lshift;
1341         for(i=0;i<n;i++) {
1342             tab[i] >>= lshift;
1343         }
1344     }
1345 }
1346
1347 /* fill the end of the frame and compute the two crcs */
1348 static int output_frame_end(AC3EncodeContext *s)
1349 {
1350     int frame_size, frame_size_58, n, crc1, crc2, crc_inv;
1351     uint8_t *frame;
1352
1353     frame_size = s->frame_size; /* frame size in words */
1354     /* align to 8 bits */
1355     flush_put_bits(&s->pb);
1356     /* add zero bytes to reach the frame size */
1357     frame = s->pb.buf;
1358     n = 2 * s->frame_size - (pbBufPtr(&s->pb) - frame) - 2;
1359     assert(n >= 0);
1360     if(n>0)
1361       memset(pbBufPtr(&s->pb), 0, n);
1362
1363     /* Now we must compute both crcs : this is not so easy for crc1
1364        because it is at the beginning of the data... */
1365     frame_size_58 = (frame_size >> 1) + (frame_size >> 3);
1366     crc1 = bswap_16(av_crc(av_crc8005, 0, frame + 4, 2 * frame_size_58 - 4));
1367     /* XXX: could precompute crc_inv */
1368     crc_inv = pow_poly((CRC16_POLY >> 1), (16 * frame_size_58) - 16, CRC16_POLY);
1369     crc1 = mul_poly(crc_inv, crc1, CRC16_POLY);
1370     frame[2] = crc1 >> 8;
1371     frame[3] = crc1;
1372
1373     crc2 = bswap_16(av_crc(av_crc8005, 0, frame + 2 * frame_size_58, (frame_size - frame_size_58) * 2 - 2));
1374     frame[2*frame_size - 2] = crc2 >> 8;
1375     frame[2*frame_size - 1] = crc2;
1376
1377     //    printf("n=%d frame_size=%d\n", n, frame_size);
1378     return frame_size * 2;
1379 }
1380
1381 static int AC3_encode_frame(AVCodecContext *avctx,
1382                             unsigned char *frame, int buf_size, void *data)
1383 {
1384     AC3EncodeContext *s = avctx->priv_data;
1385     int16_t *samples = data;
1386     int i, j, k, v, ch;
1387     int16_t input_samples[N];
1388     int32_t mdct_coef[NB_BLOCKS][AC3_MAX_CHANNELS][N/2];
1389     uint8_t exp[NB_BLOCKS][AC3_MAX_CHANNELS][N/2];
1390     uint8_t exp_strategy[NB_BLOCKS][AC3_MAX_CHANNELS];
1391     uint8_t encoded_exp[NB_BLOCKS][AC3_MAX_CHANNELS][N/2];
1392     uint8_t bap[NB_BLOCKS][AC3_MAX_CHANNELS][N/2];
1393     int8_t exp_samples[NB_BLOCKS][AC3_MAX_CHANNELS];
1394     int frame_bits;
1395
1396     frame_bits = 0;
1397     for(ch=0;ch<s->nb_all_channels;ch++) {
1398         /* fixed mdct to the six sub blocks & exponent computation */
1399         for(i=0;i<NB_BLOCKS;i++) {
1400             int16_t *sptr;
1401             int sinc;
1402
1403             /* compute input samples */
1404             memcpy(input_samples, s->last_samples[ch], N/2 * sizeof(int16_t));
1405             sinc = s->nb_all_channels;
1406             sptr = samples + (sinc * (N/2) * i) + ch;
1407             for(j=0;j<N/2;j++) {
1408                 v = *sptr;
1409                 input_samples[j + N/2] = v;
1410                 s->last_samples[ch][j] = v;
1411                 sptr += sinc;
1412             }
1413
1414             /* apply the MDCT window */
1415             for(j=0;j<N/2;j++) {
1416                 input_samples[j] = MUL16(input_samples[j],
1417                                          ac3_window[j]) >> 15;
1418                 input_samples[N-j-1] = MUL16(input_samples[N-j-1],
1419                                              ac3_window[j]) >> 15;
1420             }
1421
1422             /* Normalize the samples to use the maximum available
1423                precision */
1424             v = 14 - log2_tab(input_samples, N);
1425             if (v < 0)
1426                 v = 0;
1427             exp_samples[i][ch] = v - 9;
1428             lshift_tab(input_samples, N, v);
1429
1430             /* do the MDCT */
1431             mdct512(mdct_coef[i][ch], input_samples);
1432
1433             /* compute "exponents". We take into account the
1434                normalization there */
1435             for(j=0;j<N/2;j++) {
1436                 int e;
1437                 v = abs(mdct_coef[i][ch][j]);
1438                 if (v == 0)
1439                     e = 24;
1440                 else {
1441                     e = 23 - av_log2(v) + exp_samples[i][ch];
1442                     if (e >= 24) {
1443                         e = 24;
1444                         mdct_coef[i][ch][j] = 0;
1445                     }
1446                 }
1447                 exp[i][ch][j] = e;
1448             }
1449         }
1450
1451         compute_exp_strategy(exp_strategy, exp, ch, ch == s->lfe_channel);
1452
1453         /* compute the exponents as the decoder will see them. The
1454            EXP_REUSE case must be handled carefully : we select the
1455            min of the exponents */
1456         i = 0;
1457         while (i < NB_BLOCKS) {
1458             j = i + 1;
1459             while (j < NB_BLOCKS && exp_strategy[j][ch] == EXP_REUSE) {
1460                 exponent_min(exp[i][ch], exp[j][ch], s->nb_coefs[ch]);
1461                 j++;
1462             }
1463             frame_bits += encode_exp(encoded_exp[i][ch],
1464                                      exp[i][ch], s->nb_coefs[ch],
1465                                      exp_strategy[i][ch]);
1466             /* copy encoded exponents for reuse case */
1467             for(k=i+1;k<j;k++) {
1468                 memcpy(encoded_exp[k][ch], encoded_exp[i][ch],
1469                        s->nb_coefs[ch] * sizeof(uint8_t));
1470             }
1471             i = j;
1472         }
1473     }
1474
1475     /* adjust for fractional frame sizes */
1476     while(s->bits_written >= s->bit_rate*1000 && s->samples_written >= s->sample_rate) {
1477         s->bits_written -= s->bit_rate*1000;
1478         s->samples_written -= s->sample_rate;
1479     }
1480     s->frame_size = s->frame_size_min + (s->bits_written * s->sample_rate < s->samples_written * s->bit_rate*1000);
1481     s->bits_written += s->frame_size * 16;
1482     s->samples_written += AC3_FRAME_SIZE;
1483
1484     compute_bit_allocation(s, bap, encoded_exp, exp_strategy, frame_bits);
1485     /* everything is known... let's output the frame */
1486     output_frame_header(s, frame);
1487
1488     for(i=0;i<NB_BLOCKS;i++) {
1489         output_audio_block(s, exp_strategy[i], encoded_exp[i],
1490                            bap[i], mdct_coef[i], exp_samples[i], i);
1491     }
1492     return output_frame_end(s);
1493 }
1494
1495 static int AC3_encode_close(AVCodecContext *avctx)
1496 {
1497     av_freep(&avctx->coded_frame);
1498     return 0;
1499 }
1500
1501 #if 0
1502 /*************************************************************************/
1503 /* TEST */
1504
1505 #define FN (N/4)
1506
1507 void fft_test(void)
1508 {
1509     IComplex in[FN], in1[FN];
1510     int k, n, i;
1511     float sum_re, sum_im, a;
1512
1513     /* FFT test */
1514
1515     for(i=0;i<FN;i++) {
1516         in[i].re = random() % 65535 - 32767;
1517         in[i].im = random() % 65535 - 32767;
1518         in1[i] = in[i];
1519     }
1520     fft(in, 7);
1521
1522     /* do it by hand */
1523     for(k=0;k<FN;k++) {
1524         sum_re = 0;
1525         sum_im = 0;
1526         for(n=0;n<FN;n++) {
1527             a = -2 * M_PI * (n * k) / FN;
1528             sum_re += in1[n].re * cos(a) - in1[n].im * sin(a);
1529             sum_im += in1[n].re * sin(a) + in1[n].im * cos(a);
1530         }
1531         printf("%3d: %6d,%6d %6.0f,%6.0f\n",
1532                k, in[k].re, in[k].im, sum_re / FN, sum_im / FN);
1533     }
1534 }
1535
1536 void mdct_test(void)
1537 {
1538     int16_t input[N];
1539     int32_t output[N/2];
1540     float input1[N];
1541     float output1[N/2];
1542     float s, a, err, e, emax;
1543     int i, k, n;
1544
1545     for(i=0;i<N;i++) {
1546         input[i] = (random() % 65535 - 32767) * 9 / 10;
1547         input1[i] = input[i];
1548     }
1549
1550     mdct512(output, input);
1551
1552     /* do it by hand */
1553     for(k=0;k<N/2;k++) {
1554         s = 0;
1555         for(n=0;n<N;n++) {
1556             a = (2*M_PI*(2*n+1+N/2)*(2*k+1) / (4 * N));
1557             s += input1[n] * cos(a);
1558         }
1559         output1[k] = -2 * s / N;
1560     }
1561
1562     err = 0;
1563     emax = 0;
1564     for(i=0;i<N/2;i++) {
1565         printf("%3d: %7d %7.0f\n", i, output[i], output1[i]);
1566         e = output[i] - output1[i];
1567         if (e > emax)
1568             emax = e;
1569         err += e * e;
1570     }
1571     printf("err2=%f emax=%f\n", err / (N/2), emax);
1572 }
1573
1574 void test_ac3(void)
1575 {
1576     AC3EncodeContext ctx;
1577     unsigned char frame[AC3_MAX_CODED_FRAME_SIZE];
1578     short samples[AC3_FRAME_SIZE];
1579     int ret, i;
1580
1581     AC3_encode_init(&ctx, 44100, 64000, 1);
1582
1583     fft_test();
1584     mdct_test();
1585
1586     for(i=0;i<AC3_FRAME_SIZE;i++)
1587         samples[i] = (int)(sin(2*M_PI*i*1000.0/44100) * 10000);
1588     ret = AC3_encode_frame(&ctx, frame, samples);
1589     printf("ret=%d\n", ret);
1590 }
1591 #endif
1592
1593 AVCodec ac3_encoder = {
1594     "ac3",
1595     CODEC_TYPE_AUDIO,
1596     CODEC_ID_AC3,
1597     sizeof(AC3EncodeContext),
1598     AC3_encode_init,
1599     AC3_encode_frame,
1600     AC3_encode_close,
1601     NULL,
1602 };