hampel滤波，去除异常值

hampel滤波，去除异常值

本文源自(https://www.mathworks.com/help/signal/ref/hampel.html)

语法

y = hampel(x)

对输入向量x进行hampel滤波，检测和删除异常值。对于x的每个样本，该函数计算由样本及其周围六个样本组成的窗口的中值，每边三个。并利用中位数绝对值估计了各样本对中值的标准差。如果某个样本与中值相差超过三个标准差，则用中值替换该样本。如果x是一个矩阵，hampel将x的每一列都看成是独立的通道。

y = hampel(x,k)

可指定窗口中每个样本周围的样本数，每边k个。

y = hampel(x,k,nsigma)

可指定几倍的标准差，nsigma指一个样本的值与中值的差值

[y,j] = hampel(___)

返回一个逻辑矩阵，该矩阵在作为异常值标识的所有点的位置都为真。此语法接受来自以前语法的任何输入参数。

[y,j,xmedian,xsigma] = hampel(___)

返回x的每个元素的局部中值和估计标准差。

hampel(___)

没有输出参数的hampel()绘制经过过滤的信号，并标注被删除的异常值。

具体例子

1. 移除正弦波上的尖刺

生成100个正弦信号样本。用尖峰代替第6和第20个样品。

x = sin(2pi(0:99)/100); x(6) = 2; x(20) = -2;

使用hampel定位每个与局部中值相差超过三个标准差的样本。测量窗口由样品及其周围的六个样品组成，每边三个。

[y,i,xmedian,xsigma] = hampel(x);

绘制经过滤波的信号并标注异常值

n = 1:length(x);
plot(n,x)
hold on
plot(n,xmedian-3xsigma,n,xmedian+3xsigma)
plot(find(i),x(i),‘sk’)
hold off
legend(‘Original signal’,‘Lower limit’,‘Upper limit’,‘Outliers’)

重新计算，但现在在计算中值时，每边只取一个相邻的样本。该函数将极值视为异常值。

hampel(x,1)

2. 多通道信号的hampel滤波

生成由不同频率的正弦信号组成的双通道信号。在随机位置放置尖峰。使用NaNs随机添加缺失的样本。重置随机数生成器以获得可重复的结果。画出信号。

rng(‘default’)
n = 59;
x = sin(pi./[15 10]’(1:n)+pi/3)’;
spk = randi(2n,9,1);
x(spk) = x(spk)2;
x(randi(2n,6,1)) = NaN;
plot(x)

使用带有默认设置的hampel过滤信号。

y = hampel(x);
plot(y)

增加移动窗口的长度，并降低将样本视为异常值的阈值。

y = hampel(x,4,2); %每边4个样本，±2σ
plot(y)

输出每个通道的运行中位数。将中值覆盖在信号图上。

[y,j,xmd,xsd] = hampel(x,4,2);
plot(x)
hold on
plot(xmd,’–’)

3. 找出多通道信号中的异常值

在单位方差的高斯白噪声中产生由两个不同频率的正弦信号组成的多通道信号。

rng(‘default’)
t = 0:60;
x = sin(pi./[10;2]*t)’+randn(numel(t),2);

对信号用hampel滤波器。取那些与周围9个样本窗口的中值相差两个以上标准差的点作为异常值。输出一个逻辑矩阵，该矩阵在异常值的位置为真。

k = 4;
nsig = 2;
[y,h] = hampel(x,k,nsig);

把信号的每个通道画在它自己的一组坐标轴上。绘制原始信号、滤波信号和异常值。标注离群点位置。

for k = 1:2
hk = h(:,k);
ax = subplot(2,1,k);
plot(t,x(:,k))
hold on
plot(t,y(:,k))
plot(t(hk),x(hk,k),’*’)
hold off
ax.XTick = t(hk);
end

4. hampel滤波器返回的信号统计量

生成100个正弦信号样本。用尖峰代替第6和第20个样品。

n = 1:100;
x = sin(2pin/100);
x(6) = 2;
x(20) = -2;

使用hampel计算每个样本的局部中值和估计标准差。使用输入参数的默认值:

• 窗口大小为2×3+1=7
• 与窗口中值相差超过三个标准差的点被认为是异常值。

[y,i,xmedian,xsigma] = hampel(x);

plot(n,x)
hold on
plot(n,[1;1]xmedian+3[-1;1]*xsigma)
plot(find(i),x(i),‘sk’)
hold off
legend(‘Signal’,‘Lower’,‘Upper’,‘Outliers’)

重新计算，使用一个2×10+1=21的窗口和两个标准差作为识别异常值的标准。

sds = 2;
adj = 10;
[y,i,xmedian,xsigma] = hampel(x,adj,sds);
plot(n,x)
hold on
plot(n,[1;1]xmedian+sds[-1;1]*xsigma)
plot(find(i),x(i),‘sk’)
hold off
legend(‘Signal’,‘Lower’,‘Upper’,‘Outliers’)

输入参数

1. x——输入信号，可以是向量或矩阵
2. k——每边的相邻样本数，必须是整数，默认为3
   
3. nsigma——标准差的数量，必须是实标量，默认为3

输出参数

1. y——滤波后的信号，返回一个和x大小一样的向量或矩阵
2. j——离群值指数，返回一个和x大小一样的向量或矩阵，数据类型为逻辑运算符
3. xmedian——样本中位数，返回一个和x大小一样的向量或矩阵
4. xsigma——估计标准差，返回一个和x大小一样的向量或矩阵