第8章實(shí)驗(yàn)教程8.1實(shí)驗(yàn)一時(shí)域離散信號(hào)與系統(tǒng)分析8.2實(shí)驗(yàn)二離散傅立葉變換及其應(yīng)用8.3實(shí)驗(yàn)三IIR數(shù)字濾波器設(shè)計(jì)8.4實(shí)驗(yàn)四FIR數(shù)字濾波器設(shè)計(jì)8.5MATLAB數(shù)字信號(hào)處理基礎(chǔ)8.1實(shí)驗(yàn)一時(shí)域離散信號(hào)與系統(tǒng)分析

8.1.1實(shí)驗(yàn)?zāi)康?/p>

(1)熟悉連續(xù)信號(hào)經(jīng)理想采樣后的頻譜變化關(guān)系，加深對(duì)時(shí)域采樣定理的理解。

(2)熟悉時(shí)域離散系統(tǒng)的時(shí)域特性，利用卷積方法觀察分析系統(tǒng)的時(shí)域特性。

(3)學(xué)會(huì)離散信號(hào)及系統(tǒng)響應(yīng)的頻域分析。

(4)學(xué)會(huì)時(shí)域離散信號(hào)的MATLAB編程和繪圖。

(5)學(xué)會(huì)利用MATLAB進(jìn)行時(shí)域離散系統(tǒng)的頻率特性分析。8.1.2實(shí)驗(yàn)原理與方法

采樣是連續(xù)信號(hào)數(shù)字處理的第一個(gè)關(guān)鍵環(huán)節(jié)。對(duì)采樣過(guò)程的研究不僅可以了解采樣前后信號(hào)時(shí)域特性發(fā)生的變化及信號(hào)信息不丟失的條件，而且可以加深對(duì)傅立葉變換、Z變換和序列傅立葉變換之間關(guān)系式的理解。

我們知道，對(duì)一個(gè)連續(xù)信號(hào)xa(t)進(jìn)行理想采樣的過(guò)程可用式(8.1.1)表示：(8.1.1)其中：xa(t)為xa(t)的理想采樣，p(t)為周期沖擊脈沖，即∧(8.1.2)xa(t)的傅立葉變換xa(jΩ)為∧∧(8.1.3)式(8.1.3)表明：X

a(jΩ)為Xa(jΩ)的周期延拓，其延拓周期為采樣角頻率(ΩS=2π/T)。只有滿足采樣定理時(shí)，才不會(huì)發(fā)生頻率混疊失真。在計(jì)算機(jī)上用高級(jí)語(yǔ)言編程直接按式(8.1.3)計(jì)算理想采樣X(jué)a(t)的頻譜Xa(jΩ)很不方便。下面導(dǎo)出用序列的傅立葉變換來(lái)計(jì)算Xa(jΩ)的公式?！摹摹摹摹膶⑹?8.1.2)代入式(8.1.1)并進(jìn)行傅立葉變換：(8.1.4)式中的xa(nT)就是采樣后得到的序列x(n)，即x(n)=xa(nT)x(n)的傅立葉變換為(8.1.5)(8.1.6)比較式(8.1.5)和式(8.1.4)可知：這說(shuō)明兩者之間只在頻率度量上差一個(gè)常數(shù)因子T。實(shí)驗(yàn)過(guò)程中應(yīng)注意這一差別。

離散信號(hào)和系統(tǒng)在時(shí)域均可用序列來(lái)表示。序列圖形給人以形象直觀的印象，它可加深我們對(duì)信號(hào)和系統(tǒng)的時(shí)域特征的理解。本實(shí)驗(yàn)還將觀察分析幾種信號(hào)及系統(tǒng)的時(shí)域特性。

為了在數(shù)字計(jì)算機(jī)上觀察分析各種序列的頻域特性，通常對(duì)X(ejω)在[0，2π]上進(jìn)行M點(diǎn)采樣來(lái)觀察分析。對(duì)長(zhǎng)度為N的有限長(zhǎng)序列x(n)，有(8.1.7)其中:通常M應(yīng)取的大一些，以便觀察譜的細(xì)節(jié)變化。取模|X(ejωk)|可繪出幅頻特性曲線。

一個(gè)時(shí)域離散線性非移變系統(tǒng)的輸入/輸出關(guān)系為(8.1.8)這里，y(n)為系統(tǒng)的輸出序列，x(n)為系統(tǒng)的輸入序列。h(n)、x(n)可以是有限長(zhǎng)的。為了計(jì)算機(jī)繪圖觀察方便，主要討論有限長(zhǎng)情況。如果h(n)和x(n)的長(zhǎng)度分別為N和M，則y(n)的長(zhǎng)度為L(zhǎng)=N+M－1。這樣，式(8.1.8)所描述的卷積運(yùn)算就是序列移位、相乘和累加的過(guò)程，所以編程十分簡(jiǎn)單。 上述卷積運(yùn)算也可以在頻域使用卷積定理實(shí)現(xiàn)：

Y(ejω)=X(ejω)H(ejω)

(8.1.9)式(8.1.9)右邊的相乘是在各個(gè)頻點(diǎn){ωk}上的頻譜值相乘。8.1.3實(shí)驗(yàn)內(nèi)容

1.序列的產(chǎn)生

用MATLAB編程實(shí)現(xiàn)下列序列，并用stem語(yǔ)句繪出波形圖。

(1)單位脈沖序列：x(n)=δ(n)。

(2)矩形序列：x(n)=RN(n)，N=10。

(3)x(n)=e(0.8+3j)n，0≤n≤15。

(4)x(n)=3cos(0.3π+0.25πn)+2sin(0.2π+0.125πn)，0≤n≤15。

(5)將(3)中的序列以20為周期拓展3個(gè)周期。

(6)假設(shè)x(n)=[1，－3，2，3，－2]，編程產(chǎn)生以下序列并繪出波形圖：

x1(n)=x(n)－2x(n+2)+x(n－1)+x(n－3)

(7)編寫(xiě)MATLAB函數(shù)stepshift(n0，n1，n2)實(shí)現(xiàn)單位階躍系列的移位信號(hào)u(n-n1)，其中：n0，n2分別為階躍系列的起點(diǎn)和終點(diǎn)，n1為移位的點(diǎn)數(shù)。程序要求輸入任意n0、n1、n2都能生成u(n-n1)，并繪出波形圖。

2.采樣信號(hào)及其頻譜分析

(1)繪出時(shí)間信號(hào)x(t)=cos(50πt)sin(πt)，0≤t≤2s。

(2)對(duì)于連續(xù)信號(hào)：X(t)=Ae－atsin(Ω0t)u(t)進(jìn)行采樣，可得到采樣序列：

x(n)=x(nT)=Ae－anTsin(Ω0t)u(n)，0≤n<50

其中A為幅度因子，a為衰減因子，Ω0是模擬角頻率，T為采樣間隔。假設(shè)A=500，a=200，Ω0=50πrad/s。取采樣頻率fs=2kHz，即T=0.5ms。觀察所得采樣x(n)的幅頻特性|X(ejω)|在折疊頻率附近有無(wú)明顯差別。改變采樣頻率，fs=1000Hz,觀察|X(ejω)|的變化，并做記錄。fs=500Hz，觀察頻譜混疊是否明顯存在，說(shuō)明原因，并記錄這時(shí)的|X(ejω)|曲線。

3.系統(tǒng)的單位脈沖響應(yīng)

求以下差分方程所描述的系統(tǒng)的單位脈沖響應(yīng)h(n)(0≤n<50):

y(n)+0.2y(n－1)－y(n－2)=2x(n)－3x(n－1)

4.有限長(zhǎng)序列線性卷積

