《MATLAB輔助現(xiàn)代工程數(shù)字信號處理》課件第4章

第4章工程數(shù)字濾波器設(shè)計4.1數(shù)字濾波器的技術(shù)指標與設(shè)計方法4.2無限沖激響應(yīng)(IIR)數(shù)字濾波器的設(shè)計4.3有限沖激響應(yīng)(FIR)數(shù)字濾波器的設(shè)計4.4小結(jié)

4.1數(shù)字濾波器的技術(shù)指標與設(shè)計方法

4.1.1數(shù)字濾波器的技術(shù)指標

物理上可實現(xiàn)的線性時不變系統(tǒng)是因果的，具有頻率響應(yīng)(4.1)物理上可實現(xiàn)的低通濾波器的幅頻特性如圖4.1所示，由通帶、過渡帶和阻帶三部分組成，位于通帶頻率范圍的信號可以無損通過。圖4.1物理上可實現(xiàn)濾波器的幅度特性從圖4.1可知，數(shù)字低通濾波器的技術(shù)指標主要有:

(1)通帶截止頻率ωp;

(2)阻帶截止頻率ωs;

(3)通帶波紋δ1;

(4)阻帶波紋δ2;

(5)通帶內(nèi)允許的最大衰減αp:

αp=－20lg(1－δ1)=－20lg|H(ejωp)|

(6)阻帶內(nèi)允許的最小衰減αs:

αs=－20lg(1－δ2)=－20lg|H(ejωs)|

(7)3dB通帶截止頻率ωc:當幅度|H(ejω)|下降到0.707(即3dB)時對應(yīng)的頻率。圖4.1物理上可實現(xiàn)濾波器的幅度特性4.1.2數(shù)字濾波器的設(shè)計方法

數(shù)字濾波器的一般設(shè)計步驟如下:

(1)按照實際需要確定濾波器的性能要求;

(2)用一個因果穩(wěn)定系統(tǒng)的H(z)或h(n)去逼近該性能要求，即求h(n)的表達式，再確定系數(shù)ai、bi或零、極點ci、di，以使濾波器滿足給定的性能要求;

(3)用一個有限精度的運算去實現(xiàn)這個系統(tǒng)函數(shù)，包括選擇運算結(jié)構(gòu)(如級聯(lián)型、并聯(lián)型、卷積型、頻率采樣型、快速卷積(FFT)型等);

(4)選擇合適的字長和有效的數(shù)字處理方法，用數(shù)字電路或在計算機上編寫軟件實現(xiàn)所設(shè)計的H(z)。4.2無限沖激響應(yīng)(IIR)數(shù)字濾波器的設(shè)計

4.2.1IIR濾波器的經(jīng)典設(shè)計

從模擬濾波器設(shè)計IIR數(shù)字濾波器，就是按照一定的轉(zhuǎn)換關(guān)系將s平面上的Ha(s)轉(zhuǎn)換成z平面上的H(z)。這種映射變換應(yīng)遵循如下兩個基本原則:

(1)H(z)的頻響要能模仿Ha(s)的頻響，即s平面的虛軸應(yīng)映射到z平面的單位圓ejω上。

(2)Ha(s)的因果穩(wěn)定性映射成H(z)后保持不變，即s平面的左半平面應(yīng)映射到z平面的單位圓以內(nèi)。

IIR數(shù)字濾波器經(jīng)典設(shè)計法的步驟如下:

(1)將給定的數(shù)字濾波器的技術(shù)指標轉(zhuǎn)換為模擬濾波器的技術(shù)指標;

(2)估計模擬濾波器的最小階數(shù)和邊界條件，可利用MATLAB的工具函數(shù)buttord、cheb1ord等;

(3)根據(jù)轉(zhuǎn)換后的技術(shù)指標設(shè)計模擬低通濾波器原型，可利用MATLAB的工具函數(shù)buttap、cheb1ap、ellipap等;

(4)由頻率變換將模擬低通濾波器原型轉(zhuǎn)化為所需的模擬濾波器(低通、高通、帶通等)，可利用MATLAB

的工具函數(shù)lp2lp、lp2hp、lp2bp等;

(5)按照一定規(guī)則將模擬濾波器轉(zhuǎn)換為數(shù)字濾波器。

1.脈沖響應(yīng)不變法

脈沖響應(yīng)不變法是從濾波器的脈沖響應(yīng)出發(fā)，使數(shù)字濾波器的單位脈沖響應(yīng)序列h(n)正好等于模擬濾波器的沖激響應(yīng)ha(t)的采樣值，即

h(n)=ha(t)|t=nT

(4.2)

式中，T為采樣周期。

因此，數(shù)字濾波器的系統(tǒng)函數(shù)H(z)可由下式求得:

H(z)=Z[h(n)]=Z[ha(nT)]

(4.3)

采用脈沖響應(yīng)不變法將模擬濾波器變換為數(shù)字濾波器時，它所完成的s平面到z平面的變換，正是拉普拉斯變換到Z變換的標準變換關(guān)系z=esT。這種映射關(guān)系反映的是Ha(s)的周期延拓與H(z)的關(guān)系，而不是Ha(s)本身與H(z)的關(guān)系。脈沖響應(yīng)不變法的數(shù)字角頻率ω和模擬角頻率Ω滿足線性變換關(guān)系：ω=ΩT。

在MATLAB中，可用函數(shù)impinvar來實現(xiàn)脈沖響應(yīng)不變法，調(diào)用格式為

[bz，az]=impinvar(b，a，fs)

[bz，az]=impinvar(b，a)

其中，b和a分別為模擬濾波器的分子和分母多項式系數(shù)向量;fs為采樣頻率，缺省值為1Hz;bz和az分別為數(shù)字濾波器的分子和分母多項式系數(shù)向量。

