【Python】気象庁の強震波形をフーリエ変換でスペクトル解析

Pythonを用いて、気象庁の強震波形データをフーリエ変換し、スペクトル解析する方法について紹介します。

【Python】強震観測データのスペクトル解析(矩形窓)

気象庁では、1997年4月以降の主な地震の強震観測データをCSV形式で公開しています。
今回はこれをPythonでフーリエ変換し、スペクトル解析してみます。

【NumPy】高速フーリエ変換 (FFT)で振幅スペクトルを計算
Python言語とNumPyを用いて、信号データを高速フーリエ変換(FFT)して周波数スペクトルを求める方法をソースコード付きで解説します。

データの入手先

https://www.data.jma.go.jp/svd/eqev/data/kyoshin/jishin/110311_tohokuchiho-taiheiyouoki/index.html
※一番上にある東日本大震災「宮城県 涌谷町新町」のデータをダウンロードしました

サンプルコード


実行結果

■7.0~8.0[Hz]のバンドパスフィルタをかけた結果

時間軸データ amplitude_7.0_8.0.csv
周波数軸データ spectrum_7.0_8.0.csv

■0.0~0.1[Hz]のバンドパスフィルタをかけた結果

時間軸データ amplitude_0.0_0.1.csv
周波数軸データ spectrum_0.0_0.1.csv

測定データをそのままフーリエ変換する場合は、矩形窓を使っていることになります。

【Python】強震観測データのスペクトル解析(ハミング窓)

今度は窓関数(ハミング関数)を測定データにかけてフーリエ変換し、スペクトル解析してみます。

【Python】SciPy、NumPyで窓関数の作成(ハニング窓、ハミング窓、ブラックマン窓)
Pythonを用いて、ハニング窓、ハミング窓、ブラックマン窓などといった窓関数を作成する方法について紹介します。

サンプルコード


実行結果

■7.0~8.0[Hz]のバンドパスフィルタをかけた結果

時間軸データ amplitude_7.0_8.0.csv
周波数軸データ spectrum_7.0_8.0.csv

■0.5~35.0[Hz]のバンドパスフィルタをかけた結果

時間軸データ amplitude_7.0_8.0.csv
周波数軸データ spectrum_0.5_35.0.csv

ハミング窓を掛けると、サイドローブが小さくなり周波数領域スペクトル解析はしやすくなりますが、フィルタリング処理して時間信号に戻すと波形が大きく変わってしまいます。

コメント

  1. 機械系大学生 より:

    先週くらいから大学の研究のためpythonを勉強し始めたのですが,高速フーリエ変換したものを図示する時,それを両対数グラフとして表す方法がいまいち分からず困っております.そこで,例えばこのページにあるグラフを両対数グラフとして表示するには,何行目にどのようなものを加えれば可能でしょうか.
    ご回答いただけましたら幸いです.