張登奇， 彭 鑫， 陳海蘭

（湖南理工學(xué)院 信息與通信工程學(xué)院， 湖南 岳陽 414006）

雙線性變換法在IIR濾波器設(shè)計中的應(yīng)用

張登奇， 彭 鑫， 陳海蘭

（湖南理工學(xué)院 信息與通信工程學(xué)院， 湖南 岳陽 414006）

數(shù)字濾波器的設(shè)計是數(shù)字信號處理課程的重要教學(xué)內(nèi)容， 而雙線性變換法是間接設(shè)計IIR數(shù)字濾波器的重要方法. 本文介紹了雙線性變換法的設(shè)計原理， 并以巴特沃斯數(shù)字低通濾波器為例， 歸納了雙線性變換法的設(shè)計步驟， 對設(shè)計過程中較難理解的過渡濾波器和頻率預(yù)畸變進(jìn)行了分析， 對變換常數(shù)和采樣周期的取值問題進(jìn)行了討論， 并以實例給出了MATLAB實現(xiàn)的程序.

數(shù)字濾波器； 雙線性變換法； 預(yù)畸變； MATLAB

引言

間接設(shè)計IIR數(shù)字濾波器的常用方法有脈沖響應(yīng)不變法和雙線性變換法. 文［1］對脈沖響應(yīng)不變法進(jìn)行了介紹和分析， 但由于存在頻譜混疊失真， 故不適合設(shè)計高通和帶阻濾波器. 雙線性變換法通過頻率的非線性變換， 解決了頻譜混疊問題， 得到了廣泛應(yīng)用. 但一般教材對其原理論述簡單， 對某些參數(shù)介紹不詳， 導(dǎo)致學(xué)習(xí)理解困難. 本文從脈沖響應(yīng)不變法的頻譜混疊入手， 先分析雙線性變換法的原理思想， 再以巴特沃斯數(shù)字低通濾波器的設(shè)計為例， 歸納雙線性變換法的設(shè)計步驟， 并對設(shè)計過程中較難理解的過渡濾波器和頻率預(yù)畸變進(jìn)行深入分析， 最后以實例給出MATLAB的實現(xiàn)程序.

1　雙線性變換法

1.1設(shè)計原理

圖1　雙線性變換過程示意圖

其中C是一個正實常數(shù)， 也稱變換常數(shù). 由于

再將s1平面通過映射到z平面， 得

1.2設(shè)計步驟

雙線性變換法可由（5）式直接將）（sHa變換成）（zH，但設(shè)計）（sHa的技術(shù)指標(biāo)并不好理解.下面以巴特沃斯低通濾波器的設(shè)計為例， 結(jié)合圖2所示設(shè)計指標(biāo)及幅頻特性的變化過程， 介紹雙線性變換法的設(shè)計步驟， 并對設(shè)計指標(biāo)進(jìn)行分析說明.

圖2　設(shè)計指標(biāo)及幅頻特性的變化過程

實例 對輸入的模擬信號進(jìn)行采樣得離散時間信號，設(shè)采樣頻率Fs=20kHz， 要求保留3000Hz以下的頻率成分，幅度失真小于2dB； 濾除5000Hz以上的頻率成分， 衰減大于10dB. 試用雙線性變換法設(shè)計數(shù)字濾波器實現(xiàn)上述要求.

（1） 確定采用巴特沃斯數(shù)字低通濾波器.數(shù)字指標(biāo)為ωp=0.3π rad ，ap=2dB （即0.7943倍）， ωs=0.5π rad，as=10dB （即0.3162倍）. 在用雙線性變換法設(shè)計數(shù)字濾波器的公式中， 幅度指標(biāo)要以分貝為單位代入數(shù)據(jù)， 但用分貝為單位圖示的幅頻曲線感覺衰減不明顯. 故在圖2（a）點狀圖所示的幅頻指標(biāo)中， 幅度指標(biāo)數(shù)據(jù)的單位改用倍數(shù).

（2） 將數(shù)字低通濾波器的技術(shù)指標(biāo)轉(zhuǎn)換成虛構(gòu)的過渡模擬濾波器技術(shù)指標(biāo). （5）式描述了s平面與z平面的映射關(guān)系， 令s=jΩ， z=ejω， 代入（5）式可得過渡模擬角頻率與數(shù)字頻率之間的非線線關(guān)系

設(shè)采樣周期T=2s（采樣周期T不參與設(shè)計， 取2s計算最簡）， 轉(zhuǎn)化后的過渡模擬濾波器技術(shù)指標(biāo)為Ωp=0.5095rad/s， Ap=2dB（即0.7943倍）， Ωs=1rad/s， As=10dB （即0.3162倍）. 過渡模擬濾波器的技術(shù)指標(biāo)參見圖2（c）中的點狀圖.

（3） 按照虛構(gòu)的過渡模擬濾波器技術(shù)指標(biāo)設(shè)計過渡模擬濾波器系統(tǒng)函數(shù)［5］. 過濾模擬濾波器的最低階次及3dB頻率計算公式為：

