Hampel - Maple Help

SignalProcessing

 Hampel
 apply the Hampel filter to a 1-D signal to remove outliers

 Calling Sequence Hampel( signal ) Hampel( signal, radius ) Hampel( signal, radius, numstddevs )

Parameters

 signal - 1-D rtable or list of data. radius - (optional) positive integer which defines the radius, in number of sample points, of each window used to determine if a point is an outlier. The default is 3. numstddevs - (optional) positive numeric value which specifies the number of standard deviations a sample point must be within the median to not be considered an outlier. The default is 3.0.

Options

 • output : The type of output. The supported options are:
 – hampelsignal: 1-D Array or Vector of the Hampel-filtered signal. This is the default.
 – medians: 1-D Array or Vector containing the medians for each sample point of the original signal.
 – outliers: Vector containing the indices for the outlier sample points of the original signal.
 – standarddeviations: 1-D Array or Vector containing the standard deviations for each sample point of the original signal.
 – record: Returns a record with the previous options.
 – list of any of the above options (excluding record): Returns an expression sequence with the corresponding outputs, in the same order.

Description

 • The Hampel command takes a 1-D data set, and replaces outlier sample points with the median value with respect to neighbouring points.
 • To determine if a given sample point, say with index i, is an outlier, a window of neighbouring points is considered:
 1 Define the window to be the points signal[max(a,i-r)..min(b,i+r)], where a and b are, respectively, the lower and upper indices of signal, and r=radius.
 2 Compute the median mu and median absolute deviation mad of the window.
 3 Approximate the standard deviation of the window as sigma=rho*mad, where rho=-1/sqrt(2)/inverfc(3/2) and inverfc() is the inverse of the complementary error function. The value of rho is approximately 1.483.
 4 If abs(signal[i]-mu)>numstddevs*sigma, that is, the sample point is more than numstddevs standard deviations from the median, the point is labelled an outlier.
 • The filtered signal is formed by taking the original signal, and replacing the outliers their respective medians.
 • Maple will attempt to coerce the provided signal to a 1-D Array or Vector of float[8] datatype, and an error will be thrown if this is not possible. For this reason, it is most efficient for the input to already be a float[8] 1-D Array or Vector.
 • The input signal cannot have an indexing function, and must use rectangular storage.
 • The Hampel command is not thread safe.
 • As the underlying implementation of the SignalProcessing package is a module, it is also possible to use the form SignalProcessing:-Hampel to access the command from the package. For more information, see Module Members.

Examples

 > $\mathrm{with}\left(\mathrm{SignalProcessing}\right):$

Example 1

 > $A≔\left[5,5,10,5,5,0,5,5\right]$
 ${A}{≔}\left[{5}{,}{5}{,}{10}{,}{5}{,}{5}{,}{0}{,}{5}{,}{5}\right]$ (1)
 > $B≔\mathrm{Hampel}\left(A,1,0.5,\mathrm{output}=\left[\mathrm{hampelsignal},\mathrm{outliers}\right]\right)$
 ${B}{≔}\left[\begin{array}{c}{5.}\\ {5.}\\ {5.}\\ {5.}\\ {5.}\\ {5.}\\ {5.}\\ {5.}\end{array}\right]{,}\left[\begin{array}{c}{3}\\ {6}\end{array}\right]$ (2)

