View Raw SPL
/*****************************************************************************
*                                                                            *
*   LINMEDIAN.SPL Copyright (C) 2018 DSP Development Corporation             *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Moving mediane with zero phase                              *
*                                                                            *
*   Revisions:   23 Feb 2018  RRR  Creation                                  *
*                                                                            *
*****************************************************************************/


#if @HELP_LINMEDIAN

    LINMEDIAN

    Purpose: Moving median with zero phase shift

    Syntax:  LINMEDIAN(s, n, ramp)

              s - A series, the input data

              n - An integer, number of points to average as the series is
                  processed

           ramp - Optional, an integer. Ramp the calculation at the endpoints:

                    0: do not ramp
                    1: ramp up at the beginning and down at the end (default)

    Returns: A series

    Example:
             W1: 1..5
             W2: linmedian(w1, 3)

             W2 contains the series:

             {2, 2.25, 2.5, 2.75, 3)

             The X offset of the result is 1.0.

   Example:
             W1: integ(gnorm(1000, 1/100))
             W2: movmedian(w1, 10)
             W3: linmedian(w1, 10);overp(w1, lred);overp(w2, lgreen)

             W1 contains 1000 samples of synthesized data.
             W2 performs a 10 point standard moving median.
             W3 performs a phaseless moving median by revering the original
             data, computing a 10 point moving median, reversing the result
             and computing another 10 point moving median. Compared to
             the standard moving median, the peaks of the resulting smoothed 
             series line up with the original data.

    Remarks:
             LINMEDIAN computes a phaseless moving median by reversing the
             input series, computing an N point moving median, reversing
             the result and computing another N point moving median. The
             reversal steps help ensure that the peak locations of the
             original data are preserved.

    See Also:
             Linavg
             Linexpavg
             Movavg
             Movavg2
             Movmedian
#endif


/* compute a phaseless moving median */
ITERATE linmedian(s, n, ramp)
{
        local y;

        if (argc < 3)
        {
                if (argc < 2)
                {
                        if (argc < 1) error("linmedian - input series required");
                        
                        n = 10;
                }
                
                ramp = 1;
        }

        if (n > 1)
        {
                /* phaseless computation by time reversal */
                y = movmedian(rev(movmedian(rev(s), n, ramp)), n, ramp);

                /* remove start and end points */
                y = extract(y, n, length(s), xoffset(s));

                if (isxy(s))
                {
                        /* restore X values */
                        y = xy(xvals(s), y);
                }

                return(y);
        }
        else
        {
                return(s);
        }
}