2016年，用于脉冲星观测研究的500 m口径球面射电望远镜(Five-hundred-meter Aperture Spherical radio Telescope, FAST) 在贵州平塘落成，我国实现了依靠国内设备发现脉冲星的零的突破。脉冲星是宇宙中一类奇妙的天体，是验证强引力场、强磁场和高密度等极端环境下物理规律的天然实验室^[1]。作为脉冲星研究的利器，500 m口径球面射电望远镜采用19波束接收机每天8 h对天区进行大规模巡天，搜寻并接收脉冲星发出的电磁波信号，依据电磁脉冲信号确定天区中是否存在脉冲星。众所周知，脉冲星距地球约3 000~ 55 000光年，脉冲星发出的信号到达地球且成功被专用仪器接收的信号能量十分微弱，整个脉冲信号淹没在极强的噪声中^[2-3]。因此，去噪在脉冲星信号的检测和发现中至关重要。

小波变换作为一种多分辨分析方法能够同时在时域和频域对信号进行分析，由于具有时域局部化和多分辨特性而称为数学显微镜，非常适合处理非平稳信号^[4]。脉冲星信号是一种典型的非平稳信号，文[5]把小波分析应用于脉冲星信号去噪, 认为小波分析可以更好地对非平稳信号进行消噪处理，因而小波变换替代了传统的傅里叶变换，在脉冲星信号处理方面大显身手。此后研究者对小波变换去除脉冲星信号噪声进行了大量的研究，文[6]探索了极大极小原理、史坦(Stein) 无偏似然原理的脉冲星降噪方法；文[7]尝试改进阈值函数以提高脉冲星信号去噪效果，提出了软硬折中阈值函数、指数型阈值函数处理脉冲星小波分解系数，获得了较好的效果。

小波变换既能有效去除信号的噪声，又能较好地保持原有信号的数据，是信号处理的理想工具。但小波变换主要针对低频信息的提取和分析，经常忽略高频信息，这对整个信号分析的精确度产生较大的影响。小波变换去噪是在小波函数中选择一种合适的信号处理基函数，利用基函数进行小波分解，根据不同阈值函数的处理结果对信号进行重构，从而实现小波变换去噪。小波变换去噪因基函数和阈值函数的不同而表现出不同的去噪效果，去噪效果有一定的随机性^[8]。

随着小波理论的完善，小波包分析得到了进一步的发展。小波包分析是小波分解的一般化，具有更强的去噪能力，既能像小波分析一样对信号的低频部分进行分解，也能够对信号的高频部分进行分解，更好地提取各频段的有用信息。小波包去噪方法应用在图像、地震、心电、海波等信号领域，文[9]尝试将小波包滤波方法与X射线脉冲星噪声滤除相结合，但没有对小波包去噪进行系统的分析。本文提出将小波包阈值法应用于脉冲星信号去噪，系统分析各种阈值函数的去噪效果。实验结果显示，小波包阈值法与传统的小波变换法相比，脉冲星信号的去噪效果得到了明显的提高，为脉冲星信号去噪提供了一些思路。

1 小波变换与小波包变换 1.1 小波变换方法

观测信号可以看作为f(t)，满足

$ f(t)=s(t)+n(t), $

(1)

其中，f(t) 为接收的原始信号，即含噪信号；s(t) 为需要提取的去噪后的信号；n(t) 为方差为σ²的噪声，服从n(0, σ²) 分布。一维含噪信号f(t) 按照小波级数展开可以得到n个由多分辨分析理论得出的线性组合，具体的小波变换为

$ W_{f}\left(\psi_{j, k}\right)=2^{\frac{j}{2}} \int_{R} f \overline{\psi\left(2^{j} n-k\right)} \mathrm{d} x, $

(2)

其中，W_f(ψ_{j, k}) 为从信号变换的线性空间L²(R) 中提取的小波基ψ_{j, k}的离散小波变换；ψ(2^jn-k)为小波基函数ψ_{j, k}离散化处理后得到的函数。观测信号的小波变换系数由(2) 式确定。

1.2 小波包变换方法

小波包变换是小波变换的拓展，具备小波变换处理观测信号低频部分的能力，同时能够较好地分解高频部分，收集观测信号各个频段的有用信息，进而提高去噪精度。具体的小波包变换公式为

