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