/*****************************************************************************
*                                                                            *
*   DETREND.SPL   Copyright (C) 1996, 2011 DSP Development Corporation       *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:       Randy Race                                                 *
*                                                                            *
*   Synopsis:     Remove a linear trend from a series or table               *
*                                                                            *
*   Revisions:    18 Jan 1996  RRR  Creation                                 *
*                 11 Feb 2011  RRR  linear breakpoint option                 *
*                                                                            *
*****************************************************************************/


#if @HELP_DETREND

    DETREND

    Purpose: Removes a linear trend line or mean value from a series.

    Syntax:  DETREND(x, opt, bp)

                  x - A series, the input.

                opt - Optional, a string, the trend type.

                        "linear"   : remove a linear trend (default)
                        "constant" : remove the mean value

                 bp - Optional, a series. A list of breakpoint indices
                      that determine the locations of a piecewise
                      continuous linear trend. Defaults to a single
                      regression line over the entire input series.

    Returns: A series, the detrended output.

    Example:
             W1: gnorm(10, 1)
             W2: W1 + 100
             W3: detrend(w2, "constant")
             W4: demean(w2)

             W2 contains a 10 point random series with a mean of 10.
             Both W3 and W4 remove the mean value to recover the
             original series in W1

    Example:
             W1: gcos(1000, 1/1000, 10)
             W2: 10 * gtri(length(w1), deltax(w1), 1)
             W3: w1 + w2
             W4: detrend(w3, maxidx(w2))

             W3 contains a 1000 point cosine with a triangular linear
             trend.  W4 removes the trend and recovers W1 by setting a
             breakpoint at the index of the apex of the triangle wave.

    Remarks:
             If breakpoints are specified, each point represents the
             index shared by adjacent line segments. A piecewise
             continuous regression line is computed for each segment.

             If no breakpoints are specified, a trendline computed over
             the entire input series is removed. This is equivalent to:

                  x - linreg(x)

    See Also:
             Demean
             Linreg
             Linfit
             Polyit
#endif


/* remove a linear trend from a series */
ITERATE detrend(x, opt, bp)
{
        local y, xv, xlen, blen, slen, n, a;

        if (argc < 3) 
        {
                if (argc < 2) 
                {
                        if (argc < 1) error("detrend - input series required");
                        opt = "linear";
                }
                
                bp = {};
                
                if (isscalar(opt))
                {
                        /* detrend(x, 10) */
                        bp  = {opt};
                        opt = "linear";
                }
        }

        if (isarray(opt))
        {
                /* detrend(x, bp) */
                bp  = opt;
                opt = "linear";
        }

        if (tolower(opt) != "linear") 
        {
                /* just remove mean */
                y = demean(x);
        }
        else 
        {
                wasxy = isxy(x);

                /* remove piecewise continuous linear trend */
                if (length(bp) > 0) 
                {
                        xlen = length(x);

                        /* remove dups */
                        bp   = detrend_undup({0, bp, xlen-1});
                        blen = length(bp) - 1;

                        a = ravel(zeros(xlen, blen), ones(xlen, 1));

                        /* form regression lines */
                        loop (n = 1..blen)
                        {
                                slen = xlen - bp[n];
                                a[(1..slen) + bp[n], n] = (1..slen) / slen;
                        }

                        xv = yvals(x);

                        /* remove regression */
                        y = xv - a *^ (a \^ xv);

                        /* attributes */
                        y.deltax  = x.deltax;
                        y.xoffset = x.xoffset;
                        y.vunits  = x.vunits;
                        y.hunits  = x.hunits;
                }
                else
                {
                        /* remove full trendline */
                        y = x - linreg(x);
                }

                if (isxy(x))
                {
                        /* reconstruct XY */
                        y = xy(xvals(x), y);
                }
        }

        return(y);

}


/* sort and remove duplicates */
detrend_undup(s)
{
        local r, d;

        r = sort(s, 1);
        d = extract(r - delay(r, 1), 2, -1);
        r = delete(r, d == 0);

        return(r);
}