广东工业大学学报  2021, Vol. 38Issue (5): 48-51.  DOI: 10.12052/gdutxb.200107.
0

引用本文 

谢伟翔, 莫艳. 基于复值稀疏贝叶斯的Hilbert变换计算方法[J]. 广东工业大学学报, 2021, 38(5): 48-51. DOI: 10.12052/gdutxb.200107.
Xie Wei-xiang, Mo Yan. A Complex Sparse Bayesian Method to Compute Hilbert Transform[J]. JOURNAL OF GUANGDONG UNIVERSITY OF TECHNOLOGY, 2021, 38(5): 48-51. DOI: 10.12052/gdutxb.200107.

基金项目:

国家自然科学基金资助项目(11801095)

作者简介:

谢伟翔(1994–),男,硕士研究生,主要研究方向为计算及应用调和分析。

通信作者

莫艳(1986–),女,副教授,主要研究方向为计算及应用调和分析,E-mail:marvel-2008@163.com

文章历史

收稿日期:2020-08-24
基于复值稀疏贝叶斯的Hilbert变换计算方法
谢伟翔, 莫艳    
广东工业大学 应用数学学院,广东 广州 510520
摘要: Hilbert变换是信号分析及信号处理中的重要工具, 由于Cauchy核在原点的奇性增加了Hilbert变换计算的难度。最近, 研究者们首次提出了利用复解析的方法来计算Hilbert变换的自适应傅里叶分解(Adaptive Fourier Decomposition, AFD)方法。AFD方法通过参数化的Szegö核的线性组合来自适应逼近解析信号从而求得原始实值信号的Hilbert变换。与传统计算Hilbert变换的方法相比, AFD方法可以给出逼近的解析表达式且适用范围更广。然而AFD方法在根据最大选择原理选择参数时需要穷尽单位开圆盘的所有点, 这是非常耗时的。稀疏贝叶斯学习是近年来机器学习研究的热点, 基于Szegö核的复值稀疏贝叶斯学习算法能提供稀疏的有理逼近。本文将提出基于Szegö核的复值稀疏贝叶斯学习算法来近似计算Hilbert变换, 该方法具有AFD方法的优点且可以不需要参数控制进行迭代优化, 运算速度快。实验结果表明, 所提方法是有效的。
关键词: 复值稀疏贝叶斯    Szegö核    Hilbert变换    解析函数    
A Complex Sparse Bayesian Method to Compute Hilbert Transform
Xie Wei-xiang, Mo Yan    
School of Applied Mathematics, Guangdong University of Technology, Guangzhou 510520, China
Abstract: Hilbert transform is an important tool for signal analysis and processing. Due to the singularity of the Cauchy kernel at the origin, the calculation of the Hilbert transform becomes very difficult. Recently, researchers firstly propose the AFD (Adaptive Fourier Decomposition) method which uses complex analysis to calculate the Hilbert transform. The AFD method adaptively approximates the analytical signal through the linear combination of the parameterized Szegö kernel to obtain original signals' Hilbert transform. Compared with the traditional method, the AFD method can give approximate analytical expressions and thus has a wider application. However, when using the principle of maximum selection to select parameters, the AFD method needs to exhaust all points of the unit disk which is time-consuming. Sparse Bayesian learning has been a hot spot in machine learning research in recent years. The Szegö kernel-based complex sparse Bayesian learning algorithm can provide a sparse rational approximation. A complex sparse Bayesian learning algorithm will be proposed based on the Szegö kernel to calculate the Hilbert transform. This method has the same advantages as the AFD method, and an iterative optimization can be performed without parameter control. The calculation speed of the proposed method is fast. Experimental results show that the proposed method is effective.
Key words: complex sparse Bayesian learning    Szegö kernel    Hilbert transform    analytic function    

