compile-time fast fourier transform

cfft.cpp

#include <iostream>
#include <complex>
#include <cmath>
#include <utility>
#include <chrono>
#include <array>


// std::complex<T>::operators -- not constexpr in C++17
template <class T>
struct Complex{
    T re;
    T im;
    constexpr Complex(T r = 0, T i = 0) : re(r), im(i) {}
    constexpr T real() const {return re;}
    constexpr T imag() const {return im;}
};  

template <class T>
constexpr Complex<T> operator + (const Complex<T>& a, const Complex<T>& b){
    return { a.re + b.re, a.im + b.im};
}

template <class T>
constexpr Complex<T> operator - (const Complex<T>& a, const Complex<T>& b){
    return { a.re - b.re, a.im - b.im};
}

template <class T>
constexpr Complex<T> operator * (const Complex<T>& a, const Complex<T>& b){
    return { a.re * b.re - a.im * b.im,   a.re * b.im + a.im * b.re};
}

template <class T>
constexpr T abs(const Complex<T>& x){
    return std::hypot(x.real(), x.imag());
}


//start FFT-implementation

template <class T>
using element_type = Complex<T>;

template <class T, size_t N, int64_t I>
constexpr T circle_angle = static_cast<T>(2 * M_PI / N * I);

template <class T, size_t N, int64_t I>
constexpr element_type<T> circle_element = { static_cast<T>(cos(circle_angle<T, N, I>)),
                                             static_cast<T>(sin(circle_angle<T, N, I>)) };
template <class T, size_t N>
using element_sequence = std::array<element_type<T>, N>;



template <class T, size_t N, size_t... I>
constexpr element_sequence<T, 2*N> butterfly(const element_sequence<T, N>& first,  const element_sequence<T, N>& second,
                                             std::integer_sequence<size_t, I...> ){
    return { (first[I] + second[I] * circle_element<T, 2*N, I>)...,
             (first[I] - second[I] * circle_element<T, 2*N, I>)...  };
}


template<class T, size_t N>
constexpr auto cfft(const element_sequence<T, N>& x);

template<class T, size_t N, size_t... I>
constexpr auto cfft_impl(const element_sequence<T, N>& x, std::integer_sequence<size_t, I...> idx)
{
    static_assert (N % 2 == 0, "N must be 2^d");
    static_assert (N == 2 * sizeof...(I));
    
    return butterfly(cfft<T, N/2>({x[I * 2]...}),
                     cfft<T, N/2>({x[I * 2 + 1]...}),
                     idx);
}

template<class T, size_t N>
constexpr auto cfft(const element_sequence<T, N>& x){
    if constexpr (N == 1)
        return x;
    else
        return cfft_impl(x, std::make_integer_sequence<size_t, N/2>{});
}

//end FFT-implementation




template <class T, size_t N, size_t... I>
constexpr element_sequence<T, N> sin_sequence_helper(T freq, std::integer_sequence<size_t, I...>){
    return { element_type<T>{static_cast<T>(cos(freq * I / N)), 
                             static_cast<T>(sin(freq * I / N))} ... };
}

template <class T, size_t N>
constexpr element_sequence<T, N> sin_sequence(T freq){
    return sin_sequence_helper<T, N>(freq, std::make_integer_sequence<size_t, N>{});
}




int main()
{
    constexpr auto sin_s = sin_sequence<float, 16>(2* M_PI);
    //constexpr auto sin_s = element_sequence<float, 16>{1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0};

    constexpr auto fft = cfft<float>(sin_s);

    for (auto x : fft)
        std::cout << abs(x) << std::endl;


    return 0;
}