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);
}
}