Hilbert变换是信号分析与处理中的重要理论工具。由于复值Cauchy核在原点的奇异性,所以关于实部的Hilbert变换的计算是一项艰巨的任务[1-2]。计算Hilbert变换最常用的方法是快速傅里叶变换(Fast Fourier Transform, FFT)。最近,在文献[1-2]中给出了一些计算Hilbert变换的新方法。对于 ${L^2}$ 空间中的连续函数及 ${C^3} \cap {L^2}$ 的函数,学者们分别提出了两种样条的方法去计算Hilbert变换,从逼近精度上两种样条方法均优于FFT。随后,在文献[3]中首次使用复解析方法,提出了利用自适应傅里叶分解(Adaptive Fourier Decomposition, AFD)方法来计算Hilbert变换,充分利用解析空间中信号的实部和虚部是由Hilbert变换相关联的实值函数的特点来计算Hilbert变换。此方法适用于一般的非光滑和非连续信号,并提供显式的有理函数逼近,Takenaka-Malmquist (TM)系在有理函数逼近中起了重要的作用。虽然利用AFD计算Hilbert变换取得了较好效果,但其基于最大选择原理来自适应地选择最佳参数时,需要穷尽单位开圆盘中所有点,非常耗时。

近年来,稀疏贝叶斯学习方法成为机器学习中的一个研究热点,它是贝叶斯优化中十分重要的一类,是在贝叶斯理论的基础上发展而来,改进后的贝叶斯理论在文献[4]中也有提到。建立在稀疏贝叶斯学习基础上的学习机器称为相关向量机(Relevance Vector Machine, RVM),由Tipping[5]在2000年提出。它能够利用数据自身的特性对信号和图像实现最优的稀疏表示,已应用到智能检索、数据挖掘等领域。大多数基于稀疏贝叶斯的算法处理的是实值数据,然而在实际应用中还有许多复值信号需要处理。在文献[6]中利用同一信号实部与虚部具有相同的稀疏形式提出了复值贝叶斯压缩感知算法。本文将提出基于Szegö核的复值稀疏贝叶斯学习算法(Complex Relevance Vector Machine, CRVM),构造基于Szegö核的稀疏贝叶斯学习的Hilbert变换计算方法。基于Szegö核的稀疏贝叶斯学习算法,可以得到稀疏的有理逼近。事实上,由于TM系可通过对Szegö核做Gram-Schmidt正交化过程得到,因而基于Szegö核的稀疏贝叶斯学习算法也可以看成自适应傅里叶分解的一种,相比于AFD, 它的优势在于不需要穷尽单位开圆盘所有点去选择参数,并且最终所得到的模型是稀疏的。

对于给定能量有限信号 $s,$ 由著名的Plemelj公式[7]

$ \mathop {\lim }\limits_{y \to {0^ + }} Cs(x + {\rm{i}}y) = \frac{1}{{2{\rm{{\text{π}} i}}}}\int_{ - \infty }^\infty {\frac{{s(t)}}{{t - (x + {\rm{i}}y)}}{\rm{d}}t} =\frac{1}{2}s(t) + \frac{1}{2}{\rm{i}}Hs(t) $

其中, $Cs(z) = \dfrac{1}{{2{\rm{{\text{π}} i}}}}\displaystyle\int_{ - \infty }^\infty {\dfrac{{s(t)}}{{t - z}}} {\rm{d}}t$ $s$ 的柯西积分, $H$ 是Hilbert变换。记 ${s^ + }(t) = \dfrac{1}{2}s(t) + \dfrac{1}{2}{\rm{i}}Hs(t)$ ,则它是与 $s$ 相关的解析信号。

如果 $s(t)$ 是给定的实信号,那么有

$ \begin{split} & s = 2{\rm{Re}} {s^ + } \\& Hs = 2{\rm{Im}} {s^ + } \end{split}$

解析信号形成 ${L^2}(R)$ 的一个封闭的子空间,即Hardy空间 ${H^2}(R)$ 。相应的上半平面的解析函数构成了Hardy空间 ${H^2}({C^ + })$ 。Hardy空间 ${H^2}(R)$ ${H^2}({C^ + })$ 是等距同构的。

对于周期信号,也有类似的理论。对于一个有限能量的周期信号 $s$ ,有

${s^ + }({{\rm{e}}^{{\rm{i}}t}}) = \sum\limits_{k = - \infty }^\infty {{c_k}} {{\rm{e}}^{{\rm{i}}kt}},\sum\limits_{k = - \infty }^\infty {{{\left| {{c_k}} \right|}^2}} < \infty $

相应的解析信号可以表示为