兩個(gè)給定長(zhǎng)度的序列的卷積可以直接調(diào)用MATLAB語(yǔ)言中的卷積函數(shù)conv。conv用于兩個(gè)有限長(zhǎng)序列的卷積，它假定兩個(gè)序列都從n=0開(kāi)始。調(diào)用格式如下：

y=conv(x，h)

其中參數(shù)x、h是兩個(gè)已賦值的行向量序列。假設(shè)x=[1,-1,2，3]，試求上面所述系統(tǒng)的輸出y(n)，并繪出輸出波形圖。

5.時(shí)域離散系統(tǒng)的頻域分析

假設(shè)一因果系統(tǒng)的系統(tǒng)函數(shù)為分析該系統(tǒng)的頻率響應(yīng)，畫(huà)出幅頻和相頻特性曲線。8.1.4思考題

(1)在分析理想采樣序列特性的實(shí)驗(yàn)中，采樣頻率不同時(shí)，相應(yīng)理想采樣序列的傅立葉變換頻譜的數(shù)字頻率度量是否都相同？它們所對(duì)應(yīng)的模擬頻率是否相同？為什么？

(2)對(duì)于由兩個(gè)子系統(tǒng)級(jí)聯(lián)或者并聯(lián)的系統(tǒng)，如何用MATLAB分析它們的幅頻和相頻響應(yīng)？8.1.5部分參考程序(MATLAB語(yǔ)言)

(1)單位階躍系列的移位信號(hào)。

n0=input(′PleaseInputn0：＼n′，′s′)；

n0=str2num(n0)；

n1=input(′PleaseInputn1：(n1<n0)＼n′，′s′)；

n1=str2num(n1)；

n2=input(′PleaseInputn2：＼n′，′s′)；

n2=str2num(n2)；

ifn1<1

n1_1=n1+abs(n1)+1

n0_1=n0+abs(n1)+1

n2_1=n2+abs(n1)+1else

n1_1=n1

n0_1=n0

n2_1=n2

end

x=stepshift(n0_1，n1_1，n2_1)

n=n1：n1+length(x)-1

functiony=stepshift(n0，n1，n2)

forn=n1：n2

ifn>=n0

y(n)=1；

else

y(n)=0；

end

end

(2)采樣信號(hào)及其頻譜分析。

A=500；

a=200；

w=50*pi；

fs=500；

Xa=FF(A，a，w，fs)；%產(chǎn)生信號(hào)xa(n)

string=[′fs=′，num2str(fs)]；

DFT(Xa，50，string)；%TheFTandpaintThepicture

function[c，l]=DFT(x，N，str)

n=0：N-1；

k=-200：200；

w=(pi/100)*k；

l=w；c=x*(exp(-j*pi/100)).^(n′*k)；%計(jì)算DFT[x(n)]

magX=abs(c)；%TheAmpalitatude

angX=angle(c)；%TheAngle

subplot(2，1，1)；

t=max(x)；

n=0：N-1；

stem(n，x，′.′)；

xlabel(′n′)；

ylabel(′x(n)′)；

title(′信號(hào)的原形′)；

text((0.3*N)，(0.8*t)，str)；

holdon

n=0：N-1；

m=zeros(N)；plot(n，m)；

subplot(2，1，2)；

plot(w/pi，magX)；

xlabel(′w/pi′)；

ylabel(′|X(jw)|′)；

title(′上圖信號(hào)的傅立葉變換(|X(jw)|)′)；

%子函數(shù)：產(chǎn)生信號(hào)xa(n)

functionc=FF(A，a，w，fs)

n=0：50-1；

c=A*exp((-a)*n/fs).*sin(w*n/fs).*stepseq(0，0，49)；

8.2實(shí)驗(yàn)二離散傅立葉變換及其應(yīng)用

8.2.1實(shí)驗(yàn)?zāi)康?/p>

(1)通過(guò)實(shí)驗(yàn)進(jìn)一步加深DFT算法及FFT算法原理和基本性質(zhì)的理解。

(2)熟悉MATLAB中有關(guān)FFT算法函數(shù)及其使用方法。

(3)學(xué)會(huì)用FFT對(duì)連續(xù)信號(hào)和時(shí)域離散信號(hào)進(jìn)行譜分析的方法，了解可能出現(xiàn)的分析誤差

及原因，以便在實(shí)際中正確應(yīng)用FFT。

(4)學(xué)會(huì)應(yīng)用FFT實(shí)現(xiàn)序列的線性卷積和相關(guān)。8.2.2實(shí)驗(yàn)原理與前期準(zhǔn)備

(1)復(fù)習(xí)DFT的定義、性質(zhì)和用DFT作譜分析的有關(guān)內(nèi)容。

(2)復(fù)習(xí)FFT算法進(jìn)行譜分析的方法。對(duì)于時(shí)域離散信號(hào)，選擇的FFT點(diǎn)數(shù)應(yīng)大于等于序列的長(zhǎng)度；對(duì)于連續(xù)信號(hào)，則首先根據(jù)其最高頻率確定采樣速率fs以及由頻率分辨率選擇采樣點(diǎn)數(shù)N，然后對(duì)其進(jìn)行軟件采樣。頻率分辨率的選擇要以能分辨開(kāi)其中的關(guān)鍵頻率譜線為準(zhǔn)則。對(duì)周期序列，最好截取周期的整數(shù)倍進(jìn)行譜分析，否則有可能產(chǎn)生較大的分析誤差。

(3)編制信號(hào)產(chǎn)生子程序，產(chǎn)生以下典型信號(hào)供譜分析用：8.2.3實(shí)驗(yàn)內(nèi)容

1.序列譜分析

(1)對(duì)x1(n)分別進(jìn)行N=8和N=16的譜分析。

(2)對(duì)x2(n)、x3(n)分別進(jìn)行N=8和N=16的譜分析，畫(huà)出時(shí)域波形和頻譜曲線，比較分析結(jié)果的異同，說(shuō)明原因。

(3)對(duì)高斯序列x4(n)進(jìn)行N=16的譜分析。參數(shù)p=8時(shí)，q=2，4，8；參數(shù)q=8時(shí)，p=8，13，14，分別畫(huà)出時(shí)域波形和頻譜曲線，比較分析結(jié)果的異同，說(shuō)明是否出現(xiàn)頻率混疊以及出現(xiàn)的原因。

(4)分析衰減正弦序列x5(n)的時(shí)域波形及頻譜特性。當(dāng)α=0.1，f=0.0625時(shí)，觀察譜峰出現(xiàn)的位置，判斷該譜線對(duì)應(yīng)的頻率。改變f分別為0.4375和0.5625時(shí)，觀察譜峰出現(xiàn)的位置和頻譜形狀，判斷有無(wú)混疊和泄漏。

2.連續(xù)信號(hào)譜分析

對(duì)連續(xù)信號(hào)x6(t)，假設(shè)Δf=1/16及Δf=1/64時(shí)，試采用FFT方法，選擇合適的參數(shù)分析該信號(hào)的頻譜，畫(huà)出時(shí)域波形和頻譜曲線。

3.驗(yàn)證DFT對(duì)稱(chēng)性

(1)對(duì)x7(n)、x8(n)分別做N=8和N=16譜分析，畫(huà)出頻譜形狀。

(2)令x(n)=x7(n)+x8(n)，用FFT計(jì)算8點(diǎn)和16點(diǎn)離散傅立葉變換，X(k)=DFT[x(n)]。

根據(jù)DFT的對(duì)稱(chēng)性，由X(k)求出X7(k)=DFT[x7(n)]和X8(k)=DFT[x8(n))]，并與(1)中所得結(jié)果比較。[提示：取N=16時(shí)，x7(n)=x7(N－n)，x8(n)=－x8(N－n)]

