极点半径受限的IIR数字滤波器设计新算法

赖春露， 姚 统， 王 路

（1.哈尔滨工业大学（威海）信息科学与工程学院，山东 威海 264209；2.山东大学（威海）山东大学澳国立联合理学院，山东 威海 264209）

0 引 言

无限脉冲响应（Infinite Impulse Response，IIR）数字滤波器设计是信号处理领域的一个经典研究课题［1-3］，但仍存在值得研究的挑战性。MATLAB中的‘信号处理’及‘数字信号处理/滤波器设计’工具箱，包含了多个设计不同类型IIR数字滤波器的函数。如butter、cheby1、cheby2、ellip、iirlpnorm等，可用于巴特沃斯滤波器、切比雪夫滤波器、椭圆滤波器和具有分段线性期望幅频响应的IIR滤波器（包括微分器）的设计，得到的数字滤波器是稳定的。但现有函数大都没有用于控制最大极点半径的输入参数，无法保证任意给定的稳定裕量要求。

在对滤波器进行量化实现时，离单位圆过近的极点可能导致对应频率点的频率响应不准确，甚至导致滤波器极点跑出单位圆而变成不稳定滤波器［1］。为了提高频率响应的准确度，增大滤波器的稳定裕量，在设计时，除了要求IIR数字滤波器的极点在单位圆内，还经常要对极点半径的具体大小进行限制，使极点远离单位圆。基于此，本文考虑极点半径受限的IIR数字滤波器，提出其优化设计的新算法。

数字滤波器优化设计的一个典型准则是最大幅频响应误差最小化，即minimax设计，也称为切比雪夫设计。IIR数字滤波器的minimax设计是关于滤波器系数的一个高度非凸的优化问题。若只要求滤波器极点在单位圆内，可将该问题转化为以幅值平方响应函数的分子、分母三角多项式系数为变量的一个线性规划问题［2，4］。然后可应用谱因子分解技术［5］，并选择幅值平方响应的单位圆内极点作为滤波器极点，使滤波器稳定，但无法保证其满足最大极点半径的要求。

文献［6-12］中提出的设计算法，考虑了IIR数字滤波器的最大极点半径约束。文献［6］中将极点半径约束结合到基于Rouches定理的稳定性充分条件中，应用高斯-牛顿策略将约束最小二乘设计转化为凸二次规划问题。文献［7］中将极点半径约束结合到基于分母正实性的稳定性充分条件中，应用LSK策略将约束minimax设计转化为凸二次约束二次规划问题。文献［8-11］中则将极点半径约束与基于稳定三角形的稳定性充要条件相结合，应用约束优化、高斯-牛顿策略、序贯部分优化等方法，将非凸的约束IIR数字滤波器设计问题，转化为凸优化问题来求解。以上这些方法，可在极点半径约束下使滤波器的幅频响应和相频响应要求同时得到满足。但是，在某些应用场合，对相频响应是不做要求的。针对这些应用，按没有相频响应要求的滤波器进行设计，可带来幅频响应性能的提升［12］。

MATLAB滤波器设计工具箱中的iirlpnormc函数考虑了幅频响应的极点半径约束lp-范数最小化设计，并采用约束牛顿算法进行求解［1］。文献［13］中也考虑极点半径受限的幅频响应设计，但把问题转化为半定规划来求解。由于幅频响应逼近问题的高度非凸性，得到的滤波器经常是次优的。文献［14］中将极点半径受限的幅频响应切比雪夫设计，转化为极点半径受限的复数频率响应切比雪夫逼近问题，改进了设计结果。但该方法采用基于滤波器分母正实性的稳定性充分条件，限制了滤波器性能的提升。

本文考虑最大极点半径受限的幅频响应切比雪夫设计，采用基于单一高阶多项式分子和级联二阶多项式分母的滤波器表示，以及充要的稳定三角形条件。应用系列复数频率响应逼近方法［12］及高斯牛顿策略，将设计问题转化为二阶锥规划问题进行求解。通过设计算例，将本文方法与现有方法的设计结果进行了比较。结果表明，在满足给定最大极点半径的前提下，本文方法得到了具有更小幅频响应误差的滤波器。

