卢军志,蒋俊正

(桂林电子科技大学 信息与通信学院,广西 桂林 541004)

近年来,随着科技和社会的不断发展,信息数据呈现出海量化、多样化、非规则化等特点。传统的信号处理很难对非规则信号进行分析处理,即使可以,往往也会忽略网络本身的一些结构特性。为了更好地处理非规则信号,研究者结合图论等知识提出了图信号处理,它能结合数据本身的结构特性对数据进行分析,在处理无线传感器网络信号等非规则信号都取得了不错的效果[1-3]。作为传统信号处理的一种拓展,图信号处理同样也存在傅里叶变换、滤波器及滤波器组等概念,其中图滤波器组因其具备多分辨分析能力并能对信号进行多速率信号处理而备受研究者的关注。

关于图滤波器组的研究是先从两通道开始的。2006年一种可逆图滤波器组的设计方法[4]被提出,并用于对传感器网络上信号进行多分辨分析和处理。Narang等[5]在2012年研究了两通道临界采样正交图滤波器组的设计,通过用切比雪夫逼近的方法设计出低通分析滤波器,根据正交特性求解出其余滤波器。该方法缺点在于不能满足紧支撑特性。因此,2013年,Narang等[6]对正交性进行了松弛,变成双正交,使整个滤波器组系统满足完全重构特性和紧支撑特性,但其方法是基于对半带滤波器的分解,因此无法保证良好的频谱选择性。随后,研究者们采用不同的方法[7-8]设计滤波器组,以期达到权衡频谱选择性和完全重构特性的目的。渐渐地研究者的关注点从两通道转移到M通道上。2014年Tanaka[9]提出了M通道过采样图滤波器组的设计方法,通过构造过拉采样普拉斯矩阵实现过采样,最后实现完全重构;此后,一些基于循环图的过采样图滤波器组的设计方法[10-13]被提出。2019年,一种M通道临界采样图滤波器组的设计方法[14]被提出,通过用插值操作取代重构部分来实现重构。在目前的设计中,无论是两通道还是M通道,分析滤波器组的频带划分都是均匀的,但由于图信号的频率是离散的,频率的分布可能会出现分布不均匀且部分集中的情况,这时非均匀的频率划分会更有意义。

在图滤波器组的研究中,分布式实现[15-17]一直都是一个很重要的课题,特别是将其应用到大规模的网络数据处理时,分布式实现能够节省计算成本。目前大部分的图滤波器组都采用切比雪夫多项式逼近的方法[15]来设计,这种方法通过增大阶数获得更好的频率选择特性,在阶数小的情况下,频率选择特性较差,但具备较好的稀疏特性,能进行分布式实现。然而,这种分布式仅限于多项式形式下的滤波器,并不适用于非多项式形式的滤波器,如节点变滤波器组[18]等。2019年,Jiang等[16]提出了一种非下采样两通道图滤波器组的分布式实现方法,滤波器的形式不再局限于多项式形式,也可以是非多项式形式,更具普适性。该方法将重构问题归结为一个最小二乘问题,利用二阶信息及近似牛顿法求解,其中用矩阵近似理论和图的结构特性对海森矩阵的逆进行近似,从而实现分布式。

针对频率的分布特性,提出了一种非均匀图滤波器组的设计方法,并将其进行了分布式实现。根据图频率分布不均匀的情况,通过优化方法设计出具备良好频率特性的多项式形式的分析滤波器,用非多项式形式的滤波器对其进行近似。在已知子带信号和分析滤波器的前提下,将图滤波器组的重构问题归结成一个最小二乘问题。为了避免直接求解矩阵的逆,采用预处理梯度下降法对优化问题进行迭代求解。在迭代过程中,未利用目标函数的二阶信息,方法更简单。

1 基础知识

图信号处理中,无向图G可由三维集合

表示,其中:V为无向图的所有顶点的集合,N=|V|;E为无向图的所有边的集合;W={wij},i,j∈V为无向图的权重矩阵。根据权重矩阵可得无向图的度矩阵D,其对角线上元素

进而可定义拉普拉斯矩阵L=D-W及归一化拉普拉斯矩阵Ln=I-D-1/2WD-1/2。在无向图中,Ln是对称的,对其进行特征分解,有Ln=UΛUT,U为特征向量组成的矩阵,Λ=diag([λ1,λ2,…,λN])为对角矩阵,λ1,λ2,…,λN为Ln的特征值,Ln的特征值范围为[0,2]。无向图上信号x定义为每个节点上的信号值的集合,x=[x1,x2,…,xN]T,如图1所示。x的图傅里叶变换定义为x^=UTx,其图傅里叶逆变换为x=Ux^。根据图傅里叶变换的定义及图2可知,图频域上的频率值即为图拉普拉斯矩阵的特征值,是离散的,且λ=0.4和λ=0.8处具有明显的频率间隔。因此,需要根据频率的聚集程度对频谱进行划分,进而设计出非均匀图滤波器组。

