![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
#author("2021-04-10T22:38:34+09:00","default:LabMember","LabMember") #author("2021-04-10T22:45:53+09:00","default:LabMember","LabMember") [[PythonTutorial]] * Pythonを使ったデータ解析その3 [#xc25415b] #contents ** 準備 [#d44aa4e0] [[ex3.zip>https://granite.phys.s.u-tokyo.ac.jp/wiki/Lab/?plugin=attach&pcmd=open&file=ex1.zip&refer=PythonTutorialExample1]] [[ex3.zip>https://granite.phys.s.u-tokyo.ac.jp/takano/misc/ex3.zip]] をダウンロードして展開し，ex3.ipynbを開く *** 中身について [#y6c67512] 今回は地震計(geophone)から取得した時系列データから地面振動のスペクトルを求める． - ex3.ipynb: メインのデータ解析をするnotebook - [x/y/z]_timetrend.txt: [x/y/z]軸方向の地震計の信号 ** FFTの計算 [#e854f6f1] scipy.signalモジュールを使ってFFTをする． *** パワースペクトル密度 [#h3d4f719] パワースペクトル密度を計算するには， scipy.signal.welch を使う． frequency, PSD = welch(時系列データ, fs=サンプリングレート, window=窓関数（デフォルトはハニング）, nperseg=1セグメントの長さ, noverlap=オーバーラップする長さ（デフォルトはnperseg/2）) welch はウェルチ法によるパワースペクトルの推定を行う関数．出力は周波数とパワースペクトル密度． ウェルチ法は， + 時系列をいくつかのセグメントに分ける + それぞれのセグメントに窓関数をかける + それぞれのセグメントについてFFTをしてパワースペクトル密度を計算する + 求めたパワースペクトル密度の平均をとる という手順でパワースペクトル密度を推定する方法．特徴的なのは，セグメントをオーバーラップして分けていること． オーバーラップは50%にすることが多い．この時，例えばセグメント長全データ長の半分に取ればセグメント数は3， セグメント長全データ長の1/3に取ればセグメント数は5，といった感じになる． 窓関数はハニング，ハミング，矩形，…と色々あるが，とりあえず何も指定しないでハニング窓を使えば大体の場合はOK． なお，他にもscipy.signal.periodogram というパワースペクトル密度を計算する関数があるが，これの場合は窓関数をかけず，平均も取らないで計算をする． なのであまり使う機会はない． *** コヒーレンス関数 [#dba36609] コヒーレンス関数を計算するには scipy.signal.coherence を用いる． frequency, coherence = coherence(時系列データ1, 時系列データ2, fs=サンプリングレート, window=窓関数, nperseg=1セグメントの長さ, noverlap=オーバーラップする長さ) 基本的な使い方は welch と同じ．ただし coherence の場合は2つの時系列データを引数に入れる． ** キャリブレーション [#c01bc35f] 地震計の応答の逆数をかけて地震計の信号から地面振動のスペクトルへキャリブレーションする． *** 地震計の応答 [#ta6e2bc6] 地震計は，懸架された振り子と地面の間の相対速度が電圧として出力される．元の地面振動に換算するにはこれの逆数をかけてキャリブレーションする．~ 細かい式はnotebookを参照．