【例4.1】

用脈沖響應(yīng)不變法設(shè)計一個巴特沃斯低通數(shù)字濾波器，設(shè)計指標為ωp=0.2π,ωs=0.3π,rp=3dB,rs=15dB，采樣頻率fs=10kHz。

MATLAB程序如下:

%MATLABPROGRAM4-1

wp=0.2*pi*10000;%數(shù)字指標轉(zhuǎn)化為模擬指標

ws=0.3*pi*10000;

rp=3;

rs=15;

fs=10000;

Nn=256;

[n,wn]=buttord(wp,ws,rp,rs,′s′);%計算階數(shù)n和截止頻率wn

[z，p，k]=buttap(n);

%設(shè)計模擬低通濾波器原型

[bap，aap]=zp2tf(z，p，k);%傳遞函數(shù)形式

[b，a]=lp2lp(bap，aap，wn);%低通到低通的頻率變換

[bz，az]=impinvar(b，a，fs);%數(shù)字濾波器變換

freqz(bz，az，Nn，fs);

程序運行結(jié)果如圖4.2所示。圖4.2基于脈沖響應(yīng)不變法的數(shù)字濾波器的頻響圖

2.雙線性變換法

從s平面到z平面的雙線性變換為

(4.4)(4.5)式中，T為采樣周期。雙線性變換式(4.4)和式(4.5)具有一一對應(yīng)的映射關(guān)系，它將s平面上的一點映射為z平面上的一點，或?qū)平面上的一點映射為s平面上的一點。系統(tǒng)函數(shù)H(z)和原型模擬傳輸函數(shù)Ha(s)之間的關(guān)系為(4.6)

1)零極點增益形式

[zd，pd，kd]=bilinear(z，p，k，fs)

[zd，pd，kd]=bilinear(z，p，k，fs，fp)

其中，z、p、k和zd、pd、kd分別為s域和z域系統(tǒng)函數(shù)的零點、極點和增益，fs為采樣頻率，fp為預(yù)畸變頻率。

2)傳遞函數(shù)形式

[numd，dend]=bilinear(mun，den，fs)

[numd，dend]=bilinear(mun，den，fs，fp)

其中，numd，dend和num，den分別為數(shù)字濾波器和模擬濾波器傳遞函數(shù)的分子和分母多項式系數(shù)向量，fs和fp意義同上。

【例4.2】

設(shè)計一個巴特沃斯低通數(shù)字濾波器，其通帶截止頻率為0.2π，阻帶邊界頻率為0.3π，通帶波紋小于3dB，阻帶衰減大于15dB，采樣頻率為10kHz。

MATLAB程序如下:

%MATLABPROGRAM4-2

wp=0.2*pi;

ws=0.3*pi;

rp=3;rs=15;

fs=10000;

Nn=128;

wp=2*10000*tan(wp/2);%數(shù)字指標轉(zhuǎn)化為模擬指標

ws=2*10000*tan(ws/2);

[n，wn]=buttord(wp，ws，rp，rs，′s′);

%計算階數(shù)n和截止頻率wn

[z，p，k]=buttap(n);%設(shè)計模擬低通濾波器原型

[bap，aap]=zp2tf(z，p，k);

[b，a]=lp2lp(bap，aap，wn);

[bz，az]=bilinear(b，a，fs);

freqz(bz，az，Nn，fs);

程序運行結(jié)果如圖4.3所示。圖4.3基于雙線性變換法的數(shù)字濾波器的頻響圖

3.最小階數(shù)選擇函數(shù)

MATLAB信號處理工具箱提供了工具函數(shù)來計算各類濾波器的最小階數(shù)和截止頻率。

(1)計算巴特沃斯濾波器最小階數(shù)可選擇函數(shù)butter，調(diào)用格式為

[n，wn]=butter(wp，ws，rp，rs，′s′)

[n，wn]=butter(wp，ws，rp，rs)

(2)計算切比雪夫Ⅰ型濾波器最小階數(shù)可選擇函數(shù)cheb1ord，調(diào)用格式為

[n，wn]=cheb1ord(wp，ws，rp，rs)

[n，wn]=cheb1ord(wp，ws，rp，rs)

(3)計算切比雪夫Ⅱ型濾波器最小階數(shù)可選擇函數(shù)cheb2ord，調(diào)用格式為

[n，wn]=cheb2ord(wp，ws，rp，rs，′s′)

[n，wn]=cheb2ord(wp，ws，rp，rs)

(4)計算橢圓濾波器最小階數(shù)可選擇函數(shù)ellipord，調(diào)用格式為

[n，wn]=ellipord(wp，ws，rp，rs，′s′)

[n，wn]=ellipord(wp，ws，rp，rs)

其中，wp為通帶邊界頻率;ws為阻帶邊界頻率;rp為通帶波紋;rs為阻帶衰減;′s′表示模擬濾波器，缺省時表示數(shù)字濾波器;n為最小階數(shù);wn為截止頻率。

4.頻率變換

1)低通變換

在模擬濾波器的低通設(shè)計中，利用完全設(shè)計函數(shù)來設(shè)計模擬低通濾波器，利用雙線性變換來得到數(shù)字低通濾波器。

【例4.3】

設(shè)采樣周期T=250μs，設(shè)計一個三階巴特沃斯低通數(shù)字濾波器，其3dB截止頻率fc=1kHz。

解直接按Ωc=2πfc設(shè)計模擬濾波器Ha(s)。以截止頻率Ωc歸一化的三階巴特沃斯濾波器的傳遞函數(shù)為MATLAB程序如下:

%MATLABPROGRAM4-3

N=3;fs=4000;

