You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

1565 lines
43KB

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