基于凸优化的M通道过采样图滤波器组的优化设计方法

摘要

本发明公开一种基于凸优化的M通道过采样图滤波器组的优化设计方法，利用两步法进行目标优化，同时全面考虑图滤波器组的频谱选择性和完全重构条件，从而保证在较小的重构误差条件下改善图滤波器组的整体性能。本发明为降低图滤波器组的重构误差，实现信号的恢复重构提供了简单有效的解决方案。

说明书

技术领域

本发明涉及多速率信号处理技术领域，具体涉及一种基于凸优化的M通道过采样图滤波器组的优化设计方法。

背景技术

图信号处理在社交、生物、交通、传感器网络处理以及图像处理等方面有着重要的应用意义。在实际应用中，图数据量十分庞大，给图信号处理带来了一定的难度。多分辨率分析是一种局部化时频分析，可对信号进行多尺度细化分析，有效提取信号的特征信息，但是其对于大数据处理有一定的局限性。然而，对于大数据而言，小波和滤波器组是较为适合处理大规模的图信号的，小波和滤波器组可以对图信号进行稀疏表示，从而处理大数据更有优势。图滤波器组从频谱设计角度可以分为正交、双正交两大类，从信号采样角度分为临界采样和过采样两大类，按通道数可分为两通道和M通道两大类。两通道图滤波器组设计最早是Narang和Ortega提出的完全重构正交的临界采样的小波滤波器组设计方法，对于二分图用正交镜像滤波的方法来进行处理，此方法设计出的滤波器具有混叠消除，完全重构和正交性特性，但该设计方法是针对于二分图或可以分解为二分图的图信号的，其存在一定的局限性。于是他们又提出了对于任意无向图的紧支撑的双正交的小波滤波器组的方法，松弛正交性条件成为双正交，双正交的图小波是频域紧支撑且节点K阶局域化的，仍然满足完全重构条件，可以对任意图信号进行处理，但此方法未考虑滤波器的频谱选择性，从而有一定的局限性。图滤波器组的设计方法发生了由正交到双正交，临界采样到过采样，两通道到M通道的转变。

图信号处理从临界采样到过采样突破了原有的局限性，过采样后的信号包含了原来图信号的全部信息，更加方便图信号的进一步处理。过采样的图拉普拉斯矩阵对于图信号处理的发展和应用至关重要，过采样对于规则信号的处理有更大的设计自由。过采样对于M通道图滤波器组的设计具有极大的优势，M通道图滤波器组具有更多的子带划分，对于大规模的信号处理有重大的意义。目前，针对M通道图滤波器组的研究相对较少，进一步的研究有待提出。

发明内容

本发明所要解决的是现有M通道过采样图滤波器组整体性能较差的问题，提供一种基于凸优化的M通道过采样图滤波器组的优化设计方法。

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

基于凸优化的M通道过采样图滤波器组的优化设计方法，具体包括如下步骤：

步骤1、从频谱特性方面考虑来设计M通道图滤波器组的分析滤波器，在3dB约束条件下，以分析滤波器组的通带波纹和阻带能量为目标函数，通过半正定规划求解器来求解使得分析滤波器组的通带失真和阻带能量最小的分析滤波器；优化问题即

i＝0,1,…,M/2-1；j＝1,2,…,M/2-1

式中，E_p(h)为分析滤波器的通带波纹能量；E_s(h)为分析滤波器的阻带能量；α为分析权值；h为分析滤波器组，h＝[h₀；h₁；…；h_M/2-1]；(·)为第i个分析滤波器h_i的频率向量，L_hi为第i个分析滤波器h_i的长度，h_i为第i个分析滤波器；当i＝1,2,…,M/2-2时，h_i＝[0_i1,I_i,0_i2]h；当i＝0时，h_i＝[I_i,0_i2]h；当i＝M/2-1时，h_i＝[0_i1,I_i]h；0_i1为大小为L_hi×(L_h0+L_h1+…+L_h(i-1))的全零矩阵，0_i2为大小为L_hi×(L_h(i+1)+L_h(i+2)…+L_h(M/2-1))的全零矩阵，I_i为大小为L_hi×L_hi的单位矩阵；(·)为第j个分析滤波器h_j的频率向量，L_hj为第j个分析滤波器h_j的长度，h_j为第j个分析滤波器，当j＝1,2,…,M/2-2时，h_j＝[0_j1,I_j，0_j2]h；当j＝M/2-1时，h_j＝[0_j1,I_j]h；0_j1为大小为L_hj×(L_h0+L_h1+…+L_h(j-1))的全零矩阵，0_j2为大小为L_hj×(L_h(j+1)+L_h(j+2)…+L_h(M/2-1))的全零矩阵，I_j为大小为L_hj×L_hj的单位矩阵；T代表转置；M为通道数；

