fft.cpp

#include "benchmark_base/fft.h"

namespace comparison_benchmarks {
namespace base {
namespace fft {

complex_vector generate_input() {
  std::vector<double> known_frequencies{2, 11, 52, 88, 256};
  fft::complex_vector data(SIZE);
  // Set our input data to match a time series of the known_frequencies.
  // When applying fft to this time-series we should find these frequencies.
  for (int i = 0; i < SIZE; i++) {
    data[i] = std::complex<double>(0.0, 0.0);
    for (auto frequencie : known_frequencies) {
      data[i] += sin(2 * M_PI * frequencie * i / SIZE);
    }
  }

  return data;
}

void divide(complex_vector::iterator data, complex_vector::iterator tmp_odd_elements, int n) {
  for (int i = 0; i < n / 2; i++) {
    tmp_odd_elements[i] = data[i * 2 + 1];
  }
  for (int i = 0; i < n / 2; i++) {
    data[i] = data[i * 2];
  }
  for (int i = 0; i < n / 2; i++) {
    data[i + n / 2] = tmp_odd_elements[i];
  }
}

void combine(complex_vector::iterator data, int n) {
  for (int i = 0; i < n / 2; i++) {
    std::complex<double> even = data[i];
    std::complex<double> odd = data[i + n / 2];

    // w is the "twiddle-factor".
    // this could be cached, but we run the same 'base' algorithm parallel/serial,
    // so it won't impact the performance comparison.
    std::complex<double> w = exp(std::complex<double>(0, -2. * M_PI * i / n));

    data[i] = even + w * odd;
    data[i + n / 2] = even - w * odd;
  }
}

void conquer(complex_vector::iterator data, complex_vector::iterator swap_array, int n) {
  if (n < 2) {
    return;
  }

  divide(data, swap_array, n);
  conquer(data, swap_array, n / 2);
  conquer(data + n / 2, swap_array + n / 2, n / 2);
  combine(data, n);
}

}
}
}