www.gusucode.com > signal 案例源码程序 matlab代码 > signal/HampelFilterAlgorithmExample.m

    %% Eliminate Outliers Using Hampel Identifier
% This example shows a naive implementation of the procedure used by
% |hampel| to detect and remove outliers. The actual function is much
% faster.
% 
% Generate a random signal, |x|, containing 24 samples. Reset the random
% number generator for reproducible results.

% Copyright 2015 The MathWorks, Inc.


rng default

lx = 24;
x = randn(1,lx);

%%
% Generate an observation window around each element of |x|. Take |k| = 2
% neighbors at either side of the sample. The moving window that results
% has a length of $2\times2+1=5$ samples.

k = 2;

iLo = (1:lx)-k;
iHi = (1:lx)+k;

%%
% Truncate the window so that the function computes medians of smaller
% segments as it reaches the signal edges.

iLo(iLo<1) = 1;
iHi(iHi>lx) = lx;

%%
% Record the median of each surrounding window. Find the median of the
% absolute deviation of each element with respect to the window median.

for j = 1:lx
    w = x(iLo(j):iHi(j));
    medj = median(w);
    mmed(j) = medj;
    mmad(j) = median(abs(w-medj));
end

%%
% Scale the median absolute deviation with
% 
% $${{\displaystyle1}\over\displaystyle\sqrt2\,{\mathop{\rm
% erf}\nolimits}^{-1}({1/2})}\approx1.4826$$
% 
% to obtain an estimate of the standard deviation of a normal distribution. 

sd = mmad/(erfinv(1/2)*sqrt(2));

%% 
% Find the samples that differ from the median by more than |nd| = 2
% standard deviations. Replace each of those outliers by the value of the
% median of its surrounding window. This is the essence of the Hampel
% algorithm.

nd = 2;
ki = abs(x-mmed) > nd*sd;

yu = x;
yu(ki) = mmed(ki);

%% 
% Use the |hampel| function to compute the filtered signal and annotate 
% the outliers. Overlay the filtered values computed in this example.

hampel(x,k,nd)

hold on
plot(yu,'o')
hold off