步骤2、将步骤1所解出的分析滤波器作为已知条件，从完全重构特性方面考虑来设计M通道图滤波器组的综合滤波器，在完全重构约束条件下，以综合滤波器组的阻带能量为目标函数，通过半正定规划求解器来求解使得综合滤波器的阻带衰减最大的综合滤波器；优化问题即

s.t.|a^T(λ_k)g-b_k|≤ε_r；g₁(0)＝0

k＝0,1,…,N；l＝0,1,…,M/2-1

式中，E_s(g)为综合滤波器的阻带能量；β为综合权值；g为综合滤波器组，g＝[g₀；g₁；…；g_M/2-1]；(·)为第l个综合滤波器g_l的频率向量，L_gl为第l个综合滤波器g_l的长度，g_l为第l个综合滤波器，当l＝1,2,…,M/2-2时，g_l＝[0_gl1,I_gl,0_gl2]g；当l＝0时，g_l＝[I_gl,0_gl2]g；当l＝M/2-1时，g_l＝[0_gl1,I_gl]g；0_gl1为大小为L_gl×(L_g0+L_g1+…+L_g(l-1))的全零矩阵，0_gl2为大小为L_gl×(L_g(l+1)+L_g(l+2)+…+L_g(M/2-1))的全零矩阵，I_gl为大小为L_gl×L_gl的单位矩阵；(·)为第l个分析滤波器h_l的频率向量，L_hl为第l个分析滤波器h_l的长度，h_l为第l个分析滤波器；当l＝1,2,…,M/2-2时，h_l＝[0_hl1,I_hl,0_hl2]h；当l＝0时，h_l＝[I_hl,0_hl2]h；当l＝M/2-1时，h_l＝[0_hl1,I_hl]h；h为分析滤波器组，O_hl1为大小为L_hl×(L_h0+L_h1+…+L_h(l-1))的全零矩阵，O_hl2为大小为L_hl×(L_h(l+1)+L_h(l+2)…+L_h(M/2-1))的全零矩阵，I_hl为大小为L_hl×L_hl的单位矩阵；T代表转置；λ_k为第k个频率离散点，b_k为完全重构值；ε_r为重构误差容限；a^T(λ_k)为第k个频率离散点λ_k处的响应向量；g₁(0)是第1个综合滤波器g₁在频率零点处的取值；N+1为频率离散点数；M为通道数；

步骤3、由步骤1所得的分析滤波器h和步骤2所得的综合滤波器g即可构成M通道图滤波器组。

上述步骤1和2中，采用半正定规划求解器SDPT3有效求解分析滤波器和综合滤波器。

上述方案中，分析权值α与综合权值β相等。

上述方案中，通道数M为大于2的偶数。

上述方案中，λ_k的取值范围为[0,2]。

上述方案中，完全重构值b_k＝2。

与现有技术相比，本发明利用两步法进行目标优化，同时全面考虑图滤波器组的频谱选择性和完全重构条件，从而保证在较小的重构误差条件下改善图滤波器组的整体性能。本发明为降低图滤波器组的重构误差，实现信号的恢复重构提供了简单有效的解决方案。

附图说明

图1为M通道过采样图滤波器组的基本结构。

图2为本发明实例1中优化后所得到图滤波器的幅度响应。

图3为本发明实例2中优化后所得到图滤波器的幅度响应。

图4为本发明使用实例2中所得到的图滤波器系统与已有方法得到的图滤波器系统进行去噪实验的仿真结果。其中(a)为输入信号；(b)为噪声信号(σ＝1/2)；(c)为基于现有方法3(临界采样双正交)设计的图滤波器组去噪后的输出信号；(d)为本发明设计的图滤波器组去噪后的输出信号。

具体实施方式

为了便于理解，下面以四通道图滤波器组为例对本发明进行进一步详细说明。

将M通道图滤波器组的分析子带滤波器H_i和综合子带滤波器G_i分别表示为关于分析频谱核和综合频谱核的函数，表示出图滤波器组的完全重构条件；即