(3)令x(n)=x7(n)+jx8(n)，重復(fù)步驟(2)。

4.用FFT計(jì)算循環(huán)卷積和線性卷積

用FFT分別計(jì)算對(duì)高斯序列x4(n)(p=8，q=2)和衰減正弦序列x5(n)(α=0.1，f=0.0625)的16循環(huán)卷積和線性卷積。

5.重迭相加法和重疊保留法

產(chǎn)生一個(gè)512點(diǎn)的隨機(jī)序列x(n)，將x(n)分成4段，分別采用重疊相加法和重疊保留法實(shí)現(xiàn)x(n)與x3(n)的線性卷積，觀察卷積前后隨機(jī)序列x(n)頻譜的變化。

6.互相關(guān)函數(shù)和自相關(guān)函數(shù)

(1)分別計(jì)算x2(n)和x3(n)的自相關(guān)函數(shù)，畫(huà)出它們的自相關(guān)函數(shù)波形圖。

(2)分別計(jì)算x2(n)和x3(n)的互相關(guān)函數(shù)，畫(huà)出互相關(guān)函數(shù)波形圖。

(3)分別計(jì)算x2(n)和x2(n－5)的互相關(guān)函數(shù)，畫(huà)出相關(guān)函數(shù)波形，并指出其峰值位置及意義。8.2.4思考題

(1)在N=8時(shí)，x2(n)和x3(n)的幅頻特性會(huì)相同嗎？為什么？N=16呢？

(2)如果周期信號(hào)的周期預(yù)先不知道，如何用FFT進(jìn)行譜分析？8.2.5部分參考程序

(1)高斯序列譜分析。

%p=8，q=2

forn=1：16；

p=8；

q=2；

xa(n)=exp(-(n-p)^2/q)；

end

subplot(2，1，1)；

stem(xa)；

xlabel(′n(p=8，q=2)′)；

ylabel(′xa(n)′)；

subplot(2，1，2)；X=fft(xa，16)；

k=0：15；

plot(k，abs(X))；

xlabel(′k(p=8，q=2)′)；

ylabel(′X(k)′)；

%p=8，q=4

forn=1：16；

p=8；

q=4；

xa(n)=exp(-(n-p)^2/q)；

end

figure

subplot(2，1，1)；

stem(xa)；xlabel(′n(p=8，q=4)′)；

ylabel(′xa(n)′)；

subplot(2，1，2)；

X=fft(xa，16)；

k=0：15；

plot(k，abs(X))；

xlabel(′k(p=8，q=4)′)；

ylabel(′X(k)′)；

%p=8，q=8

forn=1：16；

p=8；

q=8；

xa(n)=exp(-(n-p)^2/q)；

end

figure；subplot(2，1，1)；

stem(xa)；

xlabel(′n(p=8，q=8)′)；

ylabel(′xa(n)′)；

subplot(2，1，2)；

X=fft(xa，16)；

k=0：15；

plot(k，abs(X))；

xlabel(′k(p=8，q=8)′)；

ylabel(′X(k)′)；

%q固定

%p=13，q=8forn=1：16；

p=13；

q=8；

xa(n)=exp(-(n-p)^2/q)；

end

figure；

subplot(2，1，1)；

stem(xa)；

xlabel(′n(p=13，q=8)′)；

ylabel(′xa(n)′)；

subplot(2，1，2)；

X=fft(xa，16)；

k=0：15；

plot(k，abs(X))；xlabel(′k(p=13，q=8)′)；

ylabel(′X(k)′)；

%p=14，q=8

forn=1：16；

p=14；

q=8；

xa(n)=exp(-(n-p)^2/q)；

end

figure

subplot(2，1，1)；

stem(xa)；

xlabel(′n(p=14，q=8)′)；

ylabel(′xa(n)′)；

subplot(2，1，2)；X=fft(xa，16)；

k=0：15；

plot(k，abs(X))；

xlabel(′k(p=14，q=8)′)；

ylabel(′X(k)′)；

%p=8，q=8

forn=1：16；

p=8；

q=8；

xa(n)=exp(-(n-p)^2/q)；

end

figure；

subplot(2，1，1)；

stem(xa)；xlabel(′n(p=8，q=8)′)；

ylabel(′xa(n)′)；

subplot(2，1，2)；

X=fft(xa，16)；

k=0：15；

plot(k，abs(X))；

xlabel(′k(p=8，q=8)′)；

ylabel(′X(k)′)；

(2)衰減正弦信號(hào)譜分析。

f=0.0625；

a=0.1；

forn=1：16

xb(n)=exp(-a*n)*sin(2*pi*f*n)；

end

subplot(2，1，1)；

n=0：15；

stem(n，xb)；

xlabel(′n(a=0.1，f=0.0625)′)；

ylabel(′xb(n)′)；

subplot(2，1，2)；

Xb=fft(xb，16)；k=0：15；

stem(k，abs(Xb))；

xlabel(′k(a=0.1，f=0.0625)′)；

ylabel(′Xb(k)′)；

%f=0.4375，a=0.1

f=0.4375；

a=0.1；

forn=1：16

xb(n)=exp(-a*n)*sin(2*pi*f*n)；

end

figure；

subplot(2，1，1)；

n=0：15；

stem(n，xb)；xlabel(′n(a=0.1，f=0.4375)′)；

ylabel(′xb(n)′)；

subplot(2，1，2)；

Xb=fft(xb，16)；

k=0：15；

stem(k，abs(Xb))；

xlabel(′k(a=0.1，f=0.4375)′)；

ylabel(′Xb(k)′)；

%f=0.5625，a=0.1

f=0.5625；

a=0.1；

forn=1：16

xb(n)=exp(-a*n)*sin(2*pi*f*n)；

endfigure；

subplot(2，1，1)；

n=0：15；

stem(n，xb)；

xlabel(′n(a=0.1，f=0.5625)′)；

ylabel(′xb(n)′)；

subplot(2，1，2)；

Xb=fft(xb，16)；

k=0：15；

stem(k，abs(Xb))；

xlabel(′k(a=0.1，f=0.5625)′)；

ylabel(′Xb(k)′)；(3)重疊相加法和重疊保留法。

N=512；

xe=rand(1，N)-0.5；

n1=0：3；

xn1=n1；

n2=4：7；

xn2=8-n2；

xc=[xn1xn2]；

figure；

subplot(2，1，1)；

n=0：length(xe)-1；

stem(n，xe)；

xlabel(′n′)；

ylabel(′xe(n)′)；subplot(2，1，2)；

Xe=fft(xe，length(xe))；

k=0：length(xe)-1；

plot(k，abs(Xe))；

xlabel(′k′)；

ylabel(′Xe(k)′)；

%%

figure；

ye=conv(xe，xc)

subplot(2，1，1)；

n=0：length(ye)-1；

stem(n，ye)；

xlabel(′n′)；

ylabel(′ye(n)′)；subplot(2，1，2)；

Ye=fft(ye，length(ye))；

k=0：length(ye)-1；

plot(k，abs(Ye))；

xlabel(′k′)；

ylabel(′Ye(k)′)；

%%xiangjiafa

%%

%ye=conv(xe，xc)

N=length(xe)/8；

M=length(xc)；

L=N+M-1

A_ye=zeros(1，M-1)；

T_ye=zeros(0)；

fori=0：7

forn=1：N

xi(n)=xe(i*N+n)；

end

yin=conv(xi，xc)；

T_add=A_ye+[yin(1：M-1)]；

A_ye=yin(N+1：L)；

yin=[T_addyin(M：N)]；

T_ye=[T_yeyin]；

end

figure；

subplot(2，1，1)；n=0：length(T_ye)-1；

