文章目录
- 功率谱密度和振幅谱密度
- gwpy实现
gwpy是一款用于引力波数据处理的Python模块，安装和加载引力波数据可见： 用gwpy处理引力波数据。
功率谱密度和振幅谱密度
为了获取谱密度，需要先将时域信号变换到频域，即做Fourier变换
h ~ ( f ) = ∫ − ∞ + ∞ h ( t ) e 2 π i f t d t \tilde h(f)=\int^{+\infty}_{-\infty}h(t)e^{2\pi ift}\mathrm dth~(f)=∫−∞+∞h(t)e2πiftdt
由于引力波数据是离散信号，即h [ n ] = h ( n Δ t ) , n = 0 , ⋯ , N − 1 h[n]=h(n\Delta t), n=0,\cdots,N-1h[n]=h(nΔt),n=0,⋯,N−1，其离散Fourier变换为
h ~ [ k ] = ∑ n = 0 N − 1 h [ n ] e − 2 π i k n / N , k = 0 , ⋯ , N − 1 \tilde{h}[k]=\sum^{N-1}_{n=0}h[n]e^{-2\pi ikn/N}, k=0,\cdots, N-1h~[k]=n=0∑N−1h[n]e−2πikn/N,k=0,⋯,N−1
对于平稳随机过程，功率谱密度定义为
S h ( f ) = 2 lim T → ∞ 1 T E [ ∣ h ~ T ( f ) ∣ 2 ] , f ⩾ 0 S_h(f)=2\lim_{T\to\infty}\frac{1}{T}\mathbb E\left[\vert \tilde{h}_T(f)\vert^2 \right], f\geqslant 0Sh(f)=2T→∞limT1E[∣h~T(f)∣2],f⩾0
式中，
- h ~ T ( f ) \tilde{h}_T(f)h~T(f)是时长T TT截断信号的Fourier变换
- E \mathbb EE表示系综平均
S h ( f ) S_h(f)Sh(f)的物理意义是，在区间[ f , f + d f ] [f, f+\mathrm df][f,f+df]之间的平均功率。
振幅谱密度则被定义为
S h 1 / 2 ( f ) = S h ( f ) S^{1/2}_h(f)=\sqrt{S_h(f)}Sh1/2(f)=Sh(f)
gwpy实现
在gwpy的TimeSeries类中，提供了生成频谱的方法，其中asd用于计算振幅谱密度(Amplitude Spectral Density, ASD)；psd用于计算功率谱密度(Power Spectral Density, PSD)。二者输入参数相同，包括
- fftlength 每段Fourier变换的时长
- overlap 相邻段重叠时长（秒）
- window 默认’hann’，时域加窗函数，抑制频谱泄漏
- method 默认 ‘welch’，为PSD估计算法
其ASD图像如下
在此图中，小于50Hz的区域噪声较大，需要滤掉，60Hz和120Hz存在一个尖峰，此为美国电网的交流电频率。绘图代码如下。
fromgwpy.timeseriesimportTimeSeries# hdata = TimeSeries.get("H1", 1126259446, 1126259478) # 获取GW150914hdata=TimeSeries.read("gw150914.hdf5")hdata.asd(8,4).plot().show()
