/
ct_dif_bi_N.c
69 lines (55 loc) · 1.46 KB
/
ct_dif_bi_N.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
/*
Written in 2019 by Alexandre Becoulet <alexandre.becoulet@free.fr>
To the extent possible under law, the author(s) have dedicated all
copyright and related and neighboring rights to this software to
the public domain worldwide. This software is distributed without
any warranty.
You should have received a copy of the CC0 Public Domain
Dedication along with this software. If not, see
http://creativecommons.org/publicdomain/zero/1.0/
*/
#include "fft.h"
#include "twiddles.h"
#include "ct_bf.h"
#include "bit_reverse.h"
/* A breadth-first implementation of the
radix 2 Cooley Tukey DIF FFT.
See ct_dit_bi_N.c for the DIT variant.
*/
static inline void
ct_fft_bi(complex double * restrict inout,
int log2_N)
{
unsigned N = 1 << log2_N;
/* stages 0 to log2_N - 1 */
for (int j = log2_N - 1; j >= 0; j--)
{
unsigned s = log2_N - j - 1;
unsigned l = 1 << j;
for (unsigned i = 0; i < (1 << s); i++)
{
complex double *t = inout + (i << (j + 1));
for (unsigned k = 0; k < l; k++)
ct_dif_bf2(l, t + k, get_twiddle(k << s, log2_N));
}
}
BIT_REVERSE_SWAP_LOOP(i, j, N, {
double complex t = inout[i];
inout[i] = inout[j];
inout[j] = t;
});
}
void
fft(unsigned N, complex double *inout)
{
unsigned log2_N = 31 - __builtin_clz(N);
ct_fft_bi(inout, log2_N);
}
void fft_init(unsigned N)
{
twiddles_init(N / 2, N);
}
void fft_cleanup()
{
twiddles_cleanup();
}