wn=2*pi*1000;

[B，A]=butter(N，wn，′s′);

[num，den]=bilinear(B，A，fs);

[h，w]=freqz(num，den);

f=w/pi*2000;%轉(zhuǎn)化為實際采樣頻率(Hz)

plot(f，20*log10(abs(h)));

xlabel(′f/Hz′);ylabel(′幅值/dB′);grid;

程序運行結(jié)果如圖4.4所示。圖4.4三階巴特沃斯濾波器的幅頻特性曲線

2)高通變換

在模擬濾波器的高通設(shè)計中，低通至高通的變換就是s變量的倒置，這一關(guān)系同樣可應(yīng)用于雙線性變換，只要將變換式中的s代之以1/s，就可得到數(shù)字高通濾波器。由于倒數(shù)關(guān)系不改變模擬濾波器的穩(wěn)定性，因此也不會影響雙線性變換后的穩(wěn)定條件，而且jΩ軸仍映射在單位圓上。

【例4.4】

設(shè)計一個切比雪夫高通濾波器，它的通帶為400～500Hz，通帶內(nèi)允許有0.5dB的波動，阻帶衰減在小于317Hz的頻帶內(nèi)至少為19dB，采樣頻率為1kHz。解

MATLAB程序如下:

%MATLABPROGRAM4-4

wp=2*1000*tan(2*pi*400/(2*1000));%數(shù)字指標轉(zhuǎn)化為模擬指標

ws=2*1000*tan(2*pi*317/(2*1000));

rp=0.5;rs=19;

[N，wn]=cheb1ord(wp，ws，rp，rs，′s′);

[B，A]=cheby1(N，0.5，wn，′high′，′s′);

[num，den]=bilinear(B，A，1000);

[h，w]=freqz(num，den);

f=w/pi*500;

plot(f，20*log10(abs(h)));

xlabel(′f/Hz′);ylabel(′幅度/dB′);grid;

程序運行結(jié)果如圖4.5所示。圖4.5切比雪夫高通濾波器的幅頻特性曲線

3)帶通變換

若數(shù)字頻域上帶通濾波器的中心頻率為ω0，則帶通變換是將s的原點映射到z=e±jω0，而將s=±j∞點映射到z=±1，滿足這一要求的雙線性變換為(4.8)帶通變換關(guān)系為(4.9)設(shè)計帶通濾波器時，一般只給出上下邊帶的截止頻率ω1、ω2作為設(shè)計要求。為了實現(xiàn)帶通變換，先需將上下邊帶參數(shù)ω1、ω2換算成中心頻率ω0及模擬低通截止頻率Ωc。為此將ω1和ω2代入變換關(guān)系式(4.9)，并由Ω1和Ω2在模擬低通濾波器中是一對鏡像頻率(Ω1=－Ω2)可得到

(4.10)又因Ω1就是模擬低通濾波器的截止頻率Ωc，故有(4.11)

【例4.5】

設(shè)計一巴特沃斯帶通濾波器，采樣頻率fs=400kHz,其3dB邊界頻率分別為f1=110kHz,f2=90kHz，在阻帶f3=120kHz處最小衰減大于10dB。

解確定數(shù)字頻域的上下邊帶的角頻率求得中心頻率求得模擬低通濾波器的通帶截止頻率Ωc與阻帶邊界頻率Ωs為

由計算結(jié)果可知，從Ωc到Ωs頻率增加了約1.05倍，衰減增加了3~10dB，故選用二階巴特沃斯濾波器可滿足指標。其對應(yīng)的歸一化系統(tǒng)函數(shù)為將Ωc代入上式，有代入變換公式，有MATLAB程序如下:

%MATLABPROGRAM4-5

w1=2*400*tan(2*pi*90/(2*400));

w2=2*400*tan(2*pi*110/(2*400));

w4=2*400*tan(2*pi*120/(2*400));

w3=2*400*tan(2*pi*(90-10)/(2*400));

wp=[w1，w2];ws=[w3，w4];

rp=3;rs=10;

[N，wn]=buttord(wp，ws，rp，rs，′s′);

[B，A]=butter(N，wn，′s′);

[num，den]=bilinear(B，A，400);

[h，w]=freqz(num，den);

f=w/pi*200;

plot(f，20*log10(abs(h)));

axis([40，160，-30，10]);xlabel(′f/kHz′);ylabel(′幅度/dB′);grid;

程序運行結(jié)果如圖4.6所示。圖4.6巴特沃斯帶通濾波器的幅頻特性曲線

4)帶阻變換

將帶通濾波器的頻率關(guān)系倒置就得到帶阻變換，即(4.12)

5.IIR數(shù)字濾波器的完全設(shè)計函數(shù)

1)函數(shù)butter

函數(shù)butter用于巴特沃斯濾波器的設(shè)計，調(diào)用格式為

[b，a]=butter(n，wn)

[b，a]=butter(n，wn，′ftype′)

[z，p，k]=butter(…)

其中，b、a為所要求的系統(tǒng)函數(shù)多項式的分子和分母系數(shù)向量;n為濾波器的階數(shù);wn為濾波器的截止頻率?！鋐type′為濾波器的類型:若缺省則為低通或帶通濾波器;′high′為高通濾波器;′stop′為帶阻濾波器，截止頻率wn=[w1，w2]。

【例4.6】

設(shè)計一巴特沃斯帶阻濾波器，采樣頻率fs=1kHz，要濾除100Hz的干擾，其3dB的邊界頻率為95Hz和105Hz，原型歸一化低通濾波器為Ha(s)=。

MATLAB程序如下:

%MATLABPROGRAM4-6

fs=1000;

w1=2*95/fs;w2=2*105/fs;

wn=[w1，w2];

