4479
Comment:
|
5400
|
Deletions are marked like this. | Additions are marked like this. |
Line 9: | Line 9: |
= TimeSeriesの計算小技・注意事項 = | = TimeSeriesの計算 = |
Line 14: | Line 14: |
四則演算や指数の演算がそのまま使える。 | 四則演算や指数の演算がそのまま使える。<<BR>> |
Line 25: | Line 25: |
== FFT length の注意 == | ----- |
Line 27: | Line 27: |
TimeSeriesからスペクトルを計算する際、fftlengthの値が小さいと下図左のように謎の冪スペクトルに埋もれてしまうことがある。(ドリフト+窓関数の影響?)<<BR>> 十分なfftlengthをとればいいのだが、平均回数を稼ぎたいときや短時間のみを見たい場合にはそうもいかない。<<BR>> そんなときは先に適当な[[https://gwpy.github.io/docs/stable/api/gwpy.timeseries.TimeSeries.html#gwpy.timeseries.TimeSeries.highpass|highpass]]フィルターをかけることで、下図右のように回避できる。<<BR>> {{attachment:wo_hp.png}} {{attachment:w_hp.png}} |
= スペクトル (FFT, ASD, PSD) の算出 = TimeSeriesのmethod [[https://gwpy.github.io/docs/stable/api/gwpy.timeseries.TimeSeries.html#gwpy.timeseries.TimeSeries.average_fft|average_fft()]], [[https://gwpy.github.io/docs/stable/api/gwpy.timeseries.TimeSeries.html#gwpy.timeseries.TimeSeries.asd|asd()]], [[https://gwpy.github.io/docs/stable/api/gwpy.timeseries.TimeSeries.html#gwpy.timeseries.TimeSeries.psd|psd()]] で計算でき、[[https://gwpy.github.io/docs/stable/api/gwpy.frequencyseries.FrequencySeries.html#gwpy.frequencyseries.FrequencySeries|FrequencySeries]]として帰ってくる。<<BR>> 引数1の数字はFFT length、引数2の数字はoverlapで、これを指定することで時系列データをfftlength分の長さ毎にぶつ切りにしてフーリエ変換をする。 {{{ TimeSeiresのデータ長さ ≧ fftlength ≧ overlap ≧ 0 }}} である必要があるので注意。 |
Line 48: | Line 54: |
== FFT length の注意 == TimeSeriesからスペクトルを計算する際、fftlengthの値が小さいと下図左のように謎の冪スペクトルに埋もれてしまうことがある。(ドリフト+窓関数の影響?)<<BR>> 十分なfftlengthをとればいいのだが、平均回数を稼ぎたいときや短時間のみを見たい場合にはそうもいかない。<<BR>> そんなときは先に適当な[[https://gwpy.github.io/docs/stable/api/gwpy.timeseries.TimeSeries.html#gwpy.timeseries.TimeSeries.highpass|highpass]]フィルターをかけることで、下図右のように回避できる。<<BR>> {{attachment:wo_hp.png}} {{attachment:w_hp.png}} |
Contents
numpy, scipyなどのメモ
工事中
TimeSeriesの計算
四則演算、指数
gwpyのTimeSeriesや FrequencySeriesは、 四則演算や指数の演算がそのまま使える。
横軸の配列(timesやfrequencies)が同一なら、TimeSeries同士の演算も可能。
スペクトル (FFT, ASD, PSD) の算出
TimeSeriesのmethod average_fft(), asd(), psd() で計算でき、FrequencySeriesとして帰ってくる。
引数1の数字はFFT length、引数2の数字はoverlapで、これを指定することで時系列データをfftlength分の長さ毎にぶつ切りにしてフーリエ変換をする。
TimeSeiresのデータ長さ ≧ fftlength ≧ overlap ≧ 0
である必要があるので注意。
FFT の窓関数について
窓関数について参考になるページ
Heinzel-Rudiger-Schilling
SciPy で使用可能な窓関数
小野精器 計測コラム
window |
周波数分解能 |
主な用途 |
推奨overlap |
hanning |
△ |
連続的な信号、一般的な信号 (とりあえずコレ) |
50% |
flattop |
× |
線スペクトルのピーク値に読み取り |
? |
boxcar |
〇 |
過渡信号、ハンマリング試験 |
0% |
gwpyのFFTでは、窓関数による減衰は補正されているっぽい。(Note参照)
試しに同じ時系列波形に対してhanningとflattopでASDを作ってband limited RMSを計算してみたら、同じ値になった。(補正なしだと2倍ほど変わるはず)
FFT length の注意
TimeSeriesからスペクトルを計算する際、fftlengthの値が小さいと下図左のように謎の冪スペクトルに埋もれてしまうことがある。(ドリフト+窓関数の影響?)
十分なfftlengthをとればいいのだが、平均回数を稼ぎたいときや短時間のみを見たい場合にはそうもいかない。
そんなときは先に適当なhighpassフィルターをかけることで、下図右のように回避できる。
TimeSeriesにフィルターをかける
引数にカットオフ周波数を入れるだけで簡単にフィルターできる。
low pass, high pass, band pass
1 from numpy import random
2 data = TimeSeries(random.random(2**15), sample_rate=1024) # white noise を生成
3 data_LP = data.lowpass(100)
4 data_HP = data.highpass(10)
5 data_BP = data.bandpass(20,80)
6
7 asd = data.asd(8,4)
8 asd_LP = data_LP.asd(8,4)
9 asd_HP = data_HP.asd(8,4)
10 asd_BP = data_BP.asd(8,4)
11
12 plot = Plot(asd, asd_LP, asd_HP, asd_BP, xscale='log', yscale='log', xlim=(1,300))
13 ax = plot.gca()
14 ax.legend(['Original', 'Low Pass', 'High Pass', 'Band Pass'])
15 plot.show()
whitening
GWpy本家のexampleを参照
1 from gwpy.timeseries import TimeSeries
2 data = TimeSeries.get('H1:ASC-Y_TR_A_NSUM_OUT_DQ', 1123084671, 1123084703)
3 white = data.whiten(4, 2)
4
5 from gwpy.plot import Plot
6 plot = Plot(data, white, separate=True, sharex=True)
7 plot.axes[0].set_ylabel('Y-arm power [counts]', fontsize=16)
8 plot.axes[1].set_ylabel('Whitened amplitude', fontsize=16)
9 plot.show()
ZPK(zeros, poles, gain) filter
GWpy本家のexampleおよびTimeSeries.zpk()を参照