式中，H_i为第i个分析子带滤波器，G_i为第i个综合子带滤波器，h_i(λ)为第i个分析子带滤波器频谱核在特征根λ处的取值，g_i(λ)为第i个综合子带滤波器频谱核在特征根λ处的取值，λ为图G的拉普拉斯矩阵的特征根，σ(G)是由图G的拉普拉斯矩阵所有特征根构成的特征空间，P_λ是特征空间的正交投影矩阵，i是子带滤波器的序号。

若g_k(λ)＝h_k+M/2(2-λ),g_k+M/2(λ)＝h_k(2-λ),则完全重构条件可以表示为：

式中，λ∈[0,2],λ是连续变量。

四通道图滤波器组的关于子带滤波器的分析频谱核和综合频谱核表示为：

式中，L_h0,L_h1,L_g0,L_g1分别表示分析滤波器h₀,h₁和综合滤波器g₀,g₁的长度。

频率的记法和表示类似传统滤波器组，λ_pd0,λ_sd0表示h₀,g₀的通带和阻带的截止频率，λ_pd1,λ_pd2,λ_sd1,λ_sd2表示h₁,g₁的通带和阻带的截止频率。低通滤波器的通带波纹可用如下积分来衡量：

当i＝0,a＝0,b＝λ_pd0时上式为E_P(h₀),即h₀的通带波纹的能量；当i＝1,a＝λ_pd1,b＝λ_pd2时上式为E_P(h₁),即h₁的通带波纹的能量。阻带衰减由阻带能量来决定：

在四通道图滤波器组中，若g₀(λ)＝h₂(2-λ)，g₁(λ)＝h₃(2-λ)，g₂(λ)＝h₀(2-λ)，g₃(λ)＝h₁(2-λ)，则图滤波器组是双正交的。完全重构条件可以表示为：

g₀(λ)h₀(λ)+g₀(2-λ)h₀(2-λ)+g₁(λ)h₁(λ)+g₁(2-λ)h₁(2-λ)＝2>

为了方便表示，记完全重构条件可表示为：

重构误差可表示为E(λ)＝l(λ)+l(2-λ)-2，设计中希望重构误差尽可能小。

充分考虑图滤波器组的频谱特性和完全重构条件，可以将图滤波器组的设计问题归结为带约束优化问题，通过两步法来设计图滤波器组。上述两个约束优化问题都是半正定规划问题，都可有效的求解。设计中将图滤波器g₀,g₁的阻带能量转换为关于滤波器g₀,g₁的函数，为了简化和方便计算记h＝[h₀；h₁]；g＝[g₀；g₁]；则上述指标都可转换为关于滤波器h,g的函数，通过半正定规划求最优解。

图1给出了一个通道数为M的过采样图滤波器组系统，基于上述结构的设计方法，即基于凸优化的M通道过采样图滤波器组的优化设计方法，包括如下步骤：

第一步：设计分析滤波器，以分析滤波器组的通带波纹和阻带能量为目标函数，在3dB约束条件下求解使得通带失真和阻带能量最小的分析滤波器，优化是半正定规划问题，可有效的求解。

3dB约束即h₀(0.5)＝1,h₁(0.5)＝1,h₁(1)＝1；α,β为权值，通常取α＝β。为了便于求解，记

h＝[h₀；h₁]；h₀＝[I₀,0]h＝B₀h；h₁＝[0,I₁]h＝B₁h>

其中，B₀是大小为L_h0×(L_h0+L_h1)的矩阵，B₁是大小为L_h1×(L_h0+L_h1)的矩阵，I₀为L_h0×L_h0的单位矩阵，I₁为L_h1×L_h1的单位矩阵，0为全零矩阵。约束求解问题可简化如下：

第二步：以解出的分析滤波器为已知条件，在完全重构约束条件下，考虑到最大化综合滤波器阻带衰减，以综合滤波器组的阻带能量为目标函数来求解综合滤波器。

其中ε_r是重构误差容限，N+1是离散点数，许多实例表明N＝100则可保证重构误差精确度。为方便计算，做如下代换：

g＝[g₀；g₁]；g₀＝[I₀,0]g＝C₀g；g₁＝[0,I₁]g＝C₁g；>

