- * @param nbits log2(MDCT size)
- */
-static av_cold int mdct_init(AVCodecContext *avctx, AC3MDCTContext *mdct,
- int nbits)
-{
- int i, n, n4, ret;
-
- n = 1 << nbits;
- n4 = n >> 2;
-
- mdct->nbits = nbits;
-
- ret = fft_init(avctx, mdct, nbits - 2);
- if (ret)
- return ret;
-
- mdct->window = ff_ac3_window;
-
- FF_ALLOC_OR_GOTO(avctx, mdct->xcos1, n4 * sizeof(*mdct->xcos1), mdct_alloc_fail);
- FF_ALLOC_OR_GOTO(avctx, mdct->xsin1, n4 * sizeof(*mdct->xsin1), mdct_alloc_fail);
- FF_ALLOC_OR_GOTO(avctx, mdct->rot_tmp, n * sizeof(*mdct->rot_tmp), mdct_alloc_fail);
- FF_ALLOC_OR_GOTO(avctx, mdct->cplx_tmp, n4 * sizeof(*mdct->cplx_tmp), mdct_alloc_fail);
-
- for (i = 0; i < n4; i++) {
- float alpha = 2.0 * M_PI * (i + 1.0 / 8.0) / n;
- mdct->xcos1[i] = FIX15(-cos(alpha));
- mdct->xsin1[i] = FIX15(-sin(alpha));
- }
-
- return 0;
-mdct_alloc_fail:
- mdct_end(mdct);
- return AVERROR(ENOMEM);
-}
-
-
-/** Butterfly 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; \
-}
-
-
-/** Complex multiply */
-#define CMUL(pre, pim, are, aim, bre, bim) \
-{ \
- pre = (MUL16(are, bre) - MUL16(aim, bim)) >> 15; \
- pim = (MUL16(are, bim) + MUL16(bre, aim)) >> 15; \
-}
-
-
-/**
- * Calculate a 2^n point complex FFT on 2^ln points.
- * @param z complex input/output samples
- * @param ln log2(FFT size)
- */
-static void fft(AC3MDCTContext *mdct, 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 = av_reverse[j] >> (8 - ln);
- if (k < j)
- FFSWAP(IComplex, z[k], z[j]);
- }
-
- /* 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);
-
- /* 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);
-
- /* 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, mdct->costab[l], -mdct->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);
-}
-
-
-/**
- * Calculate a 512-point MDCT
- * @param out 256 output frequency coefficients
- * @param in 512 windowed input audio samples
- */
-static void mdct512(AC3MDCTContext *mdct, int32_t *out, int16_t *in)
-{
- int i, re, im, n, n2, n4;
- int16_t *rot = mdct->rot_tmp;
- IComplex *x = mdct->cplx_tmp;
-
- n = 1 << mdct->nbits;
- n2 = n >> 1;
- n4 = n >> 2;
-
- /* shift to simplify computations */
- for (i = 0; i <n4; i++)
- rot[i] = -in[i + 3*n4];
- memcpy(&rot[n4], &in[0], 3*n4*sizeof(*in));
-
- /* pre rotation */
- for (i = 0; i < n4; i++) {
- re = ((int)rot[ 2*i] - (int)rot[ n-1-2*i]) >> 1;
- im = -((int)rot[n2+2*i] - (int)rot[n2-1-2*i]) >> 1;
- CMUL(x[i].re, x[i].im, re, im, -mdct->xcos1[i], mdct->xsin1[i]);
- }
-
- fft(mdct, x, mdct->nbits - 2);
-
- /* post rotation */
- for (i = 0; i < n4; i++) {
- re = x[i].re;
- im = x[i].im;
- CMUL(out[n2-1-2*i], out[2*i], re, im, mdct->xsin1[i], mdct->xcos1[i]);
- }
-}
-
-
-/**
- * Apply KBD window to input samples prior to MDCT.
- */
-static void apply_window(DSPContext *dsp, int16_t *output, const int16_t *input,
- const int16_t *window, int n)
-{
- int i;
- int n2 = n >> 1;
-
- for (i = 0; i < n2; i++) {
- output[i] = MUL16(input[i], window[i]) >> 15;
- output[n-i-1] = MUL16(input[n-i-1], window[i]) >> 15;
- }
-}
-
-
-/**
- * Calculate the log2() of the maximum absolute value in an array.
- * @param tab input array
- * @param n number of values in the array
- * @return log2(max(abs(tab[])))
- */
-static int log2_tab(int16_t *tab, int n)
-{
- int i, v;
-
- v = 0;
- for (i = 0; i < n; i++)
- v |= abs(tab[i]);
-
- return av_log2(v);
-}
-
-
-/**
- * Left-shift each value in an array by a specified amount.
- * @param tab input array
- * @param n number of values in the array
- * @param lshift left shift amount