1 滤波器表示和极点半径约束

数字滤波器的设计问题是指选择其传递函数的系数，使对应的频率响应逼近某个期望响应。一个IIR数字滤波器可用多种形式的传递函数H（z）进行表示。如分子和分母均为单一高阶多项式的传递函数：

式中，M和N分别是滤波器的分子阶数和分母阶数。

再如分子和分母均为级联二阶因子的传递函数：

式中，滤波器的阶数N不妨假设为偶数。

相比式（1）和（2）的两种滤波器表示，传递函数（1）关于其系数向量的非线性程度低，特别是，关于其分子多项式系数向量是纯线性的。传递函数（2）则适合应用充分必要的稳定三角形条件，对其最大极点半径进行限制。

本文采用式（1）和（2）相结合的滤波器表示，即传递函数的分子表示为单一多项式，分母表示为级联二阶因子，即：

式中

分别为分母和分子的系数向量。

基于式（3）表示的传递函数，滤波器的极点就是二阶分母多项式

的零点，其中an＝［an1，an2］T为An（z，an）的系数向量。设计时，要求滤波器的所有极点半径都小于r＜1，即要求所有An（z，an）的零点半径都小于r＜1。容易证明，这等价于要求：对任意n∈｛1，2，…，N/2｝成立。

在an平面上，式（5）确定的区域为图1所示的三角形，称为稳定三角形［1，8-11］。记

则IIR滤波器H（z，a，b）的所有极点半径均小于r的充分必要条件是a∈S（r）。

2 极点半径受限的幅频响应切比雪夫设计算法

2.1 约束切比雪夫设计问题

用ω表示数字角频率，D（ω）表示期望幅频响应。本文研究极点半径受限的IIR数字滤波器的幅频响应切比雪夫设计，即在最大极点半径小于r的前提下，最小化通带和阻带上幅频响应误差的最大值，数学上可描述为

式中：Ωp和Ωs代表滤波器的通带和阻带；｜H（ejω，a，b）｜表示滤波器复数频率响应H（ejω，a，b）的幅值。用δ表示通带Ωp加阻带Ωs上幅频响应误差的最大值，并注意到S（r）的表达式（6），约束切比雪夫设计问题（7）可等价地描述为

式中，x＝［aT，bT］T。

由于以下两个原因，约束切比雪夫设计问题（8）是一个高度非凸的优化问题。一方面，由式（3）可见，滤波器的传递函数H（z，x）关于其分母系数向量a是高度非线性的；另一方面，当δ给定时，由式（8b）确定的任意通带频率点ω上的幅频响应约束域是图2所示频率响应复平面上的圆环，是一个典型的非凸域。

2.2 凸优化策略

针对以上切比雪夫设计问题的非凸性，首先根据高斯-牛顿策略，应用泰勒级数展开式，在系数向量x的当前迭代点x（k）附近（其中k代表迭代步），将频率响应H（ejω，x）近似为x的线性函数，即

式中：g（z，x）≡∇xH（z，x）表示传递函数H（z，x）关于x的梯度；β＞0是某个小的正数。

其次，将幅频响应误差约束（8b）转化为如下复数频率响应误差约束：

式中，ϕ（ω）为某个相位函数。当ϕ（ω）已知时，约束（10）形成频率响应H（ejω，x）复平面上一个以D（ω）ejϕ（ω）为圆心、δ为半径的圆盘。用式（9）右端替代式（10）中的H（ejω，x）后，式（10）形成的约束域则是（x，δ）空间中的一个二阶锥。这样，高度非凸的约束切比雪夫设计问题（8），就变成一个容易求解的二阶锥规划问题。可以证明，当ϕ（ω）等于问题（7）之最优滤波器的相频响应ϕ*（ω）时，即当ϕ（ω）＝ϕ*（ω）＝∠H（ejω，x*）时，该二阶锥规划问题的解，就是原设计问题（7）的解。

