00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00029 #include "dsputil.h"
00030
00031
00032 DECLARE_ALIGNED_16(FFTSample, ff_cos_16[8]);
00033 DECLARE_ALIGNED_16(FFTSample, ff_cos_32[16]);
00034 DECLARE_ALIGNED_16(FFTSample, ff_cos_64[32]);
00035 DECLARE_ALIGNED_16(FFTSample, ff_cos_128[64]);
00036 DECLARE_ALIGNED_16(FFTSample, ff_cos_256[128]);
00037 DECLARE_ALIGNED_16(FFTSample, ff_cos_512[256]);
00038 DECLARE_ALIGNED_16(FFTSample, ff_cos_1024[512]);
00039 DECLARE_ALIGNED_16(FFTSample, ff_cos_2048[1024]);
00040 DECLARE_ALIGNED_16(FFTSample, ff_cos_4096[2048]);
00041 DECLARE_ALIGNED_16(FFTSample, ff_cos_8192[4096]);
00042 DECLARE_ALIGNED_16(FFTSample, ff_cos_16384[8192]);
00043 DECLARE_ALIGNED_16(FFTSample, ff_cos_32768[16384]);
00044 DECLARE_ALIGNED_16(FFTSample, ff_cos_65536[32768]);
00045 FFTSample *ff_cos_tabs[] = {
00046 ff_cos_16, ff_cos_32, ff_cos_64, ff_cos_128, ff_cos_256, ff_cos_512, ff_cos_1024,
00047 ff_cos_2048, ff_cos_4096, ff_cos_8192, ff_cos_16384, ff_cos_32768, ff_cos_65536,
00048 };
00049
00050 static int split_radix_permutation(int i, int n, int inverse)
00051 {
00052 int m;
00053 if(n <= 2) return i&1;
00054 m = n >> 1;
00055 if(!(i&m)) return split_radix_permutation(i, m, inverse)*2;
00056 m >>= 1;
00057 if(inverse == !(i&m)) return split_radix_permutation(i, m, inverse)*4 + 1;
00058 else return split_radix_permutation(i, m, inverse)*4 - 1;
00059 }
00060
00061 av_cold int ff_fft_init(FFTContext *s, int nbits, int inverse)
00062 {
00063 int i, j, m, n;
00064 float alpha, c1, s1, s2;
00065 int split_radix = 1;
00066 int av_unused has_vectors;
00067 int revtab_shift = 0;
00068
00069 if (nbits < 2 || nbits > 16)
00070 goto fail;
00071 s->nbits = nbits;
00072 n = 1 << nbits;
00073
00074 s->tmp_buf = NULL;
00075 s->exptab = av_malloc((n / 2) * sizeof(FFTComplex));
00076 if (!s->exptab)
00077 goto fail;
00078 s->revtab = av_malloc(n * sizeof(uint16_t));
00079 if (!s->revtab)
00080 goto fail;
00081 s->inverse = inverse;
00082
00083 s2 = inverse ? 1.0 : -1.0;
00084
00085 s->fft_permute = ff_fft_permute_c;
00086 s->fft_calc = ff_fft_calc_c;
00087 s->imdct_calc = ff_imdct_calc_c;
00088 s->imdct_half = ff_imdct_half_c;
00089 s->mdct_calc = ff_mdct_calc_c;
00090 s->exptab1 = NULL;
00091
00092 #if HAVE_MMX && HAVE_YASM
00093 has_vectors = mm_support();
00094 if (has_vectors & FF_MM_SSE && HAVE_SSE) {
00095
00096 s->imdct_calc = ff_imdct_calc_sse;
00097 s->imdct_half = ff_imdct_half_sse;
00098 s->fft_permute = ff_fft_permute_sse;
00099 s->fft_calc = ff_fft_calc_sse;
00100 } else if (has_vectors & FF_MM_3DNOWEXT && HAVE_AMD3DNOWEXT) {
00101
00102 s->imdct_calc = ff_imdct_calc_3dn2;
00103 s->imdct_half = ff_imdct_half_3dn2;
00104 s->fft_calc = ff_fft_calc_3dn2;
00105 } else if (has_vectors & FF_MM_3DNOW && HAVE_AMD3DNOW) {
00106
00107 s->imdct_calc = ff_imdct_calc_3dn;
00108 s->imdct_half = ff_imdct_half_3dn;
00109 s->fft_calc = ff_fft_calc_3dn;
00110 }
00111 #elif HAVE_ALTIVEC && !defined ALTIVEC_USE_REFERENCE_C_CODE
00112 has_vectors = mm_support();
00113 if (has_vectors & FF_MM_ALTIVEC) {
00114 s->fft_calc = ff_fft_calc_altivec;
00115 split_radix = 0;
00116 }
00117 #elif HAVE_NEON
00118 s->fft_permute = ff_fft_permute_neon;
00119 s->fft_calc = ff_fft_calc_neon;
00120 s->imdct_calc = ff_imdct_calc_neon;
00121 s->imdct_half = ff_imdct_half_neon;
00122 s->mdct_calc = ff_mdct_calc_neon;
00123 revtab_shift = 3;
00124 #endif
00125
00126 if (split_radix) {
00127 for(j=4; j<=nbits; j++) {
00128 int m = 1<<j;
00129 double freq = 2*M_PI/m;
00130 FFTSample *tab = ff_cos_tabs[j-4];
00131 for(i=0; i<=m/4; i++)
00132 tab[i] = cos(i*freq);
00133 for(i=1; i<m/4; i++)
00134 tab[m/2-i] = tab[i];
00135 }
00136 for(i=0; i<n; i++)
00137 s->revtab[-split_radix_permutation(i, n, s->inverse) & (n-1)] =
00138 i << revtab_shift;
00139 s->tmp_buf = av_malloc(n * sizeof(FFTComplex));
00140 } else {
00141 int np, nblocks, np2, l;
00142 FFTComplex *q;
00143
00144 for(i=0; i<(n/2); i++) {
00145 alpha = 2 * M_PI * (float)i / (float)n;
00146 c1 = cos(alpha);
00147 s1 = sin(alpha) * s2;
00148 s->exptab[i].re = c1;
00149 s->exptab[i].im = s1;
00150 }
00151
00152 np = 1 << nbits;
00153 nblocks = np >> 3;
00154 np2 = np >> 1;
00155 s->exptab1 = av_malloc(np * 2 * sizeof(FFTComplex));
00156 if (!s->exptab1)
00157 goto fail;
00158 q = s->exptab1;
00159 do {
00160 for(l = 0; l < np2; l += 2 * nblocks) {
00161 *q++ = s->exptab[l];
00162 *q++ = s->exptab[l + nblocks];
00163
00164 q->re = -s->exptab[l].im;
00165 q->im = s->exptab[l].re;
00166 q++;
00167 q->re = -s->exptab[l + nblocks].im;
00168 q->im = s->exptab[l + nblocks].re;
00169 q++;
00170 }
00171 nblocks = nblocks >> 1;
00172 } while (nblocks != 0);
00173 av_freep(&s->exptab);
00174
00175
00176 for(i=0;i<n;i++) {
00177 m=0;
00178 for(j=0;j<nbits;j++) {
00179 m |= ((i >> j) & 1) << (nbits-j-1);
00180 }
00181 s->revtab[i]=m;
00182 }
00183 }
00184
00185 return 0;
00186 fail:
00187 av_freep(&s->revtab);
00188 av_freep(&s->exptab);
00189 av_freep(&s->exptab1);
00190 av_freep(&s->tmp_buf);
00191 return -1;
00192 }
00193
00194 void ff_fft_permute_c(FFTContext *s, FFTComplex *z)
00195 {
00196 int j, k, np;
00197 FFTComplex tmp;
00198 const uint16_t *revtab = s->revtab;
00199 np = 1 << s->nbits;
00200
00201 if (s->tmp_buf) {
00202
00203 for(j=0;j<np;j++) s->tmp_buf[revtab[j]] = z[j];
00204 memcpy(z, s->tmp_buf, np * sizeof(FFTComplex));
00205 return;
00206 }
00207
00208
00209 for(j=0;j<np;j++) {
00210 k = revtab[j];
00211 if (k < j) {
00212 tmp = z[k];
00213 z[k] = z[j];
00214 z[j] = tmp;
00215 }
00216 }
00217 }
00218
00219 av_cold void ff_fft_end(FFTContext *s)
00220 {
00221 av_freep(&s->revtab);
00222 av_freep(&s->exptab);
00223 av_freep(&s->exptab1);
00224 av_freep(&s->tmp_buf);
00225 }
00226
00227 #define sqrthalf (float)M_SQRT1_2
00228
00229 #define BF(x,y,a,b) {\
00230 x = a - b;\
00231 y = a + b;\
00232 }
00233
00234 #define BUTTERFLIES(a0,a1,a2,a3) {\
00235 BF(t3, t5, t5, t1);\
00236 BF(a2.re, a0.re, a0.re, t5);\
00237 BF(a3.im, a1.im, a1.im, t3);\
00238 BF(t4, t6, t2, t6);\
00239 BF(a3.re, a1.re, a1.re, t4);\
00240 BF(a2.im, a0.im, a0.im, t6);\
00241 }
00242
00243
00244
00245
00246 #define BUTTERFLIES_BIG(a0,a1,a2,a3) {\
00247 FFTSample r0=a0.re, i0=a0.im, r1=a1.re, i1=a1.im;\
00248 BF(t3, t5, t5, t1);\
00249 BF(a2.re, a0.re, r0, t5);\
00250 BF(a3.im, a1.im, i1, t3);\
00251 BF(t4, t6, t2, t6);\
00252 BF(a3.re, a1.re, r1, t4);\
00253 BF(a2.im, a0.im, i0, t6);\
00254 }
00255
00256 #define TRANSFORM(a0,a1,a2,a3,wre,wim) {\
00257 t1 = a2.re * wre + a2.im * wim;\
00258 t2 = a2.im * wre - a2.re * wim;\
00259 t5 = a3.re * wre - a3.im * wim;\
00260 t6 = a3.im * wre + a3.re * wim;\
00261 BUTTERFLIES(a0,a1,a2,a3)\
00262 }
00263
00264 #define TRANSFORM_ZERO(a0,a1,a2,a3) {\
00265 t1 = a2.re;\
00266 t2 = a2.im;\
00267 t5 = a3.re;\
00268 t6 = a3.im;\
00269 BUTTERFLIES(a0,a1,a2,a3)\
00270 }
00271
00272
00273 #define PASS(name)\
00274 static void name(FFTComplex *z, const FFTSample *wre, unsigned int n)\
00275 {\
00276 FFTSample t1, t2, t3, t4, t5, t6;\
00277 int o1 = 2*n;\
00278 int o2 = 4*n;\
00279 int o3 = 6*n;\
00280 const FFTSample *wim = wre+o1;\
00281 n--;\
00282 \
00283 TRANSFORM_ZERO(z[0],z[o1],z[o2],z[o3]);\
00284 TRANSFORM(z[1],z[o1+1],z[o2+1],z[o3+1],wre[1],wim[-1]);\
00285 do {\
00286 z += 2;\
00287 wre += 2;\
00288 wim -= 2;\
00289 TRANSFORM(z[0],z[o1],z[o2],z[o3],wre[0],wim[0]);\
00290 TRANSFORM(z[1],z[o1+1],z[o2+1],z[o3+1],wre[1],wim[-1]);\
00291 } while(--n);\
00292 }
00293
00294 PASS(pass)
00295 #undef BUTTERFLIES
00296 #define BUTTERFLIES BUTTERFLIES_BIG
00297 PASS(pass_big)
00298
00299 #define DECL_FFT(n,n2,n4)\
00300 static void fft##n(FFTComplex *z)\
00301 {\
00302 fft##n2(z);\
00303 fft##n4(z+n4*2);\
00304 fft##n4(z+n4*3);\
00305 pass(z,ff_cos_##n,n4/2);\
00306 }
00307
00308 static void fft4(FFTComplex *z)
00309 {
00310 FFTSample t1, t2, t3, t4, t5, t6, t7, t8;
00311
00312 BF(t3, t1, z[0].re, z[1].re);
00313 BF(t8, t6, z[3].re, z[2].re);
00314 BF(z[2].re, z[0].re, t1, t6);
00315 BF(t4, t2, z[0].im, z[1].im);
00316 BF(t7, t5, z[2].im, z[3].im);
00317 BF(z[3].im, z[1].im, t4, t8);
00318 BF(z[3].re, z[1].re, t3, t7);
00319 BF(z[2].im, z[0].im, t2, t5);
00320 }
00321
00322 static void fft8(FFTComplex *z)
00323 {
00324 FFTSample t1, t2, t3, t4, t5, t6, t7, t8;
00325
00326 fft4(z);
00327
00328 BF(t1, z[5].re, z[4].re, -z[5].re);
00329 BF(t2, z[5].im, z[4].im, -z[5].im);
00330 BF(t3, z[7].re, z[6].re, -z[7].re);
00331 BF(t4, z[7].im, z[6].im, -z[7].im);
00332 BF(t8, t1, t3, t1);
00333 BF(t7, t2, t2, t4);
00334 BF(z[4].re, z[0].re, z[0].re, t1);
00335 BF(z[4].im, z[0].im, z[0].im, t2);
00336 BF(z[6].re, z[2].re, z[2].re, t7);
00337 BF(z[6].im, z[2].im, z[2].im, t8);
00338
00339 TRANSFORM(z[1],z[3],z[5],z[7],sqrthalf,sqrthalf);
00340 }
00341
00342 #if !CONFIG_SMALL
00343 static void fft16(FFTComplex *z)
00344 {
00345 FFTSample t1, t2, t3, t4, t5, t6;
00346
00347 fft8(z);
00348 fft4(z+8);
00349 fft4(z+12);
00350
00351 TRANSFORM_ZERO(z[0],z[4],z[8],z[12]);
00352 TRANSFORM(z[2],z[6],z[10],z[14],sqrthalf,sqrthalf);
00353 TRANSFORM(z[1],z[5],z[9],z[13],ff_cos_16[1],ff_cos_16[3]);
00354 TRANSFORM(z[3],z[7],z[11],z[15],ff_cos_16[3],ff_cos_16[1]);
00355 }
00356 #else
00357 DECL_FFT(16,8,4)
00358 #endif
00359 DECL_FFT(32,16,8)
00360 DECL_FFT(64,32,16)
00361 DECL_FFT(128,64,32)
00362 DECL_FFT(256,128,64)
00363 DECL_FFT(512,256,128)
00364 #if !CONFIG_SMALL
00365 #define pass pass_big
00366 #endif
00367 DECL_FFT(1024,512,256)
00368 DECL_FFT(2048,1024,512)
00369 DECL_FFT(4096,2048,1024)
00370 DECL_FFT(8192,4096,2048)
00371 DECL_FFT(16384,8192,4096)
00372 DECL_FFT(32768,16384,8192)
00373 DECL_FFT(65536,32768,16384)
00374
00375 static void (*fft_dispatch[])(FFTComplex*) = {
00376 fft4, fft8, fft16, fft32, fft64, fft128, fft256, fft512, fft1024,
00377 fft2048, fft4096, fft8192, fft16384, fft32768, fft65536,
00378 };
00379
00380 void ff_fft_calc_c(FFTContext *s, FFTComplex *z)
00381 {
00382 fft_dispatch[s->nbits-2](z);
00383 }
00384