hampel

햄펄(Hampel) 식별자를 사용하여 이상값 제거

구문

y = hampel(x)

y = hampel(x,k)

y = hampel(x,k,nsigma)

[y,j] = hampel(___)

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

hampel(___)

설명

y = hampel(x)는 햄펄 필터를 입력 벡터 x에 적용하여 이상값을 감지하고 제거합니다. x의 각 샘플에 대해 함수는 그 샘플과 주변 샘플 6개(측면당 3개)로 구성된 윈도우의 중앙값을 계산합니다. 또한 중앙값 절대 편차를 사용하여 각 샘플의 윈도우 중앙값에 대한 해당 샘플의 표준편차도 추정합니다. 중앙값과 3 표준편차 넘게 차이가 나면 샘플은 중앙값으로 대체됩니다. x가 행렬이면 함수는 x의 각 열을 독립적인 채널로 처리합니다.

y = hampel(x,k)는 측정 윈도우에 있는 x의 각 샘플의 양쪽 측면에 있는 이웃 수 k를 지정합니다. k의 디폴트 값은 3입니다.

y = hampel(x,k,nsigma)는 표준편차 개수 nsigma를 지정합니다. x의 샘플이 중앙값으로 대체되려면 국소 중앙값과 이 값만큼 차이가 나야 합니다. nsigma의 디폴트 값은 3입니다.

[y,j] = hampel(___)은 이상값으로 식별된 모든 지점의 위치에서 true인 논리형 행렬도 반환합니다. 이 구문은 위에 열거된 구문의 모든 입력 인수를 받습니다.

[y,j,xmedian,xsigma] = hampel(___)은 국소 중앙값과 x의 각 요소에 대한 표준편차 추정값도 반환합니다.

hampel(___)에 출력 인수를 지정하지 않으면 필터링된 신호를 플로팅하고, 제거된 이상값에 주석을 표시합니다.

예제

정현파에서 스파이크 제거하기

100개의 정현파 신호 샘플을 생성합니다. 6번째 샘플과 20번째 샘플을 스파이크로 대체합니다.

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

hampel을 사용하여, 국소 중앙값과 3 표준편차 넘게 차이가 나는 샘플을 모두 찾습니다. 측정 윈도우는 샘플과 그 주변 샘플 6개(한 쪽당 3개)로 구성됩니다.

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

필터링된 신호를 플로팅하고 이상값에 주석을 표시합니다.

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

계산을 반복하되 이번에는 중앙값을 계산할 때 각 측면에서 인접 샘플을 하나만 취합니다. 함수는 극값을 이상값으로 간주합니다.

hampel(x,1)

다중채널 신호에 햄펄 필터 적용하기

주파수가 서로 다른 정현파로 구성된 2채널 신호를 생성합니다. 임의의 위치에 스파이크를 배치합니다. NaN을 사용하여 누락된 샘플을 임의로 추가합니다. 재현 가능한 결과를 얻기 위해 난수 생성기를 재설정합니다. 신호를 플로팅합니다.

rng('default')