${s^ + }(z) = \displaystyle \sum\limits_{k = 0}^\infty {{c_k}{z^k} = \dfrac{1}{{2{\rm{{\text{π}} i}}}}\displaystyle \int_0^{2{\rm{{\text{π}} }}} {\dfrac{{s({{\rm{e}}^{{\rm{i}}u}})}}{{{{\rm{e}}^{{\rm{i}}u}} - z}}{\rm{d}}{{\rm{e}}^{{\rm{i}}u}}} } $

其边界极限为

${s^ + }({{\rm{e}}^{{\rm{i}}t}}) = \frac{1}{2}s({{\rm{e}}^{{\rm{i}}t}}) + \frac{1}{2}\tilde Hs({{\rm{e}}^{{\rm{i}}t}}) + \frac{{{c_0}}}{2}$

其中 $\tilde H$ 是循环的Hilbert变换

$ \begin{split} & \tilde Hs({{\rm{e}}^{{\rm{i}}t}}) = \sum\limits_{k \ne 0} {( - {\rm{i}}{\rm{sgn}} (k)){c_k}{{\rm{e}}^{{\rm{i}}kt}}} =\\& \;\;\;\;\;\;\;\; \frac{1}{{2{\rm{{\text{π}} }}}}{\rm{p}}{\rm{.v}}{\rm{.}}\int_0^{2{\rm{{\text{π}} }}} {\cot \left(\frac{{t - s}}{2}\right)} s({{\rm{e}}^{{\rm{i}}s}}){\rm{d}}s \end{split}$

如果 $s$ 是实值的,则有以下关系

