diff options
-rw-r--r-- | .gitignore | 1 | ||||
-rw-r--r-- | Makefile | 13 | ||||
-rw-r--r-- | debian/control | 1 | ||||
-rw-r--r-- | debian/rules | 3 | ||||
-rw-r--r-- | fft.cpp | 103 |
5 files changed, 121 insertions, 0 deletions
diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..56d2382 --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +*.o fft diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..f7d22c2 --- /dev/null +++ b/Makefile @@ -0,0 +1,13 @@ +CXX=clang++-7 + +all: fft + +fft: fft.cpp + $(CXX) -Wall -std=c++17 -o $@ $^ + +install: + +clean: + rm -f fft + +.PHONY: clean diff --git a/debian/control b/debian/control new file mode 100644 index 0000000..67ff324 --- /dev/null +++ b/debian/control @@ -0,0 +1 @@ +Depends: g++, clang diff --git a/debian/rules b/debian/rules new file mode 100644 index 0000000..cbe925d --- /dev/null +++ b/debian/rules @@ -0,0 +1,3 @@ +#!/usr/bin/make -f +%: + dh $@ @@ -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; +} + |