76e9c4139e91ab429b56447ebf34fb41b91b04f6
3 * Copyright (c) 2008 Loren Merritt
4 * Copyright (c) 2002 Fabrice Bellard
5 * Partly based on libdjbfft by D. J. Bernstein
7 * This file is part of Libav.
9 * Libav is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public
11 * License as published by the Free Software Foundation; either
12 * version 2.1 of the License, or (at your option) any later version.
14 * Libav is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with Libav; if not, write to the Free Software
21 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
26 * FFT/IFFT transforms.
31 #include "libavutil/mathematics.h"
34 /* cos(2*pi*x/n) for 0<=x<=n/4, followed by its reverse */
35 #if !CONFIG_HARDCODED_TABLES
50 COSTABLE_CONST FFTSample
* const ff_cos_tabs
[] = {
51 NULL
, NULL
, NULL
, NULL
,
52 ff_cos_16
, ff_cos_32
, ff_cos_64
, ff_cos_128
, ff_cos_256
, ff_cos_512
, ff_cos_1024
,
53 ff_cos_2048
, ff_cos_4096
, ff_cos_8192
, ff_cos_16384
, ff_cos_32768
, ff_cos_65536
,
56 static void ff_fft_permute_c(FFTContext
*s
, FFTComplex
*z
);
57 static void ff_fft_calc_c(FFTContext
*s
, FFTComplex
*z
);
59 static int split_radix_permutation(int i
, int n
, int inverse
)
62 if(n
<= 2) return i
&1;
64 if(!(i
&m
)) return split_radix_permutation(i
, m
, inverse
)*2;
66 if(inverse
== !(i
&m
)) return split_radix_permutation(i
, m
, inverse
)*4 + 1;
67 else return split_radix_permutation(i
, m
, inverse
)*4 - 1;
70 av_cold
void ff_init_ff_cos_tabs(int index
)
72 #if !CONFIG_HARDCODED_TABLES
75 double freq
= 2*M_PI
/m
;
76 FFTSample
*tab
= ff_cos_tabs
[index
];
84 av_cold
int ff_fft_init(FFTContext
*s
, int nbits
, int inverse
)
88 if (nbits
< 2 || nbits
> 16)
93 s
->revtab
= av_malloc(n
* sizeof(uint16_t));
96 s
->tmp_buf
= av_malloc(n
* sizeof(FFTComplex
));
100 s
->fft_permutation
= FF_FFT_PERM_DEFAULT
;
102 s
->fft_permute
= ff_fft_permute_c
;
103 s
->fft_calc
= ff_fft_calc_c
;
105 s
->imdct_calc
= ff_imdct_calc_c
;
106 s
->imdct_half
= ff_imdct_half_c
;
107 s
->mdct_calc
= ff_mdct_calc_c
;
110 if (ARCH_ARM
) ff_fft_init_arm(s
);
111 if (HAVE_ALTIVEC
) ff_fft_init_altivec(s
);
112 if (HAVE_MMX
) ff_fft_init_mmx(s
);
114 for(j
=4; j
<=nbits
; j
++) {
115 ff_init_ff_cos_tabs(j
);
119 if (s
->fft_permutation
== FF_FFT_PERM_SWAP_LSBS
)
120 j
= (j
&~3) | ((j
>>1)&1) | ((j
<<1)&2);
121 s
->revtab
[-split_radix_permutation(i
, n
, s
->inverse
) & (n
-1)] = j
;
126 av_freep(&s
->revtab
);
127 av_freep(&s
->tmp_buf
);
131 static void ff_fft_permute_c(FFTContext
*s
, FFTComplex
*z
)
134 const uint16_t *revtab
= s
->revtab
;
136 /* TODO: handle split-radix permute in a more optimal way, probably in-place */
137 for(j
=0;j
<np
;j
++) s
->tmp_buf
[revtab
[j
]] = z
[j
];
138 memcpy(z
, s
->tmp_buf
, np
* sizeof(FFTComplex
));
141 av_cold
void ff_fft_end(FFTContext
*s
)
143 av_freep(&s
->revtab
);
144 av_freep(&s
->tmp_buf
);
147 #define sqrthalf (float)M_SQRT1_2
149 #define BF(x,y,a,b) {\
154 #define BUTTERFLIES(a0,a1,a2,a3) {\
156 BF(a2.re, a0.re, a0.re, t5);\
157 BF(a3.im, a1.im, a1.im, t3);\
159 BF(a3.re, a1.re, a1.re, t4);\
160 BF(a2.im, a0.im, a0.im, t6);\
163 // force loading all the inputs before storing any.
164 // this is slightly slower for small data, but avoids store->load aliasing
165 // for addresses separated by large powers of 2.
166 #define BUTTERFLIES_BIG(a0,a1,a2,a3) {\
167 FFTSample r0=a0.re, i0=a0.im, r1=a1.re, i1=a1.im;\
169 BF(a2.re, a0.re, r0, t5);\
170 BF(a3.im, a1.im, i1, t3);\
172 BF(a3.re, a1.re, r1, t4);\
173 BF(a2.im, a0.im, i0, t6);\
176 #define TRANSFORM(a0,a1,a2,a3,wre,wim) {\
177 t1 = a2.re * wre + a2.im * wim;\
178 t2 = a2.im * wre - a2.re * wim;\
179 t5 = a3.re * wre - a3.im * wim;\
180 t6 = a3.im * wre + a3.re * wim;\
181 BUTTERFLIES(a0,a1,a2,a3)\
184 #define TRANSFORM_ZERO(a0,a1,a2,a3) {\
189 BUTTERFLIES(a0,a1,a2,a3)\
192 /* z[0...8n-1], w[1...2n-1] */
194 static void name(FFTComplex *z, const FFTSample *wre, unsigned int n)\
196 FFTSample t1, t2, t3, t4, t5, t6;\
200 const FFTSample *wim = wre+o1;\
203 TRANSFORM_ZERO(z[0],z[o1],z[o2],z[o3]);\
204 TRANSFORM(z[1],z[o1+1],z[o2+1],z[o3+1],wre[1],wim[-1]);\
209 TRANSFORM(z[0],z[o1],z[o2],z[o3],wre[0],wim[0]);\
210 TRANSFORM(z[1],z[o1+1],z[o2+1],z[o3+1],wre[1],wim[-1]);\
216 #define BUTTERFLIES BUTTERFLIES_BIG
219 #define DECL_FFT(n,n2,n4)\
220 static void fft##n(FFTComplex *z)\
225 pass(z,ff_cos_##n,n4/2);\
228 static void fft4(FFTComplex
*z
)
230 FFTSample t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
;
232 BF(t3
, t1
, z
[0].re
, z
[1].re
);
233 BF(t8
, t6
, z
[3].re
, z
[2].re
);
234 BF(z
[2].re
, z
[0].re
, t1
, t6
);
235 BF(t4
, t2
, z
[0].im
, z
[1].im
);
236 BF(t7
, t5
, z
[2].im
, z
[3].im
);
237 BF(z
[3].im
, z
[1].im
, t4
, t8
);
238 BF(z
[3].re
, z
[1].re
, t3
, t7
);
239 BF(z
[2].im
, z
[0].im
, t2
, t5
);
242 static void fft8(FFTComplex
*z
)
244 FFTSample t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
;
248 BF(t1
, z
[5].re
, z
[4].re
, -z
[5].re
);
249 BF(t2
, z
[5].im
, z
[4].im
, -z
[5].im
);
250 BF(t3
, z
[7].re
, z
[6].re
, -z
[7].re
);
251 BF(t4
, z
[7].im
, z
[6].im
, -z
[7].im
);
254 BF(z
[4].re
, z
[0].re
, z
[0].re
, t1
);
255 BF(z
[4].im
, z
[0].im
, z
[0].im
, t2
);
256 BF(z
[6].re
, z
[2].re
, z
[2].re
, t7
);
257 BF(z
[6].im
, z
[2].im
, z
[2].im
, t8
);
259 TRANSFORM(z
[1],z
[3],z
[5],z
[7],sqrthalf
,sqrthalf
);
263 static void fft16(FFTComplex
*z
)
265 FFTSample t1
, t2
, t3
, t4
, t5
, t6
;
271 TRANSFORM_ZERO(z
[0],z
[4],z
[8],z
[12]);
272 TRANSFORM(z
[2],z
[6],z
[10],z
[14],sqrthalf
,sqrthalf
);
273 TRANSFORM(z
[1],z
[5],z
[9],z
[13],ff_cos_16
[1],ff_cos_16
[3]);
274 TRANSFORM(z
[3],z
[7],z
[11],z
[15],ff_cos_16
[3],ff_cos_16
[1]);
283 DECL_FFT(512,256,128)
285 #define pass pass_big
287 DECL_FFT(1024,512,256)
288 DECL_FFT(2048,1024,512)
289 DECL_FFT(4096,2048,1024)
290 DECL_FFT(8192,4096,2048)
291 DECL_FFT(16384,8192,4096)
292 DECL_FFT(32768,16384,8192)
293 DECL_FFT(65536,32768,16384)
295 static void (* const fft_dispatch
[])(FFTComplex
*) = {
296 fft4
, fft8
, fft16
, fft32
, fft64
, fft128
, fft256
, fft512
, fft1024
,
297 fft2048
, fft4096
, fft8192
, fft16384
, fft32768
, fft65536
,
300 static void ff_fft_calc_c(FFTContext
*s
, FFTComplex
*z
)
302 fft_dispatch
[s
->nbits
-2](z
);