2015年9月14日,美国的激光干涉引力波观测站 LIGO 记录下这片时空里泛起的一丝波澜,即 GW150914事件([1])。后续的研究表明这正是爱因斯坦广义相对论中的引力波。LOSC上有一篇详实数据处理教程([2]),使用的是 Python语言。现在我们用 Wolfram 语言来处理这次引力波观测数据。
LIGO 在 Hanford 和 Livingston 各有一个观测站,后文以 Hanford 观测站的数据为例,Livingston 观测站数据处理方法相同。建议下载采样率为4096Hz、格式为 HDF5 的数据。
时域处理
01
初步分析
Wolfram 语言可以直接读取 HDF5 格式的文件,以 Hanford 观测站的数据为例,它的功率谱如下。
从数据的功率谱上可以看出,原始数据里含有丰富的低频分量,也存在多个很强的谐波干扰,如40Hz、60Hz。据介绍,这些谐波干扰为系统所固有的。从时域波形上也能大致看出类似现象。
02
带通滤波
根据资料,有效数据的频率范围大概在40Hz到300Hz之间,需要设计FIR 带通滤波器提取数据。通带的范围可以稍微再小一点,这里我们取40Hz到260Hz。Wolfram 语言内置的 BandpassFilter 函数采用数字域角频率,它跟模拟频率的变换关系是:\[Omega]=2\[Pi] f/Subscript[f, s]。
从上面的功率谱上看到,低频分量比我们关注的频段强 60dB 左右,我们的这个滤波器阻带衰减无法把低频分量完全滤除。为此,可以把 2 到 3 个相同的带通滤波器串联使用。
带通滤波后的数据里存在的谐波干扰功率高,带宽又非常窄,可以用陷波器剔除。从功率谱里选取若干个(比如15个)较强的功率点,各点对应的频率即陷波器候选工作频率:
可以在{35.9, 36.7, 40.97, 60.00, 120, 180}Hz等频率上进行陷波。
03
陷波
IIR陷波器可以很直观地设计出来:在共振点上设置一个零点,然后在零点附近设置一个极点。要求陷波器是因果稳定的,所以极点要在单位圆内;考虑实信号,所以把这对零极点的共轭对称点也分别设置为零极点,示意图如下:
这个IIR陷波器可以通过传递函数定义出来,其中\[Mu]表示极点和零点的接近程度。
04
大功告成
将上面的带通滤波和陷波组合起来,引力波的波形就魔术般地呈现在我们眼前:
频域处理
此类观测数据的有一个特点:大量的噪声中夹杂零星的有效数据。对于这样的数据常见的处理方法是先白化,使频谱变平缓、噪声基本变成 "白噪声",有效数据在白噪声中会变得鹤立鸡群,一目了然。
白化处理比较简单,只要先估计出一个比较平滑的低分辨率功率谱,插值后去除频谱即可 。
白化一来消除了强谐波干扰,二来低频部分也被压了下来。现在只用一个带通滤波器就可以提取引力波:
倾听宇宙的声音
LIGO 检测到的引力波频率范围在人耳的听觉范围之内,把引力波信号视为声音采样信号,我们就可以倾听来自宇宙的"啁啾"。不过这个频率很低,耳机、喇叭输出时失真很大,不妨把它整体搬移到较高的频率上,听起来更加悦耳。
在后台发送“引力波”,便可获取该文章的笔记本下载链接。
参考资料
[1]:https://losc.ligo.org/events/GW150914/
[2]:https://losc.ligo.org/tutorials/
谢谢正吉的分享!