/* Copyright (C) 2025 Joel Yliluoma - Material for https://youtu.be/C4tHxouarGc - https://iki.fi/bisqwit/ */
/* MIT license follows:
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated
documentation files (the "Software"), to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software,
and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of
the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS
OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT
OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include <vector>
#include <complex>

using complex = std::complex<float>;

// Baseclass for all DFT methods.
// Inheriting class must implement at least one of xform() or xform_many(), and also get_name().
struct DFTbase
{
    std::size_t N;
    bool        inplace;
public:
    // Inplace tells whether the method can safely use the parameter as both input and output.
    DFTbase(std::size_t n, bool is_inplace) : N(n), inplace(is_inplace) {}
    virtual ~DFTbase() {}
    /* The basic transformation: in[] into out[]. */
    virtual void xform(const complex* in, std::size_t istep, complex* out, std::size_t ostep) 
    {
        xform_many(in,istep,0, out,ostep,0, 1);
    }
    /* Many transformations in sequence. */
    virtual void xform_many(const complex* in, std::size_t istep, std::size_t istep2,
                            complex* out,      std::size_t ostep, std::size_t ostep2,
                            std::size_t num)
    {
        for(std::size_t n=0; n<num; ++n)
            xform(in + n*istep2, istep, out + n*ostep2, ostep);
    }
    /* In-place transformation: inout[] is both input and output. */
    void xform_inplace(complex* inout, std::size_t step)
    {
        xform_inplace_many(inout,step,0, 1);
    }
    /* Many in-place transformations in a sequence. */
    void xform_inplace_many(complex* inout, std::size_t step, std::size_t step2, std::size_t num)
    {
        if(inplace)
        {
            xform_many(inout,step,step2, inout,step,step2,num);
            return;
        }
        std::vector<complex> temp(N*num);
        xform_many(inout,step,step2, &temp[0],1,N, num);
        for(std::size_t n=0; n<num; ++n)
            for(std::size_t a=0; a<N; ++a)
                inout[a*step + n*step2] = temp[a + n*N];
    }
    /* A method for returning the name of the transformation method. */
    virtual std::string get_name() const = 0;

    /* Inverse transforms */
    // There are two equalities regarding inverse DFTs discussed in my thesis.
    //   IDFT(x) = (1/N) DFT(x*)*  and
    //   IDFT(x) = (1/N) DFT(x^)^
    //  where * suffix refers to the complex conjugates of the set: conj(zᵣ + izᵢ) = zᵣ - izᵢ
    //  and   ^ suffix refers to the real-imag swaps of the set:    swap(zᵣ + izᵢ) = zᵢ + izᵣ
    // Either would work. Of these, this code uses the former.
    // In the Rader and Bluestein implementations, these operations are baked-in
    // rather than calling the inverse methods, so that we can skip doing the
    // division every time the transform is performed.
    void inverse_xform(const complex* in, std::size_t istep, complex* out, std::size_t ostep)
    {
        inverse_xform_many(in,istep,0, out,ostep,0, 1);
    }
    void inverse_xform_many(const complex* in, std::size_t istep, std::size_t istep2,
                            complex* out,      std::size_t ostep, std::size_t ostep2,
                            std::size_t num)
    {
        if(inplace)
        {
            for(std::size_t n=0; n<num; ++n)
                for(std::size_t a=0; a<N; ++a)
                    out[a*ostep+n*ostep2] = std::conj(in[a*istep+n*istep2]);
            xform_many(out,ostep,ostep2, out,ostep,ostep2,num);
        }
        else
        {
            std::vector<complex> temp(N*num);
            for(std::size_t n=0; n<num; ++n)
                for(std::size_t a=0; a<N; ++a)
                    temp[a+n*N] = std::conj(in[a*istep+n*istep2]);
            xform_many(&temp[0],1,N, out,ostep,ostep2, num);
        }
        float divisor = 1.f / N;
        for(std::size_t n=0; n<num; ++n)
            for(std::size_t a=0; a<N; ++a)
                out[a*ostep+n*ostep2] = std::conj(out[a*ostep+n*ostep2]) * divisor;
    }
    void inverse_xform_inplace(complex* inout, std::size_t step)
    {
        inverse_xform_inplace_many(inout,step,0, 1);
    }
    void inverse_xform_inplace_many(complex* inout, std::size_t step, std::size_t step2, std::size_t num)
    {
        inverse_xform_many(inout,step,step2, inout,step,step2,num);
    }
};

DFTbase* FindDFT(std::size_t N, std::size_t istep,std::size_t istep2,
                                std::size_t ostep,std::size_t ostep2,
                                std::size_t num);