• Main Page
  • Related Pages
  • Modules
  • Data Structures
  • Files
  • File List
  • Globals

libavcodec/mdct.c

Go to the documentation of this file.
00001 /*
00002  * MDCT/IMDCT transforms
00003  * Copyright (c) 2002 Fabrice Bellard
00004  *
00005  * This file is part of FFmpeg.
00006  *
00007  * FFmpeg is free software; you can redistribute it and/or
00008  * modify it under the terms of the GNU Lesser General Public
00009  * License as published by the Free Software Foundation; either
00010  * version 2.1 of the License, or (at your option) any later version.
00011  *
00012  * FFmpeg is distributed in the hope that it will be useful,
00013  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00014  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00015  * Lesser General Public License for more details.
00016  *
00017  * You should have received a copy of the GNU Lesser General Public
00018  * License along with FFmpeg; if not, write to the Free Software
00019  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
00020  */
00021 #include "dsputil.h"
00022 
00028 // Generate a Kaiser-Bessel Derived Window.
00029 #define BESSEL_I0_ITER 50 // default: 50 iterations of Bessel I0 approximation
00030 av_cold void ff_kbd_window_init(float *window, float alpha, int n)
00031 {
00032    int i, j;
00033    double sum = 0.0, bessel, tmp;
00034    double local_window[n];
00035    double alpha2 = (alpha * M_PI / n) * (alpha * M_PI / n);
00036 
00037    for (i = 0; i < n; i++) {
00038        tmp = i * (n - i) * alpha2;
00039        bessel = 1.0;
00040        for (j = BESSEL_I0_ITER; j > 0; j--)
00041            bessel = bessel * tmp / (j * j) + 1;
00042        sum += bessel;
00043        local_window[i] = sum;
00044    }
00045 
00046    sum++;
00047    for (i = 0; i < n; i++)
00048        window[i] = sqrt(local_window[i] / sum);
00049 }
00050 
00051 DECLARE_ALIGNED(16, float, ff_sine_128 [ 128]);
00052 DECLARE_ALIGNED(16, float, ff_sine_256 [ 256]);
00053 DECLARE_ALIGNED(16, float, ff_sine_512 [ 512]);
00054 DECLARE_ALIGNED(16, float, ff_sine_1024[1024]);
00055 DECLARE_ALIGNED(16, float, ff_sine_2048[2048]);
00056 DECLARE_ALIGNED(16, float, ff_sine_4096[4096]);
00057 float *ff_sine_windows[6] = {
00058     ff_sine_128, ff_sine_256, ff_sine_512, ff_sine_1024, ff_sine_2048, ff_sine_4096
00059 };
00060 
00061 // Generate a sine window.
00062 av_cold void ff_sine_window_init(float *window, int n) {
00063     int i;
00064     for(i = 0; i < n; i++)
00065         window[i] = sinf((i + 0.5) * (M_PI / (2.0 * n)));
00066 }
00067 
00071 av_cold int ff_mdct_init(MDCTContext *s, int nbits, int inverse)
00072 {
00073     int n, n4, i;
00074     double alpha;
00075 
00076     memset(s, 0, sizeof(*s));
00077     n = 1 << nbits;
00078     s->nbits = nbits;
00079     s->n = n;
00080     n4 = n >> 2;
00081     s->tcos = av_malloc(n4 * sizeof(FFTSample));
00082     if (!s->tcos)
00083         goto fail;
00084     s->tsin = av_malloc(n4 * sizeof(FFTSample));
00085     if (!s->tsin)
00086         goto fail;
00087 
00088     for(i=0;i<n4;i++) {
00089         alpha = 2 * M_PI * (i + 1.0 / 8.0) / n;
00090         s->tcos[i] = -cos(alpha);
00091         s->tsin[i] = -sin(alpha);
00092     }
00093     if (ff_fft_init(&s->fft, s->nbits - 2, inverse) < 0)
00094         goto fail;
00095     return 0;
00096  fail:
00097     av_freep(&s->tcos);
00098     av_freep(&s->tsin);
00099     return -1;
00100 }
00101 
00102 /* complex multiplication: p = a * b */
00103 #define CMUL(pre, pim, are, aim, bre, bim) \
00104 {\
00105     FFTSample _are = (are);\
00106     FFTSample _aim = (aim);\
00107     FFTSample _bre = (bre);\
00108     FFTSample _bim = (bim);\
00109     (pre) = _are * _bre - _aim * _bim;\
00110     (pim) = _are * _bim + _aim * _bre;\
00111 }
00112 
00119 void ff_imdct_half_c(MDCTContext *s, FFTSample *output, const FFTSample *input)
00120 {
00121     int k, n8, n4, n2, n, j;
00122     const uint16_t *revtab = s->fft.revtab;
00123     const FFTSample *tcos = s->tcos;
00124     const FFTSample *tsin = s->tsin;
00125     const FFTSample *in1, *in2;
00126     FFTComplex *z = (FFTComplex *)output;
00127 
00128     n = 1 << s->nbits;
00129     n2 = n >> 1;
00130     n4 = n >> 2;
00131     n8 = n >> 3;
00132 
00133     /* pre rotation */
00134     in1 = input;
00135     in2 = input + n2 - 1;
00136     for(k = 0; k < n4; k++) {
00137         j=revtab[k];
00138         CMUL(z[j].re, z[j].im, *in2, *in1, tcos[k], tsin[k]);
00139         in1 += 2;
00140         in2 -= 2;
00141     }
00142     ff_fft_calc(&s->fft, z);
00143 
00144     /* post rotation + reordering */
00145     output += n4;
00146     for(k = 0; k < n8; k++) {
00147         FFTSample r0, i0, r1, i1;
00148         CMUL(r0, i1, z[n8-k-1].im, z[n8-k-1].re, tsin[n8-k-1], tcos[n8-k-1]);
00149         CMUL(r1, i0, z[n8+k  ].im, z[n8+k  ].re, tsin[n8+k  ], tcos[n8+k  ]);
00150         z[n8-k-1].re = r0;
00151         z[n8-k-1].im = i0;
00152         z[n8+k  ].re = r1;
00153         z[n8+k  ].im = i1;
00154     }
00155 }
00156 
00162 void ff_imdct_calc_c(MDCTContext *s, FFTSample *output, const FFTSample *input)
00163 {
00164     int k;
00165     int n = 1 << s->nbits;
00166     int n2 = n >> 1;
00167     int n4 = n >> 2;
00168 
00169     ff_imdct_half_c(s, output+n4, input);
00170 
00171     for(k = 0; k < n4; k++) {
00172         output[k] = -output[n2-k-1];
00173         output[n-k-1] = output[n2+k];
00174     }
00175 }
00176 
00182 void ff_mdct_calc_c(MDCTContext *s, FFTSample *out, const FFTSample *input)
00183 {
00184     int i, j, n, n8, n4, n2, n3;
00185     FFTSample re, im;
00186     const uint16_t *revtab = s->fft.revtab;
00187     const FFTSample *tcos = s->tcos;
00188     const FFTSample *tsin = s->tsin;
00189     FFTComplex *x = (FFTComplex *)out;
00190 
00191     n = 1 << s->nbits;
00192     n2 = n >> 1;
00193     n4 = n >> 2;
00194     n8 = n >> 3;
00195     n3 = 3 * n4;
00196 
00197     /* pre rotation */
00198     for(i=0;i<n8;i++) {
00199         re = -input[2*i+3*n4] - input[n3-1-2*i];
00200         im = -input[n4+2*i] + input[n4-1-2*i];
00201         j = revtab[i];
00202         CMUL(x[j].re, x[j].im, re, im, -tcos[i], tsin[i]);
00203 
00204         re = input[2*i] - input[n2-1-2*i];
00205         im = -(input[n2+2*i] + input[n-1-2*i]);
00206         j = revtab[n8 + i];
00207         CMUL(x[j].re, x[j].im, re, im, -tcos[n8 + i], tsin[n8 + i]);
00208     }
00209 
00210     ff_fft_calc(&s->fft, x);
00211 
00212     /* post rotation */
00213     for(i=0;i<n8;i++) {
00214         FFTSample r0, i0, r1, i1;
00215         CMUL(i1, r0, x[n8-i-1].re, x[n8-i-1].im, -tsin[n8-i-1], -tcos[n8-i-1]);
00216         CMUL(i0, r1, x[n8+i  ].re, x[n8+i  ].im, -tsin[n8+i  ], -tcos[n8+i  ]);
00217         x[n8-i-1].re = r0;
00218         x[n8-i-1].im = i0;
00219         x[n8+i  ].re = r1;
00220         x[n8+i  ].im = i1;
00221     }
00222 }
00223 
00224 av_cold void ff_mdct_end(MDCTContext *s)
00225 {
00226     av_freep(&s->tcos);
00227     av_freep(&s->tsin);
00228     ff_fft_end(&s->fft);
00229 }

Generated on Sat Feb 16 2013 09:23:12 for ffmpeg by  doxygen 1.7.1