一种M通道联合时间顶点非下采样滤波器组的优化设计方法与流程

本发明涉及图信号处理技术领域，具体是一种m通道联合时间顶点非下采样滤波器组的优化设计方法。

背景技术：

图信号处理是信号处理领域的一个重要部分，相比于传统信号处理，它能够处理非规则网络数据，例如传感器网络，社交网络等。作为传统信号处理的延伸，图信号处理定义了很多类似于传统信号处理的基本概念，包括图傅里叶变换，图滤波器，图滤波器组等，其中图滤波器组是图信号处理的一个重要工具，在图信号处理的多分辨分析中有着重要作用。虽然图信号处理近年来发展得很快，但是还有很多问题没有解决或完全解决。

在很多实际应用中，例如无线传感器网络，网络数据是随着时间变化而变化的，这些数据通过乘积图可以建模为时变图信号，为了处理这些时变数据，图信号处理顶点域通过乘积图模型延伸到联合时间顶点域，由此产生了时间顶点信号处理理论框架，通过笛卡尔乘积图可以定义联合傅里叶变换，其包括图傅里叶变换和离散傅里叶变换。基于联合傅里叶变换，联合时间顶点滤波器也随之提出，与之对应的应用也得到了很大发展，主要的应用有去噪、数据集分类、视频修复等。此外，对于时变图小波的研究也得以应用于捕抓联合时间顶点域的动态信息。由此可见，联合时间顶点信号处理的研究正得到快速发展，但是目前所研究的联合时间顶点滤波器和时变图小波很少能或者不能很好的对时变图信号进行多分辨分析，也就是很难对时变图信号进行更细致的频率划分，不能有效处理时变图信号，因此，具备良好的多分辨分析能力的m通道联合时间顶点滤波器组有待提出。

技术实现要素：

本发明所要解决的是现有联合时间顶点滤波器和时变图小波难以对时变图信号进行多分辨分析，不能有效处理时变图信号，提供一种m通道联合时间顶点非下采样滤波器组的优化设计方法。依据这种方法设计的滤波器组具备良好的频率划分特性和完全重构特性、能够对时变图信号进行多分辨分析，在对实测时变数据进行去噪时有良好的去噪性能。

实现本发明目的的技术方案是：

一种m通道联合时间顶点非下采样滤波器组的优化设计方法，与现有技术不同处在于，包括如下步骤：

1)设计m通道联合时间顶点非下采样分析滤波器组：联合时间顶点滤波器定义如下：

式中为普通图即顶点域上的归一化拉普拉斯矩阵，其中为单位阵，d为度矩阵，w为加权邻接矩阵，为有向循环图即时间域上的拉普拉斯矩阵，其中为单位阵，为邻接矩阵，kt和kg分别表示时域和顶点域上滤波器的长度，ak，l表示滤波器系数，表示克罗克内克积，在m通道联合时间顶点非下采样滤波器组中，假设x和分别表示输入和重构的时变图信号，{h0，…，hm-1}和{g0，…，gm-1}分别为分析滤波器组和综合滤波器组，zi＝hix为各个通道的子带信号，则滤波器组的输入输出关系如下：

其中t为传递函数，则滤波器组的完全重构条件可以表示为t＝i，即

其中i为单位矩阵，由此可得到在联合频域上的完全重构条件为：

hi(ω，λ)，gi(ω，λ)，i＝0，…，m-1分别是分析滤波器和综合滤波器的频谱核，m/2个分析滤波器的联合频谱核为：

式中ω，λ分别表示分析滤波器的时频域和图频域，联合频域范围为[-π，π]×[0，2]，m为通道数，m＞2且为偶数，kt，kg分别表示分析滤波器组中时域和顶点域上的滤波器长度，表示第i个分析滤波器的滤波器系数，另外m/2个分析滤波器频谱核可通过对m/2个分析滤波器频谱核在时频上调制得到，调制公式如下：

hi+m/2(ω，λ)＝hi(ω+π，λ)，i＝0，…，m/2-1(6)，

通过调制，可有效减少分析滤波器组的设计复杂度，不失一般性，为了避免时频域上的线性相位失真，根据滤波器线性相位的性质设定滤波器系数公式(5)可以更改为：

式中hi，r(ω，λ)＝(ci(ω，λ))^thi，为实数频率响应函数，t表示转置，其中：

上述中，假定kt为偶数，kt为奇数的情况也可以类似推出。类似于传统两维滤波器设计，联合时间顶点滤波器组的频率特性可由其通带平坦度和阻带衰减决定，因此，可以采用最小化关于滤波器系数向量hi的失真函数来控制它们，通带平坦度函数表达式如下：

其中di是与滤波器系数hi无关的常量，为分析滤波器的通带区域，以及

高阻带衰减可以通过控制阻带能量函数实现，其表达式如下：

其中为分析滤波器的阻带区域，综上，分析滤波器组的设计可以归结为一个无约束的优化问题：

其中hi为第i个分析滤波器系数向量，由公式(13)可得到m/2个分析滤波器系数向量，进而得到对应的分析滤波器频谱核，最后得到对应的分析滤波器；

