基于迭代计算的时变可分非下采样图滤波器组的设计方法与流程

本发明涉及图信号处理领域中的时变图滤波器组技术领域，具体涉及一种基于迭代计算的时变可分非下采样图滤波器组的设计方法。

背景技术：

随着目前大数据时代的到来，生活中记录的数据呈现出海量化、多样化和非规则化等特点。其中复杂的非规则的分布式网络结构，例如，疾病传播网络，谣言传播网络，全球海平面温度网络，传感器温度网络等，广泛存在于现实当中。对于复杂的非规则的网络数据，传统信号处理难以对其分析处理，因此延伸出了图信号处理。

在图信号处理领域中，图滤波器、图小波和图滤波器组成为分析图信号的有力工具。图小波和图滤波器中均含有图信号的下采样运算，对于一般的大规模和复杂图对进行采样运算时，会存在采样不精确的问题。简化的非下采样图滤波器组结构，避免了图滤波器中的采样运算，并采用分布式算法进行图滤波器组的实现。然而，目前很多图滤波器及图滤波器组是基于图的顶点域进行分析处理，而现实生活中的网络节点信号多数是随着时间变化而变化时变图信号，因此图信号的时间关联性分析是不可缺少的。

近年来，联合图傅里叶变换的提出，使其成为图信号处理最基本的工具。通过建立联合时间顶点频域的概念，联合图傅里叶变换能够在联合时间顶点频域中分析时变图信号。与联合图傅里叶变换一起，联合时间顶点图滤波器作为不可分割的部分随之出现。目前，对于时变图滤波器的研究较多，其主要应用于时变图信号的去噪，时变图信号的重构，视频的修复等，但时变图滤波器不具备多分辨分析特性，且不能稀疏表示时变图信号。而时变滤波器组却可以稀疏表示图信号，从而为大规模高维时变图信号的处理提供了有效的工具，但目前对于时变图滤波器组的研究相对较少，进一步的深入研究有待提出。

技术实现要素：

本发明所要解决的是目前图滤波器组难以处理时变图信号的的问题，提供一种基于迭代计算的时变可分非下采样图滤波器组的设计方法。

为解决上述问题，本发明是通过以下技术方案实现的：

基于迭代计算的时变可分非下采样图滤波器组的设计方法，包括如下步骤：

步骤1、设计时变非下采样图滤波器的分析滤波器组为：

其中，hll表示时域和图域都为低频的时变图滤波器，hlh表示时域低频图域高频的时变图滤波器，hhl表示时域高频图域低频的时变图滤波器，hhh表示时域和图域都为高频的时变图滤波器，it表示大小为t×t的单位矩阵，t为时域的时间节点数，ig表示大小为n×n的单位矩阵，n为图域的图节点数，lt表示有向循环图的拉普拉斯矩阵，lg表示无向图的归一化拉普拉斯矩阵，kt表示时域滤波器的长度，kg表示图域滤波器的长度；

步骤2、计算全局矩阵h和全局固定系数向量b：

其中，x为输入时变图信号；

步骤3、对于每个节点k，其中k＝1,2,…,n×t，分别执行以下迭代操作：

步骤3.1、初始化，即令迭代次数n＝1，当前输出信号和当前节点k的2r阶邻居的局部迭代系数向量

步骤3.2、计算当前节点k的2r阶邻居的局部融合系数向量

其中，表示从全局矩阵h中取出节点k的2r邻居对应的值组成的局部矩阵，表示从全局固定系数向量b中取出节点k的2r邻居对应的值组成的局部固定系数向量，表示当前节点k的2r阶邻居的局部迭代系数向量；

步骤3.3、从当前节点k的2r阶邻居的局部融合系数向量中取出节点k的r阶邻居内对应的值，组成当前节点k的r阶邻居的局部融合系数向量

步骤3.4、把当前节点k的r阶邻居的局部融合系数向量中的所有元素相加，再用相加得到的和除以元素个数得到一个平均值，并将该平均值赋给节点k的当前最小二乘解系数向量v⁽ⁿ⁾中的第k个元素；

步骤3.5、将当前最小二乘解系数向量v⁽ⁿ⁾与当前输出信号相加后，去更新当前输出信号

步骤3.6、将全局矩阵h与当前输出信号相乘后，去更新当前全局迭代系数向量θ⁽ⁿ⁾；

步骤3.7、从当前全局迭代系数向量θ⁽ⁿ⁾中取出节点k的2r阶邻居内对应的值，去更新当前节点k的2r阶邻居的局部迭代系数向量

步骤3.8、如果当前最小二乘解系数向量v⁽ⁿ⁾达到迭代终止条件，即||v⁽ⁿ⁾||∞≤ε，其中ε为设定的迭代终止阈值，则迭代终止，输出当前输出信号否则，令迭代次数n加1，返回步骤3.2继续迭代。

上述步骤中，时域滤波器的长度kt≥1，图域滤波器的长度kg≥1。

上述步骤中，节点邻居的阶数r≥1。

与现有技术相比，本发明首先基于两维可分图滤波器的性质，设计出具备频率响应的分析滤波器组。然后把综合滤波器组输出信号的重构问题转换为一个全局最小二乘问题，再将全局最小二乘问题转换为局部最小二乘问题，并采用迭代的方式进行求解。仿真结果表明，本发明迭代计算方法具有较低的迭代次数，所设计出的时变非下采样图滤波器组具备完全重构特性和具有更好的去噪性能，其分析滤波器组具备频率响应。

附图说明

图1为时变非下采样图滤波器组的基本结构。

图2为局部操作示意图。

