気象庁の強震観測データから，フーリエスペクトルの図を再現したいです．

Question

前提・実現したいこと

先月からpythonを勉強し始めた者です．今までプログラムのプの字も触れてこなかったので，おかしな点がありましたらご指摘お願いいたします．

気象庁の強震観測データから，フーリエスペクトルの図を再現したいです．

様々な方のFFTに関するの記事を拝見・拝借させていただき，少しでも再現できるように努力したのですが，図の概形は再現できても細かい値，特にフーリエスペクトルの軸がどうしても再現することができません．(エラーは発生していません)

以下に自分が再現したい気象庁の図と，様々な方の記事を参考に作成したソースコード・その実行結果を乗せるので，どこを直せば再現できるかのアドバイスをいただきたいです．

図の説明
図1：自分が再現したい図(右上の図です)　　　　　　　　　　　　　引用①
図2：一枚目の各図の説明　再現したいのは，この図でいう8です　　 引用②
図3：実行した結果

図1　自分が再現したい図(右上の図です)

図2　一枚目の各図の説明　再現したいのは，この図でいう8です

図3　実行した結果

【補足】
複数のサイトを参考に作成したものなので，変な書き方になっているかもしれないです

(引用①：
https://www.data.jma.go.jp/svd/eqev/data/kyoshin/jishin/2102132307_fukushima-oki/wave/wav20210213000157015.png
引用②：
https://www.data.jma.go.jp/svd/eqev/data/kyoshin/jishin/kaisetsu.html)

発生している問題・エラーメッセージ

エラーメッセージ

該当のソースコード

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
%matplotlib inline

# 簡単な信号の作成
N = 30000 # サンプル数
dt = 0.01 # サンプリング周期(sec):100ms =>サンプリング周波数100Hz

t = np.arange(0, N*dt, dt) # 時間軸
fq = np.linspace(0, 1.0/dt, N)  # 周波数軸 linspace(開始,終了,分割数)

# CSVのロード
df = pd.read_csv("/content/drive/MyDrive/spe-ana/quake-CSV/2021-0213-NIHONMATSU.csv",skiprows=6)

# 1列目の加速度データを参照
freq = df.iloc[:,0]

# 高速フーリエ変換(FFT)
F = np.fft.fft(freq)

# FFTの複素数結果を絶対に変換
F_abs = np.abs(F)

# 振幅をもとの信号に揃える
F_abs_amp = F_abs / N * 2 # 交流成分はデータ数で割って2倍
F_abs_amp[0] = F_abs_amp[0] / 2 # 直流成分（今回は扱わないけど）は2倍不要

#単位のつじつま合わせ 縦軸(フーリエスペクトル gal・sec)
FS = (F_abs_amp) * t

#周波数(freq.)を周期(period)に変換
peri = 1/fq

# グラフ表示
fig = plt.figure(figsize=(12, 4))

# 信号のグラフ（時間軸）
ax2 = fig.add_subplot(121)
plt.xlabel('time(sec)', fontsize=14)
plt.ylabel('Gal(cm/s^2)', fontsize=14)
plt.plot(t, freq)
plt.xlim(10,100)
plt.ylim(-650,650)
plt.grid()

# FFTのグラフ(周期軸)
ax2 = fig.add_subplot(122)
plt.xlabel('period(sec)', fontsize=14)
plt.ylabel('Fourier Spectrum(Gal sec)', fontsize=14)
plt.plot(peri[:int(N/2)+1], FS[:int(N/2)+1]) # ナイキスト定数まで表示
plt.grid()

#両対数表示　範囲指定
plt.yscale("log")
plt.xscale("log")
plt.xlim(0.01,100)
plt.ylim(0.01,1000)

試したこと

29行目
単位のつじつま合わせ 縦軸(フーリエスペクトル gal・sec)
FS = (F_abs_amp) * t

で，F_abs_ampを2乗したり(F_abs_amp　→　F_abs_amp ** 2)，そもそもtをかけないパターンも試したのですが，特に目標に近づくことはありませんでした．

補足情報（FW/ツールのバージョンなど）

pythonは先月時点で入れたので，当時の最新バージョンだと思います．
また，Google colab.で作成・実行すべてを行っています．
知り合いにおすすめされたもので，特にこだわりがあるわけではありません．

平滑化・カットオフ等のフィルター処理は省いております．
最初はどちらとも行っていたのですが，ややこしくなってとりあえず素の状態を確かめたくて無くしました．それが原因だったら申し訳ありません．

ひとえに自分の勉強不足によるところだと自覚しております．どうか初心者にもわかるご説明をお願いいたします．

Comments

[jbpb0]

データの入手手順を教えてください

[python_gekimuzu]

すみません。データとはどのデータを指しているものでしょうか。

このプログラムを書いた順番でした大まかに以下の通りです。
①「CSVデータをFFTしてグラフ化する」ソースコードの例があるサイトを参考に大枠を書く
②気象庁から強震観測データをマイドライブに保存し、①のCSV参照部分を自分用に修正
③このままグラフ化すると縦軸横軸ともに単位が違うため、つじつまを合わせるためのFSとperiを作成
④グラフが見やすいように両対数グラフにするコードと、範囲指定をするコードを追加

このような流れで作成したと思います。

質問者様が思っている回答と違いました申し訳ありません。
その時は御手数ですが改めて詳しくご質問内容を伝えていただけると幸いです。

よろしくお願い致します。

[jbpb0]

> データとはどのデータを指しているものでしょうか。

は、

> ②気象庁から強震観測データをマイドライブに保存

の「強震観測データ」のことです
そのデータをダウンロードできるURLとか、質問に掲載されてる画像の地震のデータの場所と日時とかの情報

質問に掲載のコードで、現象を再現できるようにするには、データが必要ですので

[python_gekimuzu]

「強震データ」を貼ることを失念しておりました。的外れの回答をしてしまい申し訳ありません。
それと、ご指摘ありがとうございます。

「強震データ」は以下の気象庁のリンクから
『福島県
二本松市油井＊
５強
5.4
696.9
601.5
405.7
221.1
110.0
波形
ダウンロード』
の『ダウンロード』からダウンロードいたしました。

https://www.data.jma.go.jp/svd/eqev/data/kyoshin/jishin/2102132307_fukushima-oki/index.html

[jbpb0]

FS = (F_abs_amp) * t
↓ 変更
FS = (F_abs_amp) / dt * 2
wnum = np.int(np.round(0.4 / fq[1]))
w = np.ones(wnum) / wnum
FS = np.convolve(FS, w, mode='same')

で、近くはなります

[python_gekimuzu]

早急なご回答誠にありがとうございます。
手元にPCが無いので、明日実際に打ち込んで確認致します。

[python_gekimuzu]

今確認しました。
すごいです！とても再現度の高いグラフが表示出来ました！
本当にありがとうございます。

[jbpb0]

質問に掲載のコードのように「t」を掛けるのは間違いだと思います
「t」はデータの先頭の時刻からの経過時間なので、データ取得開始がもっと遅ければ(データの先頭を少し削除する)「t」は小さくなり、それを掛け算した計算結果が小さくなりますけど、データ取得開始のタイミングが早いか遅いかでスペクトルの振幅が変わるわけないですよね
(地震が始まるよりも早くデータ取得が始まってるのが前提)

そこまでは分かるのですが、物理的な数値の意味を考えてるわけではなく、グラフの高さが合うように「t」の代わりに「2 / dt」を掛けるようにしたので、本当にそれで合ってるのかは分かりません
そこは質問者さんが考えてください

「wnum =...」からは「0.4Hzのバンド幅で平滑化を行っています。」を移動平均でやってますけど、平滑化の手法によって結果が変わるので、それでいいのかは分かりません
https://smo.kenken.go.jp/~kashima/ja/viewwave/usage/settings/calculation
とか見ると「Parzen窓」で平滑化してるみたいなので、「w =...」をParzen窓に変えたら、もっと近づくかもしれません (未確認)

[python_gekimuzu]

tをかけることについて、おっしゃる通りだと思います。
勉強して、dtをなぜかけるかについての物理的な理解ができたら、改めて報告したいと思います。

ありがとうございます。

[jbpb0]

この修正で合ってると確信したら、質問者さんが書いた回答を「自己解決」にしてください

[python_gekimuzu]

承知致しました。
回答してくださったのがjbpb0様でよかったです。
ご丁寧にありがとうございました。

Answer 1

jbpb0様が以下の解決方法を提案して下さいました．

FS = (F_abs_amp) * t
↓ 変更
FS = (F_abs_amp) / dt * 2
wnum = np.int(np.round(0.4 / fq[1]))
w = np.ones(wnum) / wnum
FS = np.convolve(FS, w, mode='same')

以上です．
これで再現度の高いグラフを表示することが出来ました．
ありがとうございます．

解決していただいたものをどこに書けばよいのかわからなかったので，自己解決ではないのですがここに書きました．
書く場所や，使い方が違っていたらすみません．

Comments

[jbpb0]

質問に掲載のコードのように「t」を掛けるのは間違いだと思います
「t」はデータの先頭の時刻からの経過時間なので、データ取得開始がもっと遅ければ(データの先頭を少し削除する)「t」は小さくなり、それを掛け算した計算結果が小さくなりますけど、データ取得開始のタイミングが早いか遅いかでスペクトルの振幅が変わるわけないですよね
(地震が始まるよりも早くデータ取得が始まってるのが前提)

そこまでは分かるのですが、物理的な数値の意味を考えてるわけではなく、グラフの高さが合うように「t」の代わりに「2 / dt」を掛けるようにしたので、本当にそれで合ってるのかは分かりません
そこは質問者さんが考えてください

「wnum =...」からは「0.4Hzのバンド幅で平滑化を行っています。」を移動平均でやってますけど、平滑化の手法によって結果が変わるので、それでいいのかは分かりません
https://smo.kenken.go.jp/~kashima/ja/viewwave/usage/settings/calculation
とか見ると「Parzen窓」で平滑化してるみたいなので、「w =...」をParzen窓に変えたら、もっと近づくかもしれません (未確認)