/*****************************************************************************
*                                                                            *
*   SGOLAYFILT.SPL Copyright (C) 2015 DSP Development Corporation            *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Savitzky-Golay smoothing filering                           *
*                                                                            *
*   Revisions:    1 Jul 2015  RRR  Creation                                  *
*                                                                            *
*****************************************************************************/


#if @HELP_SGOLAYFILT

    SGOLAYFILT  

    Purpose: 
             Filters a series with a Savitzky-Golay smoothing filter.
                                                                        
    Format:  
             SGOLAYFILT(x, d, n, w)

                x - A series, the input series to smooth.

                d - Optional. An integer or array. If an integer, the
                    Savitzky-Golay polynomial degree. If an array, the 
                    explicit Savitzky-Golay coefficients. Defaults to the
                    integer 3, design and process a 3rd order filter.

                n - Optional. An integer, the number of points for the
                    smoothing window. Defaults to d + 2 if d is odd, else
                    d + 3.

                w - Optional. A series, the weighting factors. Defaults to
                    all ones.
    Returns: 
             A series, the result of filtering the input data with the
             smoothing filter.

    Example:
             W1: gsweep(1000, 1/1000, 1, 50) + gnorm(1000, 1/1000)/5
             W2: sgolayfilt(w1)
             W3: linavg(w1, 10)
             W4: w2;overp(w3, lred)

             W1 contains 1000 samples of a noisy swept sinewave from 1 Hz
             to 50Hz.

             W2 filters the data in W1 with a 3rd order 5 point Savitzky-Golay
             filter.

             W3 filters the data in W1 with a 10 point zero phase moving 
             average

             W4 graphically compares the two results. The Savitzky-Golay
             filter more accurately preserves the peak heights of the
             high frequency oscillations of the input data.

    Example:
             W1: gsweep(1000, 1/1000, 1, 50) + gnorm(1000, 1/1000)/5
             W2: sgolayfilt(w1, 3, 11)
             W3: linavg(w1, 10)
             W4: w2;overp(w3, lred)

             Same as above except a 3rd order 11 point Savitzky-Golay filter
             is used.

    Example:
             W1: gsweep(1000, 1/1000, 1, 50) + gnorm(1000, 1/1000)/5
             W2: linavg(w1, 10)
             W3: sgolay(3, 11);
             W4: sgolayfilt(w1, w3);overp(w2, lred);

             Similar to above except the Savitzky-Golay filter coefficients
             are given as an explicit input array.

    Remarks:
             Savitzky-Golay filters perform data smoothing and preserve
             peak resolution by least squares fitting a D order polynomial
             to an N point sliding window. The data window, N, should be an
             odd value.

             If the filter coefficients are explicit, the window length N and
             weights w are ignored if provided.

    See Also:
             Conv
             Filteq
             Linavg
             Movavg
             Sgolay
#endif


static sgolay_coef = {};


/* filter data with Savitzky-Golay length N order d design */
sgolayfilt(x, d, N, w)
{
        local M, L, B, y, yit, yot, yss;

        if (argc < 4)
        {
                if (argc < 3)
                {
                        if (argc < 2)
                        {
                                if (argc < 1) error("sgolayfilt - input series required");

                                d = 3;
                        }

                        if (isarray(d))
                        {
                                N = length(d);
                        }
                        else
                        {
                                N = d + 2 + (d % 2 == 0);
                        }
                }

                w = ones(N, 1);
        }

        L = length(x);

        if (isarray(d))
        {
                /* assume coefficients already given */
                B = d;
                N = length(B);
        }
        else
        {
                /* design filter */
                B = sgolay(d, N, w);
        }

        /* local copy */
        sgolay_coef = B;

        /* filter - use static copy of coefficients for column iteration */
        y = sgolayfilt_filt(x);

        sgolay_coef = {};

        if (outargc > 1)
        {
                return(y, B);
        }
        else
        {
                return(y);
        }
}


/* Savitzky-Golay filtering, coefficients in static variable sgolay_coef */
ITERATE sgolayfilt_filt(x)
{
        local M, L, y, yit, yot, yss;

        /*
         *  filter coefficients in static sgolay_coef so we can column iterate
         *  on the input but use all columns of the coefficients
         */

        L = length(x);
        N = length(sgolay_coef);

        if (L <= N + 1) error("sgolayfilt - filter length too long");

        M = floor((N-1) / 2);

        /* input on transients */
        yit = sgolay_coef[1..M+1, ..] *^ x[1..N];

        /* input off transients */
        yot = sgolay_coef[M+1..2*M+1, ..] *^ x[L-N+1..L];

        /* steady state */
        yss = extract(filteq(sgolay_coef[.., M+1], {1}, x), N+1, length(x) - length(yit) - length(yot));

        /* combine */
        y = {yit, yss, yot};

        /* properties */
        setxoffset(y, xoffset(x));
        setdeltax(y, deltax(x));
        setvunits(y, getvunits(x));
        sethunits(y, gethunits(x));

        return(y);
}