SaF FFT

SaF_FFT.hpp

// Copyright (c) 2016 Roman Beránek. All rights reserved.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
//  This program is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.
//
//  You should have received a copy of the GNU General Public License
//  along with this program.  If not, see <http://www.gnu.org/licenses/>.

#pragma once

#include <complex>
#include <gsl/span>

#include "../common/math.hpp"

namespace reBass {
template <typename T, std::ptrdiff_t N>
class SaF_FFT
{
    using real_t = T;
    using cpx_t = std::complex<T>;

    static_assert(std::is_floating_point<T>::value);
    static_assert((N & (N - 1)) == 0, "N must be a power of 2.");

public:
    SaF_FFT()
    noexcept {
        auto const step = real_t{ -2 * math::pi<real_t> / N };
        for (auto i = 0u; i < N; ++i) {
            twiddles[i] = std::polar(real_t{1}, i * step);
        }
    }

    void
    transform_forward(
        gsl::span<cpx_t const, N> input,
        gsl::span<cpx_t, N> output
    ) const noexcept {
        step_into<false>(input, output);
    }

    void
    transform_backward(
        gsl::span<cpx_t const, N> input,
        gsl::span<cpx_t, N> output
    ) const noexcept {
        step_into<true>(input, output);
    }

private:
    template <bool Inverse, std::ptrdiff_t N_out, std::ptrdiff_t N_in>
    inline void
    step_into(
        gsl::span<cpx_t const, N_in> input,
        gsl::span<cpx_t, N_out> output
    ) const noexcept {
        auto const n_2 = N_out/2;
        auto const stride = (N_out > 0) ? N / N_out : 0;

        if (n_2 == 1) {
            output[0] = input[0];
            output[1] = input[stride];
        } else {
            auto const new_n_in = (N_in > stride) ? N_in - stride : 0;
            step_into<Inverse>(
                gsl::subspan<0, new_n_in>(input),
                gsl::subspan<0, n_2>(output)
            );
            step_into<Inverse>(
                gsl::subspan<stride, new_n_in>(input),
                gsl::subspan<n_2, n_2>(output)
            );
        }

        // special case for i = 0: the first twiddle is equal to 1,
        // therefore no multiplication is needed
        auto t = output[n_2];
        output[n_2] = output[0] - t;
        output[0] += t;

        for (auto i = 1; i < n_2; ++i) {
            t = math::fast_multiply(
                output[i + n_2],
                twiddles[Inverse ? (N - 1 - i*stride) : i*stride]
            );
            output[i + n_2] = output[i] - t;
        }
    }

    std::array<cpx_t, N> twiddles;
};
}