[B，A]=butter(1，wn，′stop′);

[h，w]=freqz(B，A);

f=w/pi*500;

plot(f，20*log10(abs(h)));

axis([50，150，-30，10]);

xlabel(′f/Hz′);ylabel(′幅度/dB′);grid;

程序運行結(jié)果如圖4.7所示。圖4.7巴特沃斯帶阻濾波器的幅頻特性曲線

2)函數(shù)cheby1

函數(shù)cheby1用于切比雪夫Ⅰ型濾波器的設(shè)計，調(diào)用格式為

[b，a]=cheby1(n，rp，wn)

[b，a]=cheby1(n，rp，wn，′ftype′)

[z，p，k]=cheby1(…)

其中，n為濾波器階數(shù)，rp為通帶波紋，其余各項同函數(shù)butter。

3)函數(shù)cheby2

函數(shù)cheby2用于切比雪夫Ⅱ型濾波器的設(shè)計，調(diào)用格式為

[b，a]=cheby2(n，rs，wn)

[b，a]=cheby2(n，rs，wn，′ftype′)

[z，p，k]=cheby2(…)

其中，n為濾波器階數(shù)，rs為阻帶衰減，其余各項同函數(shù)butter。

【例4.7】

設(shè)計一個帶通切比雪夫Ⅱ型數(shù)字濾波器，通帶頻率為150～300Hz，過渡帶寬為50Hz，通帶波紋小于1dB，阻帶衰減為20dB，采樣頻率為1kHz。

MATLAB程序如下:

%MATLABPROGRAM4-7

fs=1000;

wp=[150，300]*2/fs;

%數(shù)字指標標準化

rp=1;rs=20;

ws=[(150-50)，(300+50)]*2/fs;

[n，wn]=cheb2ord(wp，ws，rp，rs);%計算最小階數(shù)和截止頻率

[b，a]=cheby2(n，rs，wn);

[h，w]=freqz(b，a);

f=w/pi*500;

plot(f，20*log10(abs(h)));

axis([0，500，－100，2]);

xlabel(′f/Hz′);ylabel(′幅度/dB′);grid;

程序運行結(jié)果如圖4.8所示。圖4.8切比雪夫Ⅱ型帶通數(shù)字濾波器的幅頻特性曲線

4)函數(shù)ellip

函數(shù)ellip用于橢圓濾波器的設(shè)計，調(diào)用格式為

[b，a]=ellip(n，rp，rs，wn)

[b，a]=ellip(n，rp，rs，wn，′ftype′)

[z，p，k]=ellip(…)

其中，n為濾波器階數(shù)，rp為通帶波紋，rs為阻帶衰減，wn為截止頻率，其余各項同函數(shù)butter。

【例4.8】

設(shè)計一個橢圓帶通數(shù)字濾波器，通帶頻率為150～400Hz，通帶波紋小于3dB，阻帶衰減為40dB，兩邊過渡帶寬為50Hz，采樣頻率為2kHz。MATLAB程序如下:

%MATLABPROGRAM4-8

fs=2000;

wp=[150，400]*2/fs;

%數(shù)字指標標準化

rp=3;rs=40;

ws=[(150-50)，(400+50)]*2/fs;

[n，wn]=ellipord(wp，ws，rp，rs);%計算最小階數(shù)和截止頻率

[b，a]=ellip(n，rp，rs，wn);

[h，w]=freqz(b，a);

f=w/pi*500;

plot(f，20*log10(abs(h)));

xlabel(′f/Hz′);ylabel(′幅度/dB′);grid;

程序運行結(jié)果如圖4.9所示。圖4.9橢圓帶通數(shù)字濾波器的幅頻特性曲線4.2.2IIR濾波器的直接設(shè)計

MATLAB信號處理工具箱提供了函數(shù)yulewalk直接設(shè)計IIR數(shù)字濾波器，調(diào)用格式為

[b，a]=yulewalk(n，f，m)

其中，n為濾波器的階數(shù);f為給定濾波器的頻率點向量，標準化頻率取值范圍為0～1，f的第一個頻率點必須是0，最后一個頻率點必須是1，且f的頻率點必須是遞增的;m為與頻率向量f對應(yīng)的理想幅值響應(yīng)向量，m和f必須是相同維數(shù)的向量;b和a分別是濾波器系統(tǒng)函數(shù)的分子多項式和分母多項

式系數(shù)向量。函數(shù)yulewalk采用以下步驟計算分子多項式:

(1)計算與分子多項式相應(yīng)的幅值平方響應(yīng)的輔助式;

(2)由輔助分子式和分母多項式計算完全的頻率響應(yīng);

(3)計算濾波器的脈沖響應(yīng);

(4)采用最小二乘法擬合脈沖響應(yīng)，求取濾波器的多項式系數(shù)。

函數(shù)yulewalk允許自定義f和m，因此可以設(shè)計出任意形狀的幅頻響應(yīng)的濾波器，但它不能設(shè)計給定相位指標的濾波器。

【例4.9】

用直接法設(shè)計一個多頻帶數(shù)字濾波器，幅頻響應(yīng)值如下:

f=[00.10.20.30.40.50.60.70.80.91];

m=[0110011100]

MATLAB程序如下:

%MATLABPROGRAM4-9

n=10;

f=[0:0.1:1];

m=[0，1，1，0，0，1，1，1，0，0，0];

[b，a]=yulewalk(n，f，m);

[h，w]=freqz(b，a，256);

axes(′position′，[0.2，0.2，0.4，0.4]);

plot(f，m，′b-′，w/pi，abs(h)，′r:′);

xlabel(′頻率(π)′);ylabel(′幅值′);

title(′DirectIIRFilterDesign-Yulewalk′);

