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"
33 #include "fft-internal.h"
35 /* cos(2*pi*x/n) for 0<=x<=n/4, followed by its reverse */
36 #if !CONFIG_HARDCODED_TABLES
51 COSTABLE_CONST FFTSample
* const FFT_NAME(ff_cos_tabs
)[] = {
52 NULL
, NULL
, NULL
, NULL
,
59 FFT_NAME(ff_cos_1024
),
60 FFT_NAME(ff_cos_2048
),
61 FFT_NAME(ff_cos_4096
),
62 FFT_NAME(ff_cos_8192
),
63 FFT_NAME(ff_cos_16384
),
64 FFT_NAME(ff_cos_32768
),
65 FFT_NAME(ff_cos_65536
),
68 static void ff_fft_permute_c(FFTContext
*s
, FFTComplex
*z
);
69 static void ff_fft_calc_c(FFTContext
*s
, FFTComplex
*z
);
71 static int split_radix_permutation(int i
, int n
, int inverse
)
74 if(n
<= 2) return i
&1;
76 if(!(i
&m
)) return split_radix_permutation(i
, m
, inverse
)*2;
78 if(inverse
== !(i
&m
)) return split_radix_permutation(i
, m
, inverse
)*4 + 1;
79 else return split_radix_permutation(i
, m
, inverse
)*4 - 1;
82 av_cold
void ff_init_ff_cos_tabs(int index
)
84 #if !CONFIG_HARDCODED_TABLES
87 double freq
= 2*M_PI
/m
;
88 FFTSample
*tab
= FFT_NAME(ff_cos_tabs
)[index
];
90 tab
[i
] = FIX15(cos(i
*freq
));
96 static const int avx_tab
[] = {
97 0, 4, 1, 5, 8, 12, 9, 13, 2, 6, 3, 7, 10, 14, 11, 15
100 static int is_second_half_of_fft32(int i
, int n
)
105 return is_second_half_of_fft32(i
, n
/2);
107 return is_second_half_of_fft32(i
- n
/2, n
/4);
109 return is_second_half_of_fft32(i
- 3*n
/4, n
/4);
112 static av_cold
void fft_perm_avx(FFTContext
*s
)
115 int n
= 1 << s
->nbits
;
117 for (i
= 0; i
< n
; i
+= 16) {
119 if (is_second_half_of_fft32(i
, n
)) {
120 for (k
= 0; k
< 16; k
++)
121 s
->revtab
[-split_radix_permutation(i
+ k
, n
, s
->inverse
) & (n
- 1)] =
125 for (k
= 0; k
< 16; k
++) {
127 j
= (j
& ~7) | ((j
>> 1) & 3) | ((j
<< 2) & 4);
128 s
->revtab
[-split_radix_permutation(i
+ k
, n
, s
->inverse
) & (n
- 1)] = j
;
134 av_cold
int ff_fft_init(FFTContext
*s
, int nbits
, int inverse
)
138 if (nbits
< 2 || nbits
> 16)
143 s
->revtab
= av_malloc(n
* sizeof(uint16_t));
146 s
->tmp_buf
= av_malloc(n
* sizeof(FFTComplex
));
149 s
->inverse
= inverse
;
150 s
->fft_permutation
= FF_FFT_PERM_DEFAULT
;
152 s
->fft_permute
= ff_fft_permute_c
;
153 s
->fft_calc
= ff_fft_calc_c
;
155 s
->imdct_calc
= ff_imdct_calc_c
;
156 s
->imdct_half
= ff_imdct_half_c
;
157 s
->mdct_calc
= ff_mdct_calc_c
;
161 if (ARCH_ARM
) ff_fft_init_arm(s
);
162 if (HAVE_ALTIVEC
) ff_fft_init_altivec(s
);
163 if (ARCH_X86
) ff_fft_init_x86(s
);
164 if (CONFIG_MDCT
) s
->mdct_calcw
= s
->mdct_calc
;
166 if (CONFIG_MDCT
) s
->mdct_calcw
= ff_mdct_calcw_c
;
167 if (ARCH_ARM
) ff_fft_fixed_init_arm(s
);
170 for(j
=4; j
<=nbits
; j
++) {
171 ff_init_ff_cos_tabs(j
);
174 if (s
->fft_permutation
== FF_FFT_PERM_AVX
) {
179 if (s
->fft_permutation
== FF_FFT_PERM_SWAP_LSBS
)
180 j
= (j
&~3) | ((j
>>1)&1) | ((j
<<1)&2);
181 s
->revtab
[-split_radix_permutation(i
, n
, s
->inverse
) & (n
-1)] = j
;
187 av_freep(&s
->revtab
);
188 av_freep(&s
->tmp_buf
);
192 static void ff_fft_permute_c(FFTContext
*s
, FFTComplex
*z
)
195 const uint16_t *revtab
= s
->revtab
;
197 /* TODO: handle split-radix permute in a more optimal way, probably in-place */
198 for(j
=0;j
<np
;j
++) s
->tmp_buf
[revtab
[j
]] = z
[j
];
199 memcpy(z
, s
->tmp_buf
, np
* sizeof(FFTComplex
));
202 av_cold
void ff_fft_end(FFTContext
*s
)
204 av_freep(&s
->revtab
);
205 av_freep(&s
->tmp_buf
);
208 #define BUTTERFLIES(a0,a1,a2,a3) {\
210 BF(a2.re, a0.re, a0.re, t5);\
211 BF(a3.im, a1.im, a1.im, t3);\
213 BF(a3.re, a1.re, a1.re, t4);\
214 BF(a2.im, a0.im, a0.im, t6);\
217 // force loading all the inputs before storing any.
218 // this is slightly slower for small data, but avoids store->load aliasing
219 // for addresses separated by large powers of 2.
220 #define BUTTERFLIES_BIG(a0,a1,a2,a3) {\
221 FFTSample r0=a0.re, i0=a0.im, r1=a1.re, i1=a1.im;\
223 BF(a2.re, a0.re, r0, t5);\
224 BF(a3.im, a1.im, i1, t3);\
226 BF(a3.re, a1.re, r1, t4);\
227 BF(a2.im, a0.im, i0, t6);\
230 #define TRANSFORM(a0,a1,a2,a3,wre,wim) {\
231 CMUL(t1, t2, a2.re, a2.im, wre, -wim);\
232 CMUL(t5, t6, a3.re, a3.im, wre, wim);\
233 BUTTERFLIES(a0,a1,a2,a3)\
236 #define TRANSFORM_ZERO(a0,a1,a2,a3) {\
241 BUTTERFLIES(a0,a1,a2,a3)\
244 /* z[0...8n-1], w[1...2n-1] */
246 static void name(FFTComplex *z, const FFTSample *wre, unsigned int n)\
248 FFTDouble t1, t2, t3, t4, t5, t6;\
252 const FFTSample *wim = wre+o1;\
255 TRANSFORM_ZERO(z[0],z[o1],z[o2],z[o3]);\
256 TRANSFORM(z[1],z[o1+1],z[o2+1],z[o3+1],wre[1],wim[-1]);\
261 TRANSFORM(z[0],z[o1],z[o2],z[o3],wre[0],wim[0]);\
262 TRANSFORM(z[1],z[o1+1],z[o2+1],z[o3+1],wre[1],wim[-1]);\
268 #define BUTTERFLIES BUTTERFLIES_BIG
271 #define DECL_FFT(n,n2,n4)\
272 static void fft##n(FFTComplex *z)\
277 pass(z,FFT_NAME(ff_cos_##n),n4/2);\
280 static void fft4(FFTComplex
*z
)
282 FFTDouble t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
;
284 BF(t3
, t1
, z
[0].re
, z
[1].re
);
285 BF(t8
, t6
, z
[3].re
, z
[2].re
);
286 BF(z
[2].re
, z
[0].re
, t1
, t6
);
287 BF(t4
, t2
, z
[0].im
, z
[1].im
);
288 BF(t7
, t5
, z
[2].im
, z
[3].im
);
289 BF(z
[3].im
, z
[1].im
, t4
, t8
);
290 BF(z
[3].re
, z
[1].re
, t3
, t7
);
291 BF(z
[2].im
, z
[0].im
, t2
, t5
);
294 static void fft8(FFTComplex
*z
)
296 FFTDouble t1
, t2
, t3
, t4
, t5
, t6
;
300 BF(t1
, z
[5].re
, z
[4].re
, -z
[5].re
);
301 BF(t2
, z
[5].im
, z
[4].im
, -z
[5].im
);
302 BF(t5
, z
[7].re
, z
[6].re
, -z
[7].re
);
303 BF(t6
, z
[7].im
, z
[6].im
, -z
[7].im
);
305 BUTTERFLIES(z
[0],z
[2],z
[4],z
[6]);
306 TRANSFORM(z
[1],z
[3],z
[5],z
[7],sqrthalf
,sqrthalf
);
310 static void fft16(FFTComplex
*z
)
312 FFTDouble t1
, t2
, t3
, t4
, t5
, t6
;
313 FFTSample cos_16_1
= FFT_NAME(ff_cos_16
)[1];
314 FFTSample cos_16_3
= FFT_NAME(ff_cos_16
)[3];
320 TRANSFORM_ZERO(z
[0],z
[4],z
[8],z
[12]);
321 TRANSFORM(z
[2],z
[6],z
[10],z
[14],sqrthalf
,sqrthalf
);
322 TRANSFORM(z
[1],z
[5],z
[9],z
[13],cos_16_1
,cos_16_3
);
323 TRANSFORM(z
[3],z
[7],z
[11],z
[15],cos_16_3
,cos_16_1
);
332 DECL_FFT(512,256,128)
334 #define pass pass_big
336 DECL_FFT(1024,512,256)
337 DECL_FFT(2048,1024,512)
338 DECL_FFT(4096,2048,1024)
339 DECL_FFT(8192,4096,2048)
340 DECL_FFT(16384,8192,4096)
341 DECL_FFT(32768,16384,8192)
342 DECL_FFT(65536,32768,16384)
344 static void (* const fft_dispatch
[])(FFTComplex
*) = {
345 fft4
, fft8
, fft16
, fft32
, fft64
, fft128
, fft256
, fft512
, fft1024
,
346 fft2048
, fft4096
, fft8192
, fft16384
, fft32768
, fft65536
,
349 static void ff_fft_calc_c(FFTContext
*s
, FFTComplex
*z
)
351 fft_dispatch
[s
->nbits
-2](z
);