numpy, scipyなどのメモ

工事中

TimeSeriesの計算

四則演算、指数

gwpyのTimeSeriesや FrequencySeriesは、四則演算や指数の演算がそのまま使える。
横軸の配列(timesやfrequencies)が同一なら、TimeSeries同士の演算も可能。

   1 t = np.linspace(start=-10, stop=10, num=201)
   2 A = TimeSeries(np.sin(t), t0=1261872018)                                                                                                                                 
   3 plot = Plot(A, A+1, 2*A, A**2)
   4 ax = plot.gca().legend(['A', 'A+1', '2*A', 'A**2'], fontsize=18)
   5 plot.show()

max, min, peak to peak, RMS

   1 data.max()
   2 data.min()
   3 data.ptp()
   4 data.rms(2) #引数に入れた数字(秒)ごとのRMSのTimeSeries。この場合は2秒

切り取り、継ぎ接ぎ

   1 data1 = data.crop(t1, t2) # t1 ~ t2 の範囲(TPS time)を切り取ったデータになる
   2 data2 = data.crop(t2, t3)
   3 data3 = data1.append(data2) # data1とdata2はサンプリングが同じで連続的である必要がある

Sampling rateを変える (resample)

dowm sampllingもup sampllingもできるが、後者は非推奨。
エイリアシング(折り返しノイズ)は出ないように処理されているようなので、気にしなくていい。(自前でlow passをかけたりすると、その方が悪影響)

   1 import numpy as np
   2 from numpy import random
   3 from gwpy.timeseries import TimeSeries
   4 
   5 TS = TimeSeries(random.random(2**20), sample_rate=16368) /0.003 # ホワイトノイズを16kHz samplingで生成
   6 
   7 asd = TS.asd(16,8,method='median')
   8 asd2 = TS.resample(65472).asd(16,8,method='median')             # 65kHzにup sampling
   9 asd3 = TS.resample(256).asd(16,8,method='median')               # 256Hzにdown sampling
  10 asd4 = TS.lowpass(120).resample(256).asd(16,8,method='median')  # Low passしてから256Hzにdown sampling
  11 
  12 plot=asd.plot(label='16368Hz (original)',
  13               xscale='log', xlim=(0.1,2e4),
  14               yscale='log', ylim=(0.1,10))
  15 
  16 ax = plot.gca()
  17 ax.plot(asd2, label='65472Hz (resample)', alpha=0.8)
  18 ax.plot(asd3, label='256Hz (resample)', alpha=0.8)
  19 ax.plot(asd4, label='256Hz (LP + resample)', alpha=0.5)
  20 ax.set_ylim(1e-3,10)
  21 ax.legend()

スペクトル (FFT, ASD, PSD) の算出

TimeSeriesのmethod .average_fft(), .asd(), .psd() で計算でき、FrequencySeriesとして帰ってくる。
ASDとPSDは実数で、FFTは複素数。

引数1の数字はFFT length、引数2の数字はoverlapで、これを指定することで時系列データをfftlength分の長さ毎にぶつ切りにしてフーリエ変換をする。

TimeSeiresのデータ長さ ≧ fftlength  ≧ overlap ≧ 0

である必要があるので注意。特に値を指定しない場合は、fftlength = データ長, overlap = 0 として処理される。

引数で method を指定すると、ぶつ切りでフーリエ変換された結果の平均(method='welch', もしくは指定しない場合)か中央値(method='median')かを選べる。
averageがつかない.fft()は、ぶつ切りにしないでワンショットのfftを計算するもの。

引数で window を指定すると、窓関数を選べる。(指定しない場合はHann関数)

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を参照。
fftlength とoverlapを引数で指定する。

   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()を参照

KAGRA/Subgroups/PEM/PythonMemoJP/calculation (last edited 2021-07-29 15:34:38 by tatsuki.washimi)

-  ⇤ ← Revision 1 as of 2020-03-18 17:14:03 → 
  Size: 552
  Editor: tatsuki.washimi
  Comment:
+   ← Revision 17 as of 2021-07-29 15:34:38 → ⇥
  Size: 7874
  Editor: tatsuki.washimi
  Comment:
-Deletions are marked like this.
+Additions are marked like this.
 Line 1:
+<<TableOfContents(2)>>
-Line 3:
+Line 5:
-= Time Seriesなどの計算小技 =
+工事中

-----

