hample滤波器的原理及其Python实现

1. 作用及原理

功能：检测并删除异常值
  用一个一维向量 x = { x 1 , x 2 . . . x n } \boldsymbol{x}={\{x_1,x_2...x_n\}} x={x1​,x2​...xn​}表示序列，围绕每个元素生成观测窗口，假设半个窗口宽度为k，整个窗口的宽度为2k+1（包含中心的那个元素），计算该窗口中所有元素的中值，并利用中位数的绝对值估计各样本对中值的标准差.如果某个样本与中值相差超过三个标准差，则用中值替换该样本。

2. Python实现

参考了博主SineyCoder的去除异常值，使用python编写matlab的hampel(X)函数

def hampel(X，k):
    length = X.shape[0] - 1
    nsigma = 3
    iLo = np.array([i - k for i in range(0, length + 1)])
    iHi = np.array([i + k for i in range(0, length + 1)])
    iLo[iLo < 0] = 0
    iHi[iHi > length] = length
    xmad = []
    xmedian = []
    for i in range(length + 1):
        w = X[iLo[i]:iHi[i] + 1]
        medj = np.median(w)
        mad = np.median(np.abs(w - medj))
        xmad.append(mad)
        xmedian.append(medj)
    xmad = np.array(xmad)
    xmedian = np.array(xmedian)
    scale = 1.4826  # 缩放
    xsigma = scale * xmad
    xi = ~(np.abs(X - xmedian) <= nsigma * xsigma)  # 找出离群点（即超过nsigma个标准差）
 
    # 将离群点替换为中为数值
    xf = X.copy()
    xf[xi] = xmedian[xi]
    return xf
 
X = np.array([1, 2, 3, 4, 100, 4, 3, 2, 1])
res = hampel(X，3)
plt.plot(X)
plt.plot(res, '--')
plt.show()