main.cpp 2.26 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14
#include <pls/pls.h>
#include <pls/internal/helpers/profiler.h>
#include <pls/internal/helpers/mini_benchmark.h>

#include <iostream>
#include <complex>
#include <vector>

static constexpr int CUTOFF = 10;
static constexpr int NUM_ITERATIONS = 1000;
static constexpr int INPUT_SIZE = 2064;
typedef std::vector<std::complex<double>> complex_vector;

void divide(complex_vector::iterator data, int n) {
15 16 17 18 19 20 21 22 23 24
  complex_vector tmp_odd_elements(n / 2);
  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];
  }
25 26 27
}

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

32 33 34 35
    // w is the "twiddle-factor".
    // this could be cached, but we run the same 'data_structures' 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));
36

37 38 39
    data[i] = even + w * odd;
    data[i + n / 2] = even - w * odd;
  }
40 41 42
}

void fft(complex_vector::iterator data, int n) {
43 44 45
  if (n < 2) {
    return;
  }
46

47 48 49 50 51 52 53 54 55 56 57
  divide(data, n);
  if (n <= CUTOFF) {
    fft(data, n / 2);
    fft(data + n / 2, n / 2);
  } else {
    pls::invoke_parallel(
        [&] { fft(data, n / 2); },
        [&] { fft(data + n / 2, n / 2); }
    );
  }
  combine(data, n);
58 59 60
}

complex_vector prepare_input(int input_size) {
61 62
  std::vector<double> known_frequencies{2, 11, 52, 88, 256};
  complex_vector data(input_size);
63

64 65 66 67 68 69
  // 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 < input_size; i++) {
    data[i] = std::complex<double>(0.0, 0.0);
    for (auto frequencie : known_frequencies) {
      data[i] += sin(2 * M_PI * frequencie * i / input_size);
70
    }
71
  }
72

73
  return data;
74 75 76
}

int main() {
77 78
  PROFILE_ENABLE
  complex_vector initial_input = prepare_input(INPUT_SIZE);
79

80 81 82 83
  pls::internal::helpers::run_mini_benchmark([&] {
    complex_vector input = initial_input;
    fft(input.begin(), input.size());
  }, 8, 4000);
84

85
  PROFILE_SAVE("test_profile.prof")
86
}