matlab将孤立点连接起来,孤立点检测与去除

Outlier Detection and Removal [hampel]

HAMPEL(X,Y,DX,T,varargin) returns the Hampel filtered values of the

elements in Y. It was developed to detect outliers in a time series,

but it can also be used as an alternative to the standard median

filter.

References

Chapters 1.4.2, 3.2.2 and 4.3.4 in Mining Imperfect Data: Dealing with

Contamination and Incomplete Records by Ronald K. Pearson.

Acknowledgements

I would like to thank Ronald K. Pearson for the introduction to moving

window filters. Please visit his blog at:

http://exploringdatablog.blogspot.com/2012/01/moving-window-filters-and-pracma.html

X,Y are row or column vectors with an equal number of elements.

The elements in Y should be Gaussian distributed.

Input DX,T,varargin must not contain NaN values!

DX,T are optional scalar values.

DX is a scalar which defines the half width of the filter window.

It is required that DX > 0 and DX should be dimensionally equivalent to

the values in X.

T is a scalar which defines the threshold value used in the equation

|Y - Y0| > T*S0.

Standard Parameters for DX and T:

DX = 3*median(X(2:end)-X(1:end-1));

T = 3;

varargin covers addtional optional input. The optional input must be in

the form of 'PropertyName', PropertyValue.

Supported PropertyNames:

'standard': Use the standard Hampel filter.

'adaptive': Use an experimental adaptive Hampel filter. Explained under

Revision 1 details below.

Supported PropertyValues: Scalar value which defines the tolerance of

the adaptive filter. In the case of standard Hampel filter this value

is ignored.

Output YY,I,Y0,LB,UB,ADX are column vectors containing Hampel filtered

values of Y, a logical index of the replaced values, nominal data,

lower and upper bounds on the Hampel filter and the relative half size

of the local window, respectively.

NO is a scalar that specifies the Number of Outliers detected.

Examples

1. Hampel filter removal of outliers

X = 1:1000; % Pseudo Time

Y = 5000 + randn(1000, 1); % Pseudo Data

Outliers = randi(1000, 10, 1); % Index of Outliers

Y(Outliers) = Y(Outliers) + randi(1000, 10, 1); % Pseudo Outliers

[YY,I,Y0,LB,UB] = hampel(X,Y);

plot(X, Y, 'b.'); hold on; % Original Data

plot(X, YY, 'r'); % Hampel Filtered Data

plot(X, Y0, 'b--'); % Nominal Data

plot(X, LB, 'r--'); % Lower Bounds on Hampel Filter

plot(X, UB, 'r--'); % Upper Bounds on Hampel Filter

plot(X(I), Y(I), 'ks'); % Identified Outliers

2. Adaptive Hampel filter removal of outliers

DX = 1; % Window Half size

T = 3; % Threshold

Threshold = 0.1; % AdaptiveThreshold

X = 1:DX:1000; % Pseudo Time

Y = 5000 + randn(1000, 1); % Pseudo Data

Outliers = randi(1000, 10, 1); % Index of Outliers

Y(Outliers) = Y(Outliers) + randi(1000, 10, 1); % Pseudo Outliers

[YY,I,Y0,LB,UB] = hampel(X,Y,DX,T,'Adaptive',Threshold);

plot(X, Y, 'b.'); hold on; % Original Data

plot(X, YY, 'r'); % Hampel Filtered Data

plot(X, Y0, 'b--'); % Nominal Data

plot(X, LB, 'r--'); % Lower Bounds on Hampel Filter

plot(X, UB, 'r--'); % Upper Bounds on Hampel Filter

plot(X(I), Y(I), 'ks'); % Identified Outliers

3. Median Filter Based on Filter Window

DX = 3; % Filter Half Size

T = 0; % Threshold

X = 1:1000; % Pseudo Time

Y = 5000 + randn(1000, 1); % Pseudo Data

[YY,I,Y0] = hampel(X,Y,DX,T);

plot(X, Y, 'b.'); hold on; % Original Data

plot(X, Y0, 'r'); % Median Filtered Data