• 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 
00022 #include <stdlib.h>
00023 #include <string.h>
00024 #include "libavutil/common.h"
00025 #include "libavutil/mathematics.h"
00026 #include "fft.h"
00027 
00033 // Generate a Kaiser-Bessel Derived Window.
00034 #define BESSEL_I0_ITER 50 // default: 50 iterations of Bessel I0 approximation
00035 av_cold void ff_kbd_window_init(float *window, float alpha, int n)
00036 {
00037    int i, j;
00038    double sum = 0.0, bessel, tmp;
00039    double local_window[n];
00040    double alpha2 = (alpha * M_PI / n) * (alpha * M_PI / n);
00041 
00042    for (i = 0; i < n; i++) {
00043        tmp = i * (n - i) * alpha2;
00044        bessel = 1.0;
00045        for (j = BESSEL_I0_ITER; j > 0; j--)
00046            bessel = bessel * tmp / (j * j) + 1;
00047        sum += bessel;
00048        local_window[i] = sum;
00049    }
00050 
00051    sum++;
00052    for (i = 0; i < n; i++)
00053        window[i] = sqrt(local_window[i] / sum);
00054 }
00055 
00056 #include "mdct_tablegen.h"
00057 
00061 av_cold int ff_mdct_init(FFTContext *s, int nbits, int inverse, double scale)
00062 {
00063     int n, n4, i;
00064     double alpha, theta;
00065     int tstep;
00066 
00067     memset(s, 0, sizeof(*s));
00068     n = 1 << nbits;
00069     s->mdct_bits = nbits;
00070     s->mdct_size = n;
00071     n4 = n >> 2;
00072     s->permutation = FF_MDCT_PERM_NONE;
00073 
00074     if (ff_fft_init(s, s->mdct_bits - 2, inverse) < 0)
00075         goto fail;
00076 
00077     s->tcos = av_malloc(n/2 * sizeof(FFTSample));
00078     if (!s->tcos)
00079         goto fail;
00080 
00081     switch (s->permutation) {
00082     case FF_MDCT_PERM_NONE:
00083         s->tsin = s->tcos + n4;
00084         tstep = 1;
00085         break;
00086     case FF_MDCT_PERM_INTERLEAVE:
00087         s->tsin = s->tcos + 1;
00088         tstep = 2;
00089         break;
00090     default:
00091         goto fail;
00092     }
00093 
00094     theta = 1.0 / 8.0 + (scale < 0 ? n4 : 0);
00095     scale = sqrt(fabs(scale));
00096     for(i=0;i<n4;i++) {
00097         alpha = 2 * M_PI * (i + theta) / n;
00098         s->tcos[i*tstep] = -cos(alpha) * scale;
00099         s->tsin[i*tstep] = -sin(alpha) * scale;
00100     }
00101     return 0;
00102  fail:
00103     ff_mdct_end(s);
00104     return -1;
00105 }
00106 
00107 /* complex multiplication: p = a * b */
00108 #define CMUL(pre, pim, are, aim, bre, bim) \
00109 {\
00110     FFTSample _are = (are);\
00111     FFTSample _aim = (aim);\
00112     FFTSample _bre = (bre);\
00113     FFTSample _bim = (bim);\
00114     (pre) = _are * _bre - _aim * _bim;\
00115     (pim) = _are * _bim + _aim * _bre;\
00116 }
00117 
00124 void ff_imdct_half_c(FFTContext *s, FFTSample *output, const FFTSample *input)
00125 {
00126     int k, n8, n4, n2, n, j;
00127     const uint16_t *revtab = s->revtab;
00128     const FFTSample *tcos = s->tcos;
00129     const FFTSample *tsin = s->tsin;
00130     const FFTSample *in1, *in2;
00131     FFTComplex *z = (FFTComplex *)output;
00132 
00133     n = 1 << s->mdct_bits;
00134     n2 = n >> 1;
00135     n4 = n >> 2;
00136     n8 = n >> 3;
00137 
00138     /* pre rotation */
00139     in1 = input;
00140     in2 = input + n2 - 1;
00141     for(k = 0; k < n4; k++) {
00142         j=revtab[k];
00143         CMUL(z[j].re, z[j].im, *in2, *in1, tcos[k], tsin[k]);
00144         in1 += 2;
00145         in2 -= 2;
00146     }
00147     ff_fft_calc(s, z);
00148 
00149     /* post rotation + reordering */
00150     for(k = 0; k < n8; k++) {
00151         FFTSample r0, i0, r1, i1;
00152         CMUL(r0, i1, z[n8-k-1].im, z[n8-k-1].re, tsin[n8-k-1], tcos[n8-k-1]);
00153         CMUL(r1, i0, z[n8+k  ].im, z[n8+k  ].re, tsin[n8+k  ], tcos[n8+k  ]);
00154         z[n8-k-1].re = r0;
00155         z[n8-k-1].im = i0;
00156         z[n8+k  ].re = r1;
00157         z[n8+k  ].im = i1;
00158     }
00159 }
00160 
00166 void ff_imdct_calc_c(FFTContext *s, FFTSample *output, const FFTSample *input)
00167 {
00168     int k;
00169     int n = 1 << s->mdct_bits;
00170     int n2 = n >> 1;
00171     int n4 = n >> 2;
00172 
00173     ff_imdct_half_c(s, output+n4, input);
00174 
00175     for(k = 0; k < n4; k++) {
00176         output[k] = -output[n2-k-1];
00177         output[n-k-1] = output[n2+k];
00178     }
00179 }
00180 
00186 void ff_mdct_calc_c(FFTContext *s, FFTSample *out, const FFTSample *input)
00187 {
00188     int i, j, n, n8, n4, n2, n3;
00189     FFTSample re, im;
00190     const uint16_t *revtab = s->revtab;
00191     const FFTSample *tcos = s->tcos;
00192     const FFTSample *tsin = s->tsin;
00193     FFTComplex *x = (FFTComplex *)out;
00194 
00195     n = 1 << s->mdct_bits;
00196     n2 = n >> 1;
00197     n4 = n >> 2;
00198     n8 = n >> 3;
00199     n3 = 3 * n4;
00200 
00201     /* pre rotation */
00202     for(i=0;i<n8;i++) {
00203         re = -input[2*i+3*n4] - input[n3-1-2*i];
00204         im = -input[n4+2*i] + input[n4-1-2*i];
00205         j = revtab[i];
00206         CMUL(x[j].re, x[j].im, re, im, -tcos[i], tsin[i]);
00207 
00208         re = input[2*i] - input[n2-1-2*i];
00209         im = -(input[n2+2*i] + input[n-1-2*i]);
00210         j = revtab[n8 + i];
00211         CMUL(x[j].re, x[j].im, re, im, -tcos[n8 + i], tsin[n8 + i]);
00212     }
00213 
00214     ff_fft_calc(s, x);
00215 
00216     /* post rotation */
00217     for(i=0;i<n8;i++) {
00218         FFTSample r0, i0, r1, i1;
00219         CMUL(i1, r0, x[n8-i-1].re, x[n8-i-1].im, -tsin[n8-i-1], -tcos[n8-i-1]);
00220         CMUL(i0, r1, x[n8+i  ].re, x[n8+i  ].im, -tsin[n8+i  ], -tcos[n8+i  ]);
00221         x[n8-i-1].re = r0;
00222         x[n8-i-1].im = i0;
00223         x[n8+i  ].re = r1;
00224         x[n8+i  ].im = i1;
00225     }
00226 }
00227 
00228 av_cold void ff_mdct_end(FFTContext *s)
00229 {
00230     av_freep(&s->tcos);
00231     ff_fft_end(s);
00232 }

Generated on Fri Sep 16 2011 17:17:39 for FFmpeg by  doxygen 1.7.1