legend(′b－′，′理想波形′，′r:′，′實際波形′);grid;

程序運行結(jié)果如圖4.10所示。圖4.10直接法設(shè)計的多頻帶濾波器幅頻特性曲線4.2.3最大平滑IIR數(shù)字濾波器設(shè)計

MATLAB信號處理工具箱提供了函數(shù)maxflat來設(shè)計最大平滑數(shù)字濾波器，調(diào)用格式為

[b，a]=maxflat(nb，na，wn)

b=maxflat(nb，′sym′，wn)

[b，a，b1，b2]=maxflat(nb，na，wn)

[b，a，b1，b2]=maxflat(nb，na，wn，′design-flag′)

其中，nb和na分別為濾波器分子和分母多項式的系數(shù);wn為濾波器-3dB的截止頻率，取值范圍為0～1;′sym′表示對稱型FIR巴特沃斯濾波器;b為濾波器分子多項式系數(shù)向量;b1為多項式系數(shù)向量，包含全部z=－1的零點;b2為多項式系數(shù)向量，包含除－1外的其余全部零點，b=b1*b2?！鋎esign-flag′為監(jiān)測設(shè)計過程標志:trace為表格顯示設(shè)計過程;plot為繪制濾波器的幅頻圖，群延時和零、極點圖。

【例4.10】

設(shè)計通用巴特沃斯低通數(shù)字濾波器，其系統(tǒng)函數(shù)分子階數(shù)為6，分母階數(shù)為3，截止頻率為0.7π。

MATLAB程序如下:

%MATLABPROGRAM4-10

nb=6;

na=3;

wn=0.7*pi;

[b，a]=maxflat(nb，na，wn，′plot′);

程序運行結(jié)果如圖4.11所示。圖4.11最大平滑巴特沃斯低通數(shù)字濾波器特性曲線4.3有限沖激響應(yīng)(FIR)數(shù)字濾波器的設(shè)計

MATLAB信號處理工具箱提供了FIR數(shù)字濾波器的設(shè)計函數(shù)，采用的設(shè)計方法和主要函數(shù)如表4.1所示。4.3.1窗函數(shù)法設(shè)計FIR濾波器

窗函數(shù)法是根據(jù)給定的性能指標，通過對系統(tǒng)的單位脈沖響應(yīng)h(n)加窗來逼近理想的單位脈沖響應(yīng)序列hd(n)，即(4.13)一般來說，理想頻響Hd(ejω)是分段恒定的，在邊界頻率處有突變點。因此，hd(n)往往都是無限長序列，而且是非因果的。為了使系統(tǒng)變?yōu)槲锢砜蓪崿F(xiàn)的FIR，最簡單的辦法是直接截取一段hd(n)代替h(n)，這種截取可以形象地想象為h(n)是通過一個“窗口”所看到的一段hd(n)，因此，h(n)也可表達為hd(n)和一個“窗函數(shù)”w(n)的乘積，即h(n)=hd(n)w(n)

(4.14)

1.矩形窗RN(n)

矩形窗可參考第1章1.3.1節(jié)矩形序列式(1.12)的內(nèi)容。

w(n)=RN(n)

(4.15a)(4.15b)

2.漢寧(Hanning)窗

漢寧窗又稱升余弦窗。(4.16a)利用傅里葉變換的移位特性，漢寧窗頻譜的幅度函數(shù)W(ω)可用矩形窗的幅度函數(shù)表示:(4.16b)三部分矩形窗頻譜相加，使旁瓣互相抵消，能量集中在主瓣，旁瓣大大減小，主瓣寬度增加1倍，為8π/N。

3.漢明(Hamming)窗

漢明窗又稱改進的升余弦窗。(4.17a)其幅頻響應(yīng)的幅度函數(shù)為(4.17b)漢明窗是對漢寧窗的改進，在主瓣寬度(對應(yīng)第一零點的寬度)相同的情況下，旁瓣進一步減小，可使99.96%的能量集中在窗譜的主瓣內(nèi)。

4.布萊克曼(Blackman)窗(4.18a)其幅頻響應(yīng)的幅度函數(shù)為(4.18b)布萊克曼窗增加一個2次諧波余弦分量，可進一步降低旁瓣，但主瓣寬度進一步增加為12π/N。增加N可減少過渡帶。

5.凱塞(Kaiser)窗

以上四種窗函數(shù)，都是以增加主瓣寬度為代價來降低旁瓣的。凱塞窗則可自由選擇主瓣寬度和旁瓣衰減，具有較強的適應(yīng)性。(4.19)

Io(β)是零階修正貝塞爾函數(shù)，參數(shù)β可自由選擇，決定主瓣寬度與旁瓣衰減。β越大，w(n)窗越窄，其頻譜的主瓣變寬，旁瓣變小。一般取4<β<9。當β=5.44時，接近漢明窗;當β=8.5時，接近布萊克曼窗;當β=0時，為矩形窗。用窗函數(shù)法設(shè)計FIR濾波器時，要根據(jù)給定的濾波器性能指標選擇窗口寬度N和窗函數(shù)w(n)。各種窗函數(shù)的性能如表4.2所示。

【例4.11】

用凱塞窗設(shè)計一FIR低通濾波器，低通邊界頻率ωc=0.3π，阻帶邊界頻率ωr=0.5π，阻帶衰減rs小于50dB。

解首先求解hd(n)，根據(jù)指標要求，其邊界頻率應(yīng)為MATLAB程序如下:

%MATLABPROGRAM4-11

wn=kaiser(30，4.55);

nn=[0:1:29];

alfa=(30－1)/2;

hd=sin(0.4*pi*(nn-alfa))./(pi*(nn－alfa));

h=hd.*wn′;

