極點半徑受限的IIR數字濾波器設計新算法

賴春露， 姚 統， 王 路

（1.哈爾濱工業(yè)大學（威海）信息科學與工程學院，山東 威海 264209；2.山東大學（威海）山東大學澳國立聯合理學院，山東 威海 264209）

0 引 言

無限脈沖響應（Infinite Impulse Response，IIR）數字濾波器設計是信號處理領域的一個經典研究課題［1-3］，但仍存在值得研究的挑戰(zhàn)性。MATLAB中的‘信號處理’及‘數字信號處理/濾波器設計’工具箱，包含了多個設計不同類型IIR數字濾波器的函數。如butter、cheby1、cheby2、ellip、iirlpnorm等，可用于巴特沃斯濾波器、切比雪夫濾波器、橢圓濾波器和具有分段線性期望幅頻響應的IIR濾波器（包括微分器）的設計，得到的數字濾波器是穩(wěn)定的。但現有函數大都沒有用于控制最大極點半徑的輸入參數，無法保證任意給定的穩(wěn)定裕量要求。

在對濾波器進行量化實現時，離單位圓過近的極點可能導致對應頻率點的頻率響應不準確，甚至導致濾波器極點跑出單位圓而變成不穩(wěn)定濾波器［1］。為了提高頻率響應的準確度，增大濾波器的穩(wěn)定裕量，在設計時，除了要求IIR數字濾波器的極點在單位圓內，還經常要對極點半徑的具體大小進行限制，使極點遠離單位圓。基于此，本文考慮極點半徑受限的IIR數字濾波器，提出其優(yōu)化設計的新算法。

數字濾波器優(yōu)化設計的一個典型準則是最大幅頻響應誤差最小化，即minimax設計，也稱為切比雪夫設計。IIR數字濾波器的minimax設計是關于濾波器系數的一個高度非凸的優(yōu)化問題。若只要求濾波器極點在單位圓內，可將該問題轉化為以幅值平方響應函數的分子、分母三角多項式系數為變量的一個線性規(guī)劃問題［2，4］。然后可應用譜因子分解技術［5］，并選擇幅值平方響應的單位圓內極點作為濾波器極點，使濾波器穩(wěn)定，但無法保證其滿足最大極點半徑的要求。

文獻［6-12］中提出的設計算法，考慮了IIR數字濾波器的最大極點半徑約束。文獻［6］中將極點半徑約束結合到基于Rouches定理的穩(wěn)定性充分條件中，應用高斯-牛頓策略將約束最小二乘設計轉化為凸二次規(guī)劃問題。文獻［7］中將極點半徑約束結合到基于分母正實性的穩(wěn)定性充分條件中，應用LSK策略將約束minimax設計轉化為凸二次約束二次規(guī)劃問題。文獻［8-11］中則將極點半徑約束與基于穩(wěn)定三角形的穩(wěn)定性充要條件相結合，應用約束優(yōu)化、高斯-牛頓策略、序貫部分優(yōu)化等方法，將非凸的約束IIR數字濾波器設計問題，轉化為凸優(yōu)化問題來求解。以上這些方法，可在極點半徑約束下使濾波器的幅頻響應和相頻響應要求同時得到滿足。但是，在某些應用場合，對相頻響應是不做要求的。針對這些應用，按沒有相頻響應要求的濾波器進行設計，可帶來幅頻響應性能的提升［12］。

MATLAB濾波器設計工具箱中的iirlpnormc函數考慮了幅頻響應的極點半徑約束lp-范數最小化設計，并采用約束牛頓算法進行求解［1］。文獻［13］中也考慮極點半徑受限的幅頻響應設計，但把問題轉化為半定規(guī)劃來求解。由于幅頻響應逼近問題的高度非凸性，得到的濾波器經常是次優(yōu)的。文獻［14］中將極點半徑受限的幅頻響應切比雪夫設計，轉化為極點半徑受限的復數頻率響應切比雪夫逼近問題，改進了設計結果。但該方法采用基于濾波器分母正實性的穩(wěn)定性充分條件，限制了濾波器性能的提升。

