#include "DiscreteFourierTransform.hh"
#include <bits/stdc++.h>

template <class T, class E>
DiscreteFourierTransform<T, E>::DiscreteFourierTransform(vector<T> signal_, double sampleSpacing_)
{
    signal = signal_;
    sampleSpacing = sampleSpacing_;
}

template <class T, class E>
void DiscreteFourierTransform<T, E>::compute(bool computeFFT)
{
    if (computeFFT)
    {
        const int N = signal.size();
        const int powerOf2 = DiscreteFourierTransform<T, E>::highestPowerof2(N);
        // Zero padding is a simple concept; it simply refers to adding zeros to end of a time-domain signal to increase its length to reach power of 2.
        while (signal.size() < powerOf2)
            signal.push_back((T)0);

        phasors = DiscreteFourierTransform<T, E>::fft(signal);
    }
    else
        phasors = DiscreteFourierTransform<T, E>::dft(signal);
}

template <class T, class E>
std::vector<E> DiscreteFourierTransform<T, E>::createTestPoints(int N, double sampleSpacing)
{
    std::vector<E> out;
    const E t = (E)sampleSpacing;
    for (int i = 0; i < N; i++)
    {
        out.push_back((E)t * i);
    }
    return out;
}

template <class T, class E>
std::vector<E> DiscreteFourierTransform<T, E>::createTestFunction(std::vector<E> x)
{
    const int N = x.size();
    std::vector<E> out;
    for (int i = 0; i < N; i++)
    {
        const E val = (E)sin(50.0 * 2.0 * M_PI * x[i]) + 0.5 * sin(80.0 * 2.0 * M_PI * x[i]);
        out.push_back(val);
    }
    return out;
}

/**
 * @brief We use Cooley-Tukey FFT Algorithm if the size of the signal  is a power of 2 otherwise you should DFT brut. In The future we will use bluestein algorithm instead of DFT.
 * Cooley-Tukey FFT Algorithms: http://people.scs.carleton.ca/~maheshwa/courses/5703COMP/16Fall/FFT_Report.pdf
 *
 * @tparam T template type of input data
 * @tparam E template type of output data
 * @param sig the signal as an input
 * @return std::vector<std::complex<E>> the phasors (array of complex number)
 */
template <class T, class E>
std::vector<std::complex<E>> DiscreteFourierTransform<T, E>::fft(std::vector<T> sig)
{
    const int N = sig.size();
    if (N == 1)
    {
        std::vector<std::complex<E>> out;
        const std::complex<E> temp((E)sig[0], (E)0);
        out.push_back(temp);
        return out;
    }

    const std::complex<E> WN = (complex<E>)std::polar(1.0, 2 * M_PI / N);
    std::complex<E> W((E)1, (E)0);

    // divide and conquer:
    // Recurse: all even samples
    std::vector<std::complex<E>>
        x_evens = fft(getEven(sig));

    // Recurse: all odd samples
    std::vector<std::complex<E>> x_odds = fft(getOdd(sig));

    // Now, combine and perform N/2 operations!
    std::complex<E> zeroComplex((E)0, (E)0);
    std::vector<std::complex<E>> x(N, zeroComplex);
    for (int k = 0; k < N / 2; k++)
    {
        x[k] = x_evens[k] + W * x_odds[k];
        x[k + (N / 2)] = x_evens[k] - W * x_odds[k];
        W = W * WN;
    }
    return x;
}

template <class T, class E>
std::vector<std::complex<E>> DiscreteFourierTransform<T, E>::dft(std::vector<T> x)
{
    const int N = x.size();
    std::complex<E> zeroComplex((E)0, (E)0);
    std::vector<std::complex<E>> out(N, zeroComplex);
    for (int k = 0; k < N; k++)
    {
        for (int n = 0; n < N; n++)
        {
            const std::complex<E> expVal = (complex<E>)std::polar(1.0, -2 * k * n * M_PI / N);
            out[k] += ((E)x[n]) * expVal;
        }
    }
    return out;
}

template <class T, class E>
std::vector<E> DiscreteFourierTransform<T, E>::computeDSP(std::vector<std::complex<E>> x, double /*frequency_*/)
{
    const int N = x.size();
    std::vector<E> out;
    for (int k = 1; k < N / 2; k++)
    {
        const E magnitudePower2 = (E)pow(abs(x[k]), 2);
        // const E dsp = (E)(2.0 / (N * frequency_)) * magnitudePower2;
        out.push_back(magnitudePower2);
    }

    return out;
}

template <class T, class E>
std::vector<E> DiscreteFourierTransform<T, E>::getFreq(std::vector<std::complex<E>> x, double frequency_)
{
    const int N = x.size();
    std::vector<E> out;
    for (int k = 1; k < N / 2; k++)
    {
        const E freq = (E)k * frequency_ / N;
        out.push_back(freq);
    }

    return out;
}

template <class T, class E>
std::vector<E> DiscreteFourierTransform<T, E>::getPeriods(std::vector<std::complex<E>> x, double frequency_)
{
    const int N = x.size();
    std::vector<E> out;
    for (int k = 1; k < N / 2; k++)
    {
        const E period = N / ((E)k * frequency_);
        out.push_back(period);
    }

    return out;
}

template <class T, class E>
std::vector<T> DiscreteFourierTransform<T, E>::getOdd(std::vector<T> x)
{
    std::vector<T> odd;
    for (int i = 0; i < x.size(); i++)
    {
        if (i % 2 != 0)
            odd.push_back(x[i]);
    }
    return odd;
}

template <class T, class E>
std::vector<T> DiscreteFourierTransform<T, E>::getEven(std::vector<T> x)
{
    std::vector<T> even;
    for (int i = 0; i < x.size(); i++)
    {
        if (i % 2 == 0)
            even.push_back(x[i]);
    }
    return even;
}

template <class T, class E>
bool DiscreteFourierTransform<T, E>::isPowerOfTwo(int N)
{
    return (N & (N - 1)) == 0;
}

template <class T, class E>
int DiscreteFourierTransform<T, E>::highestPowerof2(int n)
{
    int p = n;
    while (!DiscreteFourierTransform<T, E>::isPowerOfTwo(p))
    {
        p += 1;
    }
    return p;
}