フーリエ変換

高速フーリエ変換

高速フーリエ変換(Fast Fourier Transformation)を利用することができます。 例えば、データとしてmydataがあるとき、

fftresult = np.fft.fft(mydata)

で、mydataの離散フーリエ変換結果が手に入ります。 fftの計算は

なので（参照：https://docs.scipy.org/doc/numpy/reference/routines.fft.html#module-numpy.fft）、結果は基本的に複素数のnumpy arrayです。プロットする際には、np.real(), np.imag(), np.abs(), np.angle()などを使って、実部, 虚部, 絶対値, 位相などを取り出すと良いでしょう。

結果を解釈するには（例えばスペクトルをプロットしたいとき）各データがどのような周波数に対応するのかという情報が必要です。これは、mydataがそもそも何秒間隔のデータなのかということに依存します。もし0.01秒間隔であれば

freqarray = np.fft.fftfreq(xdata.size, d=0.01)

とすることで、freqarrayという名前の、周波数軸に対応するarrayを作ることができます。

Power Spectral Density

ここではPower Spectral Densityの詳しい説明は省きますが、大雑把には、「フーリエ変換の絶対値の自乗を適切に規格化したもの」といえます。「適切に」の意味合いは以下のようになります。対象とするデータを $g(t)$ とするとき $(g(t))^2$ のことを、慣習的に「パワー」と呼びます（これは電気回路の解析に影響された用語です。実際に電気回路で電圧 $V(t)$ を測定している場合は、その回路で消費されるパワーは $(V(t))^2$ で書くことができます）。この瞬間的なパワーの値の統計的な平均値を、周波数ごとに分割したものがPSDです。前述のようにフーリエ変換の絶対値の自乗に比例するような量で

$$\lim_{T \to \infty} \dfrac{1}{T}\int_{-\infty}^{\infty}g^2(t) dt = \int_{-\infty}^{\infty}{S_g}(f) df= \dfrac{1}{2\pi}\int_{-\infty}^{\infty}{S_g}(\omega) d\omega$$

という性質を満たします。この $S_g(f)$ をtwo-sided power spectral densityと呼びます。 測定値 $g(t)$ とす$S_g(-f)=S_g(f)$ という性質を満たすので、新しく $S'_g(f)=2S_g(f)$ というone-sided power spectral densityを定義すると

$$\lim_{T \to \infty} \dfrac{1}{T}\int_{-\infty}^{\infty}g^2(t) dt = \int_{0}^{\infty}{S'_g}(f) df= \dfrac{1}{2\pi}\int_{0}^{\infty}{S'_g}(\omega) d\omega$$

のように正の周波数領域だけを気にすれば良くなります。

pythonではscipyモジュールを使うことでpower spectral densityを計算できます。mydataというnumpy arrayのone-sided power spectral densityを計算するには

import scipy.signal as signal

fps = 100
segmentlength=1000

mypsd=signal.welch(mydata,fps,'boxcar',segmentlength)

とします。mypsd[0]に周波数のarrayが、mypsd[1]にPSDのarrayが格納されます。 実際にやっていることは以下の図のようになっています。

segmentlength/fpsがフーリエ変換の時間範囲に対応し、周波数分解能を決めます。したがって、segmentlengthは十分な周波数分解能が得られるような値を選ぶ必要があります。窓関数やオーバーラップの具合などもオプションで設定可能です。詳細は https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.welch.html を参照してください