图3为时变非下采样图滤波器组的分析滤波器组的频率响应，其中(a)时变图滤波器hll(lt,lg)的频率响应，(b)时变图滤波器hlh(lt,lg)的频率响应，(c)时变图滤波器hhl(lt,lg)的频率响应，(d)时变图滤波器hhh(lt,lg)的频率响应。

具体实施方式

为了便于理解，下面对本发明基于迭代计算的时变可分非下采样图滤波器组的设计方法，进行进一步详细说明。

根据图1表示一个时变非下采样图滤波器组的结构，其中时变图滤波器hll(lt,lg)表示在时间频域和图频域上都为低通，hlh(lt,lg)表示在时间频域上为低通，图频域上为高通，hhl(lt,lg)表示在时间频域上为高通，图频域上为低通，hhh(lt,lg)表示在时间频域和图频域上都为高通，并与gll(lt,lg)，glh(lt,lg)，ghl(lt,lg)，ghh(lt,lg)分别组成分析图滤波器组和综合图滤波器组。x为图滤波器组的输入信号，输出信号为xll,xlh,xhl,xhh分别表示每个通道的子带系数。其中图域的归一化拉普拉斯矩阵为lg＝ig-dg^-1/2wgdg^1/2，ig、dg和wg分别表示阶数为n的单位阵、图的度矩阵和图的临界矩阵。时间域有向循环图的拉普拉斯矩阵lt＝it-at，it和at分别表示阶数为t的单位阵，有向循环图的邻接矩阵。

时变非下采样图滤波器组的输入输出关系为

其中t称为图滤波器组的传递函数。时变非下采样图滤波器组的完全重构条件为

式中i表示大小为nt×nt的单位阵。联合时间图频率响应的完全重构条件

从时变非下采样图滤波器组的结构上能观察到本文设计的时变非下采样图滤波器组对于时变信号能够从不同的频带对图信号进行划分，具备多分辨分析特性。

第一步：首先考虑分析滤波器组的设计问题，对于整个滤波器组来说，分析滤波器组的设计主要是为了对输入信号进行频率的划分。根据两维可分图滤波器的性质，两维可分时变图滤波器可以分别从图频域和时频域对期望频率响应进行设计。因此，定义时变非下采样图滤波器组的分析滤波器组为

运算表示克罗内克积，kt,kg分别表示滤波器的长度。通过图傅里叶变换和z变换可以得到联合频率响应

第二步：根据给定分析图滤波器组，把综合滤波器组的设计问题归结于一个带约束优化问题

||·||f表示矩阵的frobenius范数，通过求解带约束优化问题，得到

式中“t”代表转置。

第三步：根据求解得到的综合滤波器组中含有大规模矩阵的求逆，而时变图信号属于大规模高维数据，大规模的求逆运算不利于时变图信号处理的。为避免直接求解gll,glh,ghl,ghh，采用迭代计算的重构方式去实现综合滤波器的重构。

根据时变非下采样图滤波器组的输入输出关系，其重构信号可以通过每个通道的子带系数计算得到

其公式(17)的输出信号是下面全局最小二乘问题的唯一解

为了以迭代计算的方法求解全局问题(18)，从而引出其局部最小二乘问题

其中表示局部操作。如图2所示，t1时刻当r＝1时，表示取图上节点k的一阶邻居及前一时刻对应位置的值的局部操作。当r＝1时，表示取图上节点k的二阶邻居及节点k前两时刻和一阶邻居的前一时刻对应位置的值的局部操作。则局部最小二乘的解为

式中

第四步：采用迭代计算方法对滤波器组进行时变图信号的重构，其迭代计算方法如下表1：

表1

下面通过具体仿真实例，对本发明的性能进行说明。

仿真实例1：

本发明首先设计分析滤波器组，分析滤波器组由公式(4)-(7)构造产生，当图滤波器的长度kt＝1,kg＝1时，其联合频率响应如图3所示，以实测的一年的美国温度数据作为输入信号，采用最近距离的方式构造了温度图结构，邻接矩阵ag设为ag(i,j)＝1/(disti,j)²，如果节点i和节点j不是同一节点且有一条边相连，否则ag(i,j)＝0，disti,j表示节点i和节点j间的距离，再利用迭代计算方法重构出输出信号。表2给出了本发明设计时变非下采样图滤波器组的迭代计算方法的迭代次数和重构信噪比。实验结果表明设计所得的图滤波器组都具备完全重构特性，迭代计算方法具有较低的迭代次数。

表2

仿真实例2：

当kt＝1,kg＝1,r＝1,ε＝10-⁴时，使用时变非下采样滤波器组对时变图信号进行去噪实验，噪声为均匀随机噪声。时变图信号选取2014年1月到5月的海平面温度网络信号为例。并采用硬阈值法进行去噪实验，其中lh,hl,hh三个通道的硬阈值分别选取为τ＝0.5σ,τ＝σ,τ＝3σ，σ为噪声标准差。其中现有方法1(临界采样双正交图滤波器组)和现有方法2(基于谱域三边滤波器的图信号去噪)选取2014年3月海平面温度网络信号进行去噪。表3给出了噪声标准差σ取不同值时，本发明与已有图滤波器及图滤波器组去噪信噪比对比。对比实验仿真结果表明，与现有技术相比，本发明设计的时变图滤波器组对于实际图信号有着更好的去噪性能。

表3

需要说明的是，尽管以上本发明所述的实施例是说明性的，但这并非是对本发明的限制，因此本发明并不局限于上述具体实施方式中。在不脱离本发明原理的情况下，凡是本领域技术人员在本发明的启示下获得的其它实施方式，均视为在本发明的保护之内。