PythonでFFTとそのピーク検出をする方法

#概要
Pythonを用いて時系列データの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に限らず，このようなハイパーパラメータの設定には注意が必要である。