1
2
3
4
5
6
#include <complex>
#include <valarray>

constexpr double PI = 3.141592653589793238460;
typedef std::complex<double> Complex;
typedef std::valarray<Complex> CArray;

30.2 The DFT and FFT

DFT 를 구하기위해 FFT를 사용한다.

ref https://rosettacode.org/wiki/Fast_Fourier_transform

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
CArray RECURSIVE_FFT(CArray& a)
{

    const std::size_t n = a.size();
    if (n <= 1)
        return a;
    //w = cos(x) + i * sin(x)
    // w_n = e^{(2pi/n)*i }
    Complex wn = std::polar(1.0, 2 * PI / n);
    Complex w = 1;

    // divide
    //even
    CArray a0 = a[std::slice(0, n / 2, 2)];
    //odd
    CArray a1 = a[std::slice(1, n / 2, 2)];

    // conquer
    CArray y0 = RECURSIVE_FFT(a0);
    CArray y1 = RECURSIVE_FFT(a1);

    CArray y(Complex(0), n);
    // combine
    for (size_t k = 0; k < n / 2; ++k)
    {
        y[k] = y0[k] + w * y1[k];
        y[k + n / 2] = y0[k] - w * y1[k];
        w = w * wn;
    }
    return y;
}

1