同期の研究に首を突っ込んでいるときに、フーリエ変換の話題が出てきました。
そのときにフーリエ変換された配列の要素番号の意味がよく分からず、ググってもいまいち出てこなかったので(常識だから?)、分かったことをちょっとまとめときます。
とりあえずフーリエ変換とはなんぞや?ということで、知らない方のために個人的な見解を述べておきます(間違ってたらご指摘お願いします・・・)。
有限範囲の外でf(x)=0である関数(もしくは周期的な関数)があるとき、その関数はsinとcosの足し合わせで表すことができます。
そのsinとcosの足し合わせがどのようになっているのかを求めること、がざっくりとしたフーリエ変換の説明です。
物理学科では学部2年とかで習う基本的な物理数学ですね、その頃の私は全く理解できていませんでしたが(汗
んで本題ですが、結論から述べると、ある配列をフーリエ変換したときにその
結果のx軸は周波数に変換せねばならず、その間隔は
1 / (データ数 * サンプリングの時間間隔)
になる、ということです。
具体例として、例えば次の式を考えてみましょう。
=== 追記20160821 ===
この式はおかしいですね・・・
正しくは
f(t) = \sigma^4_0 k sin[(2k+1) \pi t / 2]
ですね・・・
=== 追記終わり ===
書き出すと
一応こんな感じの関数です。
んでこれをフーリエ変換した結果のパワースペクトルの周波数には明らかに0.25Hz、0.75Hz、1.25Hz、1.75Hz、2.25Hzが表れるはずです。
ところでPythonのscipyモジュールには"fft()"という関数があり、これを使えば与えた配列をフーリエ変換した結果を返してくれます。
そんなわけで上のf)t)をフーリエ変換した結果をそのまま描くと次のようになります。
ここでは横軸を何も指定していないので、配列の要素番号がそのまま横軸になってしまっています。
で横軸どないすればええねん?というわけですが、実はscipyにはこの周波数のスケールを作る"fftfreq()"という関数があり、それをx軸にして再度グラフを描くと次のようになります。
横軸は周波数[Hz]です。
見たいスケールが小さすぎて見えませんね・・・
0~4Hz付近を拡大しましょう。
サンプリングの点数が少ないために多少荒い図になっていますが、最初に言った周波数にパワーが表れていますね。
ただ本当は、この"fftfreq()"の中身が知りたくてこの後も色々調べたわけですが、長くなってきたので次の記事でまた見ていきましょう。。。
ちなみに、最後のフーリエ変換の結果は次のようなスクリプトになります。参考までに載せときます。
===== ex1_fftsin.py =====
import numpy as np
from scipy.fftpack import fft, fftfreq
import matplotlib.pyplot as plt
N = 1024 # データの点数
dt = 0.01 # サンプリングの時間間隔
t = np.linspace(0, (N-1)*dt, N)
y = np.sin(np.pi*t/2) + 3*np.sin(3*np.pi*t/2) + 5*np.sin(5*np.pi*t/2) + 7*np.sin(7*np.pi*t/2) 9*np.sin(9*np.pi*t/2)
y_fft = fft(y)
freq = fftfreq(N, dt)
plt.xlim(0, 4) # 周波数の見たい領域を指定
plt.plot(freq[1:N/2], abs(y_fft)[1:N/2]) # ナイキスト周波数のため、配列の半分のみ使用
plt.show()
====================