用gwpy绘制引力波时频图

文章目录

• • 时频图
  • spectrogram

时频图

一般来说，针对某一物理过程，探测器往往在时域上进行采样，获取数据之后，通过Fourier变换等方法，可将其转换到频域上。这里存在一个问题，即必须假设这种变换关系是不随时间变化的，比如正弦函数

sin

⁡

ω

t

\\sin\\omega t

sinωt，其角频率一直是

ω

\\omega

ω，而对于类似

sin

⁡

ω

t

t

\\sin\\omega_t t

sinωt​t这样频率随时间变化的函数，传统的频域分析就会失效，此时会自然地产生一个新的需求，能否在特定时刻

t

t

t附近，得到信号在普通频率上的分布，由此得到的图像就是时频图。

为了得到时频图，需要进行短时Fourier变换(STFT)，

h

~

(

t

k

,

f

m

)

=

∫

−

∞

+

∞

h

(

τ

)

w

(

τ

−

t

k

)

e

−

2

π

i

f

m

τ

d

τ

\\tilde h(t_k, f_m)=\\int^{+\\infty}_{-\\infty}h(\\tau)w(\\tau-t_k)e^{-2\\pi if_m\\tau}\\mathrm d\\tau

h~(tk​,fm​)=∫−∞+∞​h(τ)w(τ−tk​)e−2πifm​τdτ

式中，

t

k

t_k

tk​是第

k

k

k个时间窗中心，其他部分和传统Fourier变换几乎一致。其离散形式为

h

~

(

t

k

,

f

m

)

=

∑

n

=

0

M

−

1

h

[

n

+

k

L

]

w

[

n

]

e

−

2

π

i

m

n

/

M

\\tilde h(t_k, f_m)=\\sum^{M-1}_{n=0}h[n+kL]w[n]e^{-2\\pi imn/M}

h~(tk​,fm​)=n=0∑M−1​h[n+kL]w[n]e−2πimn/M

式中，

M

,

L

M,L

M,L为每段的样本数和时间步长，

h

~

\\tilde h

h~的功率谱为

P

(

t

k

,

f

m

)

=

2

f

s

S

2

1

T

∣

h

~

(

t

k

,

f

m

)

∣

2

P(t_k, f_m)=\\frac{2}{f_sS_2}\\frac{1}{T}\\vert\\tilde{h}(t_k, f_m)\\vert^2

P(tk​,fm​)=fs​S2​2​T1​∣h~(tk​,fm​)∣2

式中

• $P(t,f)$
• $\tilde{h}(t_k,f_m)$
• $T$
• $S_2$

spectrogram

【spectrogram】是gwpy中用于计算时频图的函数，其主要参数包括

• stride 为时间步长
• fftLength 为每段FFT时长
• overlap 为相邻窗重叠时长 (秒)，通常 = fftlength/2

from gwpy.timeseries import TimeSeries
# hdata = TimeSeries.get("H1", 1126259446, 1126259478) # 获取GW150914
hdata = TimeSeries.read("gw150914.hdf5")

spec = hdata.spectrogram(2, fftlength=1, overlap=.5) ** (1/2.)
plot = spec.imshow(norm="log", vmin=5e-24, vmax=1e-19)
ax = plot.gca()
ax.set_yscale("log")
ax.set_ylim(10, 2000)
ax.colorbar(label=r"Gravitational-wave amplitude [strain/$\\sqrt{\\mathrm{Hz}}$]")
plot.show()

绘图结果为

此为官方示例，地址：Plotting a Spectrogram