www.gusucode.com > signal 工具箱matlab源码程序 > signal/hampel.m

    function [xfilt, xi, xmedian, xsigma] = hampel(x, k, nsigma)
%HAMPEL Outlier removal via Hampel identifier.
%   Y = HAMPEL(X) replaces any element in vector X that is more than three
%   standard deviations from the median of itself and up to three
%   neighboring elements with that median value.  The standard deviation is
%   estimated by scaling the local median absolute deviation (MAD) by a
%   constant scale factor.  If X is a matrix, HAMPEL operates over the
%   columns of X.
% 
%   Y = HAMPEL(X,K) specifies the number of adjacent samples, K, on either
%   side of each sample in X over which to compute the Hampel identifier.
%   The default value of K is 3.
% 
%   Y = HAMPEL(X,K,NSIGMA) specifies the number of estimated standard
%   deviations above which it will replace elements in X with the local
%   median.
% 
%   [Y,I] = HAMPEL(...) returns a logical matrix, I, of the same shape as X
%   which is true when the corresponding element in X is identified as an
%   outlier.
% 
%   [Y,I,XMEDIAN,XSIGMA] = HAMPEL(...) additionally returns the local medians
%   and estimated standard deviations (scaled MAD) for each element in X.
% 
%   HAMPEL(...) without output arguments plots the filtered signal as well as
%   the outliers removed from the signal.
% 
%   % Example 1:
%   %   Remove spikes from a sinusoid
%   x = sin(2*pi*(0:99)/100);
%   x(6) = 2;
%   x(20) = -2;
%   hampel(x)
%
%   % Example 2:
%   %   plot statistics returned by hampel identifier
%   x = sin(2*pi*(0:99)/100);
%   x(6) = 2;
%   x(20) = -2;
%   [y,i,xmedian,xsigma] = hampel(x);
%   n = numel(x);
%   plot(1:n,xmedian-3*xsigma,'r', ...
%        1:n,xmedian+3*xsigma,'r', ...
%        1:n,x, ...
%        find(i),x(i),'sk ')
%   legend('lower limit','upper limit','original signal','outliers');
%
%   See also MEDFILT1, MEDIAN, FILTER, SGOLAYFILT.

%   Copyright 2015 The MathWorks, Inc.

narginchk(1,3);

if nargin<2
  k = 3;
else
  validateattributes(k,{'numeric'},{'integer','scalar','positive'});
end

% default is three standard deviations of a gaussian distributed input
if nargin<3
  nsigma=3;
else
  validateattributes(nsigma,{'numeric'},{'real','scalar','nonnegative'});
end

validateattributes(x,{'single','double'},{'real','2d','nonsparse'});

needsTranspose = isrow(x);
if needsTranspose
  x = x(:);
end

% compute the median absolute deviations and the corresponding medians over
% the size of the filter:  ignore samples that contain NaN and truncate
% at the borders of the input
[xmad,xmedian] = movmad(x,2*k+1,1,'omitnan','central');

% scale the MAD by ~1.4826 as an estimate of its standard deviation
scale = -1 /(sqrt(2)*erfcinv(3/2));
xsigma = scale*xmad;

% identify points that are either NaN or beyond the desired threshold
xi = ~(abs(x-xmedian) <= nsigma*xsigma);

% replace identified points with the corresponding median value
xf = x;
xf(xi) = xmedian(xi);

if nargout==0
  plotOutliers(x, xf, xi)
else
  xfilt = xf;
  if needsTranspose
    xfilt = xf.';
    xi = xi.';
    xmedian = xmedian.';
    xsigma = xsigma.';
  end    
end

function plotOutliers(x, xf, xi)
t = (1:size(x,1))';
if size(x,2)==1
  hLine = plot(t,x,'.-', ...
               t,xf, ...
               t(xi),x(xi),'ks');
  if numel(hLine)<3
    legend(getString(message('signal:hampel:OriginalSignal')), ...
           getString(message('signal:hampel:FilteredSignal')));
  else
    legend(getString(message('signal:hampel:OriginalSignal')), ...
           getString(message('signal:hampel:FilteredSignal')), ...
           getString(message('signal:hampel:Outliers')));
  end
else
  colors=get(0,'DefaultAxesColorOrder');
  for i=1:size(x,2)
    if i==1
      hLine = plot(t,xf(:,i),'Color',colors(1+mod(i-1,size(colors,1)),:));
      hLineOutliers = line(t(xi(:,i)),x(xi(:,i),i),'LineStyle','none','Marker','s','MarkerEdgeColor',hLine.Color);
    else
      hLine_next = line(t,xf(:,i),'Color',colors(1+mod(i-1,size(colors,1)),:));
      line(t(xi(:,i)),x(xi(:,i),i),'LineStyle','none','Marker','s','MarkerEdgeColor',hLine_next.Color);
    end
  end
  legend([hLine hLineOutliers],getString(message('signal:hampel:FilteredSignal')), ...
                                  getString(message('signal:hampel:Outliers')));

end