stem(n，T_ye)；

xlabel(′n′)；

ylabel(′T_ye(n)ADD′)；

subplot(2，1，2)；

Ye=fft(T_ye，length(T_ye))；

k=0：length(T_ye)-1；

plot(k，abs(Ye))；

xlabel(′k′)；

ylabel(′Ye(k)ADD′)；

%%baoliufa

%%%ye=conv(xe，xc)

N=length(xe)/8；

M=length(xc)；L=N+M-1

Tr_ye=zeros(0)；

forn=1：L

xi(n)=xe(n)；

end

yin=ifft((fft(xi，L).*fft(xc，L))，L)；

Tr_add=yin(M：L)；

Tr_ye=[Tr_yeTr_add]；

fori=1：7

forn=1：L

xi(n)=xe(i*N+n-M+1)；

end

yin=ifft((fft(xi，L).*fft(xc，L))，L)；

Tr_add=yin(M：L)；

Tr_ye=[Tr_yeTr_add]；

end

figure；

subplot(2，1，1)；

n=0：length(Tr_ye)-1；

stem(n，Tr_ye)；

xlabel(′n′)；

ylabel(′Tr_ye(n)′)；

subplot(2，1，2)；

Yre=fft(Tr_ye，length(Tr_ye))；

k=0：length(Tr_ye)-1；

plot(k，abs(Yre))；

xlabel(′k′)；

ylabel(′Yre(k)′)；

(4)相關(guān)運(yùn)算。

forn=1：16；

p=8；

q=2；

xa(n)=exp(-(n-p)^2/q)；

end

%f=0.0625，a=0.1

f=0.0625；

a=0.1；forn=1：16

xb(n)=exp(-a*n)*sin(2*pi*f*n)；

end

%%

k=length(xa)；

N=2*k；

rm=ifft(conj(fft(xa，N).*fft(xa，N)))；

rm=[rm(k+2：N)rm(1：k)]

m=(-k+1)：(k-1)；

subplot(2，1，1)；

stem(m，rm)；

xlabel(′m′)；

ylabel(′rma′)；

%%k=length(xb)；

N=2*k；

rm=ifft(conj(fft(xb，N).*fft(xb，N)))；

rm=[rm(k+2：N)rm(1：k)]

m=(-k+1)：(k-1)；

subplot(2，1，2)；

stem(m，rm)；

xlabel(′m′)；

ylabel(′rmb′)；

8.3實(shí)驗(yàn)三IIR數(shù)字濾波器設(shè)計(jì)

8.3.1實(shí)驗(yàn)?zāi)康?/p>

(1)掌握雙線性變換法、脈沖響應(yīng)不變法設(shè)計(jì)IIR數(shù)字濾波器的原理、方法和步驟，了解雙線性變換法、脈沖響應(yīng)不變法的特點(diǎn)。

(2)掌握模擬濾波器的設(shè)計(jì)原理與方法，熟悉巴特沃思、切比雪夫、橢圓濾波器的頻率特性。

(3)熟悉模擬低通原型進(jìn)行頻率變換設(shè)計(jì)低通、高通、帶通、帶阻濾波器的方法。

(4)學(xué)會(huì)利用MATLAB的有關(guān)函數(shù)根據(jù)給定的濾波器指標(biāo)設(shè)計(jì)各種類(lèi)型數(shù)字濾波器。8.3.2實(shí)驗(yàn)原理與方法

1.IIR數(shù)字濾波器的系統(tǒng)函數(shù)及特點(diǎn)

設(shè)IIR濾波器的輸入序列為x(n)，則IIR濾波器的輸入序列x(n)與輸出序列y(n)之間的關(guān)系可以用下面的差分方程式表示：其中，aj和bi是濾波器的系數(shù)，aj中至少有一個(gè)非零。與之相對(duì)應(yīng)的系統(tǒng)函數(shù)為由系統(tǒng)函數(shù)可以發(fā)現(xiàn)無(wú)限長(zhǎng)單位脈沖響應(yīng)濾波器有如下特點(diǎn)：

(1)單位脈沖響應(yīng)h(n)是無(wú)限長(zhǎng)的。

(2)系統(tǒng)函數(shù)H(z)在z平面上有極點(diǎn)存在。

(3)結(jié)構(gòu)上存在著輸出到輸入的反饋，也就是結(jié)構(gòu)上是遞歸型的。

2.由模擬濾波器轉(zhuǎn)換為數(shù)字濾波器的方法

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

脈沖響應(yīng)不變法是將系統(tǒng)從s平面映射到z平面的一種方法，使數(shù)字濾波器的單位脈沖響應(yīng)序列h(n)模仿模擬濾波器的沖激響應(yīng)ha(n)，其變換關(guān)系式為z=esT。由于z=esT是一個(gè)周期函數(shù)，因而s平面虛軸上每一段2π/T的線段都映射到z平面單位圓上一周。由于重疊映射，因而脈沖響應(yīng)不變法是一種多值映射關(guān)系。數(shù)字濾波器的頻率響應(yīng)是原模擬濾波器的頻率響應(yīng)的周期延拓。只有當(dāng)模擬濾波器的頻率響應(yīng)是有限帶寬的，且頻帶寬度|Ω|≤(π/T)=

(Ωs/2)時(shí)，才能避免數(shù)字濾波器的頻率響應(yīng)發(fā)生混疊現(xiàn)象。由于脈沖響應(yīng)不變法只適用于帶限的模擬濾波器，所以，在高頻區(qū)幅頻特性不等于零的高通和帶阻濾波器不能采用脈沖響應(yīng)不變法。

(2)雙線性變換法。

雙線性變換法是將整個(gè)s平面映射到整個(gè)z平面，其映射關(guān)系為雙線性變換法克服了脈沖響應(yīng)不變法從s平面到z平面的多值映射的缺點(diǎn)，消除了頻譜混疊現(xiàn)象。但其在變換過(guò)程中產(chǎn)生了非線性的畸變，在設(shè)計(jì)IIR數(shù)字濾波器的過(guò)程中需要進(jìn)行一定的頻率預(yù)畸變。頻率預(yù)畸變正公式為Ω=(2/T)tan(ω/2)。

3.設(shè)計(jì)過(guò)程

IIR數(shù)字濾波器的設(shè)計(jì)有多種方法，如頻率變換法、數(shù)字域直接設(shè)計(jì)法以及計(jì)算輔助設(shè)計(jì)等。

本實(shí)驗(yàn)只介紹頻率變換法，即由模擬低通原型濾波器到數(shù)字濾波器的轉(zhuǎn)換。其基本的設(shè)計(jì)過(guò)程如下：

(1)給定數(shù)字濾波器設(shè)計(jì)要求(低通、高通、帶通、帶阻)；

(2)轉(zhuǎn)換為模擬(低通、高通、帶通、帶阻)濾波器的技術(shù)指標(biāo)；

(3)轉(zhuǎn)換為模擬低通原型濾波器的指標(biāo)；

(4)設(shè)計(jì)得到滿足第三步要求的低通原型濾波器傳遞函數(shù)；

(5)通過(guò)頻率轉(zhuǎn)換得到模擬(低通、高通、帶通、帶阻)濾波器；

(6)變換為數(shù)字(低通、高通、帶通、帶阻)濾波器。

4.模擬原型法使用的MATLAB函數(shù)

(1)濾波器的階估計(jì)函數(shù)：Buttord，cheb1ord，cheb2ord，ellipord。

(2)低通模擬濾波器原型函數(shù)：buttap，cheb1ap，cheb2ap，ellipap。

(3)模擬原型低通濾波器到模擬濾波器(低通、高通、帶通、帶阻)變換函數(shù)：lp2lp，lp2hp，lp2bp，lp2bs。