其中， ［x］表示不超過實數(shù)x的最大整數(shù). 代入過渡模擬指標(biāo)數(shù)據(jù)， 算得N=3， Ωcp=0.5572 rad/s ， Ωcs= 0.6934 rad/s. 根據(jù)最小階次N直接得到歸一化原型系統(tǒng)函數(shù)

這里取Ωc=Ωcp， 代入（8）式算得

過渡模擬低通濾波器的幅頻曲線如圖2（c）所示. 因取Ωc=Ωcp， 所以設(shè)計的過渡模擬濾波器通帶邊界指標(biāo)精確滿足， 而阻帶邊界指標(biāo)留有余量.

（4）用雙線性變換法將設(shè)計的過渡模擬濾波器系統(tǒng)函數(shù)轉(zhuǎn)變成數(shù)字濾波器系統(tǒng)函數(shù). 取T=2s， 將（5）式代入（9）式并化簡得

用雙線性變換法設(shè)計的數(shù)字低通濾波器幅頻曲線如圖2（a）所示， 與圖2（c）對比可知幅頻曲線的形狀差別較大， 這是由于雙線性變換法中頻率的非線性關(guān)系引起的.

2　設(shè)計過程的深入分析

2.1過渡模擬濾波器

間接設(shè)計IIR數(shù)字濾波器， 要先設(shè)計一個過渡模擬濾波器， 過渡模擬濾波器可根據(jù)數(shù)字指標(biāo)自定義一個采樣周期T進(jìn)行虛構(gòu). 后面的分析可知采樣周期T不參與設(shè)計， T值不同， 算得的過渡模擬指標(biāo)會不同，由此設(shè)計的過渡模擬濾波器也肯定不同， 但按同一T值進(jìn)行轉(zhuǎn)換得到的數(shù)字濾波器是一樣的. 如果數(shù)字濾波器的技術(shù)指標(biāo)是以數(shù)字指標(biāo)形式給出， 則無需再給出采樣周期或采樣頻率； 若是以模擬指標(biāo)形式給出，則一定要注意設(shè)計指標(biāo)中用于離散連續(xù)信號的采樣周期sT和用于設(shè)計過渡模擬濾波器的采樣周期T不是同一個概念， 前者用于計算數(shù)字指標(biāo)， 后者用于虛構(gòu)模擬指標(biāo). 一些教材中混用了這兩個符號， 有的還直接用T=Ts設(shè)計過渡模擬濾波器， 理論上雖然可行， 但學(xué)生可能會認(rèn)為過渡模擬濾波器必須這樣設(shè)計， 既誤導(dǎo)了學(xué)生， 也增加了計算難度， 甚至還會引起較大誤差， 設(shè)計中應(yīng)盡量避免.

2.2頻率預(yù)畸變校正

進(jìn)行非線性變換， 則設(shè)計的Ha（s）按（4）式變換得到的H（z）一定滿足邊界頻率要求. 圖2（c）給出了通帶指標(biāo)精確滿足的模擬濾波器幅頻特性曲線， 圖2（a）給出了變換后數(shù)字濾波器的幅頻特性曲線， 其通帶指標(biāo)也精確滿足. 假如用Ωp=ωpT ， Ωs=ωsT 作為過渡模擬濾波器的邊界頻率， 則設(shè)計的Ha（s）按（4）式變換得到的H（z）就不會滿足邊界頻率要求. （11）式中的Ω與ω的非線性關(guān)系， 稱為頻率預(yù)畸變校正.

2.3變換常數(shù)C對設(shè)計的影響

2.4采樣周期T對設(shè)計的影響

在雙線性變換法的設(shè)計過程中有一個參數(shù)T， 其實設(shè)計過程中的參數(shù)T只影響變換常數(shù)C， 而C的取值不影響設(shè)計， 所以參數(shù)T也不影響設(shè)計. 當(dāng)然， 采樣周期T不影響設(shè)計還應(yīng)從圖2所示的坐標(biāo)變換中進(jìn)行理解. 對任一給定的變換常數(shù)C， 采樣周期T對從ω到Ω1的影響和從Ω1到Ω的影響可以相互抵消. 也就是說， C一定時， 采樣周期T會影響圖2（a）到圖2（b）的坐標(biāo)變換， 也會影響圖2（b）到圖2（c）的坐標(biāo)變換，這兩種影響正好使圖2（a）到圖2（c）的坐標(biāo)變換與T無關(guān). 理論上， 采樣周期T對設(shè)計結(jié)果沒有影響， 但由于在設(shè)計過程中T參與了計算， 考慮到計算機(jī)的有限字長效應(yīng)和計算簡單， 一般取T=2s.

3　MATLAB實現(xiàn)程序

用雙線性變換法設(shè)計實例中所要求的濾波器程序如下：

%技術(shù)指標(biāo)： Fs=20kHz，fp=3kHz，fs=5kHz，Ap=2dB，As=10db，

clc； close all； clear； format compact；

