Libav
|
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 }