2)依据步骤1)的数据，通过调制公式(6)对m/2个分析滤波器频谱核在时频域上调制可得到另外m/2个分析滤波器频谱核，然后得到对应的分析滤波器，最终得到m通道联合时间顶点非下采样分析滤波器组；

3)设计m通道联合时间顶点非下采样综合滤波器组：在非多项式的情况下，综合滤波器可以是一个具备稀疏特性的矩阵，也就是能够局部的重构时变图信号，因此设计非多项式形式的综合滤波器可归结为一个带约束的优化问题，如公式(14)所示：

其中约束条件为公式(3)中的完全重构条件，gi表示第i个综合滤波器，hi是由步骤1)、步骤2)求解得到的分析滤波器，i是一个单位矩阵，f表示矩阵的f范数；

4)将步骤3)的带约束优化问题采用拉格朗日乘子法求解，最终可得到综合滤波器gi的非多项式形式。

步骤1)中的无约束优化问题采用凸优化工具包cvx求解，进而得到对应的分析滤波器。

步骤4)所述的将步骤3)的带约束优化问题采用拉格朗日乘子法求解为：

拉格朗日函数的表达式为公式(15)：

其中φ是拉格朗日乘子矩阵，通过求解公式(15)的拉格朗日函数偏导数为0，则有公式(16)：

通过计算，可以得到带约束优化问题中综合滤波器组的解为公式(17)：

其中

本技术方案与现有技术相比，利用优化方法设计m通道联合时间顶点非下采样分析滤波器组，求解出具备良好频率划分特性的分析滤波器组，然后将综合滤波器组的设计问题归结为一个带约束的优化问题，并通过拉格朗日乘子法进行求解。

依据这种方法设计的滤波器组具备良好的频率划分特性和完全重构特性、能够对时变图信号进行多分辨分析，在对实测时变数据进行去噪时有良好的去噪性能。

附图说明

图1为实施例中m通道联合时间顶点非下采样滤波器组的基本结构示意图；

图2为实施例中四通道分析滤波器组在滤波器长度kt＝kg＝2条件下的频率响应示意图；

图3为实施例中四通道分析滤波器组在滤波器长度kt＝kg＝8条件下的频率响应示意图；

图4为实施例中六通道分析滤波器组在滤波器长度kt＝kg＝2条件下的频率响应示意图；

图5为实施例中六通道分析滤波器组在滤波器长度kt＝kg＝8条件下的频率响应示意图。

具体实施方式

下面结合附图和实施例对本发明的内容作进一步的阐述，但不是对本发明的限定。

实施例：

一种m通道联合时间顶点非下采样滤波器组的优化设计方法，包括如下步骤：

1)设计m通道联合时间顶点非下采样分析滤波器组：联合时间顶点滤波器定义如下：

式中为普通图即顶点域上的归一化拉普拉斯矩阵，其中为单位阵，d为度矩阵，w为加权邻接矩阵，为有向循环图即时间域上的拉普拉斯矩阵，其中为单位阵，为邻接矩阵，kt和kg分别表示时域和顶点域上滤波器的长度，ak，l表示滤波器系数，表示克罗克内克积，如图1所示的是m通道联合时间顶点非下采样滤波器组的结构，其中x和分别表示输入和重构的时变图信号，{h0，…，hm-1}和{g0，…，gm-1}分别为分析滤波器组和综合滤波器组，zi＝hix为各个通道的子带信号，滤波器组的输入输出关系如下：

其中t为传递函数，则滤波器组的完全重构条件可以表示为t＝i，即

其中i为单位矩阵，由此可得到在联合频域上的完全重构条件为：

hi(ω，λ)，gi(ω，λ)，i＝0，…，m-1分别是分析滤波器和综合滤波器的频谱核，m/2个分析滤波器的联合频谱核为：

其中ω，λ分别表示分析滤波器的时频域和图频域，联合频域范围为[-π，π]×[0，2]，m为通道数，m＞2且为偶数，kt，kg分别表示分析滤波器组中时域和顶点域上的滤波器长度，表示第i个分析滤波器的滤波器系数，另外m/2个分析滤波器频谱核可通过对m/2个分析滤波器频谱核在时频上调制得到，调制公式如下：

hi+m/2(ω，λ)＝hi(ω+π，λ)，i＝0，…，m/2-1(6)，

通过调制，可有效减少分析滤波器组的设计复杂度，不失一般性，为了避免时频域上的线性相位失真，根据滤波器线性相位的性质设定滤波器系数公式(5)可以更改为：

式中hi，r(ω，λ)＝(ci(ω，λ))^thi，为实数频率响应函数，t表示转置，其中：

上述中，假定kt为偶数，kt为奇数的情况也可以类似推出。类似于传统两维滤波器设计，联合时间顶点滤波器组的频率特性可由其通带平坦度和阻带衰减决定，因此，可以采用最小化关于滤波器系数向量hi的失真函数来控制它们，通带平坦度函数表达式如下：

