#include "Wavetable.h"

namespace serum
{

// ---------------------------------------------------------------------------
// Wavetable
// ---------------------------------------------------------------------------
void Wavetable::clear()
{
    frames.clear();
    name = {};
}

void Wavetable::buildHarmonic (const juce::String& n, HarmonicAmpFn ampFn, int numHarmonics)
{
    name = n;
    frames.assign (kFrames, std::vector<float> (kTableSize, 0.0f));

    const float invN = 1.0f / (float) kTableSize;
    constexpr double twoPi = 6.28318530717958647692;

    for (int f = 0; f < kFrames; ++f)
    {
        auto& frame = frames[(size_t) f];

        // Precompute per-harmonic rotation step (cos/sin of the angular increment).
        std::vector<double> cosStep ((size_t) numHarmonics + 1, 0.0);
        std::vector<double> sinStep ((size_t) numHarmonics + 1, 0.0);
        std::vector<double> mag    ((size_t) numHarmonics + 1, 0.0);
        std::vector<double> phase  ((size_t) numHarmonics + 1, 0.0);

        double maxAmp = 1e-9;
        for (int h = 1; h <= numHarmonics; ++h)
        {
            const float a = ampFn (f, h);
            const double m = std::abs ((double) a);
            mag[(size_t) h]   = m;
            phase[(size_t) h] = (a < 0.0f) ? twoPi * 0.5 : 0.0;   // sign -> 0 or pi/... using sine base
            const double ang = twoPi * (double) h * (double) invN;
            cosStep[(size_t) h] = std::cos (ang);
            sinStep[(size_t) h] = std::sin (ang);
            maxAmp = std::max (maxAmp, m);
        }

        double peak = 1e-9;
        for (int h = 1; h <= numHarmonics; ++h)
        {
            if (mag[(size_t) h] < 1e-9)
                continue;

            // Start the rotating phasor at the harmonic's phase offset.
            double s = std::sin (phase[(size_t) h]);
            double c = std::cos (phase[(size_t) h]);
            const double cs = cosStep[(size_t) h];
            const double ss = sinStep[(size_t) h];
            const double m = mag[(size_t) h];

            for (int i = 0; i < kTableSize; ++i)
            {
                frame[(size_t) i] += (float) (m * s);
                const double s2 = s * cs + c * ss;
                c = c * cs - s * ss;
                s = s2;
            }
        }

        // Normalise each frame to unity peak so morphing keeps a stable level.
        for (int i = 0; i < kTableSize; ++i)
            peak = std::max (peak, (double) std::abs (frame[(size_t) i]));

        if (peak > 1e-6)
        {
            const float g = (float) (1.0 / peak);
            for (int i = 0; i < kTableSize; ++i)
                frame[(size_t) i] *= g;
        }
    }
}

float Wavetable::read (float framePos, float phase) const noexcept
{
    if (frames.empty())
        return 0.0f;

    framePos = juce::jlimit (0.0f, (float) (kFrames - 1), framePos);
    int f0 = (int) framePos;
    float frac = framePos - (float) f0;
    int f1 = juce::jmin (f0 + 1, kFrames - 1);

    float p = phase * (float) kTableSize;
    int i0 = (int) p;
    if (i0 < 0) i0 = 0;
    float t = p - (float) i0;
    int i1 = (i0 + 1) & (kTableSize - 1);
    i0 &= (kTableSize - 1);

    const auto& a = frames[(size_t) f0];
    const auto& b = frames[(size_t) f1];
    const float s0 = lerp (a[(size_t) i0], a[(size_t) i1], t);
    const float s1 = lerp (b[(size_t) i0], b[(size_t) i1], t);
    return lerp (s0, s1, frac);
}

void Wavetable::copyFrame (int frameIndex, float* dest, int numSamples) const
{
    if (frames.empty() || dest == nullptr)
        return;

    frameIndex = juce::jlimit (0, kFrames - 1, frameIndex);
    const auto& frame = frames[(size_t) frameIndex];
    const int n = juce::jmin (numSamples, kTableSize);
    for (int i = 0; i < n; ++i)
        dest[i] = frame[(size_t) i];
}

// ---------------------------------------------------------------------------
// Library
// ---------------------------------------------------------------------------
void WavetableLibrary::prebuild()
{
    for (int i = 0; i < kNumWavetables; ++i)
        getTable (i);
}

const Wavetable& WavetableLibrary::getTable (int index) const
{
    index = juce::jlimit (0, kNumWavetables - 1, index);
    if (!built[(size_t) index])
        buildTable (index);
    return tables[(size_t) index];
}

void WavetableLibrary::buildTable (int index) const
{
    switch (index)
    {
        case 0: tables[(size_t) index] = makeBasic();        break;
        case 1: tables[(size_t) index] = makeSawPwm();       break;
        case 2: tables[(size_t) index] = makeSquareSync();   break;
        case 3: tables[(size_t) index] = makeTriangleFold(); break;
        case 4: tables[(size_t) index] = makeVowel();        break;
        case 5: tables[(size_t) index] = makeOrgan();        break;
        case 6: tables[(size_t) index] = makeWarmSaw();      break;
        case 7: tables[(size_t) index] = makeDigital();      break;
        case 8: tables[(size_t) index] = makeGlass();        break;
        case 9: tables[(size_t) index] = makeBass();         break;
        default: tables[(size_t) index] = makeBasic();       break;
    }
    built[(size_t) index] = true;
}

namespace
{
    inline float normT (int frame) { return (float) frame / 255.0f; }
    constexpr float kPi = 3.14159265358979323846f;
}

// sine -> saw morph
Wavetable WavetableLibrary::makeBasic()
{
    Wavetable wt;
    wt.buildHarmonic ("Basic Shapes", [] (int f, int h)
    {
        const float t = normT (f);
        const float sine = (h == 1) ? 1.0f : 0.0f;
        const float saw  = 1.0f / (float) h;
        return sine + (saw - sine) * t;
    }, 48);
    return wt;
}

// narrow pulse -> wide pulse
Wavetable WavetableLibrary::makeSawPwm()
{
    Wavetable wt;
    wt.buildHarmonic ("Saw PWM", [] (int f, int h)
    {
        const float t = normT (f);
        const float duty = 0.08f + 0.42f * t;            // 8% -> 50%
        const float a = std::sin (kPi * (float) h * duty);
        return (2.0f / ((float) h * kPi)) * a;
    }, 48);
    return wt;
}

// hard-sync style comb sweep
Wavetable WavetableLibrary::makeSquareSync()
{
    Wavetable wt;
    wt.buildHarmonic ("Square Sync", [] (int f, int h)
    {
        const float t = normT (f);
        const float ratio = 1.0f + 3.0f * t;
        const float comb = 0.5f + 0.5f * std::cos (2.0f * kPi * (float) h * ratio);
        return (1.0f / std::pow ((float) h, 0.8f)) * comb;
    }, 48);
    return wt;
}

// triangle -> folded/saw-ish
Wavetable WavetableLibrary::makeTriangleFold()
{
    Wavetable wt;
    wt.buildHarmonic ("Triangle Fold", [] (int f, int h)
    {
        const float t = normT (f);
        // Triangle: odd harmonics with alternating sign, 1/h^2.
        float tri = 0.0f;
        if ((h & 1) == 1)
        {
            const int k = (h - 1) / 2;
            const float sign = ((k & 1) == 0) ? 1.0f : -1.0f;
            tri = sign * 8.0f / (kPi * kPi * (float) (h * h));
        }
        const float saw = 1.0f / (float) h;
        return tri + (saw - tri) * t;
    }, 48);
    return wt;
}

// formant morph a -> e -> i
Wavetable WavetableLibrary::makeVowel()
{
    Wavetable wt;
    wt.buildHarmonic ("Vowel", [] (int f, int h)
    {
        const float t = normT (f);
        // Three formant sets (Hz) with f0 = 55 Hz.
        const float f0 = 55.0f;
        const float aa[3] = { 730.0f, 1090.0f, 2440.0f };
        const float ee[3] = { 530.0f, 1840.0f, 2480.0f };
        const float ii[3] = { 270.0f, 2290.0f, 3010.0f };

        float from[3], to[3];
        if (t < 0.5f)
        {
            const float u = t * 2.0f;
            for (int k = 0; k < 3; ++k) { from[k] = aa[k]; to[k] = ee[k]; }
            (void) u;
        }
        else
        {
            const float u = (t - 0.5f) * 2.0f;
            for (int k = 0; k < 3; ++k) { from[k] = ee[k]; to[k] = ii[k]; }
            (void) u;
        }

        float tt = (t < 0.5f) ? t * 2.0f : (t - 0.5f) * 2.0f;
        float sum = 0.0f;
        const float sigma = 1.8f;   // formant bandwidth in harmonic units
        for (int k = 0; k < 3; ++k)
        {
            const float fc = from[k] + (to[k] - from[k]) * tt;
            const float center = fc / f0;
            const float d = ((float) h - center) / sigma;
            sum += std::exp (-0.5f * d * d) * (k == 0 ? 1.0f : 0.6f);
        }
        return sum;
    }, 48);
    return wt;
}

// drawbar organ morph
Wavetable WavetableLibrary::makeOrgan()
{
    Wavetable wt;
    wt.buildHarmonic ("Organ", [] (int f, int h)
    {
        const float t = normT (f);
        // two drawbar registrations (16', 8', 5 1/3', 4', 2 2/3', 2', ...)
        const float regA[8] = { 0.0f, 1.0f, 0.0f, 0.35f, 0.0f, 0.2f, 0.0f, 0.1f };
        const float regB[8] = { 0.5f, 1.0f, 0.4f, 0.6f, 0.25f, 0.5f, 0.2f, 0.35f };
        if (h > 8) return 0.0f;
        const float a = regA[h - 1];
        const float b = regB[h - 1];
        return a + (b - a) * t;
    }, 8);
    return wt;
}

// warm lowpassed saw -> brighter
Wavetable WavetableLibrary::makeWarmSaw()
{
    Wavetable wt;
    wt.buildHarmonic ("Warm Saw", [] (int f, int h)
    {
        const float t = normT (f);
        const float cutoff = 10.0f + 34.0f * t;   // harmonic rolloff centre
        const float rolloff = std::exp (-((float) h / cutoff) * ((float) h / cutoff));
        const float body = 1.0f / std::pow ((float) h, 0.9f);
        return body * (0.4f + 0.6f * rolloff);
    }, 48);
    return wt;
}

// additive, brighter and slightly combed
Wavetable WavetableLibrary::makeDigital()
{
    Wavetable wt;
    wt.buildHarmonic ("Digital", [] (int f, int h)
    {
        const float t = normT (f);
        const float comb = 0.5f + 0.5f * std::cos ((float) h * 0.9f * (1.0f + t));
        return (1.0f / std::pow ((float) h, 0.6f)) * (0.5f + 0.5f * comb);
    }, 48);
    return wt;
}

// inharmonic bell-like partial clusters
Wavetable WavetableLibrary::makeGlass()
{
    Wavetable wt;
    wt.buildHarmonic ("Glass", [] (int f, int h)
    {
        const float t = normT (f);
        const float c1 = 1.0f + t * 2.5f;
        const float c2 = 3.6f + t * 3.4f;
        const float c3 = 6.2f + t * 4.0f;
        const float sigma = 0.9f;
        float sum = 0.0f;
        const float centers[3] = { c1, c2, c3 };
        const float gains[3]  = { 1.0f, 0.7f, 0.4f };
        for (int k = 0; k < 3; ++k)
        {
            const float d = ((float) h - centers[k]) / sigma;
            sum += gains[k] * std::exp (-0.5f * d * d);
        }
        return sum;
    }, 48);
    return wt;
}

// sub-heavy -> fuller bass
Wavetable WavetableLibrary::makeBass()
{
    Wavetable wt;
    wt.buildHarmonic ("Bass", [] (int f, int h)
    {
        const float t = normT (f);
        const float sub = (h == 1) ? 1.0f : ((h == 2) ? 0.4f : ((h == 3) ? 0.12f : 0.0f));
        const float full = 1.0f / std::pow ((float) h, 1.1f);
        return sub + (full - sub) * t;
    }, 48);
    return wt;
}

} // namespace serum