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 FFmpeg.
9 * FFmpeg 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 * FFmpeg 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 FFmpeg; 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"
36 #include "fft_table.h"
37 #else /* FFT_FIXED_32 */
39 /* cos(2*pi*x/n) for 0<=x<=n/4, followed by its reverse */
40 #if !CONFIG_HARDCODED_TABLES
55 COSTABLE_CONST FFTSample
* const FFT_NAME(ff_cos_tabs
)[] = {
56 NULL
, NULL
, NULL
, NULL
,
63 FFT_NAME(ff_cos_1024
),
64 FFT_NAME(ff_cos_2048
),
65 FFT_NAME(ff_cos_4096
),
66 FFT_NAME(ff_cos_8192
),
67 FFT_NAME(ff_cos_16384
),
68 FFT_NAME(ff_cos_32768
),
69 FFT_NAME(ff_cos_65536
),
72 #endif /* FFT_FIXED_32 */
74 static void fft_permute_c(FFTContext
*s
, FFTComplex
*z
);
75 static void fft_calc_c(FFTContext
*s
, FFTComplex
*z
);
77 static int split_radix_permutation(int i
, int n
, int inverse
)
80 if(n
<= 2) return i
&1;
82 if(!(i
&m
)) return split_radix_permutation(i
, m
, inverse
)*2;
84 if(inverse
== !(i
&m
)) return split_radix_permutation(i
, m
, inverse
)*4 + 1;
85 else return split_radix_permutation(i
, m
, inverse
)*4 - 1;
88 av_cold
void ff_init_ff_cos_tabs(int index
)
90 #if (!CONFIG_HARDCODED_TABLES) && (!FFT_FIXED_32)
93 double freq
= 2*M_PI
/m
;
94 FFTSample
*tab
= FFT_NAME(ff_cos_tabs
)[index
];
96 tab
[i
] = FIX15(cos(i
*freq
));
102 static const int avx_tab
[] = {
103 0, 4, 1, 5, 8, 12, 9, 13, 2, 6, 3, 7, 10, 14, 11, 15
106 static int is_second_half_of_fft32(int i
, int n
)
111 return is_second_half_of_fft32(i
, n
/2);
113 return is_second_half_of_fft32(i
- n
/2, n
/4);
115 return is_second_half_of_fft32(i
- 3*n
/4, n
/4);
118 static av_cold
void fft_perm_avx(FFTContext
*s
)
121 int n
= 1 << s
->nbits
;
123 for (i
= 0; i
< n
; i
+= 16) {
125 if (is_second_half_of_fft32(i
, n
)) {
126 for (k
= 0; k
< 16; k
++)
127 s
->revtab
[-split_radix_permutation(i
+ k
, n
, s
->inverse
) & (n
- 1)] =
131 for (k
= 0; k
< 16; k
++) {
133 j
= (j
& ~7) | ((j
>> 1) & 3) | ((j
<< 2) & 4);
134 s
->revtab
[-split_radix_permutation(i
+ k
, n
, s
->inverse
) & (n
- 1)] = j
;
140 av_cold
int ff_fft_init(FFTContext
*s
, int nbits
, int inverse
)
144 if (nbits
< 2 || nbits
> 16)
149 s
->revtab
= av_malloc(n
* sizeof(uint16_t));
152 s
->tmp_buf
= av_malloc(n
* sizeof(FFTComplex
));
155 s
->inverse
= inverse
;
156 s
->fft_permutation
= FF_FFT_PERM_DEFAULT
;
158 s
->fft_permute
= fft_permute_c
;
159 s
->fft_calc
= fft_calc_c
;
161 s
->imdct_calc
= ff_imdct_calc_c
;
162 s
->imdct_half
= ff_imdct_half_c
;
163 s
->mdct_calc
= ff_mdct_calc_c
;
169 ff_fft_lut_init(ff_fft_offsets_lut
, 0, 1 << 16, &n
);
171 #else /* FFT_FIXED_32 */
173 if (ARCH_AARCH64
) ff_fft_init_aarch64(s
);
174 if (ARCH_ARM
) ff_fft_init_arm(s
);
175 if (ARCH_PPC
) ff_fft_init_ppc(s
);
176 if (ARCH_X86
) ff_fft_init_x86(s
);
177 if (CONFIG_MDCT
) s
->mdct_calcw
= s
->mdct_calc
;
178 if (HAVE_MIPSFPU
) ff_fft_init_mips(s
);
180 if (CONFIG_MDCT
) s
->mdct_calcw
= ff_mdct_calcw_c
;
181 if (ARCH_ARM
) ff_fft_fixed_init_arm(s
);
183 for(j
=4; j
<=nbits
; j
++) {
184 ff_init_ff_cos_tabs(j
);
186 #endif /* FFT_FIXED_32 */
189 if (s
->fft_permutation
== FF_FFT_PERM_AVX
) {
194 if (s
->fft_permutation
== FF_FFT_PERM_SWAP_LSBS
)
195 j
= (j
&~3) | ((j
>>1)&1) | ((j
<<1)&2);
196 s
->revtab
[-split_radix_permutation(i
, n
, s
->inverse
) & (n
-1)] = j
;
202 av_freep(&s
->revtab
);
203 av_freep(&s
->tmp_buf
);
207 static void fft_permute_c(FFTContext
*s
, FFTComplex
*z
)
210 const uint16_t *revtab
= s
->revtab
;
212 /* TODO: handle split-radix permute in a more optimal way, probably in-place */
213 for(j
=0;j
<np
;j
++) s
->tmp_buf
[revtab
[j
]] = z
[j
];
214 memcpy(z
, s
->tmp_buf
, np
* sizeof(FFTComplex
));
217 av_cold
void ff_fft_end(FFTContext
*s
)
219 av_freep(&s
->revtab
);
220 av_freep(&s
->tmp_buf
);
225 static void fft_calc_c(FFTContext
*s
, FFTComplex
*z
) {
227 int nbits
, i
, n
, num_transforms
, offset
, step
;
229 FFTSample tmp1
, tmp2
, tmp3
, tmp4
, tmp5
, tmp6
, tmp7
, tmp8
;
231 const int fft_size
= (1 << s
->nbits
);
234 num_transforms
= (0x2aab >> (16 - s
->nbits
)) | 1;
236 for (n
=0; n
<num_transforms
; n
++){
237 offset
= ff_fft_offsets_lut
[n
] << 2;
240 tmp1
= tmpz
[0].re
+ tmpz
[1].re
;
241 tmp5
= tmpz
[2].re
+ tmpz
[3].re
;
242 tmp2
= tmpz
[0].im
+ tmpz
[1].im
;
243 tmp6
= tmpz
[2].im
+ tmpz
[3].im
;
244 tmp3
= tmpz
[0].re
- tmpz
[1].re
;
245 tmp8
= tmpz
[2].im
- tmpz
[3].im
;
246 tmp4
= tmpz
[0].im
- tmpz
[1].im
;
247 tmp7
= tmpz
[2].re
- tmpz
[3].re
;
249 tmpz
[0].re
= tmp1
+ tmp5
;
250 tmpz
[2].re
= tmp1
- tmp5
;
251 tmpz
[0].im
= tmp2
+ tmp6
;
252 tmpz
[2].im
= tmp2
- tmp6
;
253 tmpz
[1].re
= tmp3
+ tmp8
;
254 tmpz
[3].re
= tmp3
- tmp8
;
255 tmpz
[1].im
= tmp4
- tmp7
;
256 tmpz
[3].im
= tmp4
+ tmp7
;
262 num_transforms
= (num_transforms
>> 1) | 1;
264 for (n
=0; n
<num_transforms
; n
++){
265 offset
= ff_fft_offsets_lut
[n
] << 3;
268 tmp1
= tmpz
[4].re
+ tmpz
[5].re
;
269 tmp3
= tmpz
[6].re
+ tmpz
[7].re
;
270 tmp2
= tmpz
[4].im
+ tmpz
[5].im
;
271 tmp4
= tmpz
[6].im
+ tmpz
[7].im
;
277 tmp1
= tmpz
[4].re
- tmpz
[5].re
;
278 tmp2
= tmpz
[4].im
- tmpz
[5].im
;
279 tmp3
= tmpz
[6].re
- tmpz
[7].re
;
280 tmp4
= tmpz
[6].im
- tmpz
[7].im
;
282 tmpz
[4].re
= tmpz
[0].re
- tmp5
;
283 tmpz
[0].re
= tmpz
[0].re
+ tmp5
;
284 tmpz
[4].im
= tmpz
[0].im
- tmp6
;
285 tmpz
[0].im
= tmpz
[0].im
+ tmp6
;
286 tmpz
[6].re
= tmpz
[2].re
- tmp8
;
287 tmpz
[2].re
= tmpz
[2].re
+ tmp8
;
288 tmpz
[6].im
= tmpz
[2].im
+ tmp7
;
289 tmpz
[2].im
= tmpz
[2].im
- tmp7
;
291 accu
= (int64_t)Q31(M_SQRT1_2
)*(tmp1
+ tmp2
);
292 tmp5
= (int32_t)((accu
+ 0x40000000) >> 31);
293 accu
= (int64_t)Q31(M_SQRT1_2
)*(tmp3
- tmp4
);
294 tmp7
= (int32_t)((accu
+ 0x40000000) >> 31);
295 accu
= (int64_t)Q31(M_SQRT1_2
)*(tmp2
- tmp1
);
296 tmp6
= (int32_t)((accu
+ 0x40000000) >> 31);
297 accu
= (int64_t)Q31(M_SQRT1_2
)*(tmp3
+ tmp4
);
298 tmp8
= (int32_t)((accu
+ 0x40000000) >> 31);
304 tmpz
[5].re
= tmpz
[1].re
- tmp1
;
305 tmpz
[1].re
= tmpz
[1].re
+ tmp1
;
306 tmpz
[5].im
= tmpz
[1].im
- tmp2
;
307 tmpz
[1].im
= tmpz
[1].im
+ tmp2
;
308 tmpz
[7].re
= tmpz
[3].re
- tmp4
;
309 tmpz
[3].re
= tmpz
[3].re
+ tmp4
;
310 tmpz
[7].im
= tmpz
[3].im
+ tmp3
;
311 tmpz
[3].im
= tmpz
[3].im
- tmp3
;
314 step
= 1 << ((MAX_LOG2_NFFT
-4) - 4);
317 for (nbits
=4; nbits
<=s
->nbits
; nbits
++){
320 num_transforms
= (num_transforms
>> 1) | 1;
322 for (n
=0; n
<num_transforms
; n
++){
323 const FFTSample
*w_re_ptr
= ff_w_tab_sr
+ step
;
324 const FFTSample
*w_im_ptr
= ff_w_tab_sr
+ MAX_FFT_SIZE
/(4*16) - step
;
325 offset
= ff_fft_offsets_lut
[n
] << nbits
;
328 tmp5
= tmpz
[ n2
].re
+ tmpz
[n34
].re
;
329 tmp1
= tmpz
[ n2
].re
- tmpz
[n34
].re
;
330 tmp6
= tmpz
[ n2
].im
+ tmpz
[n34
].im
;
331 tmp2
= tmpz
[ n2
].im
- tmpz
[n34
].im
;
333 tmpz
[ n2
].re
= tmpz
[ 0].re
- tmp5
;
334 tmpz
[ 0].re
= tmpz
[ 0].re
+ tmp5
;
335 tmpz
[ n2
].im
= tmpz
[ 0].im
- tmp6
;
336 tmpz
[ 0].im
= tmpz
[ 0].im
+ tmp6
;
337 tmpz
[n34
].re
= tmpz
[n4
].re
- tmp2
;
338 tmpz
[ n4
].re
= tmpz
[n4
].re
+ tmp2
;
339 tmpz
[n34
].im
= tmpz
[n4
].im
+ tmp1
;
340 tmpz
[ n4
].im
= tmpz
[n4
].im
- tmp1
;
342 for (i
=1; i
<n4
; i
++){
343 FFTSample w_re
= w_re_ptr
[0];
344 FFTSample w_im
= w_im_ptr
[0];
345 accu
= (int64_t)w_re
*tmpz
[ n2
+i
].re
;
346 accu
+= (int64_t)w_im
*tmpz
[ n2
+i
].im
;
347 tmp1
= (int32_t)((accu
+ 0x40000000) >> 31);
348 accu
= (int64_t)w_re
*tmpz
[ n2
+i
].im
;
349 accu
-= (int64_t)w_im
*tmpz
[ n2
+i
].re
;
350 tmp2
= (int32_t)((accu
+ 0x40000000) >> 31);
351 accu
= (int64_t)w_re
*tmpz
[n34
+i
].re
;
352 accu
-= (int64_t)w_im
*tmpz
[n34
+i
].im
;
353 tmp3
= (int32_t)((accu
+ 0x40000000) >> 31);
354 accu
= (int64_t)w_re
*tmpz
[n34
+i
].im
;
355 accu
+= (int64_t)w_im
*tmpz
[n34
+i
].re
;
356 tmp4
= (int32_t)((accu
+ 0x40000000) >> 31);
363 tmpz
[ n2
+i
].re
= tmpz
[ i
].re
- tmp5
;
364 tmpz
[ i
].re
= tmpz
[ i
].re
+ tmp5
;
365 tmpz
[ n2
+i
].im
= tmpz
[ i
].im
- tmp6
;
366 tmpz
[ i
].im
= tmpz
[ i
].im
+ tmp6
;
367 tmpz
[n34
+i
].re
= tmpz
[n4
+i
].re
- tmp2
;
368 tmpz
[ n4
+i
].re
= tmpz
[n4
+i
].re
+ tmp2
;
369 tmpz
[n34
+i
].im
= tmpz
[n4
+i
].im
+ tmp1
;
370 tmpz
[ n4
+i
].im
= tmpz
[n4
+i
].im
- tmp1
;
381 #else /* FFT_FIXED_32 */
383 #define BUTTERFLIES(a0,a1,a2,a3) {\
385 BF(a2.re, a0.re, a0.re, t5);\
386 BF(a3.im, a1.im, a1.im, t3);\
388 BF(a3.re, a1.re, a1.re, t4);\
389 BF(a2.im, a0.im, a0.im, t6);\
392 // force loading all the inputs before storing any.
393 // this is slightly slower for small data, but avoids store->load aliasing
394 // for addresses separated by large powers of 2.
395 #define BUTTERFLIES_BIG(a0,a1,a2,a3) {\
396 FFTSample r0=a0.re, i0=a0.im, r1=a1.re, i1=a1.im;\
398 BF(a2.re, a0.re, r0, t5);\
399 BF(a3.im, a1.im, i1, t3);\
401 BF(a3.re, a1.re, r1, t4);\
402 BF(a2.im, a0.im, i0, t6);\
405 #define TRANSFORM(a0,a1,a2,a3,wre,wim) {\
406 CMUL(t1, t2, a2.re, a2.im, wre, -wim);\
407 CMUL(t5, t6, a3.re, a3.im, wre, wim);\
408 BUTTERFLIES(a0,a1,a2,a3)\
411 #define TRANSFORM_ZERO(a0,a1,a2,a3) {\
416 BUTTERFLIES(a0,a1,a2,a3)\
419 /* z[0...8n-1], w[1...2n-1] */
421 static void name(FFTComplex *z, const FFTSample *wre, unsigned int n)\
423 FFTDouble t1, t2, t3, t4, t5, t6;\
427 const FFTSample *wim = wre+o1;\
430 TRANSFORM_ZERO(z[0],z[o1],z[o2],z[o3]);\
431 TRANSFORM(z[1],z[o1+1],z[o2+1],z[o3+1],wre[1],wim[-1]);\
436 TRANSFORM(z[0],z[o1],z[o2],z[o3],wre[0],wim[0]);\
437 TRANSFORM(z[1],z[o1+1],z[o2+1],z[o3+1],wre[1],wim[-1]);\
443 #define BUTTERFLIES BUTTERFLIES_BIG
446 #define DECL_FFT(n,n2,n4)\
447 static void fft##n(FFTComplex *z)\
452 pass(z,FFT_NAME(ff_cos_##n),n4/2);\
455 static void fft4(FFTComplex
*z
)
457 FFTDouble t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
;
459 BF(t3
, t1
, z
[0].re
, z
[1].re
);
460 BF(t8
, t6
, z
[3].re
, z
[2].re
);
461 BF(z
[2].re
, z
[0].re
, t1
, t6
);
462 BF(t4
, t2
, z
[0].im
, z
[1].im
);
463 BF(t7
, t5
, z
[2].im
, z
[3].im
);
464 BF(z
[3].im
, z
[1].im
, t4
, t8
);
465 BF(z
[3].re
, z
[1].re
, t3
, t7
);
466 BF(z
[2].im
, z
[0].im
, t2
, t5
);
469 static void fft8(FFTComplex
*z
)
471 FFTDouble t1
, t2
, t3
, t4
, t5
, t6
;
475 BF(t1
, z
[5].re
, z
[4].re
, -z
[5].re
);
476 BF(t2
, z
[5].im
, z
[4].im
, -z
[5].im
);
477 BF(t5
, z
[7].re
, z
[6].re
, -z
[7].re
);
478 BF(t6
, z
[7].im
, z
[6].im
, -z
[7].im
);
480 BUTTERFLIES(z
[0],z
[2],z
[4],z
[6]);
481 TRANSFORM(z
[1],z
[3],z
[5],z
[7],sqrthalf
,sqrthalf
);
485 static void fft16(FFTComplex
*z
)
487 FFTDouble t1
, t2
, t3
, t4
, t5
, t6
;
488 FFTSample cos_16_1
= FFT_NAME(ff_cos_16
)[1];
489 FFTSample cos_16_3
= FFT_NAME(ff_cos_16
)[3];
495 TRANSFORM_ZERO(z
[0],z
[4],z
[8],z
[12]);
496 TRANSFORM(z
[2],z
[6],z
[10],z
[14],sqrthalf
,sqrthalf
);
497 TRANSFORM(z
[1],z
[5],z
[9],z
[13],cos_16_1
,cos_16_3
);
498 TRANSFORM(z
[3],z
[7],z
[11],z
[15],cos_16_3
,cos_16_1
);
507 DECL_FFT(512,256,128)
509 #define pass pass_big
511 DECL_FFT(1024,512,256)
512 DECL_FFT(2048,1024,512)
513 DECL_FFT(4096,2048,1024)
514 DECL_FFT(8192,4096,2048)
515 DECL_FFT(16384,8192,4096)
516 DECL_FFT(32768,16384,8192)
517 DECL_FFT(65536,32768,16384)
519 static void (* const fft_dispatch
[])(FFTComplex
*) = {
520 fft4
, fft8
, fft16
, fft32
, fft64
, fft128
, fft256
, fft512
, fft1024
,
521 fft2048
, fft4096
, fft8192
, fft16384
, fft32768
, fft65536
,
524 static void fft_calc_c(FFTContext
*s
, FFTComplex
*z
)
526 fft_dispatch
[s
->nbits
-2](z
);
528 #endif /* FFT_FIXED_32 */