대략적으로 다음과 같이 주어진 데이터 세트가 있다고 가정 해 봅시다.
import numpy as np
x = np.linspace(0,2*np.pi,100)
y = np.sin(x) + np.random.random(100) * 0.2
따라서 데이터 세트의 20 % 변형이 있습니다. 첫 번째 아이디어는 scipy의 UnivariateSpline 함수를 사용하는 것이었지만 문제는 이것이 작은 소음을 좋은 방식으로 고려하지 않는다는 것입니다. 주파수를 고려하면 배경이 신호보다 훨씬 작으므로 컷오프 만 스플라인으로 생각할 수 있지만 앞뒤 푸리에 변환이 수반되어 나쁜 동작이 발생할 수 있습니다. 또 다른 방법은 이동 평균이 될 수 있지만 지연의 올바른 선택이 필요합니다.
이 문제를 해결하는 방법에 대한 힌트 / 책 또는 링크가 있습니까?
답변
나는 Savitzky-Golay 필터를 선호합니다 . 최소 제곱을 사용하여 데이터의 작은 창을 다항식으로 되 돌린 다음 다항식을 사용하여 창 중앙의 점을 추정합니다. 마지막으로 창은 하나의 데이터 포인트만큼 앞으로 이동하고 프로세스가 반복됩니다. 이것은 모든 점이 이웃에 대해 최적으로 조정될 때까지 계속됩니다. 비 주기적 및 비선형 소스의 노이즈 샘플에서도 효과적입니다.
다음은 철저한 요리 책 예제 입니다. 사용하기 쉬운 방법에 대한 아이디어를 얻으려면 아래 코드를 참조하십시오. 참고 : savitzky_golay()
위에 링크 된 요리 책 예제에서 문자를 복사하여 붙여 넣을 수 있기 때문에 함수 를 정의하기위한 코드를 생략했습니다 .
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0,2*np.pi,100)
y = np.sin(x) + np.random.random(100) * 0.2
yhat = savitzky_golay(y, 51, 3) # window size 51, polynomial order 3
plt.plot(x,y)
plt.plot(x,yhat, color='red')
plt.show()
업데이트 : 내가 연결된 요리 책 예제가 중단되었습니다. 다행히도 Savitzky-Golay 필터는 @dodohjk가 지적한 것처럼 SciPy 라이브러리 에 통합 되었습니다 . SciPy 소스를 사용하여 위의 코드를 수정하려면 다음을 입력하십시오.
from scipy.signal import savgol_filter
yhat = savgol_filter(y, 51, 3) # window size 51, polynomial order 3
답변
이동 평균 상자를 기반으로 컨볼 루션을 사용하여 데이터를 부드럽게 처리하는 빠르고 더러운 방법 :
x = np.linspace(0,2*np.pi,100)
y = np.sin(x) + np.random.random(100) * 0.8
def smooth(y, box_pts):
box = np.ones(box_pts)/box_pts
y_smooth = np.convolve(y, box, mode='same')
return y_smooth
plot(x, y,'o')
plot(x, smooth(y,3), 'r-', lw=2)
plot(x, smooth(y,19), 'g-', lw=2)
답변
예를 들어주기적인 신호의 “부드러운”버전에 관심이 있다면 FFT가 올바른 방법입니다. 푸리에 변환을 수행하고 기여도가 낮은 주파수를 뺍니다.
import numpy as np
import scipy.fftpack
N = 100
x = np.linspace(0,2*np.pi,N)
y = np.sin(x) + np.random.random(N) * 0.2
w = scipy.fftpack.rfft(y)
f = scipy.fftpack.rfftfreq(N, x[1]-x[0])
spectrum = w**2
cutoff_idx = spectrum < (spectrum.max()/5)
w2 = w.copy()
w2[cutoff_idx] = 0
y2 = scipy.fftpack.irfft(w2)
신호가 완전히 주기적이 아닌 경우에도 백색 잡음을 제거하는 데 큰 도움이됩니다. 사용할 여러 유형의 필터 (고역 통과, 저역 통과 등)가 있으며 적절한 필터는 찾고자하는 것에 따라 다릅니다.
답변
이동 평균을 데이터에 맞추면 노이즈가 완화됩니다 .이를 수행하는 방법은 이 답변 을 참조하십시오 .
LOWESS 를 사용 하여 데이터에 맞추 려면 (이동 평균과 비슷하지만 더 정교한) statsmodels 라이브러리 를 사용하여 수행 할 수 있습니다 .
import numpy as np
import pylab as plt
import statsmodels.api as sm
x = np.linspace(0,2*np.pi,100)
y = np.sin(x) + np.random.random(100) * 0.2
lowess = sm.nonparametric.lowess(y, x, frac=0.1)
plt.plot(x, y, '+')
plt.plot(lowess[:, 0], lowess[:, 1])
plt.show()
마지막으로, 신호의 기능적 형태를 알고 있다면 데이터에 곡선을 맞출 수 있습니다. 아마도 최선의 방법 일 것입니다.
답변
또 다른 옵션은 statsmodels 에서 KernelReg 를 사용하는 것 입니다 .
from statsmodels.nonparametric.kernel_regression import KernelReg
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0,2*np.pi,100)
y = np.sin(x) + np.random.random(100) * 0.2
# The third parameter specifies the type of the variable x;
# 'c' stands for continuous
kr = KernelReg(y,x,'c')
plt.plot(x, y, '+')
y_pred, y_std = kr.fit(x)
plt.plot(x, y_pred)
plt.show()
답변
이것 좀 봐! 1D 신호의 평활화에 대한 명확한 정의가 있습니다.
http://scipy-cookbook.readthedocs.io/items/SignalSmooth.html
지름길:
import numpy
def smooth(x,window_len=11,window='hanning'):
"""smooth the data using a window with requested size.
This method is based on the convolution of a scaled window with the signal.
The signal is prepared by introducing reflected copies of the signal
(with the window size) in both ends so that transient parts are minimized
in the begining and end part of the output signal.
input:
x: the input signal
window_len: the dimension of the smoothing window; should be an odd integer
window: the type of window from 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'
flat window will produce a moving average smoothing.
output:
the smoothed signal
example:
t=linspace(-2,2,0.1)
x=sin(t)+randn(len(t))*0.1
y=smooth(x)
see also:
numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman, numpy.convolve
scipy.signal.lfilter
TODO: the window parameter could be the window itself if an array instead of a string
NOTE: length(output) != length(input), to correct this: return y[(window_len/2-1):-(window_len/2)] instead of just y.
"""
if x.ndim != 1:
raise ValueError, "smooth only accepts 1 dimension arrays."
if x.size < window_len:
raise ValueError, "Input vector needs to be bigger than window size."
if window_len<3:
return x
if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
raise ValueError, "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'"
s=numpy.r_[x[window_len-1:0:-1],x,x[-2:-window_len-1:-1]]
#print(len(s))
if window == 'flat': #moving average
w=numpy.ones(window_len,'d')
else:
w=eval('numpy.'+window+'(window_len)')
y=numpy.convolve(w/w.sum(),s,mode='valid')
return y
from numpy import *
from pylab import *
def smooth_demo():
t=linspace(-4,4,100)
x=sin(t)
xn=x+randn(len(t))*0.1
y=smooth(x)
ws=31
subplot(211)
plot(ones(ws))
windows=['flat', 'hanning', 'hamming', 'bartlett', 'blackman']
hold(True)
for w in windows[1:]:
eval('plot('+w+'(ws) )')
axis([0,30,0,1.1])
legend(windows)
title("The smoothing windows")
subplot(212)
plot(x)
plot(xn)
for w in windows:
plot(smooth(xn,10,w))
l=['original signal', 'signal with noise']
l.extend(windows)
legend(l)
title("Smoothing a noisy signal")
show()
if __name__=='__main__':
smooth_demo()