其中，C₀是大小为L_g0×(L_g0+L_g1)的矩阵，C₁是大小为L_g1×(L_g0+L_g1)的矩阵，I₀为L_g0×L_g0的单位矩阵，I₁为L_g1×L_g1的单位矩阵，0为全零矩阵。上述问题可等效写为：

其中b_k＝2,k＝0,…,N，记

向量a^T(λ_k)表示如下：

a^T(λ_k)＝d(λ_k)+d(2-λ_k),k＝0,…,N>

上述优化问题均为半正定规划问题，可通过半正定规划求解器SDPT3有效求解，解得的图滤波器为λ的实值函数。

上述设计过程可以推广到任意M(M>2且为偶数)通道图滤波器组。推广的基于凸优化的M通道过采样图滤波器组的优化设计方法，包括如下步骤：

步骤1、从频谱特性方面考虑来设计M通道图滤波器组的分析滤波器，在3dB约束条件下，以分析滤波器组的通带波纹和阻带能量为目标函数，通过半正定规划求解器来求解使得分析滤波器组的通带失真和阻带能量最小的分析滤波器；优化问题即

i＝0,1,…,M/2-1；j＝1,2,…,M/2-1

式中，E_p(h)为分析滤波器的通带波纹能量；E_s(h)为分析滤波器的阻带能量；α为分析权值；h为分析滤波器组，h＝[h₀；h₁；…；h_M/2-1]；(·)为第i个分析滤波器h_i的频率向量，L_hi为第i个分析滤波器h_i的长度，h_i为第i个分析滤波器；当i＝1,2,…,M/2-2时，h_i＝[0_i1,I_i,0_i2]h；当i＝0时，h_i＝[I_i,0_i2]h；当i＝M/2-1时，h_i＝[0_i1,I_i]h；0_i1为大小为L_hi×(L_h0+L_h1+…+L_h(i-1))的全零矩阵，0_i2为大小为L_hi×(L_h(i+1)+L_h(i+2)…+L_h(M/2-1))的全零矩阵，I_i为大小为L_hi×L_hi的单位矩阵；(·)为第j个分析滤波器h_j的频率向量，L_hj为第j个分析滤波器h_j的长度，h_j为第j个分析滤波器，当j＝1,2,…,M/2-2时，h_j＝[0_j1,I_j，0_j2]h；当j＝M/2-1时，h_j＝[0_j1,I_j]h；0_j1为大小为L_hj×(L_h0+L_h1+…+L_h(j-₁₎)的全零矩阵，0_j2为大小为L_hj×(L_h(j+1)+L_h(j+2)…+L_h(M/2-1))的全零矩阵，I_j为大小为L_hj×L_hj的单位矩阵；T代表转置；M为通道数；

步骤2、将步骤1所解出的分析滤波器作为已知条件，从完全重构特性方面考虑来设计M通道图滤波器组的综合滤波器，在完全重构约束条件下，以综合滤波器组的阻带能量为目标函数，通过半正定规划求解器来求解使得综合滤波器的阻带衰减最大的综合滤波器；优化问题即

s.t.|a^T(λ_k)g-b_k|≤ε_r；g₁(0)＝0

k＝0,1,…,N；l＝0,1,…,M/2-1

式中，E_s(g)为综合滤波器的阻带能量；β为综合权值；g为综合滤波器组，g＝[g₀；g₁；…；g_M/2-1]；(·)为第l个综合滤波器g_l的频率向量，L_gl为第l个综合滤波器g_l的长度，g_l为第l个综合滤波器，当l＝1,2,...,M/2-2时，g_l＝[0_gl1,I_gl,0_gl2]g；当l＝0时，g_l＝[I_gl,O_gl2]g；当l＝M/2-1时，g_l＝[O_gl1,I_gl]g；O_gl1为大小为L_gl×(L_g0+L_g1+…+L_g(l-1))的全零矩阵，O_gl2为大小为L_gl×(L_g(l+1)+L_g(l+2)+…+L_g(M/2-1))的全零矩阵，I_gl为大小为L_gl×L_gl的单位矩阵；(·)为第l个分析滤波器h_l的频率向量，L_hl为第l个分析滤波器h_l的长度，h_l为第l个分析滤波器；当l＝1,2,…,M/2-2时，h_l＝[0_hl1,I_hl,0_hl2]h；当l＝0时，h_l＝[I_hl,0_hl2]h；当l＝M/2-1时，h_l＝[0_hl1,I_hl]h；h为分析滤波器组，0_hl1为大小为L_hl×(L_h0+L_h1+…+L_h(l-1))的全零矩阵，0_hl2为大小为L_hl×(L_h(l+1)+L_h(l+2)…+L_h(M/2-1))的全零矩阵，I_hl为大小为L_hl×L_hl的单位矩阵；T代表转置；λ_k为第k个频率离散点，b_k为完全重构值；ε_r为重构误差容限；a^T(λ_k)为第k个频率离散点λ_k处的响应向量；g₁(0)是第1个综合滤波器g₁在频率零点处的取值；N+1为频率离散点数；M为通道数；

步骤3、由步骤1所得的分析滤波器h和步骤2所得的综合滤波器g即可构成M通道图滤波器组。

上述优化问题均为半正定规划问题，可通过半正定规划求解器SDPT3有效求解，解得的图滤波器为λ的实值函数。

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

仿真实例1：

设计一个图滤波器组，参数如下：L_h0＝12，L_h1＝12，L_g0＝11，L_g1＝11，λ_pd0＝0.4，λ_pd1＝0.6，λ_pd2＝0.8，λ_sd0＝0.75，λ_sd1＝0.25，λ_sd2＝1.3，α＝0.1，ε_r＝10^-9，得到的图滤波器组幅度响应如图2所示。图中PR表示l(λ)+l(2-λ),λ∈[0,2]的取值，用于度量图滤波器组的重构误差。仿真计算得到的最大重构误差和正交性的值分别为E_max＝3.9346×10^-9,Θ＝0.9353，相同的图滤波器组长度和运行环境下，表1给出了本发明与现有方法1(接近正交的M通道过采样)设计的图滤波器组的重构性能和边界比对比。

表1

R_B为边界比，仿真时取R_B的最大值进行对比，R_B等于1表明滤波器组是完全重构的，此时滤波器组的重构特性较好。本发明设计出的图滤波器组的最大边界比为1，明显优于现有方法设计的图滤波器组。通常情况下，采用正交的方法设计出的图滤波器组重构误差较双正交的大，因此，本发明在重构特性上有着显著的优势。从表1中可以看出本发明所设计的双正交图滤波器组显著降低了重构误差。从图滤波器组幅度响应可以看出两种算法设计的图滤波器组的整体频谱特性相当。在图滤波器组的设计中，为了尽可能无失真地恢复原信号，期望得到的图滤波器组的重构误差更小，本发明设计方法有着更好的重构性能，所以本发明设计所得的图滤波器组整体性能更佳。

仿真实例2：

设计一个图滤波器组，参数如下：L_h0＝8，L_h1＝8，L_g0＝7，L_g1＝7，λ_pd0＝0.3，λ_pd1＝0.6，λ_pd2＝0.8，λ_sd0＝0.75，λ_sd1＝0.25，λ_sd2＝1.3，α＝0.1，ε_r＝10^-13，得到的图滤波器组幅度响应如图3所示。相同的图滤波器组长度和运行环境下，表2给出了本发明与现有方法2(M通道过采样)设计的图滤波器组的重构性能和正交性对比。

表2

正交性Θ越接近于1表明滤波器组的正交性越好。对比可知，本发明设计的图滤波器组的重构误差明显比现有方法小，其正交性和整体频谱特性比现有方法略好。与例1相比，改变图滤波器组的参数设置，综合滤波器的频谱特性会有所变差，但是可以获得更小的重构误差，重构误差是衡量图滤波器组整体性能的主要性能指标，所以本发明设计所得的图滤波器组整体性能更好。适当松弛重构误差容限可得到频谱选择性更好的图滤波器，可根据需要灵活设计图滤波器组。

用上述例2参数设计的图滤波器组对密尼苏达交通图进行去噪实验，噪声标准差σ＝1/2时，使用现有方法3(临界采样双正交)和本发明方法设计的图滤波器组去噪的仿真结果如图4所示，其中(a)为输入信号；(b)为噪声信号(σ＝1/2)；(c)为基于现有方法3(临界采样双正交)设计的图滤波器组去噪后的输出信号；(d)为本发明设计的图滤波器组去噪后的输出信号。表3给出了噪声σ取不同值时的信噪比对比。

表3

相同的参数和运行环境下，本发明与现有方法3(临界采样双正交)相比具有较高的信噪比，与现有方法2(M通道过采样)相比，信噪比稍差1.5dB左右。