其中di是与滤波器系数hi无关的常量，为分析滤波器的通带区域，以及

高阻带衰减可以通过控制阻带能量函数实现，其表达式如下：

其中为分析滤波器的阻带区域，综上，分析滤波器组的设计可以归结为一个无约束的优化问题：

其中hi为第i个分析滤波器系数向量，由公式(13)可得到m/2个分析滤波器系数向量，进而得到对应的分析滤波器频谱核，最后得到对应的分析滤波器；

2)依据步骤1)的数据，通过调制公式(6)对m/2个分析滤波器频谱核在时频域上调制可得到另外m/2个分析滤波器频谱核，然后得到对应的分析滤波器，最终得到m通道联合时间顶点非下采样分析滤波器组；

3)设计m通道联合时间顶点非下采样综合滤波器组：在非多项式的情况下，综合滤波器可以是一个具备稀疏特性的矩阵，也就是能够局部的重构时变图信号，因此设计非多项式形式的综合滤波器可归结为一个带约束的优化问题，如公式(14)所示：

其中约束条件为公式(3)中的完全重构条件，gi表示第i个综合滤波器，hi是由步骤1)、步骤2)求解得到的分析滤波器，i是一个单位矩阵，f表示矩阵的f范数；

4)将步骤3)的带约束优化问题采用拉格朗日乘子法求解，最终可得到综合滤波器gi的非多项式形式。

步骤1)中的无约束优化问题采用凸优化工具包cvx求解，进而得到对应的分析滤波器。

步骤4)所述的将步骤3)的带约束优化问题采用拉格朗日乘子法求解为：

拉格朗日函数的表达式为公式(15)：

其中φ是拉格朗日乘子矩阵，通过求解公式(15)的拉格朗日函数偏导数为0，则有公式(16)：

通过计算，可以得到带约束优化问题中综合滤波器组的解为公式(17)：

其中

具体地：

仿真案例1：

设计m通道联合时间顶点非下采样滤波器组：本次仿真以四通道和六通道为例，分析滤波器长度分别设置为kt＝kg＝2和kt＝kg＝8，四通道中的第一个分析滤波器频谱核h0(ω，λ)的通带区域和阻带区域分别设置为[0，1.17]×[0，0.6]，[-π，π]×[1.14，2]∪[1.97，π]×[0，0.6]∪[-π，-1.97]×[0，0.6]，另外三个频谱核{hi(ω，λ)}i＝1，2，3通过对h0(ω，λ)调制可得，即h1(ω，λ)＝h0(ω，2-λ)，h2(ω，λ)＝h0(π+ω，λ)，h3(ω，λ)＝h0(π+ω，2-λ)，六通道的通带和阻带区域设置如表1所示。

表1

另外三个频谱核{hi(ω，λ)}i＝3，4，5通过对{hi(ω，λ)}i＝0，1，2在时频上调制可得，图2和图3分别表示的是四通道分析滤波器组在滤波器长度kt＝kg＝2和8的条件下的频率响应示意图，图4和图5分别表示的是六通道在kt＝kg＝2和8条件下的频率响应示意图，由此可以看出，采用本例方法设计的m通道联合时间顶点非下采样滤波器组具有良好的频率划分特性，能够对时变图信号进行更细致的频率划分，具备多分辨分析特性。综合部分，也是滤波器组的重构部分，输入时变信号为实测一年美国温度数据中的前十天温度数据，采用最近距离算法构造包含150个观测点的温度图结构，通过笛卡尔乘积对温度图和包含十个时间节点的有向循环图建模为乘积图，最后利用公式(17)求解出的综合滤波器组，对信号进行重构。表2给出了本例方法与现有方法1(可分时变非下采样图滤波器组)设计的滤波器组的重构性能对比，由表2可以看出，本例方法设计的m通道联合时间顶点非下采样滤波器组具有更好的重构性能。

表2

仿真案例2：

本次仿真案例2中去噪实验采用硬阈值法，以四通道为例对输入含噪信号进行去噪，滤波器组为仿真案例1中的四通道滤波器组，分析滤波器长度组设为kt＝kg＝2，输入原始时变图信号为2013年1月到2014年5月的十七个月的海平面温度数据，包含100个节点的海平面温度图结构由最近距离算法构造，输入的随机噪声服从均匀分布。对于四通道，另外三个通道的硬阈值{τ}0,1,2分别设置为1.3σ，2.1σ，3σ，其中σ为噪声标准差。表3给出了本例方法与现有方法1(可分时变非下采样图滤波器组)、现有方法2(非时变非下采样图滤波器组)的去噪性能的对比，其中，对于非时变非下采样图滤波器组，对含噪时变图信号处理的方式就是每一次输入滤波器组的信号为每一个时刻的含噪图信号，最后输出的每一个时刻的去噪信号组合起来得到去噪后的时变图信号。仿真结果表明，采用本例方法的去噪性能和采用现有方法1相当，基本比采用现有方法2的去噪效果要好，因此采用本例方法所设计的m通道联合时间顶点非下采样滤波器组能对实测时变信号进行良好的去噪。

表3