(4)模擬到數(shù)字濾波器變換函數(shù)：Bilinear，impinvar。8.3.3實(shí)驗(yàn)內(nèi)容

(1)基于Butterworth模擬濾波器原型，使用脈沖響應(yīng)不變法、雙線性轉(zhuǎn)換法分別設(shè)計(jì)數(shù)字低通濾波器。參數(shù)指標(biāo)為：通帶截止頻率wp＝0.2π，通帶波動(dòng)值Rp＝1dB；阻帶截止頻率ws＝0.3π，阻帶最小衰減As＝20dB；濾波器采樣頻率Fs＝2000Hz。

要求：繪出幅頻特性、相頻特性圖并比較設(shè)計(jì)結(jié)果，繪出濾波器零極圖，寫(xiě)出濾波器的系統(tǒng)函數(shù)。

(2)分別基于巴特沃思、切比雪夫1型、切比雪夫2型、橢圓濾波器原型，采用脈沖響應(yīng)不變法設(shè)計(jì)一個(gè)IIR數(shù)字帶通濾波器。其參數(shù)指標(biāo)為：下通帶截止頻率wp1＝0.3π，上通帶截止頻率wp2＝0.7π，通帶最大衰減Rp＝1dB；下阻帶截止頻率ws1＝0.2π，上阻帶截止頻率ws2＝0.8π，阻帶最小衰減As＝20dB，濾波器采樣頻率Fs＝2000Hz。

要求：繪出幅頻特性、相頻特性圖并比較設(shè)計(jì)結(jié)果，繪出濾波器零極圖，寫(xiě)出濾波器的系統(tǒng)函數(shù)。

(3)基于chebyshev2型模擬濾波器原型設(shè)計(jì)，采用雙線性轉(zhuǎn)換設(shè)計(jì)一個(gè)切比雪夫2型數(shù)字帶阻濾波器。其參數(shù)指標(biāo)為：通帶低端頻率wp1＝0.2π，通帶高端頻率wp2＝0.8π，通帶最大衰減Rp＝1dB；阻帶低端頻率ws1＝0.3π，阻帶高端頻率ws2＝0.7π，阻帶最小衰減As＝20dB,濾波器采樣頻率Fs＝2000Hz。

要求：繪出幅頻特性、相頻特性圖并驗(yàn)證設(shè)計(jì)結(jié)果，繪出濾波器零極圖，寫(xiě)出濾波器的系統(tǒng)函數(shù)。

(4)基于橢圓型模擬濾波器原型設(shè)計(jì)，采用雙線性變換方法設(shè)計(jì)數(shù)字高通濾波器。其指標(biāo)參數(shù)如下：wp＝0.45π，Rp＝1dB；ws＝0.3π，As＝40dB，濾波器采樣頻率Fs＝2000Hz。

要求：繪出幅頻特性、相頻特性圖并驗(yàn)證設(shè)計(jì)結(jié)果，繪出濾波器零極圖，寫(xiě)出濾波器的系統(tǒng)函數(shù)。8.3.4思考題

(1)總結(jié)基于Butterworth、chebyshev1、chebyshev2、橢圓型模擬濾波器原型設(shè)計(jì)IIR數(shù)字濾波器過(guò)程的異同點(diǎn)及其濾波器的特征。

(2)為什么脈沖響應(yīng)不變法不能用于設(shè)計(jì)數(shù)字高通濾波器和帶阻濾波器？數(shù)字濾波器的頻率響應(yīng)與模擬濾波器的頻率響應(yīng)有何區(qū)別？

(3)使用雙線性變換法時(shí)，模擬頻率與數(shù)字頻率有何關(guān)系？會(huì)帶來(lái)什么影響？如何解決？

(4)用FDATool重新設(shè)計(jì)實(shí)驗(yàn)內(nèi)容。8.3.5參考程序(MATLAB語(yǔ)言)

(1)程序：

%數(shù)字濾波器指標(biāo)

wp=0.2*pi；%濾波器的通帶截止頻率

ws=0.3*pi；%濾波器的阻帶截止頻率

Rp=1；As=20；%濾波器的通阻帶衰減指標(biāo)

%轉(zhuǎn)換為模擬濾波器指標(biāo)

Fs=2000；T=1/Fs；

Omgp=(2/T)*tan(wp/2)；%雙線性變換原型通帶頻率預(yù)修正

Omgs=(2/T)*tan(ws/2)；%雙線性變換原型阻帶頻率預(yù)修正

%模擬原型濾波器計(jì)算

[n，Omgc]=buttord(Omgp，Omgs，Rp，As，′s′)；%計(jì)算階數(shù)n和截止頻率[z0，p0，k0]=buttap(n)；%歸一化原型設(shè)計(jì)

b0=k0*real(poly(z0))；%求原型濾波器系數(shù)b0

a0=real(poly(p0))；%求原型濾波器系數(shù)a0

[b，a]=lp2lp(b0，a0，Omgc)；%變換為模擬低通濾波器系數(shù)b，a

%用雙線性變換法計(jì)算數(shù)字濾波器系數(shù)

[num，den]=bilinear(b，a，F(xiàn)s)%數(shù)字濾波器設(shè)計(jì)結(jié)果

%求數(shù)字濾波器的頻率特性

[H，w]=freqz(num，den)；

db_H=20*log10((abs(H)+eps)/max(abs(H)))；%轉(zhuǎn)換為分貝值

subplot(2，2，1)，plot(w/pi，abs(H))；

ylabel(′|H|′)；title(′幅度響應(yīng)′)；axis([0，1，0，1.1])；

set(gca，′XTickMode′，′manual′，′XTick′，[0，0.2，0.3，1])；ripple=10^(-Rp/20)；At=10^(-As/20)；

set(gca,′YTickMode′,′manual′,′YTick′,[0，At,ripple，1])；grid

subplot(2，2，2)，plot(w/pi，angle(H)/pi)；

ylabel(′＼phi′)；title(′相位響應(yīng)′)；axis([0，1，-1，1])；

set(gca,′XTickMode′,′manual′,′XTick′,[0，0.2，0.3，1])；

set(gca，′YTickMode′，′manual′，′YTick′，[-1，0，1])；grid

subplot(2，2，3)，plot(w/pi，db_H)；title(′幅度響應(yīng)(dB)′)；

ylabel(′dB′)；xlabel(′頻率′)；axis([0，1，-40，5])；

set(gca,′XTickMode′,′manual′,′XTick′，[0，0.2，0.3，1])；

set(gca,′YTickMode′,′manual′,′YTick′,[-50,-20,-1,0])；grid

subplot(2，2，4)，zplane(num，den)；

axis([-1.1，1.1，-1.1，1.1])；title(′零極圖′)；(2)程序：

%數(shù)字濾波器指標(biāo)

wp1=0.3*pi；wp2=0.7*pi；

ws1=0.1*pi；ws2=0.9*pi；

Rp=1；As=20；

%轉(zhuǎn)換為模擬濾波器指標(biāo)

Fs=2000；T=1/Fs；

Omgp1=wp1*Fs；Omgp2=wp2*Fs；

Omgp=[Omgp1，Omgp2]；

Omgs1=ws1*Fs；Omgs2=ws2*Fs；

Omgs=[Omgs1，Omgs2]；

bw=Omgp2-Omgp1；w0=sqrt(Omgp1*Omgp2)；%中心頻率

[n，Omgn]=cheb1ord(Omgp，Omgs，Rp，As，′s′)

[z0，p0，k0]=cheb1ap(n，Rp)；[b0,a0]=zp2tf(z0,p0,k0)；

