diff options
author | Roland Stigge <stigge@antcom.de> | 2019-01-22 22:03:06 +0100 |
---|---|---|
committer | Roland Stigge <stigge@antcom.de> | 2019-01-22 22:03:06 +0100 |
commit | 9e238820cb4df45219cf122206aa0ae153e81c60 (patch) | |
tree | c1a794e120cb1c75c9afaa95deebcb51ab04f517 | |
parent | 9cc4c3cc65df53f6c835d7473aaf3955345cd88b (diff) |
Prepared chrono
-rw-r--r-- | Makefile | 2 | ||||
-rw-r--r-- | fft.cpp | 28 |
2 files changed, 16 insertions, 14 deletions
@@ -3,7 +3,7 @@ CXX=clang++-7 all: fft fft: fft.cpp - $(CXX) -Wall -std=c++17 -o $@ $^ + $(CXX) -Wall -O2 -std=c++17 -o $@ $^ install: @@ -3,6 +3,7 @@ #include <vector> #include <complex> #include <cmath> +#include <chrono> const double PI = std::acos(-1); const std::complex<double> i(0, 1); @@ -24,13 +25,6 @@ std::vector<std::complex<double>> dft(const std::vector<std::complex<double>> &v 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) { @@ -46,8 +40,8 @@ void generate(std::vector<std::complex<double>>& v) { // 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 +void separate (std::complex<double>* a, int n) { + std::complex<double>* b = new std::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[] @@ -64,7 +58,7 @@ void separate (complex<double>* a, int n) { // 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) { +void fft2 (std::complex<double>* X, int N) { if(N < 2) { // bottom of recursion. // Do nothing here, because already X[0] = x[0] @@ -74,16 +68,24 @@ void fft2 (complex<double>* X, int N) { 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 + std::complex<double> e = X[k ]; // even + std::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) ); + std::complex<double> w = exp( std::complex<double>(0,-2.*M_PI*k/N) ); X[k ] = e + w * o; X[k+N/2] = e - w * o; } } } +// Cooley-Tukey +// TODO: Use fft2 via iterators! +std::vector<std::complex<double>> fft1(const std::vector<std::complex<double>> &v) { + std::vector<std::complex<double>> result(v.size(), 0); + + return result; +} + int main(int argc, char* argv[]) { std::vector<std::complex<double>> v(1024, 0); |