一种基于泰勒级数展开的图滤波器设计方法

葛子瑞，杨 震，2，郭振超

(1.南京邮电大学通信与信息工程学院，江苏南京 210003 2.南京邮电大学通信与网络技术国家地方联合工程研究中心，江苏南京 210003)

图信号处理（Graph Signal Processing，GSP）作为一种新兴的数据处理技术正在快速发展，其将经典的信号处理技术拓展到带有复杂结构的非规则数据集上［1－3］。目前GSP在许多领域中都有相应的应用和发展，例如在温度传感网络［4］，交通网络［5］，语音增强［6－8］， 突变点检测［9］等领域都有了大量的文献研究。在GSP中，信号位于图的顶点上，各信号点之间的关系由图的边与权重来确定［8］。正是由于GSP引入了边与权重的概念，它才会相较于传统的信号处理（Digital Signal Processing，DSP）有更多的应用场景和更好的处理效果。

与DSP的离散傅里叶变换（Discrete Fourier Transform，DFT）类似，图傅里叶变换（Graph Fourier Transform，GFT）也是GSP的一个重要信号处理方式。相较于DFT的唯一性，GFT是由图来确定，不同的图会定义不同的图傅里叶变换。例如基于图拉普拉斯矩阵的GFT，为了避免图傅里叶变换的多样性，让不同的图有着相同的滤波器设计方法，文献［10］和［11］提出了通用滤波器设计（Universal Design Approach）的概念，即让设计出的滤波器可以适用于任一图结构。本文也将遵循这一原则来设计一种基于图拉普拉斯矩阵设计相应的图滤波器。

与DSP中的滤波器设计方式类似，GSP的滤波器设计包括有限长度（Finite Impulse Response，FIR）滤波器设计和无限长度（Infinite impulse response，IIR）滤波器设计方式。在图FIR滤波器设计中，其中一种主流的方式是通过最小二乘使图频率响应拟合期望的频谱［12］，这种方式虽然能够在期望的频率点上得到较好的响应，但是频率点之间的响应往往有很大的波动。但在文献［10］中对这种方式进行改进，以此希望所有的频率点都能够满足期望的频率响应，而文献［13］是对文献［11］目标函数的另一种解法。另一种方式是通过多项式来拟合期望的频率响应，如切比雪夫多项式来设计FIR图滤波器，不过这种方式并不适用于IIR图滤波器［14］。IIR图滤波器的设计方式与FIR滤波器设计方式类似，均是通过采用优化迭代的方式来求解滤波器系数，具体的优化方程会有些不同。

目前GSP滤波器的研究经过了长时间的发展，可以设计出与DSP滤波器相对应的图滤波器。GSP滤波器的设计目前主要思路还是利用求解优化方程来计算滤波器系数。本文提出一种新的设计思路，利用现在已有的DSP滤波器来设计响应的GSP滤波器。主要设计思路为将设计好的DSP滤波器的傅里叶基函数进行泰勒展开，然后将DSP滤波器的频率替换为归一化后的GSP滤波器频率。具体的展开方式分为两种：第一种方式仅是只展开一阶DSP的傅里叶基函数，其他阶的傅里叶级数由该一阶展开表示，此方法的拟合误差较小适合于拟合高阶的DSP滤波器；第二种设计方式是展开所有阶的傅里叶基函数，但这种设计方式拟合误差较大，只适合于低阶的DSP滤波器拟合，但此方法的滤波器系数容易计算。为了改善第二种方式的拟合误差，本文将基函数的展开点由0点处转为处，经过相应的误差分析可以看出此改进可以极大减少第二种方式的拟合误差。

本文最后对比了文献［12］中提出的经典GSP滤波器设计方式以及文献［13］中提出的使用多项式平方和来拟合GSP滤波器的设计方式。实验将文献［12］和［13］的方法在设计FIR和IIR滤波器时的最优性能作为对比对象。最终的实验结果表明在设计图滤波器中，本文方法拥有更好的性能。

1 基础原理

1.1 GSP 理论基础

在本文中，用S来表示移位算子，将移位算子进行特征分解得到S＝UΛU－1，其中U表示特征向量矩阵，用来作为图傅里叶变换的基，Λ是对角矩阵，其对角线元素表示为图频率λ1，…，λN，其中图频率的高低由特征值模的大小确定［12］。使用特征向量矩阵U将信号从图域变换到图频域，可以取得图信号 x∈CN的频率表达式，分别表示信号x的图傅里叶变换（Graph Fourier Transform，GFT）及其逆变换。

1.2 图滤波器

