FFT.h

#pragma once
#include <NDEVR/Buffer.h>
#include <NDEVR/Angle.h>
#include <iostream>
#include <iomanip>
#include <complex>
 
namespace NDEVR
{
    template<class t_type>
    class FFTLogic
    {
    protected:
        static void fft(Buffer<std::complex<t_type>>& in_out, Buffer<std::complex<t_type>>& temp, uint04 offset, uint04 size)
        {
            // base case
            if (size == 1)
                return;
            if (size % 2 != 0)// Fallback to Basic Discrete Fourier Transform
            {
                for (uint04 k = 0; k < size; k++)
                    temp[k + offset] = in_out[k + offset];
                dft(temp.begin(offset), in_out.begin(offset), size, 1);
                return;
            }
            lib_assert(size % 2 == 0, "Bad FFT");
            
            uint04 even_offset = offset;
            uint04 odd_offset = offset + size / 2;
 
            // compute FFT of even terms
            for (uint04 k = 0; k < size / 2; k++)
                temp[k + even_offset] = in_out[2 * k + offset];
 
            // compute FFT of odd terms
            for (uint04 k = 0; k < size / 2; k++)
                temp[k + odd_offset] = in_out[2 * k + 1 + offset];
            
            fft(temp, in_out, even_offset, size / 2);
            fft(temp, in_out, odd_offset, size / 2);
 
            // combine
            for (uint04 k = 0; k < size / 2; k++)
            {
                fltp08 kth = -2.0 * cast<fltp08>(k) * PI<fltp08>() / size;
                std::complex<t_type> wk(std::cos(kth), std::sin(kth));
                in_out[k + even_offset] = temp[k + even_offset] + (wk * temp[k + odd_offset]);
                in_out[k + odd_offset] = temp[k + even_offset] - (wk * temp[k + odd_offset]);
            }
        }

        static void dft(const std::complex<t_type> f[], std::complex<t_type> ftilde[], uint04 N, uint04 step)
        {
            std::complex<t_type> w = std::polar(1.0, -2.0 * PI<t_type>() / N);
            std::complex<t_type> wn = 1.0;
            for (uint04 n = 0; n < N; n++)
            {
                if (n > 0) wn *= w;
                std::complex<t_type> wm = 1.0;
                ftilde[n] = 0.0;
                for (uint04 m = 0, mpos = 0; m < N; m++, mpos += step)
                {
                    if (m > 0) wm *= wn;
                    ftilde[n] += f[mpos] * wm;
                }
            }
        }

        static void ifft(Buffer<std::complex<t_type>>& in_out, Buffer<std::complex<t_type>>& temp, uint04 offset, uint04 size)
        {
            // take conjugate
            for (uint04 i = 0; i < size; i++)
                in_out[i + offset] = std::conj(in_out[i + offset]);
 
            // compute forward FFT
            fft(in_out, temp, offset, size);
 
            // take conjugate again
            t_type factor = 1.0 / size;
            for (uint04 i = 0; i < size; i++)
                in_out[i + offset] = std::conj(in_out[i + offset]) * factor;
 
        }
 
        
 
        //======================================================================

        static void iFFT(const std::complex<t_type> ftilde[], std::complex<t_type> f[], uint04 N)
        {
            Buffer<std::complex<t_type>> ftildeConjugate;
            ftildeConjugate.setSize(N);
 
            for (uint04 m = 0; m < N; m++)
                ftildeConjugate[m] = std::conj(ftilde[m]);
 
            fft(ftildeConjugate, f, N, 1);
 
            t_type factor = 1.0 / N;
            for (uint04 m = 0; m < N; m++)
                f[m] = std::conj(f[m]) * factor;
 
        }
    public:
        static void FT(Buffer<std::complex<t_type>>& in_out)
        {
            Buffer<std::complex<t_type>> temp;
            temp.setSize(in_out.size());
            fft(in_out, temp, 0, temp.size());
        }

