SciPyでFFT

一番よく使うのが、FFT。さらっと、コードを載せておく。
プロットは、模索中。

# -*- coding: utf-8

# FFT（SciPy）のテスト

import numpy as np
import scipy.fftpack as sf
import matplotlib.pyplot as plt

# 適当な信号（ランダム） 最大振幅は1.0
mu, sigma = 0, 0.1
dt = 0.01
val = np.random.normal(mu, sigma, 3000)
val /= max(abs(val))
time = np.arange(0, 30.0, dt)
# SciPyを用いてFFT
fourier_val = sf.fft(val) * dt
freq = sf.fftfreq(len(time), dt)

# 簡易プロット
fig = plt.figure(figsize=(9,6))
# 元波形
ax1 = plt.subplot(211)
plt.xlim(0.0, 30.0)
plt.ylim(-1.0, 1.0)
plt.xlabel("time(sec)")
plt.ylabel("random")
plt.plot(time, val, "-")
# フーリエスペクトル
ax2 = plt.subplot(212)
ax2.set_xscale("log")
ax2.set_yscale("log")
plt.xlim(1.0, 10.0)
plt.ylim(0.01, 1.0)
plt.xlabel("Frequency(Hz)")
plt.ylabel("amp")
n = len(freq)
plt.plot(freq[0:n/2], abs(fourier_val[0:n/2]), "-")

# 出力の方法
if 1:
  plt.savefig("fft_plot1.png")
  plt.savefig("fft_plot1.emf")

if 1: 
  plt.show()

出力される結果はこんなん。

単位は、例えば加速度（m/s2とかgal）の信号に対して、FFT後の単位が(m/sとかgal*s)とかになる感じ。
フーリエ変換後にどんな単位にしとくのがいいかは、よく分からない。

最後の出力の部分は、よくあるif文の使い方を真似てみた。
matlabとかでよく見る使い方。
昔は嫌いだったけど、今後使えるかもと、再評価してるとこ。

ランダム波形としては、本当はホワイトノイズを使いたかった。本来、こういうやり方をするべきなんだろう。 波形もスペクトルもきれいですもんね。 今回は簡単に、numpy.random.normalを参考にしました。