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


#if @HELP_SGOLAY

    SGOLAY  

    Purpose: 
             Generates Savitzky-Golay smoothing filter coefficients.
                                                                        
    Format:  
             SGOLAY(d, N, w)

                d - Optional. An integer, the Savitzky-Golay polynomial degree.
                    Defaults to 3.

                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: 
             An NxN array, the Savitzky-Golay smoothing filter coefficients.

             (B, G) = sgolay(d, N, w) returns the smoothing filter coefficients
             and derivative coefficients as two separate arrays.

    Example:
             W1: sgolay()
             W2: sgolay(4, 9)

             W1 contains a 5x5 array of 3rd order Savitzky-Golay smoothing
             filter coefficients. W2 contains a 9x9 array of 4th order
             Savitzky-Golay smoothing coefficients.

    Example:
             (b, g) = sgolay(4, 9)

             B contains a 9x9 array of 4th order Savitzky-Golay smoothing 
             filter coefficients and G contains a 9x5 array of Savitzky-Golay
             derivative coefficients.

    Example:
             W1: gsweep(1000, 1/1000, 1, 50) + gnorm(1000, 1/1000)/5
             W2: sgolay(3, 21);
             W3: sgolayfilt(w1, w2)

             W2 contains a 3rd order, 21 point Savitzky-Golay smoothing filter.
             W3 process the data in W1 with the filter.

    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.

             See SGOLAYFILT to compute the Savitzky-Golay coefficients and
             filter input data in one step.

    See Also:
             Conv
             Filteq
             Linavg
             Movavg
             Sgolayfilt
#endif


/* compute Savitzky-Golay filter coefficients */
sgolay(d, N, w)
{
        local M, B, S, q, r;

        if (argc < 3)
        {
                if (argc < 2)
                {
                        if (argc < 1)
                        {
                                d = 3;
                        }
                
                        N = d + 2 + (d % 2 == 0);
                }

                w = ones(N, 1);
        }

        d = castint(d);
        N = castint(N);

        /* weights */
        w = diag(w);

        if (length(w) != N) error("sgolay - weights must be the same length as polynomial length");

        if (d >= N) error("sgolay - degree must less than polynomial length");

        M = (N - 1) / 2;

        S = fliplr(vander(-M..M, d + 1));

        /* G = S /^ ((S' *^ w) *^ S), but use compact QR for better accuracy */

        (q, r) = qr(sqrt(w) *^ S, 0);

        r = inv(r);

        /* derivative filters */
        G = S *^ r *^ r';

        /* smoothing filters */
        B = G *^ S' *^ w;

        return(B, G);
}