图滤波器 H是移位算子 S的函数，即 H＝h（S），将图滤波器 H进行特征分解得到 H ＝Uh（Λ）U－1，其中对角矩阵 h（Λ）作为图滤波器的单位冲激响应。相应地，信号x经过图滤波器H滤波之后的输出为y＝Hx，其频率域表达式为＝， 其中分别表示输入信号和输出信号的GFT变换［8］。GSP中的滤波器与DSP滤波器相似，也有相应的FIR滤波器和IIR滤波器。其中一个K阶的FIR图滤波器［8］表示为

其中，hk表示滤波器系数。确定好滤波器的频率响应之后，需要计算式（2）来得到相关的滤波器系数。

式（2）中，H*（λ）为期望的滤波器响应［9］。 在这种情况下，与传统的图滤波器不同，这里不再限制滤波器的阶数。通过增加图滤波器的阶数可以更好地逼近期望的滤波器，但是相应的求解式（2）的复杂度会增加。

GSP中的IIR滤波器与DSP中IIR滤波器表达式相似。DSP中IIR滤波器的表达式为

与之对应的图IIR滤波器频率响应的表达式为

上述λ来自邻接矩阵的特征分解，是其特征频率；邻接矩阵需要根据不同的应用对象进行单独设计，所以不同对象的频率是变化的，这与DSP的频率有很大不同。

相应地，文献［11］给出了图IIR滤波器的图域表达式

其中，H*（λ）为期望的滤波器响应，为了计算式（6）中的系数，一般将其化成式（7）的形式来计算图IIR滤波器的系数［11］。

2 扩展DSP的滤波器设计

DSP滤波器的设计是在频率ω上进行的，最终目的是滤波器的频率响应f（ω）在［0，π］能够取得期望的特性。与此对应的图滤波器，对于归一化图拉普拉斯矩阵而言，其特征值为实数，且特征值λ的范围在0到2之间，以此移位算子设计图滤波器的最终目的也是让频率响应 h（λ）在［0，2］取得期望的特性。本文将利用这些设计方法的相似性使用已经设计好的DSP滤波器来设计GSP滤波器。

此方法的主要设计思想是利用式（1）的图滤波器，去逼近设计好的DSP滤波器，完成经典数字滤波器到图滤波器的转变。

具体步骤为：首先将归一化图拉普拉斯矩阵的特征值重新归一到［0，π］，然后对已经设计好的DSP滤波器h［n］（共N点）进行傅里叶变换

将H（ω）中的ejω进行泰勒展开得到关于ω多项式的形式，这时就可得到图滤波器相似的形式。具体的展开方式分为两种，不同的展开方式最终取得的误差不同，此外这种方式设计的图滤波器因其系数为复数，更适合处理复数图信号。

2.1 展开 ejω

终，在设计GSP滤波器时，只需要根据具体的情况设计好图滤波器的阻带和通带，然后根据此阻带和通带在DSP中设计滤波器，最后通过复数泰勒展开就可变成GSP滤波器。

由于是一种逼近设计，那么必然存在逼近误差，因此，为了估计最后泰勒展开得到的GSP滤波器频率响应曲线和原DSP频率响应曲线在每一个点的误差，以及得到图滤波器各个阶数的具体系数，下面将进行相应的系数计算和误差分析。

2.2 展开所有 ejωn

从式（21）和式（26）可以看出，虽然误差都与DSP滤波器的阶数和ejωn的展开阶数有关，但是2.1节产生的误差比本节的更小。为了达到相同的误差条件，需要将ejωn展开更多的阶数，不过在GSP滤波器系数的求解问题上本节更为简单。另一方面从式（26）可以看出，不太适合拟合阶数较高的DSP滤波器，而DSP中IIR滤波器能够以较小的阶数来得到期望的滤波器响应，所以本节方法更适合于设计GSP的IIR滤波器。

2.3 方法改进

在2.1和2.2节中，本文均是将ejωn在0处展开，这使得拟合项在ω＝0附近有较小的误差，其中就包括小于0的部分。但是对于GSP滤波器，本文将其频率限制在［0，π］内，这使得对小于0部分的拟合变得无意义（即使这一部分拟合误差很小）；另外对泰勒级数展开而言，距展开点越远，相应的拟合误差越大，此时拟合误差最大的点为ω＝π，到展开点的距离为π。为改进这些不足，本文进一步将展开点设置在区间的中点处，此时有两个好处：（1）由于展开点在ω0＝，避免了将低误差的拟合用在不需要的频率上，使得ω0＝附近都有较小的拟合误差；（2）由于展开点设置在区间中点，区间各点到中点的距离最大为，均在两端点处，相较于原来的距离π，此时距展开点更近，相应的最大误差应比原来误差更小。

3 实验对比

