Helve’s Python memo

Pythonを使った機械学習やデータ分析の備忘録

<NumPy> 高速フーリエ変換による周波数解析

NumPyのfftパッケージを使って、FFT (Fast Fourier Transform, 高速フーリエ変換) による離散信号の周波数解析を行い、信号の振幅を求める。

目次

環境

python 3.6.5
NumPy 1.14.2

元の離散信号

次の3つの信号を合成して、フーリエ変換を行う離散信号とする。

  • 周波数50[Hz], 振幅1.5の正弦波
  • 周波数120[Hz], 振幅1の正弦波
  • 定数項3

信号のデータ点数は1024, サンプリング周期は0.001[s]とする。

import numpy as np
import matplotlib.pyplot as plt

N = 1024            # サンプル数
dt = 0.001          # サンプリング周期 [s]
f1, f2 = 50, 120    # 周波数 [Hz]

t = np.arange(0, N*dt, dt) # 時間 [s]
x = 1.5*np.sin(2*np.pi*f1*t) + np.sin(2*np.pi*f2*t) + 3 # 信号

fig, ax = plt.subplots()
ax.plot(t, x)
# ax.set_xlim(0, 0.1)
ax.set_xlabel("Time [s]")
ax.set_ylabel("Signal")
ax.grid()
plt.show()

信号の波形を次のグラフに示す。
f:id:Helve:20180616224001p:plain

信号を0~0.1[s]の範囲に拡大した図は次の通り。
f:id:Helve:20180616224010p:plain

フーリエ変換

前章の信号に対して、numpy.fft.fft()により、フーリエ変換を行う。

numpy.fft.fft(a, n=None, axis=-1, norm=None)

引数の説明は以下の通り。
a: FFTを行う配列。
n: FFTを行うデータ点数。Noneとすると、aの長さに等しくなる。
axis: FFTを行う配列の軸方向。指定しなければ、配列の最大の次元の方向となる。
norm: "ortho"とすると、正規化を行う。正規化すると、変換値が1/√Nになる(Nはデータ点数)。

戻り値は、長さnの複素数配列である。


また、フーリエ変換の周波数は、numpy.fft.fftfreq()により取得する。

numpy.fft.fftfreq(n, d=1.0)

引数の説明は以下の通り。
n: FFTを行うデータ点数。
d: サンプリング周期(デフォルト値は1.0)。

戻り値は、周波数の配列となる。

先程の信号xに対してFFTを行い、変換結果の実部、虚部、周波数をプロットする。

F = np.fft.fft(x) # 変換結果
freq = np.fft.fftfreq(N, d=dt) # 周波数

fig, ax = plt.subplots(nrows=3, sharex=True, figsize=(6,6))
ax[0].plot(F.real, label="Real part")
ax[0].legend()
ax[1].plot(F.imag, label="Imaginary part")
ax[1].legend()
ax[2].plot(freq, label="Frequency")
ax[2].legend()
ax[2].set_xlabel("Number of data")
plt.show()

f:id:Helve:20180616224053p:plain

周波数freqにおいて、最初の要素freq[0]は0[Hz],
2~N/2 (=512) 番目の要素freq[1:512]は正の周波数、
(N/2+1)番目以降の要素freq[512:]は負の周波数である。

元の信号の周波数が1000[Hz]であるから、サンプリング定理(または標本化定理)より、その1/2以下の周波数領域 (500[Hz]以下) でフーリエ変換の結果は有効である。

周波数が0, 50, 120, -120, -50[Hz]のときに、実部や虚部がピーク値を持つ。


最後に、元の信号の振幅を求める。
FをN/2で割り、絶対値をとると振幅になる。
正の周波数領域 について、振幅をプロットする。

Amp = np.abs(F/(N/2)) # 振幅

fig, ax = plt.subplots()
ax.plot(freq[1:int(N/2)], Amp[1:int(N/2)])
ax.set_xlabel("Freqency [Hz]")
ax.set_ylabel("Amplitude")
ax.grid()
plt.show()

結果は以下の通り。

f:id:Helve:20180616223946p:plain

振幅は、周波数50[Hz]で約1.42, 120[Hz]で約0.98となっており、
元の値(50[Hz]で約1.5, 120[Hz]で1.0)に近い値となっている。