개요

파이썬을 사용하여 시계열 데이터의 FFT를 수행하고 그 피크 검출을 수행하는 방법을 요약한다.

데이터 준비

해석 예로 하는 시계열 데이터를 작성한다. 3개의 정현파와 노이즈를 조합한 데이터를 다음과 같이 작성했다.

import numpy as np
import matplotlib.pyplot as plt

# parameters
dt=6 #data duration
fsamp=10000 #sampling rate
f=[28, 60, 213, 355] # sin frequency
y=[0.6, 1.3, 0.3, 0.5] # sin amplitude
pi=np.pi

# data creation
# gauusian noise sigma=1
n=dt*fsamp+1 #data size![time_series.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2184372/45ed0137-bbff-fa65-3130-89ffecf07771.png)

t=np.linspace(0,dt,n) #sampling time
y=np.random.randn(n)\
+y[0]*np.sin(2*np.pi*f[0]*t)\
+y[1]*np.sin(2*np.pi*f[1]*t)\
+y[2]*np.sin(2*np.pi*f[2]*t)\
+y[3]*np.sin(2*np.pi*f[3]*t)

# plot
plt.plot(t,y)
plt.xlabel("Time [s]")
plt.ylabel("Amplitude [a.u.]")
plt.show()

결과는 아래 그림과 같습니다.

FFT

생성된 시계열 데이터를 FFT한다. FFT의 효율을 위해서는 데이터수가 2의 제곱이 되도록 하는 것이 좋지만, 이번에는 데이터수가 적기 때문에 신경쓰지 않기로 한다.
FFT 계수 y_fft는 복소수이므로 플롯하기 쉽도록 진폭 스펙트럼 amp_fft으로 변환하여 플롯합니다. 진폭 스펙트럼이 원래 신호의 진폭을 재현하도록 AC 성분에는 2/(데이터 길이)를, DC 성분에는 1/(데이터 길이)를 곱하여 규격화해 두면 플롯할 때 해석 쉬워진다.

# t: sampling time
# y: data
# fsamp: sampling rate

#FFT
y_fft=np.fft.fft(y)
y_fft=y_fft[:int(len(y_fft)/2)+1]
f_fft=fsamp/2*np.linspace(0,1,len(y)/2+1)
amp_fft=np.abs(y_fft)/len(y)*2
amp_fft[0]=amp_fft[0]/2

# plot
plt.plot(f_fft, amp_fft)
plt.xlim(1,500)
plt.xlabel("Frequency [Hz]")
plt.ylabel("Amplitude [a.u.]")
plt.grid()
plt.show()

결과는 아래 그림과 같이 되어, 설정한 4개의 정현파가 검출되고 있는 것을 알 수 있다.

피크 검출

전항에서는 진폭 스펙트럼을 플롯하여 피크를 확인했지만, 그 피크를 자동으로 검출할 수 있도록 한다. 피크 검출이 가능한 라이브러리에는 scipy.signal.argrelmax이나 scipy.signal.find_peaks가 있다. 이번에는 피크 검출 조건을 세밀하게 설정할 수 있는 후자의 find_peaks 를 사용한다.
find_peaks는 데이터 x와 피크 검출 조건 (예 : height)을 입력하고 scipy.signal.find_peaks(x, height, threshold, distance, prominence, width, wlen, rel_height, plateau_size)와 같이 사용합니다. 다음 코드에서는 피크 값이 주변에서 얼마나 돌출하는지를 나타내는 prominence를 설정하여 잡음을 피크 인정하지 않도록 조정하고 있다. 결과를 전항의 그래프에 겹치면, 확실히 피크가 검출되고 있는 것을 알 수 있다.

# peak detection
peaks,_ = signal.find_peaks(amp_fft,prominence=0.1)

#amplitude spectrum plot
plt.plot(f_fft, amp_fft)
plt.scatter(f_fft[peaks], amp_fft[peaks], color='red')
plt.xlim(1,500)
plt.xlabel("Frequency [Hz]")
plt.ylabel("Amplitude [a.u.]")
plt.grid()
plt.show()

이 때 prominence의 값이 부적절하면, 예를 들면 prominence = 1, 0.05의 경우, 피크 검출이 불충분해지거나, 잡음까지 주워 과잉으로 피크 검출해 버린다. prominence에 한정되지 않고, 이러한 하이퍼 파라미터의 설정에는 주의가 필요하다.