n = 59;
x = sin(pi./[15 10]'*(1:n)+pi/3)';

spk = randi(2*n,9,1);
x(spk) = x(spk)*2;
x(randi(2*n,6,1)) = NaN;

plot(x)

디폴트 설정으로 hampel을 사용하여 신호에 필터를 적용합니다.

y = hampel(x);
plot(y)

샘플을 이상값으로 처리하기 위해 이동 윈도우의 길이는 늘리고 임계값은 줄입니다.

y = hampel(x,4,2);
plot(y)

각 채널의 이동 중앙값을 출력합니다. 신호의 플롯에 중앙값을 겹쳐 표시합니다.

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

다중채널 신호의 이상값 찾기

단위 분산의 백색 가우스 잡음에 묻힌 서로 다른 주파수의 두 정현파로 구성된 다중 채널 신호를 생성합니다.

rng('default')

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

햄펄 필터를 신호에 적용합니다. 샘플 9개로 구성된 주변 윈도우의 중앙값에서 2 표준편차 넘게 차이가 나는 지점을 이상값으로 간주합니다. 이상값 위치에서 true인 논리형 행렬을 출력합니다.

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

햄펄 필터(Hampel Filter)에서 반환된 신호 통계량

100개의 정현파 신호 샘플을 생성합니다. 6번째 샘플과 20번째 샘플을 스파이크로 대체합니다.

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

hampel을 사용하여 각 샘플에 대한 국소 중앙값과 표준편차 추정값을 계산합니다. 입력 파라미터의 디폴트 값을 사용합니다.

• 윈도우 크기는 $$2\times 3+1=7$$입니다.

• 윈도우 중앙값과 3 표준편차 넘게 차이가 나는 지점이 이상값으로 간주됩니다.

결과를 플로팅합니다.

[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\times 10+1=21$$과 2 표준편차를 이상값 식별 기준으로 사용하여 계산을 반복합니다.

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')

입력 인수

x — 입력 신호

벡터 | 행렬

입력 신호로, 벡터나 행렬로 지정됩니다. x가 행렬이면 hampel은 x의 각 열을 독립적인 채널로 처리합니다.

예: cos(pi/4*(0:159))+randn(1,160)은 단일채널 행 벡터 신호입니다.

예: cos(pi./[4;2]*(0:159))'+randn(160,2)는 2채널 신호입니다.

데이터형: single | double

k — 양쪽 측면에 있는 이웃 수

3 (디폴트 값) | 정수 스칼라

샘플 x_s의 양쪽 측면에 있는 이웃 수로, 정수 스칼라로 지정됩니다. 신호 경계에 근접하고 한 측면의 샘플 개수가 k 미만인 샘플은 더 작은 윈도우의 중앙값과 비교됩니다.

데이터형: single | double

nsigma — 표준편차 개수

3 (디폴트 값) | 실수형 스칼라

x의 샘플이 이상값으로 간주되기 위해 국소 중앙값과 차이가 나야 하는 표준편차 개수입니다. nsigma를 실수형 스칼라로 지정하십시오. 함수는 국소 중앙값 절대 편차(MAD)를 인자 $$\kappa =\frac{1}{\sqrt{2}{\mathrm{erf}}^{-1}\left(1/2\right)}\approx 1.4826$$으로 스케일링하여 표준편차를 추정합니다.

데이터형: single | double

출력 인수

y — 필터링된 신호

벡터 | 행렬

필터링된 신호로, x와 크기가 동일한 벡터 또는 행렬로 반환됩니다.

데이터형: single | double

j — 이상값 인덱스

벡터 | 행렬

이상값 인덱스로, x와 크기가 동일한 벡터 또는 행렬로 반환됩니다.

데이터형: logical

xmedian — 국소 중앙값

벡터 | 행렬

국소 중앙값으로, x와 크기가 동일한 벡터 또는 행렬로 반환됩니다.

데이터형: single | double

xsigma — 표준편차 추정값

벡터 | 행렬

표준편차 추정값으로, x와 크기가 동일한 벡터 또는 행렬로 반환됩니다.

데이터형: single | double

세부 정보

햄펄 식별자

햄펄 식별자는 이상값의 영향을 적게 받는 로버스트 통계의 3시그마 규칙에 대한 변형입니다.

시퀀스 x₁, x₂, x₃, …, x_n과 길이가 k인 슬라이딩 윈도우가 주어진 경우, 다음을 사용하여 지점 간 중앙값과 표준편차 추정값을 정의합니다.

• 국소 중앙값 — $$m_i=\mathrm{median}\left(x_{i-k},x_{i-k+1},x_{i-k+2},\dots ,x_i,\dots ,x_{i+k-2},x_{i+k-1},x_{i+k}\right)$$

• 표준편차 — $${\sigma }_i=\kappa \mathrm{median}\left(\left|x_{i-k}-m_i\right|,\dots ,\left|x_{i+k}-m_i\right|\right),$$, 여기서 $$\kappa =\frac{1}{\sqrt{2}{\mathrm{erf}}^{-1}\left(1/2\right)}\approx 1.4826$$

수량 σ_i /κ 는 중앙값 절대 편차(MAD)라고 합니다.

샘플 x_i가 주어진 임계값 n_σ에 대해 다음과 같을 경우

$$\left|x_i-m_i\right|>n_{\sigma }{\sigma }_i$$

햄펄 식별자는 x_i를 이상값으로 선언하고 m_i로 대체합니다.

함수는 시퀀스 끝점 근처에서 m_i와 σ_i를 계산하는 데 사용된 윈도우를 자릅니다.

• i < k + 1

  $$m_i=\mathrm{median}\left(x_1,x_2,x_3,\dots ,x_i,\dots ,x_{i+k-2},x_{i+k-1},x_{i+k}\right)$$

  $${\sigma }_i=\kappa \mathrm{median}\left(\left|x_1-m_1\right|,\dots ,\left|x_{i+k}-m_i\right|\right)$$

• i > n – k

  $$m_i=\mathrm{median}\left(x_{i-k},x_{i-k+1},x_{i-k+2},\dots ,x_i,\dots ,x_{n-2},x_{n-1},x_n\right)$$

  $${\sigma }_i=\kappa \mathrm{median}\left(\left|x_{i-k}-m_i\right|,\dots ,\left|x_n-m_n\right|\right)$$

참고 문헌

[1] Liu, Hancong, Sirish Shah, and Wei Jiang. “On-line outlier detection and data cleaning.” Computers and Chemical Engineering. Vol. 28, March 2004, pp. 1635–1647.

[2] Suomela, Jukka. “Median Filtering Is Equivalent to Sorting.” 2014.

확장 기능

C/C++ 코드 생성

MATLAB® Coder™를 사용하여 C 코드나 C++ 코드를 생성할 수 있습니다.

사용법 관련 참고 및 제한 사항:

• 가변 크기 입력 행렬 또는 N차원 배열은 런타임에 행 벡터가 되지 않아야 합니다.

• 플로팅은 MATLAB^® 실행 시에만 지원됩니다.

버전 내역

R2015b에 개발됨

참고 항목

medfilt1 | median | filloutliers | filter | isoutlier | mad (Statistics and Machine Learning Toolbox) | movmad | movmedian | sgolayfilt

도움말 항목

• Eliminate Outliers Using Hampel Identifier