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