[h1，w1]=freqz(h，1);

plot(w1/pi，20*log10(abs(h1)));

axis([0，1，－80，10]);

xlabel(′歸一化頻率/pi′);ylabel(′幅度/dB′);grid;

程序運行結(jié)果如圖4.12所示。圖4.12基于凱塞窗的FIR數(shù)字濾波器幅頻特性曲線

1)標準型設(shè)計函數(shù)fir1

函數(shù)fir1用于設(shè)計標準的低通、帶通、高通和帶阻FIR濾波器。函數(shù)fir1的調(diào)用格式為

b=fir1(n，wc，′ftype′，windows)

其中，n為濾波器階數(shù);wc為截止頻率，對于帶通、帶阻濾波器，wc=[w1，w2]。′ftype′為濾波器類型:′high′為高通FIR濾波器，′stop′為帶阻FIR濾波器，缺省時為低通或帶通濾波器。windows指定窗函數(shù)類型，默認為漢明窗，可選漢寧窗、漢明窗、布萊克曼窗、三角窗(巴特利特窗)和矩形窗，每種窗都可以由MATLAB的相應(yīng)函數(shù)生成。b為FIR濾波器系數(shù)向量，長度為(n+1)。用fir1設(shè)計的FIR濾波器的群延時為n/2。

【例4.12】

用漢明窗設(shè)計一個線性相位的FIR高通濾波器，通帶邊界頻率為0.6π，阻帶邊界頻率為0.5π，阻帶衰減不小于50dB，通帶波紋不大于1dB。

MATLAB程序如下:

%MATLABPROGRAM4-12

wp=0.6*pi;ws=0.5*pi;

wd=wp-ws;

N=ceil(8*pi/wd);%計算濾波器長度

wn=(wp+ws)/2;%計算濾波器截止頻率

b=fir1(N，wn/pi，′high′，hamming(N+1));freqz(b,1,512);

程序運行結(jié)果如圖4.13所示。圖4.13基于漢明窗的FIR高通濾波器的頻率響應(yīng)曲線

2)多頻帶設(shè)計函數(shù)fir2

函數(shù)fir2用于設(shè)計任意形狀頻率響應(yīng)的FIR濾波器。函數(shù)fir2的調(diào)用格式為

b=fir2(n，f，m，npt，lap，windows)

其中，f和m分別為濾波器期望幅頻響應(yīng)和幅值向量，取值在0～1之間，f和m長度相同;npt為對頻率響應(yīng)進行內(nèi)插的點數(shù)，缺省值為512;lap定義一個區(qū)域尺寸，函數(shù)在重復(fù)頻率點周圍建立這個區(qū)域并提供光滑、陡峭的過渡頻率響應(yīng)，缺省值為25;其他項同函數(shù)fir1。

【例4.13】

設(shè)計一個30階的低通FIR濾波器，使之接近于以下理想頻率特性:

f=[00.60.61];m=[1100]

MATLAB程序如下:

%MATLABPROGRAM4-13

f=[0，0.6，0.6，1];

m=[1，1，0，0];

b=fir2(30，f，m);

[h，w]=freqz(b，1，256);

plot(f，m，w/pi，abs(h));

程序運行結(jié)果如圖4.14所示。圖4.1430階低通FIR濾波器的幅頻響應(yīng)曲線4.3.2頻域采樣法設(shè)計FIR濾波器

頻域采樣法是在頻率域?qū)硐霝V波器Hd(ejω)采樣，在采樣點上設(shè)計的濾波器H(ejω)和理想濾波器Hd(ejω)的幅度值相等，然后根據(jù)頻域的采樣值H(k)求得實際設(shè)計的濾波器的頻率特性H(ejω)。

在工程設(shè)計中，由于經(jīng)常給定頻域上的技術(shù)指標，所以采用頻域采樣法設(shè)計更直接，一般處理方法是對理想濾波器的頻率特性Hd(ejω)在[0，2π]范圍內(nèi)等間隔地取樣N個點。頻域采樣法得到的濾波器的單位脈沖響應(yīng)h(n)為(4.20)為了設(shè)計線性相位的FIR濾波器，采樣值H(k)要滿足一定的約束條件。具有線性相位的FIR濾波器的單位脈沖響應(yīng)h(n)是實序列，且滿足h(n)=±h(N－1－n)。由此得到的幅頻和相頻特性就是對H(k)的約束。

【例4.14】

用頻域采樣法設(shè)計一個帶通濾波器，滿足:低通帶邊沿ω1p=0.25π;低阻帶邊沿ω1s=0.36π;高通帶邊沿ω2p=0.68π;高阻帶邊沿ω2s=0.85π。設(shè)計過渡帶中的頻率樣本值為T1=0.108562和T2=0.58427685。MATLAB程序如下:

%MATLABPROGRAM4-14

M=40;

al=(M－1)/2;

l=[0:M－1];

T1=0.108562;

T2=0.58427685;

Hrs=[zeros(1，5)，T1，T2，ones(1，7)，T2，T1，zeros(1，9)，T1，T2，ones(1，7)，T2，T1，

zeros(1，4)];

k1=[0:floor((M－1)/2)];

k2=[floor((M－1)/2)+1:M－1];

angh=[－al*(2*pi)/M*k1，al*(2*pi)/M*(M－k2)];

H=Hrs.*exp(j*angh);

h=real(ifft(H，M));

freqz(h，1，512，1000);

程序運行結(jié)果如圖4.15所示。圖4.15用頻域采樣法設(shè)計的FIR數(shù)字濾波器頻率響應(yīng)曲線4.3.3最優(yōu)化法設(shè)計FIR濾波器

1.均方誤差最小化準則

若以E(ejω)表示逼近誤差