可惜的是，相位函数ϕ*（ω）是未知的。为此，类似文献［12］，本文采用迭代方法对ϕ*（ω）进行估计。在第l次迭代，让ϕ（ω）等于某个ϕl（ω），然后用迭代中求得的滤波器相频响应来更新ϕl（ω）。通过交替地更新ϕl（ω）和求解相应的二阶锥规划问题，可估计出ϕ*（ω）并最终求得约束切比雪夫设计问题（7）的最优滤波器。

2.3 求解算法

根据以上讨论，在给定期望幅频响应D（ω）（连同通带Ωp和阻带Ωs）及期望最大极点半径r＜1后，本文的极点半径受限的IIR数字滤波器切比雪夫设计算法可描述如下：

第1步 令ϕ0（ω）＝0、x0＝0、l＝0。

第2步 令k＝0、xl（k）＝xl。

第3步 计算H（ejω，xl（k））及g（ejω，xl（k））。

第4步 求解以下二阶锥规划问题

第5步 若‖xl（k（1）-xl（k）‖2＞ε‖xl（k）‖2（其中ε＞0是一个小的正数，如ε＝10-4），令k＝k＋1，返回第3步。否则，令xl＋1＝xl（k＋1），ϕl＋1（ω）＝∠H（ejω，xl＋1）。

第6步若l＝lmax，停止算法（其中lmax为最大外循环次数，即更新相位函数ϕl＋1（ω）的迭代次数）。否则，令l＝l＋1并返回第2步。

3 设计算例及比较

下面通过设计算例及与现有方法的比较展示本文算法的优越性。本文算法已编制成MATLAB函数iirmmc（M，N，f，d，r），其中M、N、f、d、r分别为滤波器的分子阶数、分母阶数、带边频率集、带边期望幅值集和最大极点半径。实验在MATLAB 2017中进行，频率区间［0，π］被离散化为频率点集｛ωk＝kπ/400：k＝0，1，…，400｝，算法参数ε＝10-4、β＝0.01、lmax＝400。

算例1设计一个M＝N＝12阶的低通滤波器。要求通带Ωp＝［0，0.5π］，阻带Ωs＝［0.55π，π］，极点半径不大于0.94，即f＝［0，0.5，0.55，1］，d＝［1，1，0，0］，r＝0.94。

首先用本文算法iirmmc（M，N，F，D，r）函数进行设计。图3（a）是得到的最大幅频响应误差随外循环次数l的变化曲线。由图可见，本文算法收敛很快。图3（b）是所得滤波器的零极点分布，满足预设的最大极点半径的要求。表1列出了滤波器的最大幅频响应误差、均方根幅频响应误差和最大极点半径，图4所示为滤波器的幅频响应和幅频响应误差曲线。

作为比较，用MATLAB滤波器设计工具箱中的iirlpnorm（M，N，f，f，d）函数所设计的滤波器，其最大极点半径等于0.977 4，超出了预设上界值r＝0.94。改用函数iirlpnormc（M，N，f，f，d，［1，1，1，1］，r）进行设计，得到了满足极点半径要求的滤波器，其最大和均方根幅频响应误差以及幅频响应和响应误差曲线等，分别示于表1和图4。另外，表1和图4还示出了文献［14］得到的滤波器的幅频响应和响应误差。可以看到，在极点半径满足预设要求的前提下，本文算法设计的滤波器，其幅频响应误差明显小于iirlpnormc函数及文献［14］得到的滤波器。特别是，本文算法得到的最大和均方根误差只有iirlpnormc函数的45.24%和68.72%。

表1 算例1中不同算法设计的低通滤波器的性能比较