将本文设计方法与文献［12］的最小二乘法设计方法和文献［13］中的设计方法作对比。inverse和optimal分别表示文献［13］中采用直接求解矩阵逆的方式和求解优化的方式；least_square表示文献［12］中提出的最小二乘的方式设计图滤波器；exp1x表示采取本论文的方式一设计滤波器；expnx表示采取本论文中的方式二进行展开；expnx_pi_2表示采用本论文中方式二的改进方法；文献［12］和［14］均采取实验中最好的结果。

图1为各个方法设计FIR低通滤波器的结果，截止频率为0.3π。从图1中可以看出使用本文中的方式一设计滤波器取得的效果最好，这因为使用此方法能够很好地拟合高阶DSP中的FIR滤波器；而对于本论文的方式二，此方法在拟合高阶DSP滤波器时会有很大误差。此外，从图1中还可以看出对改进的方法二，虽然效果仍旧不如方式一，但是相较于原来的方式二，改进方法的确取得了比原来好的效果。inverse方法在滤波器的通带和阻带内均有很大的波动，这主要是由于此方法在求解逆矩阵的过程中涉及到病态矩阵，导致其求得的滤波器系数最大值和最小值的模差别很大，致使滤波器有很大波动；而对于least_square方法同样存在此问题是因为该方法仅仅是对特征值处的点进行拟合，没有考虑到特征值之间的点，最终导致滤波器在通带和阻带内均有很大波动。optimal方法因为设置了过渡带以及对阻带内波动进行了抑制，所以在阻带内的波动减小，但代价是通带内波动增大。

图1 FIR滤波器频率响应

图2为各个方法设计IIR滤波器的实验结果。同样地，inverse和least_square方法在阻带有较大的波动，但是通带内的波动均小于原来FIR滤波器，这主要是因为IIR滤波器设计的阶数较少，本文inverse方法使用的阶数仅为4阶，此时病态矩阵的逆很好计算，不过使用文献［14］中的方法会使滤波器阶数无法变得很大，还有可能随着阶数提高滤波器性能变差，甚至到一定程度使得结果无法收敛。使用方法一，方法二以及方法二的改进方法设计的IIR滤波器几乎得到了相同的曲线，这主要是因为拟合的DSP滤波器的阶数较低，方法二同样可以很好地拟合DSP滤波器，不过在频率尾部有些波动。

图2 IIR滤波器频率响应

表1和表2分别为设计所得到的滤波器与理想滤波器之间的相对误差，误差的计算方式为‖H（λ）－H*（λ）‖／‖H*（λ）‖。 其中表1和表2误差较大的主要原因是所设计的滤波器截止频率较低为0.3π，由于实际得到的滤波器存在过渡带使得阻带截止频率在0.4π左右，这样造成的相对误差会比较大。从表1的实际结果来看对于FIR滤波器使用本文的方法一效果最好，改进的方法二效果次之。从表2的结果来看，本文的3个方法效果相当，均取得了较好的效果。

表1 FIR滤波器与理想滤波器的设计误差

表2 IIR滤波器与理想滤波器的设计误差

表3为各个算法的复杂度比较，least_square和optimal算法的复杂度分别根据文献［18］计算，其中｜λ｜max表示二次型矩阵特征值模的最大值，γ表示正则项系数，h0和h*分别表示迭代初始解和最优解，ε表示所要求的误差；本文计算滤波器系数的算法复杂度，可根据式（11）、（23）和（31）计算得到。 从表3中可以看出，本文算法的复杂度和GSP滤波器阶数有很大关系，且方法一的复杂度最高。对于least_square和optimal方法其迭代次数和要求误差ε成反比。对于inverse方法，其计算复杂度主要在于计算矩阵的逆，其复杂度为求逆的复杂度。

表3 算法复杂度

4 结束语

本文在现有DSP滤波器设计的基础上，通过研究GSP中拉普拉斯矩阵的性质，将DSP的滤波器设计方式引入到GSP的滤波器设计中。将DSP滤波器的傅里叶基函数进行分解得到关于角频率多项式的形式，然后将此多项式的系数作为GSP滤波器的系数。具体的分解方式有两种，一种是仅将一阶ejω进行分解，另一种是将所有ejωn都进行分解；其中第一种分解方式有较低的误差，可以拟合高阶的DSP滤波器，但是滤波器系数不易计算；第二种方式的误差相较于第一种高，只适合于低阶的DSP滤波器，由于IIR滤波器不需要很高的阶数就可以达到期望的响应，所以比较适合低阶的IIR滤波器设计，且此方法可以方便地计算出滤波器系数。同时通过与其他GSP滤波器设计方法的比较，可以看出本文提出的方法在设计FIR滤波器和IIR滤波器都具有比较好的性能。

附录：