matlab使用hampel濾波,去除異常值

原文連接：http://tecdat.cn/?p=7181

此示例顯示了Hampel用於檢測和刪除異常值的過程的 實現。 

產生一個包含24個樣本的隨機信號x。 重置隨機數生成器以得到可重複的結果。

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

圍繞x的每一個元素生成觀察窗口。 在樣本的任一邊取k = 2個鄰居。 產生的移動窗口的長度爲2×2 + 1 = 5個樣本。

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

截斷窗口，以便函數在到達信號邊緣時計算較小段的中值。

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

記錄每一個周圍窗口的中位數。 找到每一個元素相對於窗口中位數的絕對誤差的中位數。

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

縮放中位數絕對誤差

1G2erf−1(1/2)≈1.4826

以得到正態分佈標準誤差的估計值。

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

查找與中位數相差超過nd = 2個標準誤差的樣本。 將這些離羣值替換爲其周圍窗口的中間值。 這是Hampel算法的本質。

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

使用hampel 計算濾波後的信號並註釋異常值。 覆蓋在此示例中計算的過濾值。

plot(yu,'o','HandleVisibility','off')