图1 节点数为50的随机图

图2 图1中信号的图频域表示

图滤波的定义是,图信号x经过图滤波器后得到滤波后的输出f,其表达式为

其中φ(k)为滤波器第k阶的系数

根据图2的频率分布特性,主要考虑的是三通道非均匀图滤波器组。由于目前图上的采样操作并不唯一,且没有针对奇数通道的采样设计,因此不讨论图滤波器组上的采样操作。三通道非均匀非下采样图滤波器组结构框图如图3所示,x、分别为图滤波器组的输入和输出,{H0,H1,H2}为分析滤波器组,{z0,z1,z2}为滤波器的子带系数,zm=Hmx,m∈[0,1,2],{G0,G1,G2}为重构滤波器组。该框架的输入输出关系可表示为

图3 三通道非下采样图滤波器结构组框图

为满足其完全重构特性,整个系统的传递函数T应满足:

其中I为单位矩阵。

目前大部分的滤波器组设计均采用切比雪夫近似的方法同时设计分析滤波器和综合滤波器。这些方法都是设计多项式形式的,设计自由度不高,且在良好的频率选择特性与好的局部特性之间不能很好地权衡。阶数越高,频率选择性越好,但局部特性越差,越不能分布式实现。而非多项式形式的滤波器在提高设计自由度的前提下尽可能使得频率选择性更好,同时也拥有良好的局部特性。因此,选用非多项式形式中的NV图滤波器作为分析滤波器,其设计步骤为:

1)通过约束滤波器的通带平坦性和阻带衰减得到频率选择性好的多项式形式的分析滤波器[19],其优化问题可表示为

φ=[φ(0),φ(1),…,φ(K)]T,b(λ)=[1,λ,…,λK]T。式(3)的前项是对滤波器通带平坦性的约束,a、b为滤波器的通带截止频率;后项是对阻带衰减的约束,c、d为滤波器的阻带截止频率,α为可调整参数。确定好滤波器的通带、阻带截止频率后,就能通过求解式(3)的优化问题得到图滤波器的系数。

2)将设计好的多项式滤波器当作理想滤波器,再用优化的方法设计出阶数较低且具有良好频率选择性的NV 图滤波器[18]。在该方法中,通过最小化Frobenius范数

得到最优的滤波器系数,Ck为系数矩阵。NV 滤波器的优点是获得更小阶数、更好频率选择特性的滤波器,同时所需设计的系数变多,设计自由度变高,设计难度也变大。

图滤波器组的重构部分并不需要像分析滤波器一样具有良好的频谱选择特性,只需满足完全重构特性即可。因此,重构部分的实现也不一定需要通过重构滤波器来实现,且目前并无太多关于非多项式形式图滤波器组的设计方法。在给定了分析滤波器组{H0,H1,H2}及子带信号{z0,z1,z2}的前提下,将重构问题归结为一个优化问题:

此最小二乘问题的最优解为

最优解的求解需要进行矩阵求逆运算,而矩阵求逆的复杂度是O(N3)。若矩阵维度很大,且求逆矩阵无特殊性质的情况下,求逆的计算代价很高,且可能存在矩阵病态问题,因此采用一种预处理分布式梯度下降法来求解这个逆滤波问题。

2 设计方法

逆滤波问题是一个线性问题,可采用梯度法进行迭代求解,并采用预处理技术加速迭代方法的收敛速度。预处理技术的基本思想是在原本的线性问题中加入预处理算子,使之能等价转换为新的线性问题。这么做的好处是能够获得更快的收敛速度,同时在一些条件数很差的系统中,还能起到改善条件数的作用。本研究选取的预处理算子P定义为

其中:B(i,ω(H))={j∈V,ρ(i,j)≤ω(H)}为节点i的ω(H)阶以内所有邻居节点的集合;ω(H)为矩阵H的测地距离,ω(H)≥0;ρ(i,j)为从节点i到达节点j所需的最短路径。ω(H)与滤波器的阶数有关,是阶数的2倍。设置=P-1/2HP-1/2,类似归一化拉普拉斯矩阵的操作,对矩阵H进行归一化处理。本研究所选的预处理形式是右预处理形式:

将其与梯度法结合后,有

