/
cp_dit_bi_G_L_S_P.c
115 lines (97 loc) · 2.45 KB
/
cp_dit_bi_G_L_S_P.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
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
/*
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 "sr_sched_lut.h"
#include "cp_input_lut.h"
#include "cp_bf.h"
void
fft(unsigned N, complex double * restrict out,
const complex double * restrict in)
{
unsigned log2_N = 31 - clz32(N);
/* stage 0 */
for (unsigned j = 0, i = 0; i < N; i += 2)
{
unsigned i0 = cp_input_lut[i / 2];
unsigned i1 = i0 ^ (N >> 1);
if (sr_sched_off[j] * 2 > i)
{
out[i ] = in[i0];
out[i + 1] = in[i1];
}
else
{
complex double a = in[i0];
complex double b = in[i1];
out[i ] = a + b;
out[i + 1] = a - b;
j++;
}
}
/* stage 1 */
if (log2_N > 1)
{
unsigned c = sr_sched_cnt[1];
for (unsigned i = 0; i < c; i++)
{
unsigned o = sr_sched_off[i] << 2;
cp_dit_bf_0(1, out + o);
}
}
/* stage 2 */
if (log2_N > 2)
{
unsigned c = sr_sched_cnt[2];
for (unsigned i = 0; i < c; i++)
{
unsigned o = sr_sched_off[i] << 3;
cp_dit_bf_pi4(2, out + o + 1);
cp_dit_bf_0(2, out + o);
}
}
/* stage >= 3 */
if (log2_N > 3)
for (unsigned s = 3; s < log2_N; s++)
{
unsigned c = sr_sched_cnt[s];
unsigned os = 1 << (s - 1);
for (unsigned i = 0; i < c; i++)
{
unsigned o = sr_sched_off[i] << (s + 1);
cp_dit_bf_0(os, out + o);
cp_dit_bf_pi4(os, out + o + os / 2);
}
for (unsigned k = 1; k < os / 2; k++)
{
unsigned t = log2_N - s - 1;
complex double w = get_twiddle(k << t, log2_N);
for (unsigned i = 0; i < c; i++)
{
unsigned o = sr_sched_off[i] << (s + 1);
cp_dit_bf(os, out + o + k, w);
cp_dit_bf(os, out + o + os - k, -I * conj(w));
}
}
}
}
void fft_init(unsigned N)
{
twiddles_init(N / 8, N);
cp_input_lut_init(N);
sr_sched_lut_init(N);
}
void fft_cleanup()
{
sr_sched_lut_cleanup();
cp_input_lut_cleanup();
twiddles_cleanup();
}