00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
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
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[FF_KBD_WINDOW_MAX];
00040 double alpha2 = (alpha * M_PI / n) * (alpha * M_PI / n);
00041
00042 assert(n <= FF_KBD_WINDOW_MAX);
00043
00044 for (i = 0; i < n; i++) {
00045 tmp = i * (n - i) * alpha2;
00046 bessel = 1.0;
00047 for (j = BESSEL_I0_ITER; j > 0; j--)
00048 bessel = bessel * tmp / (j * j) + 1;
00049 sum += bessel;
00050 local_window[i] = sum;
00051 }
00052
00053 sum++;
00054 for (i = 0; i < n; i++)
00055 window[i] = sqrt(local_window[i] / sum);
00056 }
00057
00058 #include "mdct_tablegen.h"
00059
00063 av_cold int ff_mdct_init(FFTContext *s, int nbits, int inverse, double scale)
00064 {
00065 int n, n4, i;
00066 double alpha, theta;
00067 int tstep;
00068
00069 memset(s, 0, sizeof(*s));
00070 n = 1 << nbits;
00071 s->mdct_bits = nbits;
00072 s->mdct_size = n;
00073 n4 = n >> 2;
00074 s->mdct_permutation = FF_MDCT_PERM_NONE;
00075
00076 if (ff_fft_init(s, s->mdct_bits - 2, inverse) < 0)
00077 goto fail;
00078
00079 s->tcos = av_malloc(n/2 * sizeof(FFTSample));
00080 if (!s->tcos)
00081 goto fail;
00082
00083 switch (s->mdct_permutation) {
00084 case FF_MDCT_PERM_NONE:
00085 s->tsin = s->tcos + n4;
00086 tstep = 1;
00087 break;
00088 case FF_MDCT_PERM_INTERLEAVE:
00089 s->tsin = s->tcos + 1;
00090 tstep = 2;
00091 break;
00092 default:
00093 goto fail;
00094 }
00095
00096 theta = 1.0 / 8.0 + (scale < 0 ? n4 : 0);
00097 scale = sqrt(fabs(scale));
00098 for(i=0;i<n4;i++) {
00099 alpha = 2 * M_PI * (i + theta) / n;
00100 s->tcos[i*tstep] = -cos(alpha) * scale;
00101 s->tsin[i*tstep] = -sin(alpha) * scale;
00102 }
00103 return 0;
00104 fail:
00105 ff_mdct_end(s);
00106 return -1;
00107 }
00108
00109
00110 #define CMUL(pre, pim, are, aim, bre, bim) \
00111 {\
00112 FFTSample _are = (are);\
00113 FFTSample _aim = (aim);\
00114 FFTSample _bre = (bre);\
00115 FFTSample _bim = (bim);\
00116 (pre) = _are * _bre - _aim * _bim;\
00117 (pim) = _are * _bim + _aim * _bre;\
00118 }
00119
00126 void ff_imdct_half_c(FFTContext *s, FFTSample *output, const FFTSample *input)
00127 {
00128 int k, n8, n4, n2, n, j;
00129 const uint16_t *revtab = s->revtab;
00130 const FFTSample *tcos = s->tcos;
00131 const FFTSample *tsin = s->tsin;
00132 const FFTSample *in1, *in2;
00133 FFTComplex *z = (FFTComplex *)output;
00134
00135 n = 1 << s->mdct_bits;
00136 n2 = n >> 1;
00137 n4 = n >> 2;
00138 n8 = n >> 3;
00139
00140
00141 in1 = input;
00142 in2 = input + n2 - 1;
00143 for(k = 0; k < n4; k++) {
00144 j=revtab[k];
00145 CMUL(z[j].re, z[j].im, *in2, *in1, tcos[k], tsin[k]);
00146 in1 += 2;
00147 in2 -= 2;
00148 }
00149 ff_fft_calc(s, z);
00150
00151
00152 for(k = 0; k < n8; k++) {
00153 FFTSample r0, i0, r1, i1;
00154 CMUL(r0, i1, z[n8-k-1].im, z[n8-k-1].re, tsin[n8-k-1], tcos[n8-k-1]);
00155 CMUL(r1, i0, z[n8+k ].im, z[n8+k ].re, tsin[n8+k ], tcos[n8+k ]);
00156 z[n8-k-1].re = r0;
00157 z[n8-k-1].im = i0;
00158 z[n8+k ].re = r1;
00159 z[n8+k ].im = i1;
00160 }
00161 }
00162
00168 void ff_imdct_calc_c(FFTContext *s, FFTSample *output, const FFTSample *input)
00169 {
00170 int k;
00171 int n = 1 << s->mdct_bits;
00172 int n2 = n >> 1;
00173 int n4 = n >> 2;
00174
00175 ff_imdct_half_c(s, output+n4, input);
00176
00177 for(k = 0; k < n4; k++) {
00178 output[k] = -output[n2-k-1];
00179 output[n-k-1] = output[n2+k];
00180 }
00181 }
00182
00188 void ff_mdct_calc_c(FFTContext *s, FFTSample *out, const FFTSample *input)
00189 {
00190 int i, j, n, n8, n4, n2, n3;
00191 FFTSample re, im;
00192 const uint16_t *revtab = s->revtab;
00193 const FFTSample *tcos = s->tcos;
00194 const FFTSample *tsin = s->tsin;
00195 FFTComplex *x = (FFTComplex *)out;
00196
00197 n = 1 << s->mdct_bits;
00198 n2 = n >> 1;
00199 n4 = n >> 2;
00200 n8 = n >> 3;
00201 n3 = 3 * n4;
00202
00203
00204 for(i=0;i<n8;i++) {
00205 re = -input[2*i+n3] - input[n3-1-2*i];
00206 im = -input[n4+2*i] + input[n4-1-2*i];
00207 j = revtab[i];
00208 CMUL(x[j].re, x[j].im, re, im, -tcos[i], tsin[i]);
00209
00210 re = input[2*i] - input[n2-1-2*i];
00211 im = -(input[n2+2*i] + input[n-1-2*i]);
00212 j = revtab[n8 + i];
00213 CMUL(x[j].re, x[j].im, re, im, -tcos[n8 + i], tsin[n8 + i]);
00214 }
00215
00216 ff_fft_calc(s, x);
00217
00218
00219 for(i=0;i<n8;i++) {
00220 FFTSample r0, i0, r1, i1;
00221 CMUL(i1, r0, x[n8-i-1].re, x[n8-i-1].im, -tsin[n8-i-1], -tcos[n8-i-1]);
00222 CMUL(i0, r1, x[n8+i ].re, x[n8+i ].im, -tsin[n8+i ], -tcos[n8+i ]);
00223 x[n8-i-1].re = r0;
00224 x[n8-i-1].im = i0;
00225 x[n8+i ].re = r1;
00226 x[n8+i ].im = i1;
00227 }
00228 }
00229
00230 av_cold void ff_mdct_end(FFTContext *s)
00231 {
00232 av_freep(&s->tcos);
00233 ff_fft_end(s);
00234 }