再经过变换后,可得

在整个方法中,避免了复杂矩阵的求逆问题。矩阵的测地距离ω(H)是滤波器阶数的2倍,而由于所设计的滤波器阶数小,故ω(H)很小。因此,在每次迭代过程中,节点只与其少量的邻居进行信息交换,整个过程可分布式实现。该方法对每个节点i上信号进行分布式重构的算法流程如下:

将所有的(M)(i)进行组合,可得到所有节点最终的重构信号值。由算法流程可看出,在整个迭代过程中,计算复杂度最高的操作就是矩阵与向量的相乘,此计算复杂度比矩阵求逆的计算复杂度低了一个数量级。因此,该方法的计算复杂度小于直接求解的复杂度,且在图规模很大的情况,这种优势会更明显。

3 仿真实验

首先给出关于三通道非均匀图滤波器组的分析滤波器组的设计及分析,然后给出该方法与现有的分布式迭代方法的对比。仿真所用的图信号如图1所示,所有仿真均在相同环境下进行。

实验1根据图2的频率分布特性,通过优化非均匀分析滤波器组的系数,设计出如图4所示的非均匀频谱划分的分析滤波器组。h0、h1、h2分别表示低通、带通和高通非均匀分析滤波器的频率响应,3个滤波器的阶数都是16。从图4可看出,低通、高通分析滤波器的通带都比较平坦,频率选择特性都很好,高阶多项式滤波器虽然频率选择特性好,但局部特性很差,并不适合应用于分布式方法中,因此采用NV滤波器[18]对设计好的多项式滤波器进行近似。在实验中,设计了阶数依次为1、2和3的NV 滤波器,分别用K1、K2、K3表示。但由于NV 滤波器并不存在频率响应的说法,为了说明所得到的NV滤波器的滤波性能与图4的性能相近,将高频子带信号进行比较来说明近似效果,观察NV滤波器经过高频滤波后的得到的高频子带信号的频谱图是否与多项式形式滤波器一致,从而判断NV形式的滤波器是否达到了较好的近似效果。图5为输入信号经过不同形式的高通滤波器后得到的子带信号频域表示。从图5可看出,NV形式滤波器滤波得到的子带信号频谱尽管在少许频率分量上会有幅度上的小差别,但整体上滤波效果与多项式形式的效果接近。因此,采用NV形式滤波器来近似高阶多项式形式滤波器的结果是可行的,能取得同样的效果。

图4 设计得到的三通道非均匀分析滤波器组的幅度响应

图5 相同信号经过不同滤波器得到的高频子带信号的频域表示

实验2对图信号进行分布式重构的实验仿真,总的迭代次数为200,分析滤波器基于例1的多项式分析滤波器,是通过优化设计出的非多项式形式的分析滤波器。当滤波器分别为K1、K2和K3时,不同迭代方法的收敛曲线分别如图6、7和图8所示。相对误差为‖(m)-x‖2/‖x‖2。在仿真中发现,当步长取0.9时,分布式梯度法[20]能达到最优的效果。因此,在不同的阶数下,该方法均采用相同的步长。根据图6～8可知,本方法的迭代收敛次数比分布式梯度法[20]要少很多,与分布式牛顿法[17]相差不大,在不使用二阶信息的前提下,取得了比一般分布式梯度法更快的收敛速度,表明预处理技术能够有效加快迭代收敛的速度;在阶数不断增大时,所有算法的收敛次数都在减少,原因在于阶数增大,节点与周围邻居节点的通信次数变多,获得的有用信息增多,通过不断的信息交换,就能更快地达到收敛。表1为迭代完成后本方法在不同阶数条件下的重构误差。从表1可看出,本方法能够进行分布式重构,且重构误差比分布式梯度法[20]小,与分布式牛顿法[17]差距不大。

表1 迭代完成后本算法在不同滤波器阶数下的重构误差

图6 采用K1时,不同方法的收敛曲线

图7 采用K2时,不同方法的收敛曲线

图8 采用K3时,不同方法的收敛曲线

4 结束语

针对频率分布不均匀的情况,提出了一种非均匀非下采样图滤波器组的设计方法,并进行了分布式实现。根据图频率的分布特性,设计了具备良好频率划分特性及稀疏特性的低阶非多项式分析滤波器。在已知分析滤波器和子带信号的前提下,将图滤波器组的重构问题归结为一个最小二乘问题,并采用预处理梯度法对其进行迭代求解。仿真结果表明,加入预处理算子后,能够加快收敛速度,且比利用二阶信息的迭代方法拥有更低的计算成本。