[ba,aa]=lp2bp(b0,a0,w0,bw)%低通原型轉(zhuǎn)換為帶通

[bd,ad]=impinvar(ba,aa,Fs)

[H,w]=freqz(bd,ad)；

dbH=20*log10(abs(H)/max(abs(H)))；

%繪圖

subplot(2,2,1),plot(w/pi,abs(H))；

ylabel(′幅度′)；xlabel(′頻率′)；axis([0,1,0,1.1])；

set(gca,′XTickMode′,′manual′,′XTick′,[0.2,0.3,0.7,0.8])；

ripple=10^(-Rp/20)；

At=10^(-As/20)；

set(gca,′YTickMode′,′manual′,′YTick′,[0,At,ripple,1])；grid

subplot(2,2,2)，plot(w/pi,angle(H)/pi)；ylabel(′相位′)；xlabel(′頻率′)；axis([0，1，-1，1])；

set(gca,′XTickMode′,′manual′,′XTick′,[0.2,0.3,0.7,0.8])；

set(gca,′YTickMode′,′manual′,′YTick′，[-1,0,1])；grid

subplot(2,2，3)，plot(w/pi，dbH)；

ylabel(′幅度(dB)′)；xlabel(′頻率′)；axis([0，1，-40，5])；

set(gca，′XTickMode′，′manual′，′XTick′，[0.2，0.3，0.7，0.8])；

set(gca，′YTickMode′，′manual′，′YTick′，[-50，-20，-1，0])；grid

subplot(2，2，4)，zplane(bd，ad)；

axis([-1.1，1.1，-1.1，1.1])；ylabel(′零極圖′)；(3)程序：

%數(shù)字濾波器指標(biāo)

ws1=0.3*pi；ws2=0.7*pi；

wp1=0.2*pi；wp2=0.8*pi；

Rp=1；As=20；

%轉(zhuǎn)換為模擬濾波器指標(biāo)

Fs=2000；T=1/Fs；

Omgp1=(2/T)*tan(wp1/2)；Omgp2=(2/T)*tan(wp2/2)；

Omgp=[Omgp1，Omgp2]；

Omgs1=(2/T)*tan(ws1/2)；Omgs2=(2/T)*tan(ws2/2)；

Omgs=[Omgs1，Omgs2]；

bw=Omgp2-Omgp1；w0=sqrt(Omgp1*Omgp2)；%阻帶中心頻率

%低通原型設(shè)計(jì)

[n，Omgn]=cheb2ord(Omgp，Omgs，Rp，As，′s′)；[z0，p0，k0]=cheb2ap(n，As)；

[b0，a0]=zp2tf(z0，p0，k0)；

[ba，aa]=lp2bs(b0,a0,w0,bw)；%模擬低通原型轉(zhuǎn)換為模擬帶阻濾波器

[bd，ad]=bilinear(ba，aa，F(xiàn)s)

[H，w]=freqz(bd，ad)；

dbH=20*log10(abs(H)/max(abs(H)))；

%繪圖

subplot(2，2，1)，plot(w/pi，abs(H)，′k′)；

ylabel(′幅度′)；xlabel(′頻率′)；axis([0，1，0，1.1])；

set(gca，′XTickMode′，′manual′，′XTick′，[0.2，0.3，0.7，0.8])；

ripple=10^(-Rp/20)；

Attn=10^(-As/20)；

set(gca,′YTickMode′,′manual′,′YTick′,[0,Attn,ripple,1])；gridsubplot(2，2，2)，plot(w/pi，angle(H)/pi*180，′k′)；

ylabel(′相位′)；xlabel(′頻率′)；axis([0，1，-180，180])；

set(gca，′XTickMode′，′manual′，′XTick′，[0.2，0.3，0.7，0.8])；

set(gca，′YTickMode′，′manual′，′YTick′，[-180，-90，0，90，180])；grid

subplot(2，2，3)，plot(w/pi，dbH，′k′)；

ylabel(′幅度(dB)′)；xlabel(′頻率′)；axis([0，1，-60，5])；

set(gca，′XTickMode′，′manual′，′XTick′，[0.2，0.3，0.7，0.8])；

set(gca，′YTickMode′，′manual′，′YTick′，[-80，-20，-1，0])；grid

subplot(2，2，4)，zplane(bd，ad)；

axis([-1.1，1.1，-1.1，1.1])；ylabel(′零極圖′)；(4)程序：

%數(shù)字濾波器指標(biāo)

wp=0.45*pi；

ws=0.3*pi；

Rp=1；As=40；

%轉(zhuǎn)換為模擬濾波器指標(biāo)

Fs=2000；T=1/Fs；

Omgp=(2/T)*tan(wp/2)；

Omgs=(2/T)*tan(ws/2)；

%模擬低通原型設(shè)計(jì)

[n，Omgc]=ellipord(Omgp，Omgs，Rp，As，′s′)

[z0，p0，k0]=ellipap(n，Rp，As)；

[b0，a0]=zp2tf(z0，p0，k0)；[ba，aa]=lp2hp(b0，a0，Omgc)；%模擬低通原型轉(zhuǎn)換為模擬高通濾波器

[bd，ad]=bilinear(ba，aa，F(xiàn)s)；

freqz(bd，ad)

figure；zplane(bd，ad)；

8.4實(shí)驗(yàn)四FIR數(shù)字濾波器設(shè)計(jì)

8.4.1實(shí)驗(yàn)?zāi)康?/p>

(1)掌握窗函數(shù)法、頻率采樣法及優(yōu)化設(shè)計(jì)法設(shè)計(jì)FIR數(shù)字濾波器的方法，學(xué)會(huì)相應(yīng)的MATLAB編程；

(2)熟悉線性相位FIR濾波器的幅頻、相位特性；

(3)了解各種不同窗函數(shù)對(duì)濾波器性能的影響。8.4.2實(shí)驗(yàn)原理與方法

1.窗函數(shù)設(shè)計(jì)法

1)設(shè)計(jì)原理

如果系統(tǒng)的單位脈沖響應(yīng)hd(n)為已知，則系統(tǒng)的輸入/輸出關(guān)系為

y(n)=x(n)*hd(n)

對(duì)于低通濾波器，只要設(shè)計(jì)出低通濾波器的單位脈沖響應(yīng)，就可以由上式得到系統(tǒng)的輸出了。

假設(shè)所希望的數(shù)字濾波器的頻率響應(yīng)為Hd(ejω)，它是頻域的周期函數(shù)，周期為2π，那么Hd(ejω)的傅立葉系數(shù)為以hd(n)為單位脈沖響應(yīng)的數(shù)字濾波器將具有頻域響應(yīng)Hd(ejω)。

但是將hd(n)作為FIR濾波器脈沖響應(yīng)有兩個(gè)問(wèn)題：

(1)它是無(wú)限長(zhǎng)的，與FIR濾波器脈沖響應(yīng)有限長(zhǎng)這一前提不一致；

(2)它是非因果的，hd(n)≠0，n<0。

常用的窗函數(shù)有矩形窗、哈明窗和布萊克曼窗等。

窗函數(shù)設(shè)計(jì)線性相位FIR濾波器的步驟如下：

(1)確定所要設(shè)計(jì)的濾波器的類(lèi)型和技術(shù)指標(biāo)。

(2)確定窗函數(shù)。查表5.2.1和表5.2.2，根據(jù)阻帶衰減指標(biāo)選擇窗函數(shù)，選擇原則是：在符合指標(biāo)要求情況下，選擇最簡(jiǎn)單的窗函數(shù)。

(3)確定目標(biāo)濾波器的頻率特性和hd(n)。以理想濾波器模型并考慮線性相位確定目標(biāo)濾波器的頻率特性函數(shù)Hd(ejω)，通過(guò)求反離散傅立葉變換得到目標(biāo)濾波器的單位脈沖響應(yīng)hd(n)。

