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
|
#include <iostream>
#include <string>
#include <vector>
#include <complex>
#include <cmath>
const double PI = std::acos(-1);
const std::complex<double> i(0, 1);
using namespace std::complex_literals;
// (Slow) DFT
std::vector<std::complex<double>> dft(const std::vector<std::complex<double>> &v) {
std::vector<std::complex<double>> result(v.size(), 0);
const double N = v.size();
for (int k = 0; k < v.size(); k++) {
for (int j = 0; j < result.size(); j++) {
result[k] += v[j] * std::exp(-i * 2. * PI * double(k * j) / N);
}
}
return result;
}
// Cooley-Tukey
std::vector<std::complex<double>> fft1(const std::vector<std::complex<double>> &v) {
std::vector<std::complex<double>> result(v.size(), 0);
return result;
}
std::vector<double> freqs{ 2, 5, 11, 17, 29 }; // known freqs for testing
void generate(std::vector<std::complex<double>>& v) {
for (int i = 0; i < v.size(); i++) {
v[i] = 0.;
// sum several known sinusoids into v[]
for(int j = 0; j < freqs.size(); j++)
v[i] += sin(2 * M_PI * freqs[j] * i / v.size() );
}
}
// separate even/odd elements to lower/upper halves of array respectively.
// Due to Butterfly combinations, this turns out to be the simplest way
// to get the job done without clobbering the wrong elements.
void separate (complex<double>* a, int n) {
complex<double>* b = new complex<double>[n/2]; // get temp heap storage
for(int i=0; i<n/2; i++) // copy all odd elements to heap storage
b[i] = a[i*2+1];
for(int i=0; i<n/2; i++) // copy all even elements to lower-half of a[]
a[i] = a[i*2];
for(int i=0; i<n/2; i++) // copy all odd (from heap) to upper-half of a[]
a[i+n/2] = b[i];
delete[] b; // delete heap storage
}
// N must be a power-of-2, or bad things will happen.
// Currently no check for this condition.
//
// N input samples in X[] are FFT'd and results left in X[].
// Because of Nyquist theorem, N samples means
// only first N/2 FFT results in X[] are the answer.
// (upper half of X[] is a reflection with no new information).
void fft2 (complex<double>* X, int N) {
if(N < 2) {
// bottom of recursion.
// Do nothing here, because already X[0] = x[0]
} else {
separate(X,N); // all evens to lower half, all odds to upper half
fft2(X, N/2); // recurse even items
fft2(X+N/2, N/2); // recurse odd items
// combine results of two half recursions
for(int k=0; k<N/2; k++) {
complex<double> e = X[k ]; // even
complex<double> o = X[k+N/2]; // odd
// w is the "twiddle-factor"
complex<double> w = exp( complex<double>(0,-2.*M_PI*k/N) );
X[k ] = e + w * o;
X[k+N/2] = e - w * o;
}
}
}
int main(int argc, char* argv[]) {
std::vector<std::complex<double>> v(1024, 0);
generate(v);
std::vector<std::complex<double>> v_dft = dft(v);
std::vector<std::complex<double>> v_fft = fft1(v);
if (v_dft != v_fft) {
std::cout << "Error: Results differ!\n";
return 1;
}
return 0;
}
|