        static void IFT(Buffer<std::complex<t_type>>& in_out)
        {
            Buffer<std::complex<t_type>> temp;
            temp.setSize(in_out.size());
            ifft(in_out, temp, 0, temp.size());
        }
    };
        
        /*static Buffer<std::complex<t_type>> fft(const Buffer<std::complex<t_type>>& f)           // Fast Fourier Transform
        {
            // base case
            if (f.size() == 1)
            {
                return { f[0] };
            }
            lib_assert(f.size() % 2 == 0, "Bad FFT");
            
            Buffer<std::complex<t_type>> buffer;
            buffer.setSize(f.size() / 2);
            // compute FFT of even terms
            for (uint04 k = 0; k < f.size() / 2; k++)
                buffer[k] = f[2 * k];
            Buffer<std::complex<t_type>> even_fft = fft(buffer);
 
            // compute FFT of odd terms
            for (uint04 k = 0; k < f.size() / 2; k++)
                buffer[k] = f[2 * k + 1];
            Buffer<std::complex<t_type>> odd_fft = fft(buffer);
 
            // combine
            Buffer<std::complex<t_type>> y;
            y.setSize(f.size());
            for (uint04 k = 0; k < f.size() / 2; k++) 
            {
                fltp08 kth = -2.0 * cast<fltp08>(k) * Angle::PI<fltp08>() / f.size();
                std::complex<t_type> wk(std::cos(kth), std::sin(kth));
                y[k               ] = even_fft[k] + (wk * odd_fft[k]);
                y[k + f.size() / 2] = even_fft[k] - (wk * odd_fft[k]);
            }
            return y;
        }
 
        // compute the inverse FFT of x[], assuming its length n is a power of 2
        static Buffer<std::complex<t_type>> ifft(const Buffer<std::complex<t_type>>& x)
        {
            Buffer<std::complex<t_type>> y;
            y.setSize(x.size());
 
            // take conjugate
            for (uint04 i = 0; i < x.size(); i++)
                y[i] = std::conj(x[i]);
 
            // compute forward FFT
            y = fft(y);
 
            // take conjugate again
            t_type factor = 1.0 / x.size();
            for (uint04 i = 0; i < x.size(); i++)
                y[i] = std::conj(y[i]) * factor;
            return y;
 
        }
 
        //======================================================================
 
        static void DFT(const std::complex<t_type> f[], std::complex<t_type> ftilde[], uint04 N, uint04 step)           // Basic Discrete Fourier Transform
        {
            std::complex<t_type> w = std::polar(1.0, -2.0 * Angle::PI<t_type>() / N);
            std::complex<t_type> wn = 1.0;
            for (uint04 n = 0; n < N; n++)
            {
                if (n > 0) wn *= w;
                std::complex<t_type> wm = 1.0;
                ftilde[n] = 0.0;
                for (uint04 m = 0, mpos = 0; m < N; m++, mpos += step)
                {
                    if (m > 0) wm *= wn;
                    ftilde[n] += f[mpos] * wm;
                }
            }
        }
 
        //======================================================================
 
        static void iFFT(const std::complex<t_type> ftilde[], std::complex<t_type> f[], uint04 N)                    // Inverse Fast Fourier Transform
        {
            Buffer<std::complex<t_type>> ftildeConjugate;
            ftildeConjugate.setSize(N);
 
            for (uint04 m = 0; m < N; m++) 
                ftildeConjugate[m] = std::conj(ftilde[m]);
 
            FFT(ftildeConjugate, f, N, 1);
 
            t_type factor = 1.0 / N;
            for (uint04 m = 0; m < N; m++)
                f[m] = std::conj(f[m]) * factor;
 
        }
    public:
        static Buffer<std::complex<t_type>> FFT(const Buffer<std::complex<t_type>>& complex, uint04 step = 1)
        {
            Buffer<std::complex<t_type>> output;
            output.setSize(complex.size());
            FFT(complex.ptr(), output.ptr(), complex.size(), step);
            return output;
        }
        static Buffer<std::complex<t_type>> IFFT(const Buffer<std::complex<t_type>>& complex)
        {
            return ifft(complex);
        }
    };*/
}