SignalProcessing
Hampel
apply the Hampel filter to a 1-D signal to remove outliers
Calling Sequence
Parameters
Options
Description
Examples
Compatibility
Hampel( signal )
Hampel( signal, radius )
Hampel( signal, radius, numstddevs )
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.
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.
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:
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.
Compute the median mu and median absolute deviation mad of the window.
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.
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.
withSignalProcessing:
Example 1
A≔5,5,10,5,5,0,5,5
B≔HampelA,1,0.5,output=hampelsignal,outliers
B≔5.5.5.5.5.5.5.5.,36
Example 2
Consider the following signal:
n≔51:
X≔Vectorn,i↦5+cos4⋅π⋅i−1n−1,datatype=float8
X≔6.5.968583161128635.876306680043865.728968627421415.535826794979005.309016994374955.062790519529314.812618685414284.574220708434934.36257601025131⋮51 element Vector[column]
To this, add outliers:
Y≔copyX:
Y3≔X3+4.0:
Yfloorn2≔Xfloorn2+2.5:
Y−2≔X−2−3.0:
Y=Y
Y=6.5.968583161128639.876306680043865.728968627421415.535826794979005.309016994374955.062790519529314.812618685414284.574220708434934.36257601025131⋮51 element Vector[column]
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≔HampelY,3,2.0,output=record
H≔Recordhampelsignal=6.5.968583161128635.848775894275025.728968627421415.535826794979005.309016994374955.062790519529314.812618685414284.574220708434934.36257601025131⋮51 element Vector[column],medians=5.984291580564325.968583161128635.848775894275025.728968627421415.535826794979005.309016994374955.062790519529314.812618685414284.574220708434934.36257601025131⋮51 element Vector[column],outliers=32550,standarddeviations=0.2009158571348160.3552530392605080.3440921117674970.4018317142696330.6416055485258700.6226212228197380.7013246314153260.6672342686188060.5681890684009100.433442459362165⋮51 element Vector[column]
For the plot, we require only the filtered signal:
Z≔ArrayTools:-AliasHhampelsignal:
Finally, plot all three signals:
dataplotX,Y,Z,view=1..n,0..10,labels=time,amplitude,title=Signals,font=Verdana,15,legend=original signal,original signal with outliers,filtered signal,legendstyle=font=Verdana,15,symbolsize=5,color=green,red,blue,thickness=6,3,3
The SignalProcessing[Hampel] command was introduced in Maple 2021.
For more information on Maple 2021 changes, see Updates in Maple 2021.
See Also
ArrayTools
dataplot
erfc
Download Help Document
What kind of issue would you like to report? (Optional)
