[SciPy] 사비츠키-골레이 필터 (Savitzky-Golay Filter)

사비츠키-골레이 필터 (Savitzky-Golay Filter)

 

사비츠키-골레이 (Savitzky-Golay) 평활화 필터는 잡음이 일부 섞여 있지만 잡음 없는 영역이 주파수 범위의 큰 부분을 차지하는 신호를 "평활화"하는 데 사용된다. 이 필터는 디지털 평활화 다항식 필터 또는 최소제곱 평활화 필터라고도 한다. 사비츠키-골레이 평활화 필터는 표준 평균 FIR 필터보다 신호의 고주파 성분을 적게 제거하는 경향이 있다. 일부 응용 사례에서는 표준 평균 FIR 필터보다 사비츠키-골레이 필터가 성능이 더 좋다. 표준 평균 FIR 필터는 잡음과 함께 고주파 성분도 필터링하는 경향이 있다. 그러나 잡음 수준이 특히 높은 경우 잡음 제거 측면에서는 표준 평균 FIR 필터보다 효율이 떨어진다. 사비츠키-골레이 필터는 고주파 신호 성분을 유지하는 데는 더 효과적이지만 잡음 제거 측면에서는 덜 효과적이다.

 

필터링은 인근 점에서 측정한 값들의 기본적인 수준이 거의 동일하다는 가정하에, 신호의 각 점을 그 점을 중심으로 하는 이동 윈도우에 포함된 신호 값들의 조합으로 대체하여 이루어진다. 예를 들어, 이동 평균 필터는 각 데이터 점을 주변 데이터 점의 국소 평균으로 대체한다. 주어진 데이터 점의 왼쪽에 점이 k개 있고 오른쪽에 점이 k개 있다면, 이동평균 필터는 총 L = 2k + 1의 윈도우 길이에 대해 다음과 같이 대체한다.

 

 

사비츠키-골레이 필터는 윈도우 내의 신호 값 전체에 대해 n차 다항식을 최소제곱 피팅하고 피팅된 다항식 곡선의 계산된 중심 점을 새로운 평활화된 데이터 점으로 택하여 이 개념을 일반화한다. 이는 주어진 점 xs에 대해 다음과 같다.

 

 

행렬로 표현하면 다음과 같다.

 

 

사비츠키-골레이 추정값을 구하려면 H의 의사 역행렬을 사용하여 a를 계산하고 결과 앞에 H를 곱한다.

 

 

sgolay는 조건이 나빠지지 않도록 qr 함수를 사용하여 B = QQT가 되는 H의 효율적 크기의 분해 H = QR을 계산한다. B는 한 번만 계산하면 된다. 대부분의 신호 점에 대해 사비츠키-골레이 추정값은 신호를 B의 중앙 행과 컨벌루션하여 얻게 된다. 필터링된 신호의 정상 상태 부분이 결과로 반환된다. B의 처음 k개 행은 초기 과도 상태의 값을 생성하고 B의 마지막 k개 행은 최종 과도 상태의 값을 생성한다. 윈도우 길이를 늘려서 잡음 억제를 개선할 수 있지만, 이렇게 하면 모든 과도 상태 근처에서 깁스 현상과 비슷한 링잉 현상이 발생한다.

 

scipy.signal.savgol_filter

 

1-D의 필터로 x의 치수가 1보다 크면, 축이 필터를 적용할 축을 결정한다. mode는 'mirror', 'constant', 'nearest', 'wrap' 또는 'interp'이어야 한다. 이는 필터가 적용되는 패딩 신호에 사용할 확장의 유형을 결정한다. mode가 'constant'이면, 패딩 값은 cval로 주어진다. 'interp' 모드를 선택한 경우(기본값), 확장이 사용되지 않는다. 대신 에지의 마지막 window_length 값에 차수 다항식이 적합하며, 이 다항식은 마지막 window_length / 2 출력 값을 평가하는 데 사용된다.

 

예를 들어 입력이 [1, 2, 3, 4, 5, 6, 7, 8]이고 window_length가 7인 경우 다양한 mode 옵션에 대한 확장 데이터를 보여준다 (cval을 0으로 가정).

 

mode       |   Ext   |         Input          |   Ext
-----------+---------+------------------------+---------
'mirror'   | 4  3  2 | 1  2  3  4  5  6  7  8 | 7  6  5
'nearest'  | 1  1  1 | 1  2  3  4  5  6  7  8 | 8  8  8
'constant' | 0  0  0 | 1  2  3  4  5  6  7  8 | 0  0  0
'wrap'     | 6  7  8 | 1  2  3  4  5  6  7  8 | 1  2  3

import numpy as np
from scipy.signal import savgol_filter

np.set_printoptions(precision=2)  # For compact display

x = np.array([2, 2, 5, 2, 1, 0, 1, 4, 9])

savgol_filter(x, 5, 2)

array([1.66, 3.17, 3.54, 2.86, 0.66, 0.17, 1.  , 4.  , 9.  ])

savgol_filter(x, 5, 2, mode='nearest')

array([1.74, 3.03, 3.54, 2.86, 0.66, 0.17, 1.  , 4.6 , 7.97])