Merge pull request #276 from antonte/cooley_tukey_cpp

Minor clean up for cooley tukey in C++

Author: Gathros

contents/cooley_tukey/code/c++/fft.cpp

@@ -17,10 +17,19 @@ using std::swap;
 
 using std::size_t;
 
-using c64 = std::complex<double>;
-template <typename T>
-constexpr T pi() {
-  return 3.14159265358979323846264338327950288419716;
+using complex = std::complex<double>;
+static const double pi = 3.14159265358979323846264338327950288419716;
+
+template <typename Iter>
+void dft(Iter X, Iter last) {
+  const auto N = last - X;
+  std::vector<complex> tmp(N);
+  for (auto i = 0; i < N; ++i) {
+    for (auto j = 0; j < N; ++j) {
+      tmp[i] += X[j] * exp(complex(0, -2.0 * M_PI * i * j / N));
+    }
+  }
+  std::copy(std::begin(tmp), std::end(tmp), X);
 }
 
 // `cooley_tukey` does the cooley-tukey algorithm, recursively
@@ -30,12 +39,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<c64>(size / 2);
-    for (size_t i = 0; i < size / 2; ++i) {
+    auto temp = std::vector<complex>(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 +54,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<double>() * 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 +92,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<double>() / 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 +112,7 @@ int main() {
   std::mt19937 rng(random_device());
   std::uniform_real_distribution<double> distribution(0.0, 1.0);
 
-  std::array<c64, 64> initial;
+  std::array<complex, 64> initial;
   std::generate(
       begin(initial), end(initial), [&] { return distribution(rng); });
 
@@ -117,7 +126,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)

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:35-66, lang:"c_cpp"](code/c++/fft.cpp)
 {% sample lang="hs" %}
 [import:6-19, lang:"haskell"](code/haskell/fft.hs)
 {% sample lang="py" %}