/**
Various window types.

Copyright: Guillaume Piolat 2015.
License:   $(LINK2 http://www.boost.org/LICENSE_1_0.txt, Boost License 1.0)
*/
module dplug.dsp.window;

import std.math: PI, sqrt, cos, isClose;

import dplug.core.vec;

// Useful reference:
// https://en.wikipedia.org/wiki/Window_function
// These window generation functions are explicit with the end-points, 
// else it would be very easy to get confused and ignorant of these issues.

enum WindowType
{
    /// Constant window.
    rect,
    
    /// Triangular window. You probably want something better.
    bartlett,

    /// Good default choice if you lack ideas.
    hann,

    hamming,
    
    nuttall,
    
    blackmannNutall, 
    
    blackmannHarris,

    flatTopSR785,    // Flat top window

    /// Kaiser-Bessel window, this one need a parameter alpha (typical values: 1.0 to 4.0)
    /// Note: Kaiser windows are talked about in terms of parameters alpha or beta.
    kaiserBessel,

    // Older names

    /// Please use the new name instead.
    RECT = rect,

    /// Please use the new name instead.
    BARTLETT = bartlett,

    /// Please use the new name instead.
    HANN = hann,

    /// Please use the new name instead.
    HAMMING = hamming,

    /// Please use the new name instead.
    NUTTALL = nuttall,

    /// Please use the new name instead.
    BLACKMANN_NUTTALL = blackmannNutall,

    /// Please use the new name instead. 
    BLACKMANN_HARRIS = blackmannHarris,

    /// Please use the new name instead.
    FLATTOP_SR785 = flatTopSR785,

    /// Please use the new name instead.
    KAISER_BESSEL = kaiserBessel,
}


/// How "aligned" the window is in its support.
/// Very important, see Issue #236.
/// When choosing a window, you probably want to take a hard look about `WindowSymmetry` 
/// because it has implications for latency and correctness.
enum WindowAlignment
{
    /// The window is asymmetric and if it goes to zero (eg: HANN), it will have one zero 
    /// at the FIRST coefficient output[N-1].
    /// Its center is exactly at output[(N/2)-1] which is an integer delay for even window
    /// lengths.
    /// This loose one sample of latency, however it has the easiest to compute latency
    /// in most usage.
    right,

    /// The window is asymmetric and if it goes to zero (eg: HANN), it will have one zero 
    /// at the LAST coefficient output[N-1].
    /// Its center is exactly at output[(N/2)-1] which is an integer delay for even window
    /// lengths.
    /// This has the best latency characteristics, however this might make the latency
    /// computation itself trickier.
    left,    

    /// The window is symmetric and if it goes to zero (eg: HANN), it will have two zero 
    /// coefficients. 
    /// Its center is exactly between samples at position output[(N - 1)/2] which is NOT 
    /// an integer delay for even window lengths.
    /// Such a window might also break derivativeness.
    /// However FOR HISTORICAL REASONS THIS WAS THE DEFAULT SETTINGS.
    /// IT IS ADVISED NOT TO USE IT, FOR THE CONSEQUENCES ARE SUBTLE AND OFTEN HARMFUL.
    symmetric,
}


struct WindowDesc
{
public:
nothrow:
@nogc:

    /// Construct a window description, with a foolproof constructor.
    this(WindowType type, WindowAlignment alignment, float param = float.nan)
    {
        this.type = type;
        this.alignment = alignment;
        this.param = param;
    }

private:
    WindowType type;

    // TODO: this is a bad default! You'll probably want WindowAlignment.right instead.
    // Make sure it isn't used implicitely anymore.
    // Then change that default
    WindowAlignment alignment = WindowAlignment.symmetric; 

    float param;
}

/// Generates a window described by `windowDesc`, with periodicity of 
/// `outputWindow.length`.
void generateWindow(T)(WindowDesc desc, T[] outputWindow) pure nothrow @nogc
{
    int N = cast(int)(outputWindow.length);
    for (int i = 0; i < N; ++i)
    {
        outputWindow[i] = cast(T)(evalWindow(desc, i, N));
    }
}

/// Multiplies the given slice in-place by a window described by `windowDesc`,
/// whose periodicity is `inoutImpulse.length`.
void multiplyByWindow(T)(T[] inoutImpulse, WindowDesc windowDesc) pure nothrow @nogc
{
    int N = cast(int)(inoutImpulse.length);
    for (int i = 0; i < N; ++i)
    {
        inoutImpulse[i] *= evalWindow(windowDesc, i, N);
    }
}

// Note: N is the number of output samples, not necessarily the periodicity
double evalWindow(WindowDesc desc, int n, int N) pure nothrow @nogc
{
    final switch(desc.alignment) with (WindowAlignment)
    {
        case right:
            return evalWindowInternal(desc, n, N);

        case left:
            // left are just a rotation of the right-aligned window
            if (n == 0)
                return evalWindowInternal(desc, N - 1, N);
            else
                return evalWindowInternal(desc, n + 1, N);

        case symmetric:
            // Symmetric is just a shorter window
            return evalWindowInternal(desc, n, N - 1);
    }
}


struct Window(T) if (is(T == float) || is(T == double))
{
    void initialize(WindowDesc desc, int lengthInSamples) nothrow @nogc
    {
        _lengthInSamples = lengthInSamples;
        _window.reallocBuffer(lengthInSamples);
        generateWindow!T(desc, _window);
    }