= TimeSeriesの計算 =
-Line 6:
+Line 12:
-gwpyのTime SeriesやFrequency Seriesは、四則演算や指数の演算がそのまま使える。
+gwpyの[[https://gwpy.github.io/docs/stable/api/gwpy.timeseries.TimeSeries.html#gwpy.timeseries.TimeSeries.psd|TimeSeries]]や
[[https://gwpy.github.io/docs/stable/api/gwpy.frequencyseries.FrequencySeries.html#gwpy.frequencyseries.FrequencySeries|FrequencySeries]]は、
四則演算や指数の演算がそのまま使える。<<BR>>
横軸の配列(timesやfrequencies)が同一なら、TimeSeries同士の演算も可能。
-Line 14:
+Line 23:
+{{attachment:sample_calic.png}}

== max, min, peak to peak, RMS ==
{{{#!python
data.max()
data.min()
data.ptp()
data.rms(2) #引数に入れた数字(秒)ごとのRMSのTimeSeries。この場合は2秒
}}}


== 切り取り、継ぎ接ぎ ==
{{{#!python
data1 = data.crop(t1, t2) # t1 ~ t2 の範囲(TPS time)を切り取ったデータになる
data2 = data.crop(t2, t3)
data3 = data1.append(data2) # data1とdata2はサンプリングが同じで連続的である必要がある
}}}


== Sampling rateを変える (resample) ==
dowm sampllingもup sampllingもできるが、後者は非推奨。<<BR>>
エイリアシング(折り返しノイズ)は出ないように処理されているようなので、気にしなくていい。(自前でlow passをかけたりすると、その方が悪影響)
{{{#!python
import numpy as np
from numpy import random
from gwpy.timeseries import TimeSeries

TS = TimeSeries(random.random(2**20), sample_rate=16368) /0.003 # ホワイトノイズを16kHz samplingで生成

asd = TS.asd(16,8,method='median')
asd2 = TS.resample(65472).asd(16,8,method='median')             # 65kHzにup sampling
asd3 = TS.resample(256).asd(16,8,method='median')               # 256Hzにdown sampling
asd4 = TS.lowpass(120).resample(256).asd(16,8,method='median')  # Low passしてから256Hzにdown sampling

plot=asd.plot(label='16368Hz (original)',
              xscale='log', xlim=(0.1,2e4),
              yscale='log', ylim=(0.1,10))

ax = plot.gca()
ax.plot(asd2, label='65472Hz (resample)', alpha=0.8)
ax.plot(asd3, label='256Hz (resample)', alpha=0.8)
ax.plot(asd4, label='256Hz (LP + resample)', alpha=0.5)
ax.set_ylim(1e-3,10)
ax.legend()
}}}
{{attachment:resample2.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>>
ASDとPSDは実数で、FFTは複素数。

引数1の数字はFFT length、引数2の数字はoverlapで、これを指定することで時系列データをfftlength分の長さ毎にぶつ切りにしてフーリエ変換をする。
{{{
TimeSeiresのデータ長さ ≧ fftlength  ≧ overlap ≧ 0
}}}
である必要があるので注意。特に値を指定しない場合は、fftlength = データ長, overlap = 0 として処理される。

引数で method を指定すると、ぶつ切りでフーリエ変換された結果の平均(method='welch', もしくは指定しない場合)か中央値(method='median')かを選べる。<<BR>>
averageがつかない[[https://gwpy.github.io/docs/stable/api/gwpy.timeseries.TimeSeries.html#gwpy.timeseries.TimeSeries.fft|.fft()]]は、ぶつ切りにしないでワンショットのfftを計算するもの。

引数で window を指定すると、窓関数を選べる。(指定しない場合はHann関数)

== FFT の窓関数について ==
窓関数について参考になるページ<<BR>>
[[https://holometer.fnal.gov/GH_FFT.pdf|Heinzel-Rudiger-Schilling]] <<BR>>
[[https://org-technology.com/posts/scipy-window-function.html|SciPy で使用可能な窓関数]] <<BR>>
[[https://www.onosokki.co.jp/HP-WK/eMM_back/emm150.pdf|小野精器 計測コラム]]

||'''window''' || '''周波数分解能''' || '''主な用途'''                              || '''推奨overlap''' ||
|| hanning     || △                 || 連続的な信号、一般的な信号 (とりあえずコレ) || 50%               ||
|| flattop     || ×                 || 線スペクトルのピーク値に読み取り            ||                  ?||
|| boxcar      || 〇                 || 過渡信号、ハンマリング試験                  ||  0%               ||


gwpyのFFTでは、窓関数による減衰は補正されているっぽい。([[https://gwpy.github.io/docs/stable/api/gwpy.timeseries.TimeSeries.html#gwpy.timeseries.TimeSeries.fft|Note参照]])<<BR>>
試しに同じ時系列波形に対してhanningとflattopでASDを作ってband limited RMSを計算してみたら、同じ値になった。(補正なしだと2倍ほど変わるはず)


== 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}}

-----

= TimeSeriesにフィルターをかける =
== low pass, high pass, band pass ==
引数にカットオフ周波数を入れるだけで簡単にフィルターできる。
{{{#!python
from numpy import random
data = TimeSeries(random.random(2**15), sample_rate=1024) # white noise を生成
data_LP = data.lowpass(100)
data_HP = data.highpass(10)
data_BP = data.bandpass(20,80)

asd = data.asd(8,4)
asd_LP = data_LP.asd(8,4)
asd_HP = data_HP.asd(8,4)
asd_BP = data_BP.asd(8,4)

plot = Plot(asd, asd_LP, asd_HP, asd_BP, xscale='log', yscale='log', xlim=(1,300))
ax = plot.gca()
ax.legend(['Original', 'Low Pass', 'High Pass', 'Band Pass'])
plot.show()
}}}
{{attachment:LP-HP-BP.png}}

== whitening ==
GWpy本家の[[https://gwpy.github.io/docs/stable/examples/timeseries/whiten.html|example]]を参照。<<BR>>
fftlength とoverlapを引数で指定する。 
{{{#!python
from gwpy.timeseries import TimeSeries
data = TimeSeries.get('H1:ASC-Y_TR_A_NSUM_OUT_DQ', 1123084671, 1123084703)
white = data.whiten(4, 2)

from gwpy.plot import Plot
plot = Plot(data, white, separate=True, sharex=True)
plot.axes[0].set_ylabel('Y-arm power [counts]', fontsize=16)
plot.axes[1].set_ylabel('Whitened amplitude', fontsize=16)
plot.show()
}}}
{{https://gwpy.github.io/docs/stable/examples/timeseries/whiten-3.png}}

== ZPK(zeros, poles, gain) filter ==
GWpy本家の[[https://gwpy.github.io/docs/stable/examples/timeseries/filter.html|example]]および[[https://gwpy.github.io/docs/stable/api/gwpy.timeseries.TimeSeries.html#gwpy.timeseries.TimeSeries.zpk|TimeSeries.zpk()]]を参照