E(ejω)=Hd(ejω)－H(ejω)

(4.21)

則均方誤差為(4.22)均方誤差最小化準則就是選擇一組時域采樣值，使均方誤差ε2最小。這一方法注重的是在整個(－π～π)頻率區(qū)間內(nèi)總誤差的全局最小，但不能保證局部頻率點的性能，有些頻率點可能會有較大的誤差。窗函數(shù)法設(shè)計FIR濾波器，實際上是采用有限項的h(n)逼近理想的hd(n)。逼近誤差為(4.23)

MATLAB信號處理模塊提供了函數(shù)firls來設(shè)計基于均方誤差最小化的FIR濾波器。函數(shù)firls是函數(shù)fir1和fir2的擴展，調(diào)用格式為

b=firls(n，f，a，w，′d′)其中，n為濾波器階數(shù);f為濾波器期望頻率特性的頻率向量標準化頻率，取值為0~1，遞增向量，允許定義重復(fù)頻率點;a為濾波器所期望的幅值向量，f與a同長度且為偶數(shù);w為權(quán)向量，為f和a向量的一半，一個頻帶必須對應(yīng)一個權(quán)值，可缺省;′d′表示選擇項，若所設(shè)計的濾波器為微分濾波器，則可缺省;b為函數(shù)返回的濾波器系數(shù)，長度為n＋1，且具有對偶關(guān)系。

【例4.15】

設(shè)計一個25階的高通濾波器，通帶邊界頻率為0.55π，幅值為1，阻帶邊界頻率為0.45π，幅值為0。

MATLAB程序如下:

%MATLABPROGRAM4-15

n=25;

f=[0，0.45，0.55，1];

a=[0，0，1，1];

b=firls(n，f，a);

[h，w]=freqz(b);

axes(′position′，[0.2，0.2，0.5，0.5，]);

plot(w/pi，abs(h));

xlabel(′歸一化頻率′);ylabel(′幅值′);grid;

程序運行結(jié)果如圖4.16所示。圖4.16最優(yōu)FIR高通濾波器幅頻曲線

2.最大誤差最小化準則

最大誤差最小化準則也稱切比雪夫最佳一致逼近準則，表示為

max|E(ejω)|=minω∈F

(4.24)

式中，F(xiàn)是根據(jù)要求預(yù)先給定的一個頻率取值范圍，可以是通帶或阻帶。

雷米茲(Remez)算法給出了求解切比雪夫最佳一致逼近問題的方法：

(1)在頻率取值范圍F上均勻等間隔地選取M+2個頻率值ω0，ω1，…，ωM+1作為初值，并計算ρ:

(4.25)式中

(2)由ωi求H(ω)和E(ω)，利用重心形式的拉格朗日插值公式，可得(4.26)式中，(4.27)由此得到(4.28)若在F范圍內(nèi)，對所有頻率都有|E(ω)|≤ρ，則ρ為所求，ω0，ω1，…，ωM+1可視為極值點頻率。

(3)對上次確定的極值點頻率ω0，ω1，…，ωM+1中的每一點，在其附近檢查是否在某一頻率處有|E(ω)|>ρ，若有，則以該頻率點作為新的局部極值點。對M+2個極值點頻率依次進行檢查，得到一組新的極值點頻率。重復(fù)步驟(1)、(2)，求出ρ、H(ω)、E(ω)，完成一次迭代。重復(fù)上述步驟，直到ρ的值改變很小，則迭代結(jié)束，這個ρ即為所求的h(n)最大誤差最小值。由最后一組極值點頻率求出H(ω)，反變換得到ωc和ωs，完成設(shè)計。

MATLAB信號處理工具箱中的函數(shù)remez可實現(xiàn)Parks-McClellan算法，這種算法利用雷米茲交換算法和切比雪夫近似理論來設(shè)計濾波器，使實際頻率響應(yīng)擬合期望頻率響應(yīng)達到最優(yōu)。函數(shù)調(diào)用格式為

b=remez(n，f，m)

b=remez(n，f，m，w，′h′)

b=remez(n，f，m，w，′d′)

其中，′h′為選擇項，表示設(shè)計的濾波器是奇對稱線性相位濾波器，濾波器可實現(xiàn)信號的赫爾伯特(Hilbert)變換。其他參數(shù)與函數(shù)firls相同。

函數(shù)調(diào)用根據(jù)所設(shè)計濾波器的最優(yōu)形式的不同略有不同，有基本形式的最優(yōu)濾波器、加權(quán)最優(yōu)濾波器、反對稱(赫爾伯特)FIR濾波器及微分濾波器。

【例4.16】

采用Parks-McClellan算法，設(shè)計一個17階的帶通濾波器，并畫出期望的幅頻特性曲線和實際的幅頻特性曲線。其中，f=[00.30.40.60.71];m=[001100]。

MATLAB程序如下:

%MATLABPROGRAM4-16

f=[0，0.3，0.4，0.6，0.7，1];

m=[0，0，1，1，0，0];

b=remez(17，f，m);

[h，w]=freqz(b，1，512);

plot(f，m，′b-′，w/pi，abs(h)，′b:′);

xlabel(′歸一化頻率′);ylabel(′幅值′);grid;

legend(′desired′，′remez′);

程序運行結(jié)果如圖4.17所示。圖4.17最優(yōu)帶通FIR濾波器的幅頻特性曲線

【例4.17】

利用函數(shù)firls和remez設(shè)計一個23階的高通反對稱線性相位濾波器，并繪制其幅頻特性圖。其中，f=[00.20.31];m=[0011]。

MATLAB程序如下:

%MATLABPROGRAM4-17

clf;

n=23;

f=[0，0.2，0.3，1];