wp=0.3*pi， ws=0.5*pi， ap=2， as=10， Ts=1/20000；%輸入數(shù)字技術(shù)指標(biāo)

T=2，%改變T值不影響設(shè)計結(jié)果

Wp=（2/T）*tan（wp/2）， Ws=（2/T）*tan（ws/2）， Ap=ap， As=as， %雙線性變換法指標(biāo)轉(zhuǎn)換

Mn=log10（（10^（0.1*Ap）-1）/（10^（0.1*As）-1）），Md=（2*log10（Wp/Ws）），%先算分子分母

N=ceil（Mn/Md），%計算過渡模擬濾波器最小階次

Wcpd=（10^（0.1*Ap）-1）^（1/（2*N）），Wcsd=（10^（0.1*As）-1）^（1/（2*N））， %先算分母

Wcp=Wp/Wcpd； Wcs=Ws/Wcsd； %計算截止頻率

Wc=Wcp，%通帶精確滿足，取Wc=Wcs時阻帶精確滿足

［bp1，ap1］ = butter（N，1，'s'）； %%計算歸一化低通原型模擬濾波器的系統(tǒng)函數(shù)Ga（p）

tf（bp1，ap1， 'variable'， 'p'）， %%顯示歸一化低通原型模擬濾波器的系統(tǒng)函數(shù)

［bs，as］ = lp2lp（bp1，ap1，Wc）； %%計算過渡模擬低通濾波器的系統(tǒng)函數(shù)Ha（s）

tf（bs，as）， %%顯示過渡模擬濾波器的系統(tǒng)函數(shù)

［bz，az］ = bilinear （bs，as，1/T），%將過渡模擬濾波器轉(zhuǎn)換為數(shù)字濾波器

tf（bz，az，Ts，'variable'，'z^-1'），%顯示設(shè)計的數(shù)字濾波器系統(tǒng)函數(shù)

freqz（bz，az，200），%繪圖并驗證設(shè)計技術(shù)指標(biāo)

該程序的運行結(jié)果在MATLAB的命令窗或圖形窗中直接顯示. 從圖形窗中顯示的幅頻曲線可以直接測出ω=0.3π時衰減為2dB， ω=0.5π時衰減為15.37dB， 通帶指標(biāo)精確滿足， 阻帶指標(biāo)富有余量. 如果改用Ωc=Ωcs重新設(shè)計， 測出ω=0.3π時衰減0.6352dB， ω=0.5π時衰減10dB， 通帶指標(biāo)富有余量， 阻帶指標(biāo)精確滿足.改變采樣周期T， 設(shè)計結(jié)果相同.

4　結(jié)束語

［1］ 張登奇， 彭 鑫. 脈沖響應(yīng)不變法在IIR濾波器設(shè)計中的應(yīng)用［J］. 湖南理工學(xué)院學(xué)報（自然科學(xué)版）. 2015（04）： 23～26

［2］ 吳鎮(zhèn)揚(yáng). 數(shù)字信號處理的原理與實現(xiàn)［M］. 南京： 東南大學(xué)出版社， 2001： 77

［3］ 李芬華. 雙線性變換中的映射［J］. 河北大學(xué)學(xué)報（自然科學(xué)版）. 1991（02）： 65～71

［4］ 吳湘淇. 信號、系統(tǒng)與信號處理（下）［M］. 北京： 電子工業(yè)出版社，1996： 60～67

［5］ 劉順蘭， 吳 杰. 數(shù)字信號處理［M］. 西安： 西安電子科技大學(xué)出版社， 2003： 158

［6］ 高西全， 丁玉美. 數(shù)字信號處理［M］. 第3版. 西安： 西安電子科技大學(xué)出版社， 2008： 184～189

Application of the Bilinear Transformation in IIR Digital Filter Design

ZHANG Deng-qi， PENG Xin， CHEN Hai-lan
（College of Information and Communication Engineering， Hunan Institute of Science and Technology， Yueyang 414006， China）

Digital filter design is the most important component in teaching of digital signal processing. Bilinear Transformation Method is an important method in the indirect design of IIR filter. We deduce the procedure of IIR filter design based on the fundamentals of Bilinear Transformation. And taking Butterworth low pass filter for instance， we introduce the design steps of filter. Furthermore， we deeply analyze that the unintelligible steps in the design procedure including transitional filters and frequency prewarping. And then， the value of transform constant and sampling period is discussed. Finally， the MATLAB source code of the filter design is given.

digital filter； bilinear Transformation； frequency prewarping transformation； MATLAB

TN713.7

A

1672-5298（2016）03-0021-05

2016-04-28

湖南省自然科學(xué)基金項目（14JJ3130）； 湖南省教育廳重點項目（12A057）； 湖南理工學(xué)院教學(xué)改革研究項目（2014B20）

張登奇（1968- ）， 男， 湖南臨湘人， 碩士， 湖南理工學(xué)院信息與通信工程學(xué)院副教授. 主要研究方向： 信號與信息處理