【NumPy】高速フーリエ変換でパワースペクトル解析

スポンサーリンク
ビッグバナー(上2)

この記事では、Python言語とNumPyを用いて、高速フーリエ変換(FFT)でパワースペクトルを計算する方法をソースコード付きで解説します。

パワースペクトル解析

パワースペクトルとは、信号の振幅と周波数の関係を示す指標です。
フーリエ変換(Fourier Transform)によりパワースペクトルを求めることができます。

今回は、PythonモジュールNumPyのnumpy.fft.fftを用いて高速フーリエ変換を行い、周波数スケールで振幅と位相をグラフ表示してみました。

【書式】
numpy.fft.fft(ndarray)

■引数
ndarray:1次元のNumpy配列

■返り値
フーリエ変換の結果(複素数)

■使用したデータ:nikkei.csv
(2016年の日経平均株価のデータ)

ソースコード

サンプルプログラムのソースコードは下記の通りです。

【補足】
離散フーリエ変換の定義から、FFT後のデータに1/Nを掛けて正規化しています。
また、FFT後のデータは、対称な共役複素数です。
ナイキスト周波数を中心とした対称な周波数に対応おり、対称同士の周波数の波は区別できません。
これをエイリアシング現象といいます。。
そのため、ナイキスト周波数より高い周波数帯のデータは除外しています。
そして、直流成分を1/2倍(交流成分を2倍)して揃えています。

実行結果

サンプルプログラム実行結果は下記の通りです。
2016年の1年分の日経平均株価の終値をサンプルとしています。

左上のグラフは、サンプルデータ(日経平均株価の終値)です。(単位は千円)
右上のグラフは、FFT処理されたサンプルデータのうち、振幅を周波数スケールで表示したものです。
左下のグラフは、FFT処理されたサンプルデータのうち、位相を周波数スケールで表示したものです。

右上のグラフを見てやると、どの周波数帯の信号が強いかわかります。
※元のサンプルデータには、周波数が0付近の信号が多く含まれていることがわかります。

【関連記事】
Python入門 基本文法
NumPy入門 サンプルプログラム集

スポンサーリンク
レクタングル(下2)
レクタングル(下2)

シェア&フォローお願いします!