当 Wolfram 语言遇上引力波

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/

谢谢正吉的分享！