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.

1492 lines
40KB

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