X-Git-Url: https://git.sesse.net/?a=blobdiff_plain;f=libavcodec%2Fmdct.c;h=6f645342732a018a5157fa4e0ac2a02df739770d;hb=8400b126acb153329c418433c377d96afd1d1e02;hp=a0f5671771b3446df1d76deedf5f73780a41f942;hpb=983e3246b7528e6ffc0cc347d54b86e161140fe6;p=ffmpeg diff --git a/libavcodec/mdct.c b/libavcodec/mdct.c index a0f5671771b..6f645342732 100644 --- a/libavcodec/mdct.c +++ b/libavcodec/mdct.c @@ -1,90 +1,107 @@ /* * MDCT/IMDCT transforms - * Copyright (c) 2002 Fabrice Bellard. + * Copyright (c) 2002 Fabrice Bellard * - * This library is free software; you can redistribute it and/or + * This file is part of Libav. + * + * Libav is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public * License as published by the Free Software Foundation; either - * version 2 of the License, or (at your option) any later version. + * version 2.1 of the License, or (at your option) any later version. * - * This library is distributed in the hope that it will be useful, + * Libav 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 * Lesser General Public License for more details. * * You should have received a copy of the GNU Lesser General Public - * License along with this library; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * License along with Libav; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA */ -#include "dsputil.h" + +#include +#include +#include "libavutil/common.h" +#include "libavutil/mathematics.h" +#include "fft.h" +#include "fft-internal.h" /** - * @file mdct.c + * @file * MDCT/IMDCT transforms. */ +#if CONFIG_FFT_FLOAT +# define RSCALE(x) (x) +#else +# define RSCALE(x) ((x) >> 1) +#endif + /** * init MDCT or IMDCT computation. */ -int ff_mdct_init(MDCTContext *s, int nbits, int inverse) +av_cold int ff_mdct_init(FFTContext *s, int nbits, int inverse, double scale) { int n, n4, i; - float alpha; + double alpha, theta; + int tstep; memset(s, 0, sizeof(*s)); n = 1 << nbits; - s->nbits = nbits; - s->n = n; + s->mdct_bits = nbits; + s->mdct_size = n; n4 = n >> 2; - s->tcos = av_malloc(n4 * sizeof(FFTSample)); + s->mdct_permutation = FF_MDCT_PERM_NONE; + + if (ff_fft_init(s, s->mdct_bits - 2, inverse) < 0) + goto fail; + + s->tcos = av_malloc(n/2 * sizeof(FFTSample)); if (!s->tcos) goto fail; - s->tsin = av_malloc(n4 * sizeof(FFTSample)); - if (!s->tsin) + + switch (s->mdct_permutation) { + case FF_MDCT_PERM_NONE: + s->tsin = s->tcos + n4; + tstep = 1; + break; + case FF_MDCT_PERM_INTERLEAVE: + s->tsin = s->tcos + 1; + tstep = 2; + break; + default: goto fail; + } + theta = 1.0 / 8.0 + (scale < 0 ? n4 : 0); + scale = sqrt(fabs(scale)); for(i=0;itcos[i] = -cos(alpha); - s->tsin[i] = -sin(alpha); + alpha = 2 * M_PI * (i + theta) / n; + s->tcos[i*tstep] = FIX15(-cos(alpha) * scale); + s->tsin[i*tstep] = FIX15(-sin(alpha) * scale); } - if (fft_init(&s->fft, s->nbits - 2, inverse) < 0) - goto fail; return 0; fail: - av_freep(&s->tcos); - av_freep(&s->tsin); + ff_mdct_end(s); return -1; } -/* complex multiplication: p = a * b */ -#define CMUL(pre, pim, are, aim, bre, bim) \ -{\ - float _are = (are);\ - float _aim = (aim);\ - float _bre = (bre);\ - float _bim = (bim);\ - (pre) = _are * _bre - _aim * _bim;\ - (pim) = _are * _bim + _aim * _bre;\ -} - /** - * Compute inverse MDCT of size N = 2^nbits - * @param output N samples + * Compute the middle half of the inverse MDCT of size N = 2^nbits, + * thus excluding the parts that can be derived by symmetry + * @param output N/2 samples * @param input N/2 samples - * @param tmp N/2 samples */ -void ff_imdct_calc(MDCTContext *s, FFTSample *output, - const FFTSample *input, FFTSample *tmp) +void ff_imdct_half_c(FFTContext *s, FFTSample *output, const FFTSample *input) { int k, n8, n4, n2, n, j; - const uint16_t *revtab = s->fft.revtab; + const uint16_t *revtab = s->revtab; const FFTSample *tcos = s->tcos; const FFTSample *tsin = s->tsin; const FFTSample *in1, *in2; - FFTComplex *z = (FFTComplex *)tmp; + FFTComplex *z = (FFTComplex *)output; - n = 1 << s->nbits; + n = 1 << s->mdct_bits; n2 = n >> 1; n4 = n >> 2; n8 = n >> 3; @@ -98,25 +115,37 @@ void ff_imdct_calc(MDCTContext *s, FFTSample *output, in1 += 2; in2 -= 2; } - fft_calc(&s->fft, z); + s->fft_calc(s, z); /* post rotation + reordering */ - /* XXX: optimize */ - for(k = 0; k < n4; k++) { - CMUL(z[k].re, z[k].im, z[k].re, z[k].im, tcos[k], tsin[k]); - } for(k = 0; k < n8; k++) { - output[2*k] = -z[n8 + k].im; - output[n2-1-2*k] = z[n8 + k].im; + FFTSample r0, i0, r1, i1; + CMUL(r0, i1, z[n8-k-1].im, z[n8-k-1].re, tsin[n8-k-1], tcos[n8-k-1]); + CMUL(r1, i0, z[n8+k ].im, z[n8+k ].re, tsin[n8+k ], tcos[n8+k ]); + z[n8-k-1].re = r0; + z[n8-k-1].im = i0; + z[n8+k ].re = r1; + z[n8+k ].im = i1; + } +} - output[2*k+1] = z[n8-1-k].re; - output[n2-1-2*k-1] = -z[n8-1-k].re; +/** + * Compute inverse MDCT of size N = 2^nbits + * @param output N samples + * @param input N/2 samples + */ +void ff_imdct_calc_c(FFTContext *s, FFTSample *output, const FFTSample *input) +{ + int k; + int n = 1 << s->mdct_bits; + int n2 = n >> 1; + int n4 = n >> 2; - output[n2 + 2*k]=-z[k+n8].re; - output[n-1- 2*k]=-z[k+n8].re; + ff_imdct_half_c(s, output+n4, input); - output[n2 + 2*k+1]=z[n8-k-1].im; - output[n-2 - 2 * k] = z[n8-k-1].im; + for(k = 0; k < n4; k++) { + output[k] = -output[n2-k-1]; + output[n-k-1] = output[n2+k]; } } @@ -124,19 +153,17 @@ void ff_imdct_calc(MDCTContext *s, FFTSample *output, * Compute MDCT of size N = 2^nbits * @param input N samples * @param out N/2 samples - * @param tmp temporary storage of N/2 samples */ -void ff_mdct_calc(MDCTContext *s, FFTSample *out, - const FFTSample *input, FFTSample *tmp) +void ff_mdct_calc_c(FFTContext *s, FFTSample *out, const FFTSample *input) { int i, j, n, n8, n4, n2, n3; - FFTSample re, im, re1, im1; - const uint16_t *revtab = s->fft.revtab; + FFTDouble re, im; + const uint16_t *revtab = s->revtab; const FFTSample *tcos = s->tcos; const FFTSample *tsin = s->tsin; - FFTComplex *x = (FFTComplex *)tmp; + FFTComplex *x = (FFTComplex *)out; - n = 1 << s->nbits; + n = 1 << s->mdct_bits; n2 = n >> 1; n4 = n >> 2; n8 = n >> 3; @@ -144,32 +171,33 @@ void ff_mdct_calc(MDCTContext *s, FFTSample *out, /* pre rotation */ for(i=0;ifft, x); - + s->fft_calc(s, x); + /* post rotation */ - for(i=0;itcos); - av_freep(&s->tsin); - fft_end(&s->fft); + ff_fft_end(s); }