开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 多媒体 > 流媒体/Mpeg4/MP4 > 查看源码
ac3enc.c
资源名称：ffmpeg-0.3.4.tar.gz　[点击查看]
上传用户：jxp0626
上传日期：2007-01-08
资源大小：102k
文件大小：39k
源码类别：

流媒体/Mpeg4/MP4

开发平台：

Unix_Linux

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