(4)確定濾波器的最小階數(shù)。查表5.2.1和表5.2.2，根據(jù)技術(shù)指標(biāo)要求的過(guò)渡帶和所選用的窗函數(shù)的過(guò)渡帶寬度，計(jì)算出窗的寬度N，這也是濾波器的最小階數(shù)。

(5)確定濾波器階數(shù)。根據(jù)濾波器的類(lèi)型和線性相位的約束條件，選擇線性相位類(lèi)型，并以此確定濾波器的階數(shù)和群延時(shí)。

(6)計(jì)算h(n)=hd(n)WR(n)，即為所要設(shè)計(jì)濾波器的單位脈沖響應(yīng)。如果要求還可進(jìn)一步求出濾波器的系統(tǒng)函數(shù)。

(7)求H(ejw)。分析其幅頻特性，若不滿足要求，可適當(dāng)改變窗函數(shù)形式或長(zhǎng)度N，重復(fù)上述設(shè)計(jì)過(guò)程，以得到滿意的結(jié)果。

2)MATLAB函數(shù)

(1)矩形窗：w=boxcar(n)，可產(chǎn)生長(zhǎng)度為n點(diǎn)的矩形窗。

(2)哈明窗：w=hamming(n)，可產(chǎn)生長(zhǎng)度為n點(diǎn)的hamming窗。

(3)布萊克曼窗：w=blackman(n)，可產(chǎn)生長(zhǎng)度為n點(diǎn)的blackman窗。

(4)窗函數(shù)設(shè)計(jì)FIR濾波器函數(shù)：b＝fir1(n，Wn，′ftype′，Window)。使用的窗類(lèi)型由window

參數(shù)指定，省略時(shí)使用hamming窗；濾波器類(lèi)型由參數(shù)ftype指定：

當(dāng)省略ftype參數(shù)時(shí)，設(shè)計(jì)低通帶通FIR濾波器；當(dāng)ftype＝high時(shí)，設(shè)計(jì)高通FIR濾波器；

當(dāng)ftype＝stop時(shí)，設(shè)計(jì)帶阻FIR濾波器。

在設(shè)計(jì)高通和帶阻濾波器時(shí)，fir1函數(shù)總是使用偶對(duì)稱(chēng)N為奇數(shù)(即第一類(lèi)線性相位FIR濾波器)的結(jié)構(gòu)，因此當(dāng)輸入的階次為偶數(shù)時(shí)，fir1函數(shù)會(huì)自動(dòng)加1。

2.頻率采樣法

1)設(shè)計(jì)原理

設(shè)計(jì)濾波器時(shí)，通常給出的是幅頻特性的技術(shù)指標(biāo)要求，可直接在頻域進(jìn)行處理。按照理想的頻率特性H(ejω)，在ω＝0到2π之間等間隔采樣N點(diǎn)，得到然后用H(k)的傅里葉逆變換作為濾波器的系數(shù)：

b(n)＝h(n)＝IDFT[H(k)]

2)MATLAB函數(shù)

b＝fir2(n，f，m)；設(shè)計(jì)一個(gè)n階的FIR濾波器，其濾波器的頻率特性由矢量f和m決定。

參數(shù)f為頻率點(diǎn)矢量，且f∈[0，1]，當(dāng)f＝1時(shí)，相應(yīng)于0.5Fs。矢量f中按升序排列，且第一個(gè)必須為0，最后一個(gè)必須為1，并允許出現(xiàn)相同的頻率值。

參數(shù)m為與f相對(duì)應(yīng)的濾波器幅度響應(yīng)期望值。

注意：矢量f和m的長(zhǎng)度必須相同。

3.均方誤差最小化設(shè)計(jì)

在MATLAB中，采樣均方誤差最小化設(shè)計(jì)的函數(shù)是firls，使用格式為

h=Firls(M，F(xiàn)，A)

其中M為FIR數(shù)字濾波器的階數(shù)(濾波器的長(zhǎng)度為N＝M＋1)，數(shù)組f和A給出它的預(yù)期頻率響應(yīng)。算出的濾波器系數(shù)為h，其長(zhǎng)度為N。矢量f包含各邊緣頻率，單位為π，即

0.0≤f≤1.0。矢量A為各指定頻率上預(yù)期的幅度響應(yīng)。8.4.3實(shí)驗(yàn)內(nèi)容

(1)用窗函數(shù)法設(shè)計(jì)一個(gè)FIR數(shù)字低通濾波器，并驗(yàn)證設(shè)計(jì)結(jié)果的幅頻特性和相頻特性。要求：通帶截止頻率wp=0.3π，Rp=0.5dB，阻帶截止頻率ws=0.45π，Rr=50dB。

(2)用窗函數(shù)法設(shè)計(jì)一個(gè)FIR數(shù)字高通濾波器，并驗(yàn)證設(shè)計(jì)結(jié)果的幅頻特性和相頻特性。要求：通帶截止頻率wp=0.45π，Rp=0.5dB，阻帶截止頻率ws=0.3π，Rr=50dB。

(3)用窗函數(shù)法設(shè)計(jì)一個(gè)FIR數(shù)字帶通濾波器，并驗(yàn)證設(shè)計(jì)結(jié)果的幅頻特性和相頻特性。要求：低端阻帶截止頻率為0.2π，最小衰減為65dB；低端通帶截止頻率為0.3π，最大衰減為0.05dB。高端通帶截止頻率為0.7π，最大衰減為0.05dB；高端阻帶截止頻率為0.8π，最小衰減為65dB。

(4)用窗函數(shù)法設(shè)計(jì)一個(gè)FIR數(shù)字帶阻濾波器，并驗(yàn)證設(shè)計(jì)結(jié)果的幅頻特性和相頻特性。要求：低端通帶截止頻率為0.2π，最大衰減為0.05dB；低端阻帶截止頻率為0.3π，最小衰減為50dB；高端阻帶截止頻率為0.7π，最小衰減為50dB；高端通帶截止頻率為0.8π，最大衰減為0.05dB。

(5)用Matlab的fir1函數(shù)(窗函數(shù)設(shè)計(jì)法)設(shè)計(jì)實(shí)驗(yàn)(1)～(4)指標(biāo)的FIR數(shù)字濾波器。

(6)用fir2函數(shù)(頻率采樣法)設(shè)計(jì)實(shí)驗(yàn)(1)～(4)指標(biāo)的FIR數(shù)字濾波器。

(7)用Matlab的firls函數(shù)(均方誤差最小優(yōu)化法)設(shè)計(jì)實(shí)驗(yàn)(1)～(4)指標(biāo)的FIR數(shù)字濾波器。

(8)用凱塞窗設(shè)計(jì)一專(zhuān)用線性相位濾波器，濾波器的頻率特性如圖8.4.1所示。當(dāng)N=45，β=4、6、10時(shí)，分別設(shè)計(jì)并比較結(jié)果。

(9)用頻率采樣法設(shè)計(jì)如圖8.4.1所示的特性濾波器，采用過(guò)渡帶分別不插入過(guò)度點(diǎn)及插入一個(gè)過(guò)渡點(diǎn)0.3904，比較兩種設(shè)計(jì)結(jié)果。如要進(jìn)一步提高阻帶衰減，試優(yōu)化濾波器設(shè)計(jì)結(jié)果。

(10)采用雷米茲法設(shè)計(jì)如圖8.4.1所示的特性濾波器，并將設(shè)計(jì)結(jié)果與(8)、(9)的最好設(shè)計(jì)結(jié)果比較。