Example 2

 • Consider the following signal:
 > $n≔51:$
 > $X≔\mathrm{Vector}\left(n,i↦5+\mathrm{cos}\left(\frac{4\cdot \mathrm{\pi }\cdot \left(i-1\right)}{n-1}\right),\mathrm{datatype}=\mathrm{float}\left[8\right]\right)$
 ${X}{≔}\begin{array}{c}\left[\begin{array}{c}{6.}\\ {5.96858316112863}\\ {5.87630668004386}\\ {5.72896862742141}\\ {5.53582679497900}\\ {5.30901699437495}\\ {5.06279051952931}\\ {4.81261868541428}\\ {4.57422070843493}\\ {4.36257601025131}\\ {⋮}\end{array}\right]\\ \hfill {\text{51 element Vector[column]}}\end{array}$ (3)
 > $Y≔\mathrm{copy}\left(X\right):$
 > $Y\left[3\right]≔X\left[3\right]+4.0:$
 > $Y\left[\mathrm{floor}\left(\frac{n}{2}\right)\right]≔X\left[\mathrm{floor}\left(\frac{n}{2}\right)\right]+2.5:$
 > $Y\left[-2\right]≔X\left[-2\right]-3.0:$
 > $'Y'=Y$
 ${Y}{=}\begin{array}{c}\left[\begin{array}{c}{6.}\\ {5.96858316112863}\\ {9.87630668004386}\\ {5.72896862742141}\\ {5.53582679497900}\\ {5.30901699437495}\\ {5.06279051952931}\\ {4.81261868541428}\\ {4.57422070843493}\\ {4.36257601025131}\\ {⋮}\end{array}\right]\\ \hfill {\text{51 element Vector[column]}}\end{array}$ (4)
 • Now, compute the Hampel-filtered signal Z of X with, say, radius of 3 sample points and 2.0 standard deviations. Here, we return the full record:
 > $H≔\mathrm{Hampel}\left(Y,3,2.0,\mathrm{output}=\mathrm{record}\right)$
 ${H}{≔}{\mathrm{Record}}{}\left({\mathrm{hampelsignal}}{=}\begin{array}{c}\left[\begin{array}{c}{6.}\\ {5.96858316112863}\\ {5.84877589427502}\\ {5.72896862742141}\\ {5.53582679497900}\\ {5.30901699437495}\\ {5.06279051952931}\\ {4.81261868541428}\\ {4.57422070843493}\\ {4.36257601025131}\\ {⋮}\end{array}\right]\\ \hfill {\text{51 element Vector[column]}}\end{array}{,}{\mathrm{medians}}{=}\begin{array}{c}\left[\begin{array}{c}{5.98429158056432}\\ {5.96858316112863}\\ {5.84877589427502}\\ {5.72896862742141}\\ {5.53582679497900}\\ {5.30901699437495}\\ {5.06279051952931}\\ {4.81261868541428}\\ {4.57422070843493}\\ {4.36257601025131}\\ {⋮}\end{array}\right]\\ \hfill {\text{51 element Vector[column]}}\end{array}{,}{\mathrm{outliers}}{=}\left[\begin{array}{c}{3}\\ {25}\\ {50}\end{array}\right]{,}{\mathrm{standarddeviations}}{=}\begin{array}{c}\left[\begin{array}{c}{0.200915857134816}\\ {0.355253039260508}\\ {0.344092111767497}\\ {0.401831714269633}\\ {0.641605548525870}\\ {0.622621222819738}\\ {0.701324631415326}\\ {0.667234268618806}\\ {0.568189068400910}\\ {0.433442459362165}\\ {⋮}\end{array}\right]\\ \hfill {\text{51 element Vector[column]}}\end{array}\right)$ (5)
 • For the plot, we require only the filtered signal:
 > $Z≔\mathrm{ArrayTools}:-\mathrm{Alias}\left(H\left[\mathrm{hampelsignal}\right]\right):$
 • Finally, plot all three signals:
 > $\mathrm{dataplot}\left(\left[X,Y,Z\right],\mathrm{view}=\left[1..n,0..10\right],\mathrm{labels}=\left[\mathrm{time},\mathrm{amplitude}\right],\mathrm{title}="Signals",\mathrm{font}=\left[\mathrm{Verdana},15\right],\mathrm{legend}=\left["original signal","original signal with outliers","filtered signal"\right],\mathrm{legendstyle}=\left[\mathrm{font}=\left[\mathrm{Verdana},15\right]\right],\mathrm{symbolsize}=5,\mathrm{color}=\left["green","red","blue"\right],\mathrm{thickness}=\left[6,3,3\right]\right)$

Compatibility

 • The SignalProcessing[Hampel] command was introduced in Maple 2021.