m=[0，0，1，1];

b1=firls(n，f，m，′h′);%用firls設(shè)計FIR濾波器

[h1，w1]=freqz(b1，1，512);

b2=remez(n，f，m，′h′);%用remez設(shè)計FIR濾波器

[h2，w2]=freqz(b2，1，512);

plot(f，m，′b－′，w1/pi，abs(h1)，′k:′，w2/pi，abs(h2)，′r-.′);

xlabel(′歸一化頻率′);ylabel(′幅值′);grid;

legend(′b－′，′desired′，′k:′，′firls′，′r－.′，′remez′);

程序運行結(jié)果如圖4.18所示。圖4.18反對稱高通FIR濾波器的幅頻特性曲線4.3.4約束最小二乘FIR濾波器

1.函數(shù)fircls

函數(shù)fircls用于實現(xiàn)基于約束最小二乘法的線性相位多頻帶、分段常數(shù)FIR濾波器，調(diào)用格式為

b=fircls(n，f，a，up，lo，′flag′)

其中，n為濾波器的階數(shù);f和a分別為給定濾波器的期望幅頻特性向量和幅值向量，f為標準化頻率，在0～1之間，a的長度為length(f)－1;up和lo分別為每個頻帶的上邊界和下邊界頻率，均為向量，且長度等于a的長度;b為返回的FIR濾波器系數(shù)向量，長度為(n＋1);′flag′為選擇項，用于監(jiān)視濾波器設(shè)計，trace表示文字跟蹤，plot表示繪制濾波器的幅頻圖、群延時及零極點圖。

【例4.18】

利用CLS法設(shè)計一個127階的多頻帶濾波器，滿足以下要求:

當頻率為0～0.25π時，幅值為0，允許變化范圍為[－0.005，0.005];

當頻率為0.25π～0.45π時，幅值為0.5，允許變化范圍為[0.48，0.52];

當頻率為0.45π～0.65π時，幅值為0，允許變化范圍為

[－0.035，0.035];

當頻率為0.65π～0.85π時，幅值為1，允許變化范圍為[0.97，1.03];

當頻率為0.85π～1.0π時，幅值為0，允許變化范圍為

[－0.04，0.04]。MATLAB程序如下:

%MATLABPROGRAM4-18

n=127;

f=[0，0.25，0.45，0.65，0.85，1];

a=[0，0.5，0，1，0];

f1=[0，0.25，0.25，0.45，0.45，0.65，0.65，0.85，0.85，1];

a1=[0，0，0.5，0.5，0，0，1，1，0，0];

up=[0.005，0.52，0.035，1.03，0.04];

lo=[-0.005，0.48，-0.035，0.97，-0.04];

b=fircls(n，f，a，up，lo);

[h，w]=freqz(b，1，512);

axes(′position′，[0.2，0.2，0.5，0.5]);

plot(f1，a1，′b-′，w/pi，abs(h)，′m:′);

xlabel(′歸一化頻率′);ylabel(′幅值′);grid;

legend(′desired′，′fircls′，2);

程序運行結(jié)果如圖4.19所示。圖4.19多頻帶濾波器的幅頻特性曲線

2.函數(shù)fircls1

函數(shù)fircls1用于基本線性相位的低通和高通濾波器的設(shè)計，低通濾波器的調(diào)用格式為

b=fircls1(n，w0，dp，ds)

b=fircls1(n，w0，dp，ds，wt，′flag′)

其中，n為濾波器的階數(shù);w0為濾波器的截止頻率，為標準化頻率;dp為通帶距幅值為1的最大偏差;ds為阻帶距幅值為0的最大偏差;′flag′為設(shè)計監(jiān)測標志;b為濾波器系數(shù)向量;wt為一個定義頻率，以確保設(shè)計的濾波器滿足通帶或阻帶的邊界要求。高通濾波器的調(diào)用格式為

b=fircls1(n，w0，dp，ds，′high′)

b=fircls1(n，w0，dp，ds，wt，′high′，′flag′)

其中，′high′表示設(shè)計高通濾波器，其余各項意義同低通濾波器。

【例4.19】

用CLS法設(shè)計一個55階的高通濾波器，截止頻率為0.4π，通帶允許最大波紋為0.03dB，

阻帶允許最大波紋為0.009dB。

MATLAB程序如下:

%MATLABPROGRAM4-19

n=55;

w0=0.4;

dp=0.03;

ds=0.009

b=fircls1(n，w0，dp，ds，′high′，′plot′);

程序運行結(jié)果如圖4.20所示。圖4.20CLS高通濾波器的幅頻及偏差圖4.3.5任意響應(yīng)法設(shè)計FIR濾波器

1.多頻帶復(fù)響應(yīng)濾波器的設(shè)計

函數(shù)cremez最基本的調(diào)用格式為

b=cremez(n，f，′frep′)

b=cremez(n，f，{′frep′，p1，p2…}，w)

b=cremez(n，f，a，w)

其中，n為濾波器的階數(shù);f為濾波器的多頻帶邊界頻率向量，取值可在－1～1之間;w為權(quán)向量，對每個頻段上幅值擬合度加權(quán)，為非負值，w向量長度為頻率向量f的一半;′frep′為濾波器預(yù)先定義的頻率響應(yīng)函數(shù):

(1)lowpass(低通)、highpass(高通)、bandpass(帶通)、bandstop(帶阻)等用于設(shè)計標準型濾波器。例如:

b=cremez(n，f，′lowpass′，…)

(2)multiband用于設(shè)計多頻帶線性相位任意幅頻響應(yīng)濾波器。例如:

b=cremez(n，f，{′multiband′，a}，…)

其中，a為邊界頻帶向量f點上的期望幅值。

(3)differenti