Python FFTを実施する例を作成してみた。
そのコードは以下のようになった。
フーリエ変換については，係数の部分が定義によって異なるので，この例のnp.fft.fft()の使用方法ではMATLABやGNU Octaveのように配列要素数分係数がかかるので物理的な意味を持つ振幅に直す場合にはlen(f)で割る必要がある。さらに，下の例では，折り返す分を半分しか表示していないので振幅としては半分になってしまう。そこで2をかけて調整している。
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
N = 1000 # Number of Samples
dt = 0.01 # Sampling Interval
f1, f2 = 23, 36 # Frequency
t = np.arange(0, N * dt, dt) # Time range
freq = np.linspace(0, 1.0 / dt, N) # Frequecny range
# Signal (Sinusoidal waves, frequencies of f1 and f2 + noise)
f = np.sin(2 * np.pi * f1 * t) \
+ np.sin(2 * np.pi * f2 * t) \
+ 0.8 * np.random.randn(N)
F = np.fft.fft(f) / len(f) * 2
F[0] = F[0] / 2.0
# 振幅スペクトルを計算
Amp = np.abs(F)
# plot
plt.figure()
plt.subplots_adjust(wspace=0.4, \
hspace=0.0)
plt.subplot(121)
plt.plot(t, f, label='Raw signal')
plt.xlabel("time")
plt.ylabel("signal")
plt.grid()
plt.ylim([-7, 7])
leg = plt.legend(loc=1)
plt.subplot(122)
plt.plot(freq[0:int(len(F)/2)], \
Amp[0:int(len(F)/2)], \
label='Amplitude')
plt.xlabel('frequency')
plt.ylabel('amplitude')
plt.grid()
plt.ylim([0, 2])
leg = plt.legend(loc=1)
plt.show()