    ~this() nothrow @nogc
    {
        _window.reallocBuffer(0);
    }

    double sumOfWindowSamples() pure const nothrow @nogc
    {
        double result = 0;
        foreach(windowValue; _window)
            result += windowValue;
        return result;
    }

    @disable this(this);

    T[] _window = null;
    int _lengthInSamples;
    alias _window this;
}


private:

// Generates WindowAlignment.right window types.
// N is the periodicity.
// Does NOT look at desc.alignment, which is taken care of in `evalWindow`.
double evalWindowInternal(WindowDesc desc, int n, int N) pure nothrow @nogc
{
    static double computeKaiserFunction(double alpha, int n, int N) pure nothrow @nogc
    {
        static double I0(double x) pure nothrow @nogc
        {
            double sum = 1;
            double mx = x * 0.5;
            double denom = 1;
            double numer = 1;
            for (int m = 1; m <= 32; m++) 
            {
                numer *= mx;
                denom *= m;
                double term = numer / denom;
                sum += term * term;
            }
            return sum;
        }

        double piAlpha = PI * alpha;
        double C = 2.0 * n / cast(double)N - 1.0; 
        double result = I0(piAlpha * sqrt(1.0 - C * C)) / I0(piAlpha);
        return result;
    }

    final switch(desc.type)
    {
        case WindowType.RECT:
            return 1.0;

        case WindowType.BARTLETT:
        {
            double nm1 = cast(double)N / 2;
            double t = n - nm1;
            if (t < 0)
                t = -t;
            return 1 - t / nm1;
        }

        case WindowType.HANN:
        {
            double phi = (2 * PI * n ) / N;
            return 0.5 - 0.5 * cos(phi);
        }

        case WindowType.HAMMING:
        {
            double phi = (2 * PI * n ) / N;
            return 0.54 - 0.46 * cos(phi);
        }

        case WindowType.NUTTALL:
        {
            double phi = (2 * PI * n ) / N;
            return 0.355768 - 0.487396 * cos(phi) + 0.144232 * cos(2 * phi) - 0.012604 * cos(3 * phi);            
        }

        case WindowType.BLACKMANN_HARRIS:
        {
            double phi = (2 * PI * n ) / N;
            return 0.35875 - 0.48829 * cos(phi) + 0.14128 * cos(2 * phi) - 0.01168 * cos(3 * phi);
        }

        case WindowType.BLACKMANN_NUTTALL:
        {
            double phi = (2 * PI * n ) / N;
            return 0.3635819 - 0.4891775 * cos(phi) + 0.1365995 * cos(2 * phi) - 0.0106411 * cos(3 * phi);
        }

        case WindowType.FLATTOP_SR785:
        {
            double phi = (2 * PI * n ) / cast(double)N;
            return 1 - 1.93 * cos(phi) 
                     + 1.29 * cos(2 * phi) 
                     - 0.388 * cos(3 * phi) 
                     + 0.028 * cos(4 * phi);
        }

        case WindowType.KAISER_BESSEL:
        {
            return computeKaiserFunction(desc.param, n, N);
        }
    }
}


// Tests the general shape of windows.
unittest
{
    void checkIsSymmetric(T)(T[] inp)
    {
        for (int i = 0; i < inp.length/2; ++i)
        {
            // Window generation should be precise
            assert(isClose(inp[i], inp[$-1-i], 1e-10));
        }
    }

    void checkMaxInCenter(T)(T[] inp)
    {
        bool odd = (inp.length & 1) != 0;
        if (odd)
        {
            for (int i = 0; i < inp.length/2; ++i)
            {
                assert(inp[i] <= inp[inp.length/2]);
                assert(inp[$-1-i] <= inp[inp.length/2]);
            }
        }
        else
        {
            for (int i = 0; i < inp.length/2-1; ++i)
            {
                // take either of the centers
                assert(inp[i] <= inp[inp.length/2-1]);
                assert(inp[$-1-i] <= inp[inp.length/2-1]);
            }
        }
    }

    void testAllWindows(T)(int size)
    {
        Window!T window;

        foreach(type; WindowType.min..WindowType.max+1)
        {
            foreach(alignment; WindowAlignment.min..WindowAlignment.max+1)
            {
                // Note: only Kaiser-Bessel have a parameter for now, so take 2
                WindowDesc desc = WindowDesc(cast(WindowType)type, 
                                             cast(WindowAlignment)alignment, 
                                             2.0);

                window.initialize(desc, size);

                final switch(alignment) with (WindowAlignment)
                {
                    case WindowAlignment.symmetric:
                        checkIsSymmetric(window[0..$]);
                        checkMaxInCenter(window[0..$]);
                        break;

                    case WindowAlignment.right:
                        checkIsSymmetric(window[1..$]);
                        checkMaxInCenter(window[1..$]);
                        break;

                    case WindowAlignment.left:
                        checkIsSymmetric(window[0..$-1]);
                        checkMaxInCenter(window[0..$-1]);
                        break;
                }
            }
        }
    }

    testAllWindows!float(64);
    testAllWindows!double(32);

    // It should be possible to generate odd-sized window, 
    // though in practice with WindowAlignement you don't need them.
    testAllWindows!float(5);
}