本文考慮最大極點半徑受限的幅頻響應切比雪夫設計，采用基于單一高階多項式分子和級聯二階多項式分母的濾波器表示，以及充要的穩(wěn)定三角形條件。應用系列復數頻率響應逼近方法［12］及高斯牛頓策略，將設計問題轉化為二階錐規(guī)劃問題進行求解。通過設計算例，將本文方法與現有方法的設計結果進行了比較。結果表明，在滿足給定最大極點半徑的前提下，本文方法得到了具有更小幅頻響應誤差的濾波器。

1 濾波器表示和極點半徑約束

數字濾波器的設計問題是指選擇其傳遞函數的系數，使對應的頻率響應逼近某個期望響應。一個IIR數字濾波器可用多種形式的傳遞函數H（z）進行表示。如分子和分母均為單一高階多項式的傳遞函數：

式中，M和N分別是濾波器的分子階數和分母階數。

再如分子和分母均為級聯二階因子的傳遞函數：

式中，濾波器的階數N不妨假設為偶數。

相比式（1）和（2）的兩種濾波器表示，傳遞函數（1）關于其系數向量的非線性程度低，特別是，關于其分子多項式系數向量是純線性的。傳遞函數（2）則適合應用充分必要的穩(wěn)定三角形條件，對其最大極點半徑進行限制。

本文采用式（1）和（2）相結合的濾波器表示，即傳遞函數的分子表示為單一多項式，分母表示為級聯二階因子，即：

式中

分別為分母和分子的系數向量。

基于式（3）表示的傳遞函數，濾波器的極點就是二階分母多項式

的零點，其中an＝［an1，an2］T為An（z，an）的系數向量。設計時，要求濾波器的所有極點半徑都小于r＜1，即要求所有An（z，an）的零點半徑都小于r＜1。容易證明，這等價于要求：對任意n∈｛1，2，…，N/2｝成立。

在an平面上，式（5）確定的區(qū)域為圖1所示的三角形，稱為穩(wěn)定三角形［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）是一個高度非凸的優(yōu)化問題。一方面，由式（3）可見，濾波器的傳遞函數H（z，x）關于其分母系數向量a是高度非線性的；另一方面，當δ給定時，由式（8b）確定的任意通帶頻率點ω上的幅頻響應約束域是圖2所示頻率響應復平面上的圓環(huán)，是一個典型的非凸域。

2.2 凸優(yōu)化策略

針對以上切比雪夫設計問題的非凸性，首先根據高斯-牛頓策略，應用泰勒級數展開式，在系數向量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），就變成一個容易求解的二階錐規(guī)劃問題。可以證明，當?（ω）等于問題（7）之最優(yōu)濾波器的相頻響應?*（ω）時，即當?（ω）＝?*（ω）＝∠H（ejω，x*）時，該二階錐規(guī)劃問題的解，就是原設計問題（7）的解。

可惜的是，相位函數?*（ω）是未知的。為此，類似文獻［12］，本文采用迭代方法對?*（ω）進行估計。在第l次迭代，讓?（ω）等于某個?l（ω），然后用迭代中求得的濾波器相頻響應來更新?l（ω）。通過交替地更新?l（ω）和求解相應的二階錐規(guī)劃問題，可估計出?*（ω）并最終求得約束切比雪夫設計問題（7）的最優(yōu)濾波器。

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步 求解以下二階錐規(guī)劃問題

第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，停止算法（其中l(wèi)max為最大外循環(huán)次數，即更新相位函數?l＋1（ω）的迭代次數）。否則，令l＝l＋1并返回第2步。

3 設計算例及比較

下面通過設計算例及與現有方法的比較展示本文算法的優(yōu)越性。本文算法已編制成MATLAB函數iirmmc（M，N，f，d，r），其中M、N、f、d、r分別為濾波器的分子階數、分母階數、帶邊頻率集、帶邊期望幅值集和最大極點半徑。實驗在MATLAB 2017中進行，頻率區(qū)間［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）是得到的最大幅頻響應誤差隨外循環(huán)次數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）進行的設計，最大幅頻響應誤差隨算法外循環(huán)次數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階全帶微分器。文獻方法包括約束優(yōu)化方法（Constrained Optimization Method，COM）［8］、序貫部分優(yōu)化（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階全帶微分器，逼近精度高、群延遲低、全頻帶上的群延遲偏差也不大，可在今后的教學和科研中加以應用。