$ \left\{ \begin{aligned} & {s = 2{\rm{Re}}{s^ + } - {c_0}}\\& {Hs = 2{\rm{Im}}{s^ + }} \end{aligned}\right.$ (1)

有限能量的周期解析信号构成了 ${L^2}(\partial D)$ 的一个闭子空间,称其为单位圆周的Hardy空间,记为 ${H^2}(\partial D)$ ,单位圆周的Hardy空间 ${H^2}(\partial D)$ 与单位圆盘的Hardy空间是 ${H^2}(D)$ 同构的。本文将讨论单位圆周的情况。

1 基于Szegö核的CRVM

有理正交系,也称为TM系,在自适应傅里叶分解中起了十分重要的作用,在单位开圆盘空间 $D$ 中,TM系定义为

$ \begin{split} & {B_k}(z) = \frac{{\sqrt {1 - {{\left| {{a_k}} \right|}^2}} }}{{1 - {{\bar a}_k}z}}{\varPhi _{k - 1}}(z)\\& {\varPhi _k}(z) = \prod\limits_{j = 1}^k {\frac{{z - {a_j}}}{{1 - {{\bar a}_j}z}}} \end{split}$

其中 ${\varPhi _0}(z) = 1$ $\left| {{a_k}} \right| < 1$ $k \in {\bf{N}}$ $\left\{ {{B_k}} \right\}_{k = 1}^\infty $ 构成 ${H^2}(D)$ 的一个基当且仅当 $\displaystyle \sum\limits_{k = 1}^\infty {(1 - \left| {{a_k}} \right|) = \infty }$ 。当 $k = 1$ 时, ${B_1}(z) \!=\! \dfrac{{\sqrt {1 \!-\! {{\left| {{a_1}} \right|}^2}} }}{{1 \!-\! {{\bar a}_1}z}}$ 是单位化的带参数 ${a_1}$ 的Szegö核。众所周知的是TM系可以由一列带参数的Szegö核经过Gram-Schmidt正交化过程得来。

CRVM将基于贝叶斯公式进行推导,它研究的是将一个函数运用贝叶斯模型概率表示为多个基本函数的线性叠加。稀疏贝叶斯学习与支持向量机都具有共同的核函数形式,然而稀疏贝叶斯学习可通过迭代优化删除大量的训练样本和核函数且不需要进行参数控制,仅根据训练样本自动调整,它能够充分挖掘和利用数据的先验信息,不仅可以减少异常值的影响,而且可以“剔除”大量不必要的基。

$K(z,a) \!=\! \dfrac{{\sqrt {1 \!-\! |a{|^2}} }}{{1\! -\! \bar az}}$ ,对于一个实值信号 $s \in {L^2}(\partial D)$ 或者 $s \in {L^2}(R)$ ,本文对信号 $s$ 进行分解 $s = {s^ + } + {s^ - }$ 。对于给定的点集 $\{ ({z_n},{s^ + }({z_n}))\} _{n = 1}^N$ ,其中,N为点集个数,基于Szegö核的稀疏贝叶斯模型具有以下形式:

${s^ + }(z) = \sum\limits_{m = 1}^M {{w_m}} K(z,{a_m}) + v$

其向量形式为

${{{s}}^ + } = {{\varPhi w}} + {{v}}$ (2)

其中, ${{w}} = {[{w_1}, \cdots ,{w_M}]^{\rm{T}}} \in {C^{M \times 1}}$ 是待求的权值向量,M为待定模型的阶数, ${{{s}}^ + } = {[{s^ + }({z_1}), \cdots ,{s^ + }({z_N})]^{\rm{T}}} \in {C^{N \times {\rm{1}}}}$ 是给定的测量数据, ${{v}}$ 是未知的噪声矢量,

$ {{\varPhi}} = \left[ {\begin{array}{*{20}{c}} {K({z_1},{a_1})}& \cdots &{K({z_1},{a_M})} \\ \vdots & \; & \vdots \\ {K({z_N},{a_1})}& \cdots &{K({z_N},{a_M})} \end{array}} \right] \in {C^{N \times M}} $

${a_i} \in D,i = 1, \cdots ,M$ ,并且 ${{\varPhi}} $ 中各列均线性独立。此处本文假设测量所得的解析信号 ${s^ + }({z_i})$ 之间相互独立。将复数模型(2)经过实部、虚部分离变换为实数模型:

${{S}} = {{\varPsi W}} + {{V}}$ (3)

其中,

$ \begin{split} & {{S}} = \left[ {\begin{array}{*{20}{c}} {{\rm{Re}} ({{{s}}^ + })} \\ {{\rm{Im}} ({{{s}}^ + })} \end{array}} \right],{{W}} = \left[ {\begin{array}{*{20}{c}} {{\rm{Re}} ({{w}})} \\ {{\rm{Im}} ({{w}})} \end{array}} \right],{{V}}{\rm{ = }}\left[ {\begin{array}{*{20}{c}} {{\rm{Re}} ({{v}})} \\ {{\rm{Im}} ({{v}})} \end{array}} \right],\\& {{\varPsi}} = \left[ {\begin{array}{*{20}{c}} {{\rm{Re}} ({{\varPhi}} )}&{ - {\rm{Im}} ({{\varPhi}} )} \\ {{\rm{Im}} ({{\varPhi}} )}&{{\rm{Re}} ({{\varPhi}} )} \end{array}} \right] \end{split} $

假设 ${{W}}$ 服从参数化的均值为0的多维高斯分布 $p({{W}}|{{\alpha}} ) = N({{W}}|{{{\textit{0}}}_{2M}},{{{A}}^{ - 1}})$ ,其中 ${{A}} = {\rm{diag}}({{\alpha}} )$ ${{\alpha}} {\rm{ = (}}{\alpha _1}, \cdots , $ $ {\alpha _{2M}}{{\rm{)}}^{\rm{T}}}$ ${\alpha _i}$ 为未知超参数,控制着 ${{w}}$ 的稀疏程度。由于解析信号的实部和虚部具有相同的稀疏性,因此两者的控制参数也相同,表示为 ${\alpha _{l + M}} = {\alpha _l},l = 1, \cdots ,M$

假设噪声 ${{V}}$ 为高斯白噪声向量,即

$ {{V}} \sim N({{{\textit{0}}}_{2N}},{\delta ^2}{{{I}}_{2N}}) $

${{{\textit{0}}}_{2N}} = {(0, \cdots ,0)^{\rm{T}}} \in {R^{2N \times 1}}$ ${{{I}}_{2N}}$ $2N$ 阶单位矩阵, ${\delta ^2}$ 为噪声方差,根据以上假设,可以得到样本集的似然函数为

$p({{S}}|{{W}},{\delta ^2}) = N({{S}}|{{\varPsi W}},{\delta ^2}{I_{2N}}) $ (4)

为了推导的方便,此处省去了先验信息 ${{Z}} = $ $ \left[ {\begin{array}{*{20}{c}} {{\rm{Re}} ({{z}})} \\ {{\rm{Im}} ({{z}})} \end{array}} \right]$ 的书写。根据先验概率分布和似然分布,再利用贝叶斯公式计算权值的后验概率分布,即

$ p({{W}}|{{S}},{{\alpha}} ,{\delta ^{\rm{2}}}){\rm{ = }}\frac{{p({{S}}|{{W}},{\delta ^2})p({{W}}|{{\alpha}} )}}{{p({{S}}|{{\alpha}} ,{\delta ^2})}} $

其后验分布属于多变量的高斯分布,即

$p({{W}}|{{S}},{{\alpha}} ,{\delta ^{\rm{2}}})= N({{\mu}} ,{{\varSigma}} ) $ (5)

其中,

${{\varSigma}} = {({\delta ^{ - 2}}{{{\varPsi}} ^{\rm{T}}}{{\varPsi}} + {{A}})^{ - 1}} $ (6)
${{\mu}} {\rm{ = }}{\delta ^{{\rm{ - 2}}}}{{\varSigma}} {{{\varPsi}} ^{\rm{T}}}{{S}} $ (7)

训练目标值的似然分布通过对权值变量进行积分,即

$p({{S}}|{{\alpha}} ,{\delta ^2}){\rm{ = }}\int {p({{S}}|{{W}},{\delta ^2})} p({{W}}|{{\alpha}} ){\rm{d}}{{W}} $

从而求得超参数的边缘似然函数分布为

$p({{S}}|{{\alpha}} ,{\delta ^2}) = N({{{\textit{0}}}_{2N}},{\delta ^2}{{{I}}_{2N}} + {{\varPsi}} {{{A}}^{ - 1}}{{{\varPsi}} ^{\rm{T}}}) $

超参数 ${{\alpha}} $ 和噪声方差 ${{{\delta}} ^2}$ 的估计值分别为 ${{{\alpha }}_{{\rm{MP}}}}$ $\delta _{{\rm{MP}}}^2$ ,它们可以通过最大期望算法(Expectationmaximization Algorithm, EM算法)或第二类最大似然估计迭代得到。在迭代学习中,绝大部分的 ${\alpha _m}$ 将会趋于无穷大。当 ${\alpha _m}$ 趋于无穷大时,由式(7)可得到对应的 ${w_m}$ 趋于零。这样得到的模型是稀疏的,只有部分样本点在起作用,这些样本点称为相关向量。当超参数 ${{\alpha}} $ 和噪声方差 ${\delta ^2}$ 的值被估计出来后,权值 ${{W}}$ 的最大后验估计(Maximum a Posteriori, MAP)由后验概率分布式(5)的均值 ${{\mu}} $ 给出,进而可以得出 ${{w}}$ 的估计值。

在此框架下,

$ {s^ + }(z) \approx \sum\limits_{m = 1}^p {{w_m}} K(z,{z_m}),{\text{其中}},p \ll M。$

由式(1),对于单位圆的情况,可以得到实值信号 $s$ 的展开式为

$ s \approx 2{\rm{Re}} \left\{ {\sum\limits_{m = 1}^p {{w_m}K(z,{z_m})} } \right\} - {c_0}= \sum\limits_{m = 1}^p {{\rho _m}(t)\cos {\theta _m}\left( t \right)} - {c_0} $

其Hilbert变换为

$ Hs \approx 2{\rm{Im}} \left\{ {\sum\limits_{m = 1}^p {{w_m}K(z,{z_m})} } \right\}= \sum\limits_{m = 1}^p {{\rho _m}(t)\sin {\theta _m}\left( t \right)} $

其中, $ {\rho _m}(t){{\rm{e}}^{{\rm{i}}{\theta _m}(t)}} = 2{w_m}K(z,{z_m})$

2 收敛性

基于Szegö核的稀疏贝叶斯学习采用最大期望算法(EM算法)或第二类最大似然估计进行参数的迭代,由算法的性质可以保证每次迭代都降低成本函数,直到达到一个稳定点。因而基于Szegö核的稀疏贝叶斯算法是全局收敛的。

3 实验

下面将以一个例子来说明该方法的有效性。本文研究了文献[3]中的单位圆盘内的奇异函数:

${s^ + }(z) = \frac{{1 + 2{z^2}}}{{(z - 2)(z - 3)}} - \frac{1}{{(z + 2)(z + 3)}}{{\rm{e}}^{\tfrac{{z + 1}}{{z - 1}} + \tfrac{{z + {\rm{i}}}}{{z - {\rm{i}}}} + \tfrac{{z - {\rm{i}}}}{{z + {\rm{i}}}}}} $

显然 ${s^ + }$ 是一个解析信号或说是在 ${H^2}(D)$ 内解析的信号。由于解析信号的虚部是其实部的Hilbert变换,当有理解析函数序列逼近一个解析信号时,该序列中函数的实部和虚部分别逼近原始的实值信号及实值信号的Hilbert变换。图1图2图3分别给出了利用CRVM、FFT及AFD求Hilbert 变换的结果。

图 1 ${s^ + }(z)$ 的虚部与CRVM得到的Hilbert变换比较 Figure 1 Comparison of the imaginary part of ${s^ + }(z)$ and the Hilbert transform given by CRVM
图 2 ${s^ + }(z)$ 的虚部与FFT得到的Hilbert变换比较 Figure 2 Comparison of the imaginary part of ${s^ + }(z)$ and the Hilbert transform given by FFT

近似结果显示,所提方法CRVM效果优于AFD和FFT。在解析信号 ${s^ + }(z)$ 虚部的高振荡部分,CRVM和AFD逼近效果优于FFT, 而且CRVM和AFD能给出逼近函数的解析表达式,便于后续应用,而FFT只能得到一组近似的离散数据。比较AFD和CRVM,显然CRVM所用基的个数少于AFD,而且其逼近效果优于AFD。另外,CRVM由于不需要每一步在单位圆盘内根据最大选择原理选择参数,所用计算时间及复杂度远远小于AFD。

图 3 ${s^ + }(z)$ 的虚部与用AFD得到的Hilbert变换比较 Figure 3 Comparison of the imaginary part of ${s^ + }(z)$ and the Hilbert Transform given by AFD
4 总结

本文提出了基于复值稀疏贝叶斯定理近似计算Hilbert变换的方法。该方法是一种复解析方法,适用于所有能量有限的信号,能够在计算Hilbert变换的同时给出解析信号的稀疏有理逼近。与AFD方法相比,所提方法不需要穷尽单位圆盘的点去寻找参数,运算速度快且具有稀疏性。实验结果验证了所提方法的有效性。

参考文献
[1]
ZHOU C, YANG L, LIU Y, et al. A novel method for computing the Hilbert transform with Haar multiresolution approximation[J]. Journal of Computational and Applied Mathematics, 2009, 223: 585-579. DOI: 10.1016/j.cam.2008.02.006.
[2]
MICCHELLI C A, XU Y S, YU B. On computing with the Hilbert spline transform[J]. Advances in Computational Mathematics, 2013, 38(3): 623-646. DOI: 10.1007/s10444-011-9252-x.
[3]
MO Y, QIAN T, MAI W X, et al. The AFD methods to compute Hilbert transform[J]. Applied Mathematics Letters, 2015, 45: 18-24. DOI: 10.1016/j.aml.2014.12.017.
[4]
杨婷, 滕少华. 改进的贝叶斯分类方法在电信客户流失中的研究与应用[J]. 广东工业大学学报, 2015, 32(3): 67-72.
YANG T, TENG S H. Research and application of improved Bayes algorithm for the telecommunication customer churn[J]. Journal of Guangdong University of Technology, 2015, 32(3): 67-72. DOI: 10.3969/j.issn.1007-7162.2015.03.013.
[5]
TIPPING M E. The relevance vector machine[C]//Advances in Neural Information Processing Systems. MIT: DBLP, 2000: 652-658.
[6]
WU Q S, ZHANG Y D, AMIN M G, et al. Complex multitask Bayesian compressive sensing[C]//IEEE International Conference on Acoustics. Florence: SSP, 2014: 3375-3379.
[7]
GARNETT J B. Bounded analytic functions[M]. New York: Academic Press, 1981: 27-45.