(11)使用雷米茲法設(shè)計(jì)一線性相位FIR帶通數(shù)字濾波器，并滿足如下技術(shù)指標(biāo)要求：對(duì)模擬信號(hào)進(jìn)行采樣的周期T=0.001s，在Ωr=[300π800π]弧度/秒處的衰減大于50dB，在Ωp=[500π，600π]弧度/秒處的衰減小于3dB。8.4.4思考題

(1)總結(jié)窗函數(shù)設(shè)計(jì)法的設(shè)計(jì)步驟，比較不同窗函數(shù)的設(shè)計(jì)結(jié)果。

(2)比較基于窗函數(shù)原理設(shè)計(jì)法、fir1、fir2、firls的設(shè)計(jì)結(jié)果及特點(diǎn)，說(shuō)明其應(yīng)用場(chǎng)合。

(3)采用窗函數(shù)法、頻率采樣法、雷米茲法設(shè)計(jì)時(shí)，濾波器階數(shù)如何確定？

(4)利用課外時(shí)間，用FDATool實(shí)現(xiàn)實(shí)驗(yàn)內(nèi)容所列濾波器設(shè)計(jì)，并與實(shí)驗(yàn)結(jié)果比較。8.4.5實(shí)驗(yàn)參考程序(MATLAB語(yǔ)言)

(1)參考程序：

wp=0.45*pi；ws=0.3*pi；

deltaw=wp-ws；

N0=ceil(6.6*pi/deltaw)；

N=N0+mod(N0+1，2)

windows=(hamming(N))′；

wc=(ws+wp)/2；

alfa=(N-1)/2；

n=0：N-1；

hd=sin(wc*(n-alfa+eps))./(pi*(n-alfa+eps))；%理想低通濾波器

b=hd.*windows；

freqz(b，1)；(2)參考程序：

wp=0.45*pi；ws=0.3*pi；

deltaw=wp-ws；

N0=ceil(6.6*pi/deltaw)；

N=N0+mod(N0+1，2)

windows=(hamming(N))′；

wc=(ws+wp)/2；

alfa=(N-1)/2；

n=0：N-1；

hd=sin(pi*(n-alfa+eps))./(pi*(n-alfa+eps))-sin(wc*(n-alfa+eps)).

/(pi*(n-alfa+eps))；%建立理想高通濾波器

b=hd.*windows；

freqz(b，1)；(3)參考程序：

wp1=0.3*pi；wp2=0.7*pi；

ws1=0.2*pi；ws2=0.8*pi；

deltaw=wp1-ws1；

N0=ceil(11*pi/deltaw)；

N=N0+mod(N0+1，2)

windows=(blackman(N))′；

wc1=(ws1+wp1)/2；wc2=(ws2+wp2)/2；

alfa=(N-1)/2；

n=0：N-1；

hd=sin(wc2*(n-alfa+eps))./(pi*(n-alfa+eps))-sin(wc1*(n-alfa+eps))./(pi*(n-alfa+eps))；%建立理想帶通濾波器

b=hd.*windows；

freqz(b，1)；(4)參考程序：

wp1=0.2*pi；wp2=0.8*pi；

ws1=0.3*pi；ws2=0.7*pi；

deltaw=ws1-wp1；

N0=ceil(6.6*pi/deltaw)；

N=N0+mod(N0+1，2)

windows=(hamming(N))′；

wc1=(ws1+wp1)/2；wc2=(ws2+wp2)/2；

alfa=(N-1)/2；

n=0：N-1；

hd=sin(wc1*(n-alfa+eps))./(pi*(n-alfa+eps))+sin(pi*(n-alfa+eps))./(pi*(n-alfa+eps))-sin(wc2*(n-alfa+eps))./(pi*(n-alfa+eps))；%建立理想帶阻濾波器

b=hd.*windows；

freqz(b，1)；(5)

5-1參考程序：

wp=0.3*pi；ws=0.45*pi；

deltaw=ws-wp；

N0=ceil(6.6*pi/deltaw)；

N=N0+mod(N0+1，2)

windows=(hamming(N+1))′；

wc=(ws+wp)/2；

b=fir1(N，wc/pi，windows)；

freqz(b，1)；5-2參考程序：

wp=0.45*pi；ws=0.3*pi；

deltaw=wp-ws；

N0=ceil(6.6*pi/deltaw)；

N=N0+mod(N0+1，2)

wc=(ws+wp)/2；

b=fir1(N，wc/pi，′high′)；

freqz(b，1)；5-3參考程序：

wp1=0.3*pi；wp2=0.7*pi；

ws1=0.2*pi；ws2=0.8*pi；

deltaw=wp1-ws1；

N0=ceil(11*pi/deltaw)；

N=N0+mod(N0+1，2)

wc1=(ws1+wp1)/2；wc2=(ws2+wp2)/2；

b=fir1(N，[wc1，wc2]/pi，′bandpass′，blackman(N+1))；

freqz(b，1)；5-4參考程序：

wp1=0.2*pi；wp2=0.8*pi；

ws1=0.3*pi；ws2=0.7*pi；

deltaw=ws1-wp1；

N0=ceil(6.6*pi/deltaw)；

N=N0+mod(N0+1，2)

wc1=(ws1+wp1)/2；wc2=(ws2+wp2)/2；

b=fir1(N，[wc1，wc2]/pi，′stop′)

freqz(b，1)；(6)參考程序

6-1參考程序：

wp=0.3*pi；ws=0.45*pi；

deltaw=ws-wp；

rp=0.5；rs=50；

N0=ceil(6.6*pi/deltaw)；

N=N0+mod(N0+1，2)

f=[0，wp/pi，ws/pi，1]；

A=[1，10^(-rp/20)，10^(-rs/20)，0]；

b=fir2(N-1，f，A)；

freqz(b，1)；6-2參考程序：

ws=0.3*pi；wp=0.45*pi；

deltaw=wp-ws；

rp=0.5；rs=50；

N0=ceil(6.6*pi/deltaw)；

N=N0+mod(N0+1，2)

f=[0，ws/pi，wp/pi，1]；

A=[0，10^(-rs/20)，10^(-rp/20)，1]

b=fir2(N-1，f，A)；

freqz(b，1)；6-3參考程序：

wp1=0.3*pi；wp2=0.7*pi；

ws1=0.2*pi；ws2=0.8*pi；

deltaw=wp1-ws1；

N0=ceil(11*pi/deltaw)；

N=N0+mod(N0+1，2)

rs=65；rp=0.05；

f=[0，ws1/pi，wp1/pi，wp2/pi，ws2/pi，1]；

A=[0，10^(-rs/20)，10^(-rp/20)，10^(-rp/20)，10^(-rs/20)，0]

b=fir2(N-1，f，A)；

freqz(b，1)；6-4參考程序：

wp1=0.2*pi；wp2=0.8*pi；

ws1=0.3*pi；ws2=0.7*pi；

deltaw=ws1-wp1；

N0=ceil(6.6*pi/deltaw)；

N=N0+mod(N0+1，2)

rs=65；rp=0.05；

f=[0，wp1/pi，ws1/pi，ws2/pi，wp2/pi，1]；

A=[1，10^(-rp/20)，10^(-rs/20)，10^(-rs/20)，10^(-rp/20)，1]

b=fir2(N-1，f，A)；

freqz(b，1)；(7)參考程序

7-1參考程序：

wp=0.3*pi；ws=0.45*pi；

deltaw=ws-wp；

rp=0.5；rs=50；

N0=ceil(6.6*pi/deltaw)；

N=N0+mod(N0+1，2)

f=[0，wp/pi，ws/pi，1]；

A=[1，10^(-rp/20)，10^(-rs/20)，0]；

b=firls(N-1，f，A)；

freqz(b