$ \left\{\begin{array}{l} \varphi(x)=\sqrt{2} \sum\limits_{n \in {\bf{Z}}} h_{n} \varphi(2 x-n) \\ \psi(x)=\sqrt{2} \sum\limits_{n \in {\bf{Z}}} g_{n} \varphi(2 x-n) \end{array},\right. $

(3)

其中，φ(x)和ψ(x) 分别为信号在L²(R) 上多分辨分析中的尺度函数和小波函数；h_n和g_n为滤波器系数，满足

$ \sum\limits_{n \in {\bf{Z}}} h_{n} \varphi(2 x-n)=\delta_{k l}, $

(4)

$ \sum\limits_{n \in {\bf{Z}}} h_{n}=\sqrt{2}, $

(5)

$ g_{n}=(-1)^{2} h_{n-1} . $

(6)

2 小波包阈值去噪法 2.1 小波包分解与重构

小波包理论中的小波包变换具有分解与重组的性质^[10]。小波包分解算法的数学表达式为

$ w_{j+1}^{2 n}(k)=\sum\limits_{l \in {\bf{Z}}} h_{l-2 k} w_{j}^{n}(l), $

(7)

$ w_{j+1}^{2 n+1}(k)=\sum\limits_{l \in {\bf{Z}}} g_{l-2 k} w_{j}^{n}(l), $

(8)

其中，w_jⁿ(k) 表示在分解尺度为j的第n条分叉树上第k个小波包的分解系数；w_j+1²ⁿ(k)和w_j+1²ⁿ⁺¹(k)为w_jⁿ(k) 的两个分叉。相应的小波包重构算法的数学表达式为

$ w_{j}^{n}(k)=\sum\limits_{l \in {\bf{Z}}} h_{k-2 l} w_{j+1}^{2 n}(l)+\sum\limits_{l \in {\bf{Z}}} g_{k-2 l} w_{j+1}^{2 n+1}(l), $

(9)

其中，w为信号经过小波包分解后得到的分解系数；h和g为小波包分解的滤波器系数；j和n为小波包分解中各分解节点的编号；l和k为信号进行小波包分解所在的分解层数。

2.2 最优小波包基的确定

小波包分析处理观测信号时，首先对信号进行小波包分解。小波包分解与小波分解有所不同，小波包分解是按照二叉树的形式，信号的低频和高频部分都进行相应的完全二叉树形式的分解。信号分解的层数直接影响信号的去噪效果，处理的信号在逐层分解过程中必须确定一个最优分解层数。分解层次过低造成频段中混合有用信号，分解层次过高导致结果有误差^[11]。在小波分析去噪过程中，分解层数过多造成信号信息量的丢失，从而导致信噪比下降，影响最终的去噪效果。在小波分析的实际应用中，分解层数一般取3~5^[12]，本文取5进行对比分析。

小波包变换中存在多种小波构成的基底库，因此需要根据信号的特性，结合代价函数的评价标准进行选择。选择小波函数时需要充分考虑脉冲信号的特点，选择能够满足给定区间的紧支性和足够的消失矩，以便能够更加有效地去除噪声，提取真正的目标脉冲信号。为了选出符合脉冲星信号特点的小波基函数，本节选取Symlets, Coiflets, Daubechies和Biorthogonal系列小波以及Haar小波进行比较分析。为了便于比较分析，本节统一采用固定阈值规则，利用不同的小波基进行5层分解，分析结果如表 1、表 2、表 3和表 4。其中N表示小波阶数；Sym N, Coif N, Db N和Bior N分别表示Symlets, Coiflets, Daubechies和Biorthogonal系列小波。从表 1~表 4统计的信噪比和均方根误差可以看出，利用Sym N，Coif N，Db N，Bior N和Haar小波进行脉冲星信号去噪整体上效果较好，其中，Sym N小波的效果最好，Coif N小波次之，Bior N小波最差。此外，小波基阶数不同，去噪效果也不同，就单个小波基而言，Sym5小波去噪效果最好。因此，本文选取近似对称的紧支集正交小波Sym5作为小波基函数对脉冲信号进行小波包分析，分解层数为5层，得到最优小波包基，也就是小波包分解的最优化树。

2.3 阈值函数的选择

利用小波包变换处理信号的过程中经常出现翻转现象，也就是在信号经过高通滤波器时的一种频率交叠现象，导致分解的小波包系数出现频率大小顺序错位，即信号的低频部分从小到大排列，信号的高频部分则从大到小排列。因此，在进行小波包阈值处理之前需要先对变换后的小波包系数进行排序，再根据排序后的小波包系数选取对应的信号处理阈值。

小波包分析常用的阈值选取准则有4种，分别为固定形式(Sqtwolog) 阈值准则、自适应(Rigrsure) 阈值准则、启发式(Heursure)阈值准则和极大极小值(Minimaxi) 阈值准则，本文采用固定形式阈值准则对信号进行处理。信号分解时，每层的系数选定阈值以后，需要确定对信号处理的阈值函数。常用的阈值函数有软阈值函数和硬阈值函数，其中，硬阈值函数采用将绝对值不大于阈值λ的所有元素用0取代，这样容易造成处理后的信号不连续。硬阈值函数的表达式为

$ \overline{w_{j}^{n}}= \begin{cases}w_{j}^{n} & \left|w_{j}^{n}\right| \geqslant \lambda \\ 0 & \left|w_{j}^{n}\right|<\lambda\end{cases}. $

(10)

软阈值函数克服了硬阈值函数的不足，有效避免信号间断，具有较好的连续性。然而当小波系数较大时，经过阈值函数处理后容易造成信号的高频部分损失，进而影响信号重构。软阈值函数的表达式为

$ \overline{w_{j}^{n}}=\left\{\begin{array}{cc} \operatorname{sign}\left(w_{j}^{n}\right)\left(\left|w_{j}^{n}\right|-\lambda\right) & \left|w_{j}^{n}\right| \geqslant \lambda \\ 0 & \left|w_{j}^{n}\right|<\lambda \end{array}\right.. $

(11)

其中，w_jⁿ为信号在分解尺度为j的第n条分叉树上的分解系数；λ为选定的阈值；w_jⁿ为经过阈值函数处理后得到的分解系数。鉴于软、硬阈值函数处理信号时存在各自的不足，都有一定的局限性，研究人员提出了软硬阈值折中法，具体的表达式为

$ \overline{w_{j}^{n}}=\left\{\begin{array}{cc} \operatorname{sign}\left(w_{j}^{n}\right)\left(\left|w_{j}^{n}\right|-\alpha \lambda\right) & \left|w_{j}^{n}\right| \geqslant \lambda \\ 0 & \left|w_{j}^{n}\right|<\lambda \end{array}\right., $

(12)

其中，α取值范围为0≤α≤1。当α取0时相当于硬阈值函数，当α取1时相当于软阈值函数，适当调整α的大小可以得到更优的信号去噪效果。

3 脉冲星信号去噪分析

实验采用的脉冲信号来源于帕克斯(Parkes) 脉冲星观测数据，该数据经过PRESTO去干扰、消色散后得到.dat数据文件，以二进制形式记录，采样频率为1 024 Hz。

图 1为原始脉冲信号消色散后得到的图形，图 2为原始脉冲信号消色散后的功率谱。本文对目标脉冲星信号文件按照Sym5进行小波包分解，分解层数为5层，分别选择小波分析的硬阈值、软阈值和软硬折中阈值函数，采用固定形式阈值准则进行处理。

为了比较采用不同去噪方法获得信号的去噪效果，本文参照文[13]采用信噪比、峰值信噪比和平滑度作为评价指标比较最终的去噪效果，信噪比越大、峰值信噪比越大、平滑度越小，去噪效果越好。具体公式定义如下：

(1) 信噪比的计算公式为

$ S N R=10 \log _{10}\left\{\left[\sum\limits_{i=1}^{m} f(i)^{2}\right] /\left[\sum\limits_{i=1}^{m} s(i)^{2}-f(i)^{2}\right]\right\}; $

(13)

(2) 均方根误差的计算公式为

$ R M S E=\sqrt{\frac{1}{m} \sum\limits_{i=1}^{m}[s(i)-f(i)]^{2}}; $

(14)

(3) 峰值信噪比是信号的最大功率和噪声功率的比值，用PSNR表示，具体公式为

$ P S N R=10 \log _{10}\left\{\max \left[f(i)^{2}\right] / \frac{1}{m}\left[\sum\limits_{i=1}^{m} s(i)^{2}-f(i)^{2}\right]\right\}; $

(15)

(4) 平滑度是去噪后信号差分数的方差根和信号之间差分数的方差根之比，常用r表示，具体公式为

$ r=\sum\limits_{i=1}^{m-1}\left\{[f(i+1)-f(i)]^{2}\right\} / \sum\limits_{i=1}^{m-1}\left\{[s(i+1)-s(i)]^{2}\right\} . $

(16)

以上各式中，s(i) 为含噪信号；f(i) 为去噪后的信号。

小波变换和小波包变换去噪法采用不同阈值处理的结果图和功率谱图见图 3~图 10，去噪后的各评价指标见表 5。由表 5可知，采用小波包阈值法去噪相比小波变换能够取得更好的效果，尽管小波包变换各类阈值函数处理之后，得到的均方根误差、功率谱变化并不明显，但是在信噪比、峰值信噪比以及平滑度方面都有相应的提高。如实验显示，采用小波包硬阈值法与小波变换硬阈值法相比，平滑度降低至0.000 068，信噪比提高至128.574 093 dB，峰值信噪比提高至50.216 641 dB，小波包软阈值、软硬折中阈值法去噪效果的各项评价指标都有相应的提高。

4 总结与展望

本文针对小波变换在去噪方面存在的不足，提出了使用小波包阈值法对脉冲星信号去噪。小波包阈值法首先需要对小波包分解后的系数按从小到大的顺序排列，然后采用固定阈值准则求出阈值，利用常见的软阈值函数、硬阈值函数和软硬折中阈值函数对小波包变换后的系数进行处理，并对处理后的信号进行重构，得到去噪后的信号。本文首先从理论方面对小波包阈值法去噪进行了相关论述，之后通过实验验证了小波包阈值法去除脉冲星信号噪声的效果优于小波分解法。为了使实验结果得到量化标准，本文采用信噪比、峰值信噪比、平滑度以及均方根误差进行量化，比较不同方法去噪的效果。实验结果表明，小波包阈值法较小波变换有更好的去噪效果。由于目前本文只使用Sym5小波基和固定阈值准则对选用的帕克斯望远镜脉冲星信号进行了测试，是否适用于其他脉冲星信号，是否其他小波基或者阈值准则也能够获得同样的去噪效果有待进一步探索研究。

致谢: 感谢贵州省信息与计算重点实验室提供数据支持。