From b9acf8c4700337a12758e5c784af25965a643d2a Mon Sep 17 00:00:00 2001 From: Anton Te Date: Thu, 19 Jul 2018 00:17:35 -0700 Subject: [PATCH 1/3] Minor clean up for cooley tukey in C++ - rename c64 to complex - make pi a constant instead of function - fix complilation warnings about comparing signed int with unsigned int --- contents/cooley_tukey/code/c++/fft.cpp | 31 ++++++++++++-------------- 1 file changed, 14 insertions(+), 17 deletions(-) diff --git a/contents/cooley_tukey/code/c++/fft.cpp b/contents/cooley_tukey/code/c++/fft.cpp index 56e2554bb..b58e3f5ab 100644 --- a/contents/cooley_tukey/code/c++/fft.cpp +++ b/contents/cooley_tukey/code/c++/fft.cpp @@ -17,11 +17,8 @@ using std::swap; using std::size_t; -using c64 = std::complex; -template -constexpr T pi() { - return 3.14159265358979323846264338327950288419716; -} +using complex = std::complex; +static const double pi = 3.14159265358979323846264338327950288419716; // `cooley_tukey` does the cooley-tukey algorithm, recursively template @@ -30,12 +27,12 @@ void cooley_tukey(Iter first, Iter last) { if (size >= 2) { // split the range, with even indices going in the first half, // and odd indices going in the last half. - auto temp = std::vector(size / 2); - for (size_t i = 0; i < size / 2; ++i) { + auto temp = std::vector(size / 2); + for (int i = 0; i < size / 2; ++i) { temp[i] = first[i * 2 + 1]; first[i] = first[i * 2]; } - for (size_t i = 0; i < size / 2; ++i) { + for (int i = 0; i < size / 2; ++i) { first[i + size / 2] = temp[i]; } @@ -45,8 +42,8 @@ void cooley_tukey(Iter first, Iter last) { cooley_tukey(split, last); // now combine each of those halves with the butterflies - for (size_t k = 0; k < size / 2; ++k) { - auto w = std::exp(c64(0, -2.0 * pi() * k / size)); + for (int k = 0; k < size / 2; ++k) { + auto w = std::exp(complex(0, -2.0 * pi * k / size)); auto& bottom = first[k]; auto& top = first[k + size / 2]; @@ -83,11 +80,11 @@ void iterative_cooley_tukey(Iter first, Iter last) { // perform the butterfly on the range auto size = last - first; - for (size_t stride = 2; stride <= size; stride *= 2) { - auto w = exp(c64(0, -2.0 * pi() / stride)); - for (size_t j = 0; j < size; j += stride) { - auto v = c64(1.0); - for (size_t k = 0; k < stride / 2; k++) { + for (int stride = 2; stride <= size; stride *= 2) { + auto w = exp(complex(0, -2.0 * pi / stride)); + for (int j = 0; j < size; j += stride) { + auto v = complex(1.0); + for (int k = 0; k < stride / 2; k++) { first[k + j + stride / 2] = first[k + j] - v * first[k + j + stride / 2]; first[k + j] -= (first[k + j + stride / 2] - first[k + j]); @@ -103,7 +100,7 @@ int main() { std::mt19937 rng(random_device()); std::uniform_real_distribution distribution(0.0, 1.0); - std::array initial; + std::array initial; std::generate( begin(initial), end(initial), [&] { return distribution(rng); }); @@ -117,7 +114,7 @@ int main() { // Check if the arrays are approximately equivalent std::cout << std::right << std::setw(16) << "idx" << std::setw(16) << "rec" << std::setw(16) << "it" << std::setw(16) << "subtracted" << '\n'; - for (int i = 0; i < initial.size(); ++i) { + for (size_t i = 0; i < initial.size(); ++i) { auto rec = recursive[i]; auto it = iterative[i]; std::cout << std::setw(16) << i << std::setw(16) << std::abs(rec) From 3227558a72dcee31ac9ef0b180f135a3cdec5b50 Mon Sep 17 00:00:00 2001 From: Anton Te Date: Sat, 4 Aug 2018 21:48:36 -0700 Subject: [PATCH 2/3] Implement dft --- contents/cooley_tukey/code/c++/fft.cpp | 13 +++++++++++++ contents/cooley_tukey/cooley_tukey.md | 2 +- 2 files changed, 14 insertions(+), 1 deletion(-) diff --git a/contents/cooley_tukey/code/c++/fft.cpp b/contents/cooley_tukey/code/c++/fft.cpp index b58e3f5ab..ade544442 100644 --- a/contents/cooley_tukey/code/c++/fft.cpp +++ b/contents/cooley_tukey/code/c++/fft.cpp @@ -20,6 +20,19 @@ using std::size_t; using complex = std::complex; static const double pi = 3.14159265358979323846264338327950288419716; +template +void dft(Iter X, Iter last) { + const auto N = last - X; + std::vector tmp(N); + for (auto i = 0; i < N; ++i) { + for (auto j = 0; j < N; ++j) { + using namespace std::literals::complex_literals; + tmp[i] += X[j] * exp(-2.0 * M_PI * i * j / N * 1i); + } + } + std::copy(std::begin(tmp), std::end(tmp), X); +} + // `cooley_tukey` does the cooley-tukey algorithm, recursively template void cooley_tukey(Iter first, Iter last) { diff --git a/contents/cooley_tukey/cooley_tukey.md b/contents/cooley_tukey/cooley_tukey.md index 41b00170a..91d24eaff 100644 --- a/contents/cooley_tukey/cooley_tukey.md +++ b/contents/cooley_tukey/cooley_tukey.md @@ -119,7 +119,7 @@ In the end, the code looks like: {% sample lang="c" %} [import:20-39, lang:"c_cpp"](code/c/fft.c) {% sample lang="cpp" %} -[import:27-57, lang:"c_cpp"](code/c++/fft.cpp) +[import:36-67, lang:"c_cpp"](code/c++/fft.cpp) {% sample lang="hs" %} [import:6-19, lang:"haskell"](code/haskell/fft.hs) {% sample lang="py" %} From fe358f8f444744edaa46afec6a43e00f523cf4b0 Mon Sep 17 00:00:00 2001 From: Anton Te Date: Mon, 6 Aug 2018 21:16:39 -0700 Subject: [PATCH 3/3] Shorter version of formula, one line less --- contents/cooley_tukey/code/c++/fft.cpp | 3 +-- contents/cooley_tukey/cooley_tukey.md | 2 +- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/contents/cooley_tukey/code/c++/fft.cpp b/contents/cooley_tukey/code/c++/fft.cpp index ade544442..5d4772f10 100644 --- a/contents/cooley_tukey/code/c++/fft.cpp +++ b/contents/cooley_tukey/code/c++/fft.cpp @@ -26,8 +26,7 @@ void dft(Iter X, Iter last) { std::vector tmp(N); for (auto i = 0; i < N; ++i) { for (auto j = 0; j < N; ++j) { - using namespace std::literals::complex_literals; - tmp[i] += X[j] * exp(-2.0 * M_PI * i * j / N * 1i); + tmp[i] += X[j] * exp(complex(0, -2.0 * M_PI * i * j / N)); } } std::copy(std::begin(tmp), std::end(tmp), X); diff --git a/contents/cooley_tukey/cooley_tukey.md b/contents/cooley_tukey/cooley_tukey.md index 91d24eaff..e414e8c16 100644 --- a/contents/cooley_tukey/cooley_tukey.md +++ b/contents/cooley_tukey/cooley_tukey.md @@ -119,7 +119,7 @@ In the end, the code looks like: {% sample lang="c" %} [import:20-39, lang:"c_cpp"](code/c/fft.c) {% sample lang="cpp" %} -[import:36-67, lang:"c_cpp"](code/c++/fft.cpp) +[import:35-66, lang:"c_cpp"](code/c++/fft.cpp) {% sample lang="hs" %} [import:6-19, lang:"haskell"](code/haskell/fft.hs) {% sample lang="py" %}