Prepared FFT test program (WIP)

1 files changed, 103 insertions, 0 deletions

diff --git a/fft.cpp b/fft.cpp
new file mode 100644
index 0000000..6b2f447
--- /dev/null
+++ b/fft.cpp
@@ -0,0 +1,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;
+}
+