소개

파이썬에서 Numpy를 사용하여 시계열 데이터를 FFT (Fast Fourier Transform : 고속 푸리에 변환)하는 방법과 시계열 데이터의 추세를 제거하는 방법에 대해 소개하려고합니다. FFT는 DFT (Discrete Fourier Transform : 이산 푸리에 변환)를 고속 처리하는 계산 방법입니다. 이 기사에서는 이론을 다루지 않고 FFT를 실행하는 최소 코드를 보여줍니다. 참고로 한 문헌은 「스펙트럼 해석 저:히노 미츠오(아사쿠라 서점)」. 푸리에 해석의 기초부터 FFT의 이론까지, 이 책 1권으로 충분합니다.

목차

1. 시계열 데이터
2.FFT 실행
3. 트렌드 제거
4. 푸리에 성분을 주파수 평활화(스무딩)
5.Appendix

주제

1. 시계열 데이터

샘플링 주파수 10Hz의 30분간 데이터. 오렌지 선은 이동 평균입니다. 트렌드가 있는 것을 알 수 있군요.

그림 1. 원래 시계열 데이터

이 데이터를 FFT합니다.

2. FFT 실행

N =len(X)      #データ長
fs=10          #サンプリング周波数
dt =1/fs       # サンプリング間隔
t = np.arange(0.0, N*dt, dt) #時間軸
freq = np.linspace(0, fs,N) #周波数軸
fn=1/dt/2     #ナイキスト周波数

FFT는 데이터 길이가 2의 제곱인 경우에 가장 계산 속도가 빨라지는 계산 방법입니다만, 그렇지 않은 경우에도 계산은 가능합니다(처리 시간이 다소 길어집니다만). 다만, 데이터 길이가 소수인 경우에는, 2의 멱승시에 비해 상대적으로 상당히 처리 시간이 걸리므로, 2의 멱승이 되도록 0패딩 하는 것이 좋을 것 같습니다. Appendix에 그것을 확인하는 코드를 올렸으므로 꼭 확인해보십시오.

F=np.fft.fft(X)/(N/2)
F[(freq>fn)]=0 #ナイキスト周波数以降をカット

plt.plot(freq,np.abs(F))#
plt.xlabel("[Hz]")
plt.ylabel("Amp")
plt.xlim(-0.01,0.5)
plt.grid()
plt.show()

np.fft.fft에서는 복소 푸리에 성분으로 주어지므로, 그 절대치를 취한 것을 플롯한 것이 아래 그림입니다. 트렌드 성분이 0Hz 부근에 있네요.
다음 섹션에서는 추세를 제거해 보겠습니다.

그림 2. 원래 시계열 데이터에 FFT 실행

3. 트렌드 제거

그림 1의 이동 평균(오렌지선)을 보면 시계열 데이터의 축이 x축과 괴리하고 있어 데이터에 트렌드가 있는 것을 알 수 있습니다. 다음과 같이 0.03Hz 이하를 0으로 설정하여 추세를 제거합시다.

F=np.fft.fft(X)/(N/2)
F[(freq>fn)]=0
F[(freq<=0.03)]=0 #0.03HZ以下を除去
X_1=np.real(np.fft.ifft(F))*N

plt.xlabel("Time [s]")
plt.ylabel("Signal")
plt.xlim(-50,1850)
plt.grid()
plt.show()

그림 3. 추세 제거 후 시계열 데이터

x축을 기점으로 주기 함수가 되어 있는 것을 알 수 있습니다.

4. 푸리에 성분을 주파수 평활화(스무딩)

그림 3의 데이터를 FFT하면 그림 4가됩니다.

그림 4. 추세 제거 후 시계열 데이터에 대한 FFT

덜컹 거리고 있기 때문에, 평활화 윈도우를 작용시켜 스무딩을 해 봅니다.

window=np.ones(5)/5 #平滑化ウィンドー
F3=np.convolve(F2,window,mode='same') #畳み込み
F3=np.convolve(F3,window,mode='same') #畳み込み
F3=np.convolve(F3,window,mode='same') #畳み込み

plt.plot(freq,np.abs(F3))
plt.xlabel("[Hz]")
plt.ylabel("Amp")
plt.xlim(-0.01,0.5)
plt.grid()
plt.show()

그림 5. 매끄러운

5. Appendix

데이터 길이를 3 종류 준비하고 계산 시간을 비교해 보겠습니다.
①2^19(2의 멱승) ②2^(19)-1 (소수) ③2^(19)-2

import time

if __name__ == '__main__':
    start = time.time()
    x2 = np.random.uniform(size=2**19-2)#2**19 , 2**19-1
    print(np.fft.fft(x2))
    elapsed_time = time.time() - start
    print ("elapsed_time:{0}".format(elapsed_time) + "[sec]")

계산 결과
①0.04197[sec]
②0.1679[sec]
③0.05796[sec]

데이터 길이가 소수인 경우에는 0패딩을 실시해 2의 거듭제곱으로 하는 것이 좋을 것 같네요.