Mercurial > libavcodec.hg
view ac3enc.c @ 0:986e461dc072 libavcodec
Initial revision
| author | glantau |
|---|---|
| date | Sun, 22 Jul 2001 14:18:56 +0000 |
| parents | |
| children | 2f603eb2f23d |
line wrap: on
line source
/* * The simplest AC3 encoder * Copyright (c) 2000 Gerard Lantau. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */ #include <stdlib.h> #include <stdio.h> #include <string.h> #include <math.h> #include "avcodec.h" #include "ac3enc.h" #include "ac3tab.h" //#define DEBUG //#define DEBUG_BITALLOC #define NDEBUG #include <assert.h> #define MDCT_NBITS 9 #define N (1 << MDCT_NBITS) #define NB_BLOCKS 6 /* number of PCM blocks inside an AC3 frame */ /* new exponents are sent if their Norm 1 exceed this number */ #define EXP_DIFF_THRESHOLD 1000 /* exponent encoding strategy */ #define EXP_REUSE 0 #define EXP_NEW 1 #define EXP_D15 1 #define EXP_D25 2 #define EXP_D45 3 static void fft_init(int ln); static void ac3_crc_init(void); static inline INT16 fix15(float a) { int v; v = (int)(a * (float)(1 << 15)); if (v < -32767) v = -32767; else if (v > 32767) v = 32767; return v; } static inline int calc_lowcomp1(int a, int b0, int b1) { if ((b0 + 256) == b1) { a = 384 ; } else if (b0 > b1) { a = a - 64; if (a < 0) a=0; } return a; } static inline int calc_lowcomp(int a, int b0, int b1, int bin) { if (bin < 7) { if ((b0 + 256) == b1) { a = 384 ; } else if (b0 > b1) { a = a - 64; if (a < 0) a=0; } } else if (bin < 20) { if ((b0 + 256) == b1) { a = 320 ; } else if (b0 > b1) { a= a - 64; if (a < 0) a=0; } } else { a = a - 128; if (a < 0) a=0; } return a; } /* AC3 bit allocation. The algorithm is the one described in the AC3 spec with some optimizations because of our simplified encoding assumptions. */ void parametric_bit_allocation(AC3EncodeContext *s, UINT8 *bap, INT8 *exp, int start, int end, int snroffset, int fgain) { int bin,i,j,k,end1,v,v1,bndstrt,bndend,lowcomp,begin; int fastleak,slowleak,address,tmp; INT16 psd[256]; /* scaled exponents */ INT16 bndpsd[50]; /* interpolated exponents */ INT16 excite[50]; /* excitation */ INT16 mask[50]; /* masking value */ /* exponent mapping to PSD */ for(bin=start;bin<end;bin++) { psd[bin]=(3072 - (exp[bin] << 7)); } /* PSD integration */ j=start; k=masktab[start]; do { v=psd[j]; j++; end1=bndtab[k+1]; if (end1 > end) end1=end; for(i=j;i<end1;i++) { int c,adr; /* logadd */ v1=psd[j]; c=v-v1; if (c >= 0) { adr=c >> 1; if (adr > 255) adr=255; v=v + latab[adr]; } else { adr=(-c) >> 1; if (adr > 255) adr=255; v=v1 + latab[adr]; } j++; } bndpsd[k]=v; k++; } while (end > bndtab[k]); /* excitation function */ bndstrt = masktab[start]; bndend = masktab[end-1] + 1; lowcomp = 0; lowcomp = calc_lowcomp1(lowcomp, bndpsd[0], bndpsd[1]) ; excite[0] = bndpsd[0] - fgain - lowcomp ; lowcomp = calc_lowcomp1(lowcomp, bndpsd[1], bndpsd[2]) ; excite[1] = bndpsd[1] - fgain - lowcomp ; begin = 7 ; for (bin = 2; bin < 7; bin++) { lowcomp = calc_lowcomp1(lowcomp, bndpsd[bin], bndpsd[bin+1]) ; fastleak = bndpsd[bin] - fgain ; slowleak = bndpsd[bin] - s->sgain ; excite[bin] = fastleak - lowcomp ; if (bndpsd[bin] <= bndpsd[bin+1]) { begin = bin + 1 ; break ; } } end1=bndend; if (end1 > 22) end1=22; for (bin = begin; bin < end1; bin++) { lowcomp = calc_lowcomp(lowcomp, bndpsd[bin], bndpsd[bin+1], bin) ; fastleak -= s->fdecay ; v = bndpsd[bin] - fgain; if (fastleak < v) fastleak = v; slowleak -= s->sdecay ; v = bndpsd[bin] - s->sgain; if (slowleak < v) slowleak = v; v=fastleak - lowcomp; if (slowleak > v) v=slowleak; excite[bin] = v; } for (bin = 22; bin < bndend; bin++) { fastleak -= s->fdecay ; v = bndpsd[bin] - fgain; if (fastleak < v) fastleak = v; slowleak -= s->sdecay ; v = bndpsd[bin] - s->sgain; if (slowleak < v) slowleak = v; v=fastleak; if (slowleak > v) v = slowleak; excite[bin] = v; } /* compute masking curve */ for (bin = bndstrt; bin < bndend; bin++) { v1 = excite[bin]; tmp = s->dbknee - bndpsd[bin]; if (tmp > 0) { v1 += tmp >> 2; } v=hth[bin >> s->halfratecod][s->fscod]; if (v1 > v) v=v1; mask[bin] = v; } /* compute bit allocation */ i = start ; j = masktab[start] ; do { v=mask[j]; v -= snroffset ; v -= s->floor ; if (v < 0) v = 0; v &= 0x1fe0 ; v += s->floor ; end1=bndtab[j] + bndsz[j]; if (end1 > end) end1=end; for (k = i; k < end1; k++) { address = (psd[i] - v) >> 5 ; if (address < 0) address=0; else if (address > 63) address=63; bap[i] = baptab[address]; i++; } } while (end > bndtab[j++]) ; } typedef struct IComplex { short re,im; } IComplex; static void fft_init(int ln) { int i, j, m, n; float alpha; n = 1 << ln; for(i=0;i<(n/2);i++) { alpha = 2 * M_PI * (float)i / (float)n; costab[i] = fix15(cos(alpha)); sintab[i] = fix15(sin(alpha)); } for(i=0;i<n;i++) { m=0; for(j=0;j<ln;j++) { m |= ((i >> j) & 1) << (ln-j-1); } fft_rev[i]=m; } } /* butter fly op */ #define BF(pre, pim, qre, qim, pre1, pim1, qre1, qim1) \ {\ int ax, ay, bx, by;\ bx=pre1;\ by=pim1;\ ax=qre1;\ ay=qim1;\ pre = (bx + ax) >> 1;\ pim = (by + ay) >> 1;\ qre = (bx - ax) >> 1;\ qim = (by - ay) >> 1;\ } #define MUL16(a,b) ((a) * (b)) #define CMUL(pre, pim, are, aim, bre, bim) \ {\ pre = (MUL16(are, bre) - MUL16(aim, bim)) >> 15;\ pim = (MUL16(are, bim) + MUL16(bre, aim)) >> 15;\ } /* do a 2^n point complex fft on 2^ln points. */ static void fft(IComplex *z, int ln) { int j, l, np, np2; int nblocks, nloops; register IComplex *p,*q; int tmp_re, tmp_im; np = 1 << ln; /* reverse */ for(j=0;j<np;j++) { int k; IComplex tmp; k = fft_rev[j]; if (k < j) { tmp = z[k]; z[k] = z[j]; z[j] = tmp; } } /* pass 0 */ p=&z[0]; j=(np >> 1); do { BF(p[0].re, p[0].im, p[1].re, p[1].im, p[0].re, p[0].im, p[1].re, p[1].im); p+=2; } while (--j != 0); /* pass 1 */ p=&z[0]; j=np >> 2; do { BF(p[0].re, p[0].im, p[2].re, p[2].im, p[0].re, p[0].im, p[2].re, p[2].im); BF(p[1].re, p[1].im, p[3].re, p[3].im, p[1].re, p[1].im, p[3].im, -p[3].re); p+=4; } while (--j != 0); /* pass 2 .. ln-1 */ nblocks = np >> 3; nloops = 1 << 2; np2 = np >> 1; do { p = z; q = z + nloops; for (j = 0; j < nblocks; ++j) { BF(p->re, p->im, q->re, q->im, p->re, p->im, q->re, q->im); p++; q++; for(l = nblocks; l < np2; l += nblocks) { CMUL(tmp_re, tmp_im, costab[l], -sintab[l], q->re, q->im); BF(p->re, p->im, q->re, q->im, p->re, p->im, tmp_re, tmp_im); p++; q++; } p += nloops; q += nloops; } nblocks = nblocks >> 1; nloops = nloops << 1; } while (nblocks != 0); } /* do a 512 point mdct */ static void mdct512(INT32 *out, INT16 *in) { int i, re, im, re1, im1; INT16 rot[N]; IComplex x[N/4]; /* shift to simplify computations */ for(i=0;i<N/4;i++) rot[i] = -in[i + 3*N/4]; for(i=N/4;i<N;i++) rot[i] = in[i - N/4]; /* pre rotation */ for(i=0;i<N/4;i++) { re = ((int)rot[2*i] - (int)rot[N-1-2*i]) >> 1; im = -((int)rot[N/2+2*i] - (int)rot[N/2-1-2*i]) >> 1; CMUL(x[i].re, x[i].im, re, im, -xcos1[i], xsin1[i]); } fft(x, MDCT_NBITS - 2); /* post rotation */ for(i=0;i<N/4;i++) { re = x[i].re; im = x[i].im; CMUL(re1, im1, re, im, xsin1[i], xcos1[i]); out[2*i] = im1; out[N/2-1-2*i] = re1; } } /* XXX: use another norm ? */ static int calc_exp_diff(UINT8 *exp1, UINT8 *exp2, int n) { int sum, i; sum = 0; for(i=0;i<n;i++) { sum += abs(exp1[i] - exp2[i]); } return sum; } static void compute_exp_strategy(UINT8 exp_strategy[NB_BLOCKS][AC3_MAX_CHANNELS], UINT8 exp[NB_BLOCKS][AC3_MAX_CHANNELS][N/2], int ch) { int i, j; int exp_diff; /* estimate if the exponent variation & decide if they should be reused in the next frame */ exp_strategy[0][ch] = EXP_NEW; for(i=1;i<NB_BLOCKS;i++) { exp_diff = calc_exp_diff(exp[i][ch], exp[i-1][ch], N/2); #ifdef DEBUG printf("exp_diff=%d\n", exp_diff); #endif if (exp_diff > EXP_DIFF_THRESHOLD) exp_strategy[i][ch] = EXP_NEW; else exp_strategy[i][ch] = EXP_REUSE; } /* now select the encoding strategy type : if exponents are often recoded, we use a coarse encoding */ i = 0; while (i < NB_BLOCKS) { j = i + 1; while (j < NB_BLOCKS && exp_strategy[j][ch] == EXP_REUSE) j++; switch(j - i) { case 1: exp_strategy[i][ch] = EXP_D45; break; case 2: case 3: exp_strategy[i][ch] = EXP_D25; break; default: exp_strategy[i][ch] = EXP_D15; break; } i = j; } } /* set exp[i] to min(exp[i], exp1[i]) */ static void exponent_min(UINT8 exp[N/2], UINT8 exp1[N/2], int n) { int i; for(i=0;i<n;i++) { if (exp1[i] < exp[i]) exp[i] = exp1[i]; } } /* update the exponents so that they are the ones the decoder will decode. Return the number of bits used to code the exponents */ static int encode_exp(UINT8 encoded_exp[N/2], UINT8 exp[N/2], int nb_exps, int exp_strategy) { int group_size, nb_groups, i, j, k, recurse, exp_min, delta; UINT8 exp1[N/2]; switch(exp_strategy) { case EXP_D15: group_size = 1; break; case EXP_D25: group_size = 2; break; default: case EXP_D45: group_size = 4; break; } nb_groups = ((nb_exps + (group_size * 3) - 4) / (3 * group_size)) * 3; /* for each group, compute the minimum exponent */ exp1[0] = exp[0]; /* DC exponent is handled separately */ k = 1; for(i=1;i<=nb_groups;i++) { exp_min = exp[k]; assert(exp_min >= 0 && exp_min <= 24); for(j=1;j<group_size;j++) { if (exp[k+j] < exp_min) exp_min = exp[k+j]; } exp1[i] = exp_min; k += group_size; } /* constraint for DC exponent */ if (exp1[0] > 15) exp1[0] = 15; /* Iterate until the delta constraints between each groups are satisfyed. I'm sure it is possible to find a better algorithm, but I am lazy */ do { recurse = 0; for(i=1;i<=nb_groups;i++) { delta = exp1[i] - exp1[i-1]; if (delta > 2) { /* if delta too big, we encode a smaller exponent */ exp1[i] = exp1[i-1] + 2; } else if (delta < -2) { /* if delta is too small, we must decrease the previous exponent, which means we must recurse */ recurse = 1; exp1[i-1] = exp1[i] + 2; } } } while (recurse); /* now we have the exponent values the decoder will see */ encoded_exp[0] = exp1[0]; k = 1; for(i=1;i<=nb_groups;i++) { for(j=0;j<group_size;j++) { encoded_exp[k+j] = exp1[i]; } k += group_size; } #if defined(DEBUG) printf("exponents: strategy=%d\n", exp_strategy); for(i=0;i<=nb_groups * group_size;i++) { printf("%d ", encoded_exp[i]); } printf("\n"); #endif return 4 + (nb_groups / 3) * 7; } /* return the size in bits taken by the mantissa */ int compute_mantissa_size(AC3EncodeContext *s, UINT8 *m, int nb_coefs) { int bits, mant, i; bits = 0; for(i=0;i<nb_coefs;i++) { mant = m[i]; switch(mant) { case 0: /* nothing */ break; case 1: /* 3 mantissa in 5 bits */ if (s->mant1_cnt == 0) bits += 5; if (++s->mant1_cnt == 3) s->mant1_cnt = 0; break; case 2: /* 3 mantissa in 7 bits */ if (s->mant2_cnt == 0) bits += 7; if (++s->mant2_cnt == 3) s->mant2_cnt = 0; break; case 3: bits += 3; break; case 4: /* 2 mantissa in 7 bits */ if (s->mant4_cnt == 0) bits += 7; if (++s->mant4_cnt == 2) s->mant4_cnt = 0; break; case 14: bits += 14; break; case 15: bits += 16; break; default: bits += mant - 1; break; } } return bits; } static int bit_alloc(AC3EncodeContext *s, UINT8 bap[NB_BLOCKS][AC3_MAX_CHANNELS][N/2], UINT8 encoded_exp[NB_BLOCKS][AC3_MAX_CHANNELS][N/2], UINT8 exp_strategy[NB_BLOCKS][AC3_MAX_CHANNELS], int frame_bits, int csnroffst, int fsnroffst) { int i, ch; /* compute size */ for(i=0;i<NB_BLOCKS;i++) { s->mant1_cnt = 0; s->mant2_cnt = 0; s->mant4_cnt = 0; for(ch=0;ch<s->nb_channels;ch++) { parametric_bit_allocation(s, bap[i][ch], encoded_exp[i][ch], 0, s->nb_coefs[ch], (((csnroffst-15) << 4) + fsnroffst) << 2, fgaintab[s->fgaincod[ch]]); frame_bits += compute_mantissa_size(s, bap[i][ch], s->nb_coefs[ch]); } } #if 0 printf("csnr=%d fsnr=%d frame_bits=%d diff=%d\n", csnroffst, fsnroffst, frame_bits, 16 * s->frame_size - ((frame_bits + 7) & ~7)); #endif return 16 * s->frame_size - frame_bits; } #define SNR_INC1 4 static int compute_bit_allocation(AC3EncodeContext *s, UINT8 bap[NB_BLOCKS][AC3_MAX_CHANNELS][N/2], UINT8 encoded_exp[NB_BLOCKS][AC3_MAX_CHANNELS][N/2], UINT8 exp_strategy[NB_BLOCKS][AC3_MAX_CHANNELS], int frame_bits) { int i, ch; int csnroffst, fsnroffst; UINT8 bap1[NB_BLOCKS][AC3_MAX_CHANNELS][N/2]; /* init default parameters */ s->sdecaycod = 2; s->fdecaycod = 1; s->sgaincod = 1; s->dbkneecod = 2; s->floorcod = 4; for(ch=0;ch<s->nb_channels;ch++) s->fgaincod[ch] = 4; /* compute real values */ s->sdecay = sdecaytab[s->sdecaycod] >> s->halfratecod; s->fdecay = fdecaytab[s->fdecaycod] >> s->halfratecod; s->sgain = sgaintab[s->sgaincod]; s->dbknee = dbkneetab[s->dbkneecod]; s->floor = floortab[s->floorcod]; /* header size */ frame_bits += 65; if (s->acmod == 2) frame_bits += 2; /* audio blocks */ for(i=0;i<NB_BLOCKS;i++) { frame_bits += s->nb_channels * 2 + 2; if (s->acmod == 2) frame_bits++; frame_bits += 2 * s->nb_channels; for(ch=0;ch<s->nb_channels;ch++) { if (exp_strategy[i][ch] != EXP_REUSE) frame_bits += 6 + 2; } frame_bits++; /* baie */ frame_bits++; /* snr */ frame_bits += 2; /* delta / skip */ } frame_bits++; /* cplinu for block 0 */ /* bit alloc info */ frame_bits += 2*4 + 3 + 6 + s->nb_channels * (4 + 3); /* CRC */ frame_bits += 16; /* now the big work begins : do the bit allocation. Modify the snr offset until we can pack everything in the requested frame size */ csnroffst = s->csnroffst; while (csnroffst >= 0 && bit_alloc(s, bap, encoded_exp, exp_strategy, frame_bits, csnroffst, 0) < 0) csnroffst -= SNR_INC1; if (csnroffst < 0) { fprintf(stderr, "Error !!!\n"); return -1; } while ((csnroffst + SNR_INC1) <= 63 && bit_alloc(s, bap1, encoded_exp, exp_strategy, frame_bits, csnroffst + SNR_INC1, 0) >= 0) { csnroffst += SNR_INC1; memcpy(bap, bap1, sizeof(bap1)); } while ((csnroffst + 1) <= 63 && bit_alloc(s, bap1, encoded_exp, exp_strategy, frame_bits, csnroffst + 1, 0) >= 0) { csnroffst++; memcpy(bap, bap1, sizeof(bap1)); } fsnroffst = 0; while ((fsnroffst + SNR_INC1) <= 15 && bit_alloc(s, bap1, encoded_exp, exp_strategy, frame_bits, csnroffst, fsnroffst + SNR_INC1) >= 0) { fsnroffst += SNR_INC1; memcpy(bap, bap1, sizeof(bap1)); } while ((fsnroffst + 1) <= 15 && bit_alloc(s, bap1, encoded_exp, exp_strategy, frame_bits, csnroffst, fsnroffst + 1) >= 0) { fsnroffst++; memcpy(bap, bap1, sizeof(bap1)); } s->csnroffst = csnroffst; for(ch=0;ch<s->nb_channels;ch++) s->fsnroffst[ch] = fsnroffst; #if defined(DEBUG_BITALLOC) { int j; for(i=0;i<6;i++) { for(ch=0;ch<s->nb_channels;ch++) { printf("Block #%d Ch%d:\n", i, ch); printf("bap="); for(j=0;j<s->nb_coefs[ch];j++) { printf("%d ",bap[i][ch][j]); } printf("\n"); } } } #endif return 0; } static int AC3_encode_init(AVCodecContext *avctx) { int freq = avctx->sample_rate; int bitrate = avctx->bit_rate; int channels = avctx->channels; AC3EncodeContext *s = avctx->priv_data; int i, j, k, l, ch, v; float alpha; static unsigned short freqs[3] = { 48000, 44100, 32000 }; avctx->frame_size = AC3_FRAME_SIZE; avctx->key_frame = 1; /* always key frame */ /* number of channels */ if (channels == 1) s->acmod = 1; else if (channels == 2) s->acmod = 2; else return -1; s->nb_channels = channels; /* frequency */ for(i=0;i<3;i++) { for(j=0;j<3;j++) if ((freqs[j] >> i) == freq) goto found; } return -1; found: s->sample_rate = freq; s->halfratecod = i; s->fscod = j; s->bsid = 8 + s->halfratecod; s->bsmod = 0; /* complete main audio service */ /* bitrate & frame size */ bitrate /= 1000; for(i=0;i<19;i++) { if ((bitratetab[i] >> s->halfratecod) == bitrate) break; } if (i == 19) return -1; s->bit_rate = bitrate; s->frmsizecod = i << 1; s->frame_size_min = (bitrate * 1000 * AC3_FRAME_SIZE) / (freq * 16); /* for now we do not handle fractional sizes */ s->frame_size = s->frame_size_min; /* bit allocation init */ for(ch=0;ch<s->nb_channels;ch++) { /* bandwidth for each channel */ /* XXX: should compute the bandwidth according to the frame size, so that we avoid anoying high freq artefacts */ s->chbwcod[ch] = 50; /* sample bandwidth as mpeg audio layer 2 table 0 */ s->nb_coefs[ch] = ((s->chbwcod[ch] + 12) * 3) + 37; } /* initial snr offset */ s->csnroffst = 40; /* compute bndtab and masktab from bandsz */ k = 0; l = 0; for(i=0;i<50;i++) { bndtab[i] = l; v = bndsz[i]; for(j=0;j<v;j++) masktab[k++]=i; l += v; } bndtab[50] = 0; /* mdct init */ fft_init(MDCT_NBITS - 2); for(i=0;i<N/4;i++) { alpha = 2 * M_PI * (i + 1.0 / 8.0) / (float)N; xcos1[i] = fix15(-cos(alpha)); xsin1[i] = fix15(-sin(alpha)); } ac3_crc_init(); return 0; } /* output the AC3 frame header */ static void output_frame_header(AC3EncodeContext *s, unsigned char *frame) { init_put_bits(&s->pb, frame, AC3_MAX_CODED_FRAME_SIZE, NULL, NULL); put_bits(&s->pb, 16, 0x0b77); /* frame header */ put_bits(&s->pb, 16, 0); /* crc1: will be filled later */ put_bits(&s->pb, 2, s->fscod); put_bits(&s->pb, 6, s->frmsizecod + (s->frame_size - s->frame_size_min)); put_bits(&s->pb, 5, s->bsid); put_bits(&s->pb, 3, s->bsmod); put_bits(&s->pb, 3, s->acmod); if (s->acmod == 2) { put_bits(&s->pb, 2, 0); /* surround not indicated */ } put_bits(&s->pb, 1, 0); /* no LFE */ put_bits(&s->pb, 5, 31); /* dialog norm: -31 db */ put_bits(&s->pb, 1, 0); /* no compression control word */ put_bits(&s->pb, 1, 0); /* no lang code */ put_bits(&s->pb, 1, 0); /* no audio production info */ put_bits(&s->pb, 1, 0); /* no copyright */ put_bits(&s->pb, 1, 1); /* original bitstream */ put_bits(&s->pb, 1, 0); /* no time code 1 */ put_bits(&s->pb, 1, 0); /* no time code 2 */ put_bits(&s->pb, 1, 0); /* no addtional bit stream info */ } /* symetric quantization on 'levels' levels */ static inline int sym_quant(int c, int e, int levels) { int v; if (c >= 0) { v = (levels * (c << e)) >> 25; v = (levels >> 1) + v; } else { v = (levels * ((-c) << e)) >> 25; v = (levels >> 1) - v; } assert (v >= 0 && v < levels); return v; } /* asymetric quantization on 2^qbits levels */ static inline int asym_quant(int c, int e, int qbits) { int lshift, m, v; lshift = e + qbits - 24; if (lshift >= 0) v = c << lshift; else v = c >> (-lshift); /* rounding */ v = (v + 1) >> 1; m = (1 << (qbits-1)); if (v >= m) v = m - 1; assert(v >= -m); return v & ((1 << qbits)-1); } /* Output one audio block. There are NB_BLOCKS audio blocks in one AC3 frame */ static void output_audio_block(AC3EncodeContext *s, UINT8 exp_strategy[AC3_MAX_CHANNELS], UINT8 encoded_exp[AC3_MAX_CHANNELS][N/2], UINT8 bap[AC3_MAX_CHANNELS][N/2], INT32 mdct_coefs[AC3_MAX_CHANNELS][N/2], INT8 global_exp[AC3_MAX_CHANNELS], int block_num) { int ch, nb_groups, group_size, i, baie; UINT8 *p; UINT16 qmant[AC3_MAX_CHANNELS][N/2]; int exp0, exp1; int mant1_cnt, mant2_cnt, mant4_cnt; UINT16 *qmant1_ptr, *qmant2_ptr, *qmant4_ptr; int delta0, delta1, delta2; for(ch=0;ch<s->nb_channels;ch++) put_bits(&s->pb, 1, 0); /* 512 point MDCT */ for(ch=0;ch<s->nb_channels;ch++) put_bits(&s->pb, 1, 1); /* no dither */ put_bits(&s->pb, 1, 0); /* no dynamic range */ if (block_num == 0) { /* for block 0, even if no coupling, we must say it. This is a waste of bit :-) */ put_bits(&s->pb, 1, 1); /* coupling strategy present */ put_bits(&s->pb, 1, 0); /* no coupling strategy */ } else { put_bits(&s->pb, 1, 0); /* no new coupling strategy */ } if (s->acmod == 2) { put_bits(&s->pb, 1, 0); /* no matrixing (but should be used in the future) */ } #if defined(DEBUG) { static int count = 0; printf("Block #%d (%d)\n", block_num, count++); } #endif /* exponent strategy */ for(ch=0;ch<s->nb_channels;ch++) { put_bits(&s->pb, 2, exp_strategy[ch]); } for(ch=0;ch<s->nb_channels;ch++) { if (exp_strategy[ch] != EXP_REUSE) put_bits(&s->pb, 6, s->chbwcod[ch]); } /* exponents */ for (ch = 0; ch < s->nb_channels; ch++) { switch(exp_strategy[ch]) { case EXP_REUSE: continue; case EXP_D15: group_size = 1; break; case EXP_D25: group_size = 2; break; default: case EXP_D45: group_size = 4; break; } nb_groups = (s->nb_coefs[ch] + (group_size * 3) - 4) / (3 * group_size); p = encoded_exp[ch]; /* first exponent */ exp1 = *p++; put_bits(&s->pb, 4, exp1); /* next ones are delta encoded */ for(i=0;i<nb_groups;i++) { /* merge three delta in one code */ exp0 = exp1; exp1 = p[0]; p += group_size; delta0 = exp1 - exp0 + 2; exp0 = exp1; exp1 = p[0]; p += group_size; delta1 = exp1 - exp0 + 2; exp0 = exp1; exp1 = p[0]; p += group_size; delta2 = exp1 - exp0 + 2; put_bits(&s->pb, 7, ((delta0 * 5 + delta1) * 5) + delta2); } put_bits(&s->pb, 2, 0); /* no gain range info */ } /* bit allocation info */ baie = (block_num == 0); put_bits(&s->pb, 1, baie); if (baie) { put_bits(&s->pb, 2, s->sdecaycod); put_bits(&s->pb, 2, s->fdecaycod); put_bits(&s->pb, 2, s->sgaincod); put_bits(&s->pb, 2, s->dbkneecod); put_bits(&s->pb, 3, s->floorcod); } /* snr offset */ put_bits(&s->pb, 1, baie); /* always present with bai */ if (baie) { put_bits(&s->pb, 6, s->csnroffst); for(ch=0;ch<s->nb_channels;ch++) { put_bits(&s->pb, 4, s->fsnroffst[ch]); put_bits(&s->pb, 3, s->fgaincod[ch]); } } put_bits(&s->pb, 1, 0); /* no delta bit allocation */ put_bits(&s->pb, 1, 0); /* no data to skip */ /* mantissa encoding : we use two passes to handle the grouping. A one pass method may be faster, but it would necessitate to modify the output stream. */ /* first pass: quantize */ mant1_cnt = mant2_cnt = mant4_cnt = 0; qmant1_ptr = qmant2_ptr = qmant4_ptr = NULL; for (ch = 0; ch < s->nb_channels; ch++) { int b, c, e, v; for(i=0;i<s->nb_coefs[ch];i++) { c = mdct_coefs[ch][i]; e = encoded_exp[ch][i] - global_exp[ch]; b = bap[ch][i]; switch(b) { case 0: v = 0; break; case 1: v = sym_quant(c, e, 3); switch(mant1_cnt) { case 0: qmant1_ptr = &qmant[ch][i]; v = 9 * v; mant1_cnt = 1; break; case 1: *qmant1_ptr += 3 * v; mant1_cnt = 2; v = 128; break; default: *qmant1_ptr += v; mant1_cnt = 0; v = 128; break; } break; case 2: v = sym_quant(c, e, 5); switch(mant2_cnt) { case 0: qmant2_ptr = &qmant[ch][i]; v = 25 * v; mant2_cnt = 1; break; case 1: *qmant2_ptr += 5 * v; mant2_cnt = 2; v = 128; break; default: *qmant2_ptr += v; mant2_cnt = 0; v = 128; break; } break; case 3: v = sym_quant(c, e, 7); break; case 4: v = sym_quant(c, e, 11); switch(mant4_cnt) { case 0: qmant4_ptr = &qmant[ch][i]; v = 11 * v; mant4_cnt = 1; break; default: *qmant4_ptr += v; mant4_cnt = 0; v = 128; break; } break; case 5: v = sym_quant(c, e, 15); break; case 14: v = asym_quant(c, e, 14); break; case 15: v = asym_quant(c, e, 16); break; default: v = asym_quant(c, e, b - 1); break; } qmant[ch][i] = v; } } /* second pass : output the values */ for (ch = 0; ch < s->nb_channels; ch++) { int b, q; for(i=0;i<s->nb_coefs[ch];i++) { q = qmant[ch][i]; b = bap[ch][i]; switch(b) { case 0: break; case 1: if (q != 128) put_bits(&s->pb, 5, q); break; case 2: if (q != 128) put_bits(&s->pb, 7, q); break; case 3: put_bits(&s->pb, 3, q); break; case 4: if (q != 128) put_bits(&s->pb, 7, q); break; case 14: put_bits(&s->pb, 14, q); break; case 15: put_bits(&s->pb, 16, q); break; default: put_bits(&s->pb, b - 1, q); break; } } } } /* compute the ac3 crc */ #define CRC16_POLY ((1 << 0) | (1 << 2) | (1 << 15) | (1 << 16)) static void ac3_crc_init(void) { unsigned int c, n, k; for(n=0;n<256;n++) { c = n << 8; for (k = 0; k < 8; k++) { if (c & (1 << 15)) c = ((c << 1) & 0xffff) ^ (CRC16_POLY & 0xffff); else c = c << 1; } crc_table[n] = c; } } static unsigned int ac3_crc(UINT8 *data, int n, unsigned int crc) { int i; for(i=0;i<n;i++) { crc = (crc_table[data[i] ^ (crc >> 8)] ^ (crc << 8)) & 0xffff; } return crc; } static unsigned int mul_poly(unsigned int a, unsigned int b, unsigned int poly) { unsigned int c; c = 0; while (a) { if (a & 1) c ^= b; a = a >> 1; b = b << 1; if (b & (1 << 16)) b ^= poly; } return c; } static unsigned int pow_poly(unsigned int a, unsigned int n, unsigned int poly) { unsigned int r; r = 1; while (n) { if (n & 1) r = mul_poly(r, a, poly); a = mul_poly(a, a, poly); n >>= 1; } return r; } /* compute log2(max(abs(tab[]))) */ static int log2_tab(INT16 *tab, int n) { int i, v; v = 0; for(i=0;i<n;i++) { v |= abs(tab[i]); } return log2(v); } static void lshift_tab(INT16 *tab, int n, int lshift) { int i; if (lshift > 0) { for(i=0;i<n;i++) { tab[i] <<= lshift; } } else if (lshift < 0) { lshift = -lshift; for(i=0;i<n;i++) { tab[i] >>= lshift; } } } /* fill the end of the frame and compute the two crcs */ static int output_frame_end(AC3EncodeContext *s) { int frame_size, frame_size_58, n, crc1, crc2, crc_inv; UINT8 *frame; frame_size = s->frame_size; /* frame size in words */ /* align to 8 bits */ flush_put_bits(&s->pb); /* add zero bytes to reach the frame size */ frame = s->pb.buf; n = 2 * s->frame_size - (s->pb.buf_ptr - frame) - 2; assert(n >= 0); memset(s->pb.buf_ptr, 0, n); /* Now we must compute both crcs : this is not so easy for crc1 because it is at the beginning of the data... */ frame_size_58 = (frame_size >> 1) + (frame_size >> 3); crc1 = ac3_crc(frame + 4, (2 * frame_size_58) - 4, 0); /* XXX: could precompute crc_inv */ crc_inv = pow_poly((CRC16_POLY >> 1), (16 * frame_size_58) - 16, CRC16_POLY); crc1 = mul_poly(crc_inv, crc1, CRC16_POLY); frame[2] = crc1 >> 8; frame[3] = crc1; crc2 = ac3_crc(frame + 2 * frame_size_58, (frame_size - frame_size_58) * 2 - 2, 0); frame[2*frame_size - 2] = crc2 >> 8; frame[2*frame_size - 1] = crc2; // printf("n=%d frame_size=%d\n", n, frame_size); return frame_size * 2; } int AC3_encode_frame(AVCodecContext *avctx, unsigned char *frame, int buf_size, void *data) { AC3EncodeContext *s = avctx->priv_data; short *samples = data; int i, j, k, v, ch; INT16 input_samples[N]; INT32 mdct_coef[NB_BLOCKS][AC3_MAX_CHANNELS][N/2]; UINT8 exp[NB_BLOCKS][AC3_MAX_CHANNELS][N/2]; UINT8 exp_strategy[NB_BLOCKS][AC3_MAX_CHANNELS]; UINT8 encoded_exp[NB_BLOCKS][AC3_MAX_CHANNELS][N/2]; UINT8 bap[NB_BLOCKS][AC3_MAX_CHANNELS][N/2]; INT8 exp_samples[NB_BLOCKS][AC3_MAX_CHANNELS]; int frame_bits; frame_bits = 0; for(ch=0;ch<s->nb_channels;ch++) { /* fixed mdct to the six sub blocks & exponent computation */ for(i=0;i<NB_BLOCKS;i++) { INT16 *sptr; int sinc; /* compute input samples */ memcpy(input_samples, s->last_samples[ch], N/2 * sizeof(INT16)); sinc = s->nb_channels; sptr = samples + (sinc * (N/2) * i) + ch; for(j=0;j<N/2;j++) { v = *sptr; input_samples[j + N/2] = v; s->last_samples[ch][j] = v; sptr += sinc; } /* apply the MDCT window */ for(j=0;j<N/2;j++) { input_samples[j] = MUL16(input_samples[j], ac3_window[j]) >> 15; input_samples[N-j-1] = MUL16(input_samples[N-j-1], ac3_window[j]) >> 15; } /* Normalize the samples to use the maximum available precision */ v = 14 - log2_tab(input_samples, N); if (v < 0) v = 0; exp_samples[i][ch] = v - 8; lshift_tab(input_samples, N, v); /* do the MDCT */ mdct512(mdct_coef[i][ch], input_samples); /* compute "exponents". We take into account the normalization there */ for(j=0;j<N/2;j++) { int e; v = abs(mdct_coef[i][ch][j]); if (v == 0) e = 24; else { e = 23 - log2(v) + exp_samples[i][ch]; if (e >= 24) { e = 24; mdct_coef[i][ch][j] = 0; } } exp[i][ch][j] = e; } } compute_exp_strategy(exp_strategy, exp, ch); /* compute the exponents as the decoder will see them. The EXP_REUSE case must be handled carefully : we select the min of the exponents */ i = 0; while (i < NB_BLOCKS) { j = i + 1; while (j < NB_BLOCKS && exp_strategy[j][ch] == EXP_REUSE) { exponent_min(exp[i][ch], exp[j][ch], s->nb_coefs[ch]); j++; } frame_bits += encode_exp(encoded_exp[i][ch], exp[i][ch], s->nb_coefs[ch], exp_strategy[i][ch]); /* copy encoded exponents for reuse case */ for(k=i+1;k<j;k++) { memcpy(encoded_exp[k][ch], encoded_exp[i][ch], s->nb_coefs[ch] * sizeof(UINT8)); } i = j; } } compute_bit_allocation(s, bap, encoded_exp, exp_strategy, frame_bits); /* everything is known... let's output the frame */ output_frame_header(s, frame); for(i=0;i<NB_BLOCKS;i++) { output_audio_block(s, exp_strategy[i], encoded_exp[i], bap[i], mdct_coef[i], exp_samples[i], i); } return output_frame_end(s); } #if 0 /*************************************************************************/ /* TEST */ #define FN (N/4) void fft_test(void) { IComplex in[FN], in1[FN]; int k, n, i; float sum_re, sum_im, a; /* FFT test */ for(i=0;i<FN;i++) { in[i].re = random() % 65535 - 32767; in[i].im = random() % 65535 - 32767; in1[i] = in[i]; } fft(in, 7); /* do it by hand */ for(k=0;k<FN;k++) { sum_re = 0; sum_im = 0; for(n=0;n<FN;n++) { a = -2 * M_PI * (n * k) / FN; sum_re += in1[n].re * cos(a) - in1[n].im * sin(a); sum_im += in1[n].re * sin(a) + in1[n].im * cos(a); } printf("%3d: %6d,%6d %6.0f,%6.0f\n", k, in[k].re, in[k].im, sum_re / FN, sum_im / FN); } } void mdct_test(void) { INT16 input[N]; INT32 output[N/2]; float input1[N]; float output1[N/2]; float s, a, err, e, emax; int i, k, n; for(i=0;i<N;i++) { input[i] = (random() % 65535 - 32767) * 9 / 10; input1[i] = input[i]; } mdct512(output, input); /* do it by hand */ for(k=0;k<N/2;k++) { s = 0; for(n=0;n<N;n++) { a = (2*M_PI*(2*n+1+N/2)*(2*k+1) / (4 * N)); s += input1[n] * cos(a); } output1[k] = -2 * s / N; } err = 0; emax = 0; for(i=0;i<N/2;i++) { printf("%3d: %7d %7.0f\n", i, output[i], output1[i]); e = output[i] - output1[i]; if (e > emax) emax = e; err += e * e; } printf("err2=%f emax=%f\n", err / (N/2), emax); } void test_ac3(void) { AC3EncodeContext ctx; unsigned char frame[AC3_MAX_CODED_FRAME_SIZE]; short samples[AC3_FRAME_SIZE]; int ret, i; AC3_encode_init(&ctx, 44100, 64000, 1); fft_test(); mdct_test(); for(i=0;i<AC3_FRAME_SIZE;i++) samples[i] = (int)(sin(2*M_PI*i*1000.0/44100) * 10000); ret = AC3_encode_frame(&ctx, frame, samples); printf("ret=%d\n", ret); } #endif AVCodec ac3_encoder = { "ac3", CODEC_TYPE_AUDIO, CODEC_ID_AC3, sizeof(AC3EncodeContext), AC3_encode_init, AC3_encode_frame, NULL, };