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