算例2设计一个M＝N＝8阶的带通滤波器，要求通带Ωp＝［0，0.3π］∪［0.85π，π］，阻带Ωs＝［0.4π，0.8π］，极点半径不大于0.925，即f＝［0，0.3，0.4，0.8，0.85，1］，d＝［0，0，1，1，0，0］，r＝0.925。

用本文算法iirmmc（M，N，f，d，r）进行的设计，最大幅频响应误差随算法外循环次数l的变化曲线及算法收敛后得到的滤波器零极点分布如图5所示，最大和均方根幅频响应误差及最大极点半径列于表2，幅频响应和幅频响应误差曲线如图6所示。函数iirlpnorm（M，N，f，f，d，［1，1，1，1，1，1］）设计的滤波器，最大极点半径为0.951 1，不满足预设值r＝0.925的要求。函数iirlpnormc（M，N，f，f，d，［1，1，1，1，1，1］，r）设计的滤波器，其幅频响应误差比本文方法的结果大，也见表2和图6。本文得到的最大和均方根幅频响应误差分别是iirlpnormc方法的79.08%和89.86%。

表2 算例2中不同方法设计的带通滤波器性能比较

算例3设计一个极点半径不大于0.92、M＝N＝4阶的全带数字微分器（D（ω）＝ω，ω∈［0，π］）。此设计中，f＝［0，1］，d＝［0，1］，r＝0.92。为使得到的滤波器在ω＝0附近的低频段有好的逼近性能，设计时给滤波器的传递函数预设一个z＝1（即ω＝0）的零点。这可通过在算法第4步中给问题（11）增加一个滤波器分母系数之和等于0（即H（1）＝0）的等式约束来实现。

本文算法iirmmc（M，N，f，d，r）得到的微分器，分子多项式及分母二阶因子的系数如表3所示，幅频响应及幅频响应误差曲线示于图7。该微分器的最大极点半径等于0.92，满足预设值要求；其最大和均方根幅频响应误差分别等于6.457×10-3和3.237×10-3。

表3 本方法得到的4阶全带数字微分器

本文也用iirlpnorm（M，N，f，f，d，［1，1］）及iirlpnormc（M，N，f，f，d，［1，1］，r）函数设计了本例的4阶全带微分器，所得微分器的最大极点半径都小于r＝0.92。最大幅频响应误差分别为5.915×10-2和7.061×10-2，比本文的结果大很多。均方根误差分别为4.092×10-2和4.670×10-2，也比本文的结果大很多。表明iirlpnorm和iirlpnormc函数不适合数字微分器的设计。

表4和图8比较了本文方法及四种新近文献方法得到的4阶全带微分器。文献方法包括约束优化方法（Constrained Optimization Method，COM）［8］、序贯部分优化（Sequential Partial Optimization，SPO）［11］方法、蝙蝠算法（Bat Algorithm，BA）［15］及Salp群算法（Salp Swarm Algorithm，SSA）［16］。这些文献方法得到的数字微分器都满足最大极点半径小于r＝0.92的要求。但从表4可见，本文算法比文献方法得到了小很多的最大幅频响应误差，均方根幅频响应误差也比COM方法［8］、SPO方法［11］和BA方法［15］的结果小。由图8可见，本文所得微分器的幅频响应误差，在低频段、中频段和高频段，都比文献［8，11］中所得微分器的幅频响应误差小，群延迟也小。文献［15-16］中的微分器，低、中频段上幅频响应误差都很小，但在高频段特别是ω＝π附近，幅频响应误差明显比本文滤波器的大。

表4 本文方法与文献方法所得滤波器的性能比较

4 结 语

本文针对最大极点半径受限的IIR数字滤波器和微分器的minimax设计问题，提出一个新算法，可以在最大极点半径不超过预设值的前提下，最小化幅频响应误差的最大值。设计算例表明，该算法收敛快，比现有算法得到了更小的幅频响应误差。特别是，本文算法设计的4阶全带微分器，逼近精度高、群延迟低、全频带上的群延迟偏差也不大，可在今后的教学和科研中加以应用。