1、本科生實驗報告實驗課程 數(shù)字信號處理 學院名稱 信息科學與技術學院 專業(yè)名稱 物聯(lián)網(wǎng)工程 學生姓名 曹林鑫 學生學號 201413060301 指導教師 羅耀耀 實驗地點 6B607 實驗成績 二一六年 十月 二一六年 十二月實驗一 熟悉Matlab環(huán)境一、實驗目的1. 熟悉MATLAB的主要操作命令。2. 學會簡單的矩陣輸入和數(shù)據(jù)讀寫。3. 掌握簡單的繪圖命令。4. 用MATLAB編程并學會創(chuàng)建函數(shù)。5. 觀察離散系統(tǒng)的頻率響應。二、實驗內(nèi)容（1）數(shù)組的加、減、乘、除和乘方運算。輸入A=1 2 3 4，B=3 4 5 6，求C=A+B，D=A-B，E=A.*B，F(xiàn)=A./B，G=A.B并用s

2、tem語句畫出A、B、C、D、E、F、G。clear all;a=1 2 3 4; b=3 4 5 6;c=a+b;d=a-b;e=a.*b;f=a./b;g=a.b;n=1:4;subplot(4,2,1);stem(n,a);xlabel('n');xlim(0 5);ylabel('A');subplot(4,2,2);stem(n,b);xlabel('n');xlim(0 5);ylabel('B');subplot(4,2,3);stem(n,c);xlabel('n');xlim(0 5);ylabe

3、l('C');subplot(4,2,4);stem(n,d);xlabel('n');xlim(0 5);ylabel('D');subplot(4,2,5);stem(n,e);xlabel('n');xlim(0 5);ylabel('E');subplot(4,2,6);stem(n,f);xlabel('n');xlim(0 5);ylabel('F');subplot(4,2,7);stem(n,g);xlabel('n');xlim(0 5);ylabe

4、l('G');（2）用MATLAB實現(xiàn)下列序列：1) x(n)=0.8n 0n15 2) x(n)=e(0.2+3j)n 0n15 3) x(n)=3cos(0.125n+0.2)+2sin(0.25n+0.1) 0n15 4) 將3)中的x(n)擴展為以16為周期的函數(shù)x16(n)=x(n+16),繪出四個周期。5) 將3)中的x(n)擴展為以10為周期的函數(shù)x10(n)=x(n+10),繪出四個周期。解：1) x(n)=0.8n 0n15clear all;xa=0.8.N;figure;subplot(2,1,1);stem(N,xa);xlabel('n'

5、;);xlim(0 16);ylabel('xa'); 2) x(n)=e(0.2+3j)n 0n15xb=exp(0.2+3*j)*N);subplot(2,1,2);stem(N,xb);xlabel('n');xlim(0 16);ylabel('xb');figure; 3) x(n)=3cos(0.125n+0.2)+2sin(0.25n+0.1) 0n15xc=3*cos(0.125*pi*N+0.2*pi)+2*sin(0.25*pi*N+0.1*pi);subplot(3,1,1);stem(N,xc);xlabel('n

6、');xlim(0 16);ylabel('xc');4) 將3)中的x(n)擴展為以16為周期的函數(shù)x16(n)=x(n+16),繪出四個周期。 k=0:3;m=0;for i=1:4 for j=1:16 m=m+1; n(m)=N(j)+16*k(i); x16(m)=3*cos(0.125*pi*n(m)+0.2*pi)+2*sin(0.25*pi*n(m)+0.1*pi); endend subplot(3,1,2);stem(n,x16);xlabel('n');ylabel('x16');5) 將3)中的x(n)擴展為以10

7、為周期的函數(shù)x10(n)=x(n+10),繪出四個周期。 for j=1:10 x10(j)=x16(j); end for i=1:3 for m=1:10 x10(i*10+m)=x10(m); end endn=1:40;subplot(3,1,3);stem(n,x10);xlabel('n');ylabel('x10');（3）x(n)=1,-1,3,5,產(chǎn)生并繪出下列序列的樣本: 1) x1(n)=2x(n+2)-x(n-1)-2x(n)2) clear alln=1:4;T=4;x=1 -1 3 5;x(5:8)=x(1:4);subplot(2,

8、1,1);stem(1:8,x);grid;for i=1:4 if i-1<0 x1(i)=2*x(i+2)-x(i-1)-2*x(i); else x1(i)=2*x(i+2)-x(i-1+T)-2*x(i); endendx1(5:8)=x1(1:4);subplot(2,1,2);stem(1:8,x1);grid; （4）繪出下列時間函數(shù)的圖形，對x軸、y軸以及圖形上方均須加上適當?shù)臉俗ⅲ?a) x(t)=sin(2t) 0t10sb) x(t)=cos(100t)sin(t) 0t4sta=0:0.05:10;xa=sin(2*pi*ta);subplot(2,1,1);pl

9、ot(ta,xa);xlabel('t');ylabel('幅度');tb=0:0.01:4;xb=cos(100*pi*tb).*sin(pi*tb);subplot(2,1,2);plot(tb,xb);xlabel('t');ylabel('幅度'); （5）編寫函數(shù)stepshift(n0,n1,n2)實現(xiàn)u(n-n0),n1<n0<n2,繪出該函數(shù)的圖形，起點為n1，終點為n2。n0=5;ns=1;nf=10;%ns為起點;nf為終點；在=n=n0處生成單位階躍序列n=ns:nf;x=(n-n0)>=0

10、;stem(n,x); （6）給一定因果系統(tǒng)求出并繪制H(z)的幅頻響應與相頻響應。clear all;b=1,sqrt(2),1;a=1,-0.67,0.9;h,w=freqz(b,a);am=20*log10(abs(h);subplot(2,1,1);plot(w,am);ph=angle(h);subplot(2,1,2);plot(w,ph); （7）計算序列8 -2 -1 2 3和序列2 3 -1 -3的離散卷積，并作圖表示卷積結果。clear all;a=8 -2 -1 2 3;b=2 3 -1 -3;c=conv(a,b); %計算卷積M=length(c)-1;n=0:1:M

11、;stem(n,c);xlabel('n');ylabel('幅度'); （8）求以下差分方程所描述系統(tǒng)的單位脈沖響應h(n),0n50 y(n)+0.1y(n-1)-0.06y(n-2)=x(n)-2x(n-1)clear all;N=50;a=1 -2;b=1 0.1 -0.06;x=1 zeros(1,N-1);k=0:1:N-1;y=filter(a,b,x);stem(k,y);xlabel('n');ylabel('幅度 '); 3、 思考題1. 比較實驗內(nèi)容第2題中（4）和（5）兩小題的結果，試說明對于周期性信號，應

12、當如何采樣，才能保證周期擴展后與原信號保持一致？ 答：對于周期性信號，在進行采樣時，其采樣周期必須滿足采樣定理，即采樣頻率應該大于信號最高頻率的兩倍，這樣才能避免迭混，以便采樣后仍能準確的恢復原信號。2. 對于有限長序列，如何用MATLAB計算其DTFT？ 答：用函數(shù)Freqz可以計算序列在給定的離散頻率點上的DTFT，函數(shù)freqz返回的頻率響應值為向量H。在H = freqz(num,den,w)中，0到&#960；之間指定的頻率集由向量w給出。freqz函數(shù)的自變量k就是頻率點的總數(shù)。3. 對于由兩個子系統(tǒng)級聯(lián)或并聯(lián)的系統(tǒng)，如何用MATLAB計算它們的幅頻響應與相頻響應？ 系統(tǒng)的

13、級聯(lián)或并聯(lián)實現(xiàn)涉及到了因式分解。在MATLAB中，我們可以用函數(shù)roots來實現(xiàn)多項式的因式分解。 實驗二 快速Fourier變換(FFT)及其應用一、實驗目的1. 在理論學習的基礎上，通過本實驗，加深對FFT的理解，熟悉MATLAB中的有關函數(shù)。2. 熟悉應用FFT對典型信號進行頻譜分析的方法。3. 了解應用FFT進行信號頻譜分析過程中可能出現(xiàn)的問題，以便在實際中正確應用FFT。4. 熟悉應用FFT實現(xiàn)兩個序列的線性卷積和相關的方法。二、實驗原理與方法  在各種信號序列中，有限長序列信號處理占有很重要地位，對有限長序列，我們可以使用離散Fouier變換(DFT)。這一變換不但可以很

14、好的反映序列的頻譜特性，而且易于用快速算法在計算機上實現(xiàn)，當序列x(n)的長度為N時，它的DFT定義為：反變換為：    有限長序列的DFT是其Z變換在單位圓上的等距采樣，或者說是序列Fourier變換的等距采樣，因此可以用于序列的譜分析。    FFT并不是與DFT不同的另一種變換，而是為了減少DFT運算次數(shù)的一種快速算法。它是對變換式進行一次次分解，使其成為若干小點數(shù)的組合，從而減少運算量。常用的FFT是以2為基數(shù)的，其長度 。它的效率高，程序簡單，使用非常方便，當要變換的序列長度不等于2的整數(shù)次方時，為了使用以2為基數(shù)的FFT，

15、可以用末位補零的方法，使其長度延長至2的整數(shù)次方。     (一)、在運用DFT進行頻譜分析的過程中可能產(chǎn)生三種誤差：    (1)    混疊    序列的頻譜時被采樣信號的周期延拓，當采樣速率不滿足Nyquist定理時，就會發(fā)生頻譜混疊，使得采樣后的信號序列頻譜不能真實的反映原信號的頻譜。避免混疊現(xiàn)象的唯一方法是保證采樣速率足夠高，使頻譜混疊現(xiàn)象不致出現(xiàn)，即在確定采樣頻率之前，必須對頻譜的性質(zhì)有所了解，在一般情況下，為了保證高于折疊頻率的分量不會出現(xiàn)，在采樣前，先用低通模

16、擬濾波器對信號進行濾波。    (2)    泄漏    實際中我們往往用截短的序列來近似很長的甚至是無限長的序列，這樣可以使用較短的DFT來對信號進行頻譜分析，這種截短等價于給原信號序列乘以一個矩形窗函數(shù)，也相當于在頻域將信號的頻譜和矩形窗函數(shù)的頻譜卷積，所得的頻譜是原序列頻譜的擴展。泄漏不能與混疊完全分開，因為泄漏導致頻譜的擴展，從而造成混疊。為了減少泄漏的影響，可以選擇適當?shù)拇昂瘮?shù)使頻譜的擴散減至最小。    (3)    柵欄效應

17、0;   DFT是對單位圓上Z變換的均勻采樣，所以它不可能將頻譜視為一個連續(xù)函數(shù)，就一定意義上看，用DFT來觀察頻譜就好像通過一個柵欄來觀看一個圖景一樣，只能在離散點上看到真實的頻譜，這樣就有可能發(fā)生一些頻譜的峰點或谷點被“尖樁的柵欄”所攔住，不能別我們觀察到。    減小柵欄效應的一個方法就是借助于在原序列的末端填補一些零值，從而變動DFT的點數(shù)，這一方法實際上是人為地改變了對真實頻譜采樣的點數(shù)和位置，相當于搬動了每一根“尖樁柵欄”的位置，從而使得頻譜的峰點或谷點暴露出來。    (二)、用FFT計算線性卷積&#

18、160;   用FFT可以實現(xiàn)兩個序列的圓周卷積。在一定的條件下，可以使圓周卷積等于線性卷積。一般情況，設兩個序列的長度分別為N1和N2，要使圓周卷積等于線性卷積的充要條件是FFT的長度NN1N2對于長度不足N的兩個序列，分別將他們補零延長到N。    當兩個序列中有一個序列比較長的時候，我們可以采用分段卷積的方法。有兩種方法：重疊相加法。將長序列分成與短序列相仿的片段，分別用FFT對它們作線性卷積，再將分段卷積各段重疊的部分相加構成總的卷積輸出。  重疊保留法。這種方法在長序列分段時，段與段之間保留有互相重疊的部分，在構成總的卷積輸

19、出時只需將各段線性卷積部分直接連接起來，省掉了輸出段的直接相加。     (三)、用周期圖法(平滑周期圖的平均法)對隨機信號作譜分析    實際中許多信號往往既不具有有限能量，由非周期性的。無限能量信號的基本概念是隨機過程，也就是說無限能量信號是一隨機信號。周期圖法是隨機信號作譜分析的一種方法，它特別適用于用FFT直接計算功率譜的估值。    將長度為N的實平穩(wěn)隨機序列的樣本x(n)再次分割成K段，每段長度為L,即L=N/K。每段序列仍可表示為：xi(n)=x(n+(i-1)L)，0nL-1，1iK 但是

20、這里在計算周期圖之前，先用窗函數(shù)w(n)給每段序列xi(n)加權，K個修正的周期圖定義為其中U表示窗口序列的能量，它等于在此情況下，功率譜估計量可表示為 三、實驗內(nèi)容及步驟    實驗中用到的信號序列：a)       Gaussian序列b)      衰減正弦序列c)      三角波序列d)      反三角波序列  &#

21、160; 上機實驗內(nèi)容：(1)、觀察高斯序列的時域和幅頻特性，固定信號xa(n)中參數(shù)p=8，改變q的值，使q分別等于2，4，8，觀察它們的時域和幅頻特性，了解當q取不同值時，對信號序列的時域幅頻特性的影響；固定q=8，改變p，使p分別等于8，13，14，觀察參數(shù)p變化對信號序列的時域及幅頻特性的影響，觀察p等于多少時，會發(fā)生明顯的泄漏現(xiàn)象，混疊是否也隨之出現(xiàn)？記錄實驗中觀察到的現(xiàn)象，繪出相應的時域序列和幅頻特性曲線。% %定義高斯序列% function Xa,Fa =gauss(p,q)% n=0:15;% Xa(n+1)=exp(-(n+1-p).2./q);% F=fft(Xa);%

22、Fa=abs(F);clear all;% p=8,q=2 %Xa1,Fa1= gauss(8,2);k=0:15;subplot(5,2,1);plot(k,Xa1);xlabel('n');ylabel('時域特性');text(10,0.5,'p=8,q=2');subplot(5,2,2);plot(k,Fa1);xlabel('n');ylabel('幅頻特性');text(8,3,'p=8,q=2');% p=8,q=4 %Xa2,Fa2= gauss(8,4);subplot(5,2,

23、3);plot(k,Xa2);xlabel('n');ylabel('時域特性');text(10,0.5,'p=8,q=4');subplot(5,2,4);plot(k,Fa2);xlabel('n');ylabel('幅頻特性');text(8,3,'p=8,q=4');% p=8,q=8 %Xa3,Fa3= gauss(8,8);subplot(5,2,5);plot(k,Xa3);xlabel('n');ylabel('時域特性');text(10,0.5,

24、'p=8,q=8');subplot(5,2,6);plot(k,Fa3);xlabel('n');ylabel('幅頻特性');text(8,3,'p=8,q=8');% p=13,q=8 %Xa4,Fa4= gauss(13,8);subplot(5,2,7);plot(k,Xa4);xlabel('n');ylabel('時域特性');text(10,0.5,'p=13,q=8');subplot(5,2,8);plot(k,Fa4);xlabel('n');y

25、label('幅頻特性');text(8,3,'p=13,q=8');% p=14,q=8 %Xa5,Fa5= gauss(14,8);subplot(5,2,9);plot(k,Xa5);xlabel('n');ylabel('時域特性');text(10,0.5,'p=14,q=8');subplot(5,2,10);plot(k,Fa5);xlabel('n');ylabel('幅頻特性');text(8,3,'p=14,q=8');(2)、觀察衰減正弦序列xb

26、(n)的時域和幅頻特性，a=0.1，f=0.0625，檢查譜峰出現(xiàn)位置是否正確，注意頻譜的形狀，繪出幅頻特性曲線，改變f，使f分別等于0.4375和0.5625，觀察這兩種情況下，頻譜的形狀和譜峰出現(xiàn)位置，有無混疊和泄漏現(xiàn)象？說明產(chǎn)生現(xiàn)象的原因。% 定義衰減正弦序列% function Xb,Fb = downsin(a,f)% n=0:15;% Xb(n+1)=exp(-a.*n).*sin(2*pi*f.*n); %自然對數(shù)的底:e=:2.71828 18284 59045 23536 % F = fft(Xb);% Fb=abs(F);clear all;k=0:15;% a=0.1,f

27、=0.0.0625 %Xb,Fb=downsin(0.1,0.0625);subplot(3,2,1); plot(k,Xb);xlabel('n');ylabel('時域特性');text(8,0.5,'a=0.1,f=0.0625');subplot(3,2,2); plot(k,Fb);xlabel('n');ylabel('幅值特性');text(10,3,'a=0.1,f=0.0625');% a=0.1,f=0.4375 %Xb1,Fb1=downsin(0.1,0.4375);subp

28、lot(3,2,3); plot(k,Xb1);xlabel('n');ylabel('時域特性');text(8,0.5,'a=0.1,f=0.4375');subplot(3,2,4); plot(k,Fb1);xlabel('n');ylabel('幅值特性');text(10,3,'a=0.1,f=0.4375');% a=0.1,f= 0.5625 %Xb2,Fb2=downsin(0.1,0.5625);subplot(3,2,5); plot(k,Xb2);xlabel('n&

29、#39;);ylabel('時域特性');text(8,0.5,'a=0.1,f=0.5625');subplot(3,2,6); plot(k,Fb2);xlabel('n');ylabel('幅值特性');text(10,3,'a=0.1,f=0.5625');    (3)、觀察三角波和反三角波序列的時域和幅頻特性，用N=8點FFT分析信號序列xc(n)和xd(n)的幅頻特性，觀察兩者的序列形狀和頻譜曲線有什么異同？繪出兩序列及其幅頻特性曲線。在xc(n)和xd(n)末尾補零，用

30、N=16點FFT分析這兩個信號的幅頻特性，觀察幅頻特性發(fā)生了什么變化？兩情況的FFT頻譜還有相同之處嗎？這些變化說明了什么？clear all;n=0:3;k=1:8;%定義三角波序列Xc(n+1) = n;Xc(n+5) =4-n;%定義反三角波序列Xd(n+1) = 4-n;Xd(n+5) =n;% 三角波特性 %subplot(2,2,1);plot(k-1,Xc);xlabel('n');ylabel('時域特性');text(1,3,'三角波');subplot(2,2,2);plot(k-1,abs(fft(Xc);xlabel(&#

31、39;k');ylabel('幅頻特性');text(4,10,'三角波');% 反三角波特性 %subplot(2,2,3);plot(k-1,Xd);xlabel('n');ylabel('時域特性');text(3,3,'反三角波');subplot(2,2,4);plot(k-1,abs(fft(Xd);xlabel('k');ylabel('幅頻特性');text(4,10,'反三角波');%末尾補0，計算32點FFTXc(9:32)=0;Xd(9:

32、32)=0;k=1:32;figure;% 三角波特性 %subplot(2,2,1);plot(k-1,Xc);xlabel('n');ylabel('時域特性');text(1,3,'三角波');subplot(2,2,2);plot(k-1,abs(fft(Xc);xlabel('k');ylabel('幅頻特性');text(4,10,'三角波');% fan三角波特性 %subplot(2,2,3);plot(k-1,Xd);xlabel('n');ylabel('

33、時域特性');text(3,3,'反三角波');subplot(2,2,4);plot(k-1,abs(fft(Xd);xlabel('k');ylabel('幅頻特性');text(4,10,'反三角波');    (4)、一個連續(xù)信號含兩個頻率分量，經(jīng)采樣得x(n)=sin2*0.125n+cos2*(0.125+f)n n=0,1,N-1已知N=16，f分別為1/16和1/64，觀察其頻譜；當N=128時，f不變，其結果有何不同，為什么？clear all;% N = 16 %N=16;

34、detf=1/16;n=0:N-1;x1(n+1)=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detf).*n);detf = 1/64;x2(n+1)=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detf).*n);% N = 16,detf = 1/16 %subplot(2,2,1);stem(n,x1);hold;plot(n,x1);xlabel('n');ylabel('時域特性');text(6,1,'N=16,detf=1/16');subplot(2,2,2);stem(n,

35、abs(fft(x1);xlabel('n');ylabel('幅值特性');text(6,4,'N=16,detf=1/16');% N = 16,detf = 1/64 %subplot(2,2,3);stem(n,x2);xlabel('n');ylabel('時域特性');text(6,1,'N=16,detf=1/64');subplot(2,2,4);stem(n,abs(fft(x2);xlabel('n');ylabel('幅值特性');text(6,

36、4,'N=16,detf=1/64');% N = 128 %N=128;detf=1/16;n=0:N-1;x3(n+1)=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detf).*n);detf = 1/64;x4(n+1)=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detf).*n);% N = 128,detf = 1/16 %figure;subplot(2,2,1);stem(n,x3);xlabel('n');ylabel('時域特性');axis(0 128 -2 2);tex

37、t(6,1.5,'N=128,detf=1/16');subplot(2,2,2);stem(n,abs(fft(x3);xlabel('n');ylabel('幅值特性');axis(0 128 -10 70);text(40,60,'N=128,detf=1/16');% N = 128,detf = 1/64 %subplot(2,2,3);stem(n,x3);xlabel('n');ylabel('時域特性');axis(0 128 -2 2);text(6,1.5,'N=128,

38、detf=1/16');subplot(2,2,4);stem(n,abs(fft(x4);xlabel('n');ylabel('幅值特性');axis(0 128 -10 70);text(40,60,'N=128,detf=1/16');(5)、用FFT分別實現(xiàn)xa(n)（p8，q2）和 xb(n)（a0.1，f0.0625）的16點圓周卷積和線性卷積。clear all;N=16;n=0:N-1;p=8;q=2;Xa(n+1)=exp(-(n-p).2./q);a=0.1;f=0.0625;Xb(n+1)=exp(-a.*n).*

39、sin(2*pi*f.*n);%16點循環(huán)卷積Fa=fft(Xa); Fb=fft(Xb);Fx=Fa.*Fb;X51=ifft(Fx);stem(n,X51);%16點線性卷積Xa(N+1:2*N-1)=0;Xb(N+1:2*N-1)=0;Fa=fft(Xa); Fb=fft(Xb);Fc=Fa.*Fb;X52=ifft(Fc);figure;stem(1:2*N-1,X52);（6）用FFT分別計算xa(n)(p=8,q=2)和xb(n)(a=0.1,f=0.0625)的16點循環(huán)相關和線性相關，問一共有多少種結果，他們之間有何異同點。clear all;N=16;n=0:N-1;p=8;

40、q=2;Xa(n+1)=exp(-(n-p).2./q);a=0.1;f=0.0625;Xb(n+1)=exp(-a.*n).*sin(2*pi*f.*n); N=length(Xa); %16點循環(huán)相關Fa=fft(Xa,2*N); Fb=fft(Xb,2*N);Fx=conj(Fa).*Fb;X71=real(ifft(Fx);X71=X71(N+2:2*N) X71(1:N);n=(-N+1):(N-1); stem(n,X71);%16點線性相關Xa(N+1:2*N-1)=0;Xb(N+1:2*N-1)=0;Fa=fft(Xa); Fb=fft(Xb);Fc=conj(Fa).*Fb;

41、X72=real(ifft(Fc);figure;stem(1:2*N-1,X72);（7）用FFT分別計算xa(n)(p=8,q=2)和xb(n)(a=0.1,f=0.0625)的自相關函數(shù)。clear all;N=16;n=0:N-1;p=8;q=2;Xa(n+1)=exp(-(n-p).2./q);a=0.1;f=0.0625;Xb(n+1)=exp(-a.*n).*sin(2*pi*f.*n); %自然對數(shù)的底:e=:2.71828 18284 59045 23536 N=length(Xa); % Xa(n) 16點自相關Fa=fft(Xa,2*N); Fb=fft(Xb,2*N);

42、F1=conj(Fa).*Fa;X81=real(ifft(F1);X81=X81(N+2:2*N) X81(1:N);n=(-N+1):(N-1); subplot(2,1,1);stem(n,X81);xlabel('n'); ylabel('幅度');% Xb(n) 16點自相關Fb=fft(Xb,2*N);F2=conj(Fb).*Fb;X82=real(ifft(F2);X82=X82(N+2:2*N) X82(1:N);% n=(-N+1):(N-1); subplot(2,1,2);stem(n,X82);xlabel('n');

43、ylabel('幅度');四、思考題 實驗中的信號序列和，在單位圓上的Z變換頻譜和會相同嗎？如果不同，說出哪一個低頻分量更多一些，為什么？答：xc和xd的幅頻和相頻完全一樣。兩者為同一個信號序列。只是采樣周期的起始位置不同，相差相位180度而已。但由FFT對其所得的圖像去不同。 因此，兩者在單位圓上的z變換頻譜完全一樣；并且直流分量大小完全一樣。實驗三 IIR數(shù)字濾波器的設計一、實驗目的 1. 掌握雙線性變換法及脈沖響應不變法設計IIR數(shù)字濾波器的具體設計方法及其原理，熟悉用雙線性變換法及脈沖響應不變法設計低通、高通和帶通IIR數(shù)字濾波器的MATLAB編程。 2. 觀

44、察雙線性變換及脈沖響應不變法設計的濾波器的頻域特性，了解雙線性變換法及脈沖響應不變法的特點。 3. 熟悉Butterworth濾波器、Chebyshev濾波器和橢圓濾波器的頻率特性。 2、 實驗原理 1 脈沖響應不變法用數(shù)字濾波器的單位脈沖響應序列模仿模擬濾波器的沖激響應,讓正好等于的采樣值，即，其中為采樣間隔，如果以及分別表示的拉式變換及的Z變換，則2雙線性變換法    S平面與z平面之間滿足以下映射關系：s平面的虛軸單值地映射于z平面的單位圓上，s平面的左半平面完全映射到z平面的單位圓內(nèi)。雙線性變換不存在混疊問題。    雙線性變

45、換是一種非線性變換 ，這種非線性引起的幅頻特性畸變可通過預畸而得到校正。    IIR低通、高通、帶通數(shù)字濾波器設計采用雙線性原型變換公式：變換類型           變換關系式     備   注 低通 高通   帶通 ：帶通的上下邊帶臨界頻率 以低通數(shù)字濾波器為例，將設計步驟歸納如下：1. 確定數(shù)字濾波器的性能指標：通帶臨界頻率、阻帶臨界頻率、通帶波動、阻帶內(nèi)的最小衰減、采樣周期、采樣頻率 ；2.  

46、確定相應的數(shù)字角頻率，； 3.  計算經(jīng)過預畸的相應模擬低通原型的頻率，；4. 根據(jù)c和r計算模擬低通原型濾波器的階數(shù)N，并求得低通原型的傳遞函數(shù)； 5. 用上面的雙線性變換公式代入，求出所設計的傳遞函數(shù)； 6.  分析濾波器特性，檢查其指標是否滿足要求。 三、實驗內(nèi)容：(1)，,;設計一切比雪夫高通濾波器，觀察其通帶損耗和阻帶衰減是否滿足要求。wc=2*1000*tan(2*pi*300/(2*1000);wt=2*1000*tan(2*pi*200/(2*1000); %對臨界頻率進行預畸變；N,wn=cheb1ord(wc,wt,0.8,20,'s')

47、; %給定通帶邊界頻率wc和阻帶邊界頻率wr、通帶波動0.8dB和最小阻帶衰減20dB，求滿足指標的模擬濾波器的最低階數(shù)N和通帶邊界頻率wn；B,A=cheby1(N,0.8,wn,'high','s'); %給定模擬濾波器的最低階數(shù)N、通帶邊界頻率wn和通帶波動0.8dB，設計N階模擬巴特沃斯高通濾波器，B和A分別表示系統(tǒng)函數(shù)的分子和分母多項式的系數(shù)；num,den=bilinear(B,A,1000); %給定模擬濾波器系統(tǒng)函數(shù)H(s)=B(s)/A(s)和采樣頻率1000Hz，根據(jù)雙線性變換法求出數(shù)字濾波器的系統(tǒng)函數(shù)H(z)=B(z)/A(z)，num和d

48、en分別為數(shù)字濾波器的分子分母多項式系數(shù)；h,w=freqz(num,den); %計算以(num,den)為參數(shù)的濾波器的頻率向量w和復頻率響應向量h；f=w/pi*500;plot(f,20*log10(abs(h);axis(0,500,-80,10);grid;xlabel('頻率/Hz');ylabel('幅度/dB');分析：f=200Hz時阻帶衰減大于30dB，通過修改axis(0,fs/2,-80,10)為axis(200,fs/2,-1,1)發(fā)現(xiàn)通帶波動rs滿足<0.8。（2），,;分別用脈沖響應不變法及雙線性變換法設計一巴特沃思數(shù)字低通

49、濾波器，觀察所設計數(shù)字濾波器的幅頻特性曲線，記錄帶寬和衰減量，檢查是否滿足要求。比較這兩種方法的優(yōu)缺點。Wc=2*pi*200;Wr=2*pi*300;Rp=1;Rs=25;N,Wn=buttord(Wc,Wr,Rp,Rs,s); B,A=butter(N,Wn,s); num1,den1=impinvar(B,A,1000); h1,w=freqz(num1,den1);B,A=butter(N,Wn,s);num2,den2=bilinear(B,A,1000); h2,w=freqz(num2,den2);f=w/pi*500;plot(f,abs(h1),-.,f,abs(h2),-)

50、;grid;xlabel(頻率/Hz);ylabel(幅度);bz1 =0.0000 0.0002 0.0153 0.0995 0.1444 0.0611 0.0075 0.0002 0.0000 0az1 =1.0000 -1.9199 2.5324 -2.2053 1.3868 -0.6309 0.2045 -0.0450 0.0060 -0.0004因此脈沖響應不變法的系統(tǒng)函數(shù)為:bz2 =0.0179 0.1072 0.2681 0.3575 0.2681 0.1072 0.0179az2 =1.0000 -0.6019 0.9130 -0.2989 0.1501 -0.0208 0.

51、0025因此雙線性變換法的系統(tǒng)函數(shù)為:分析：脈沖響應不變法的N=9，雙線性變換法的N=6，由圖知它們都滿足要求，但脈沖響應的衰減較快，雙線性變換的過渡帶窄一些，且階數(shù)比脈沖小，容易實現(xiàn)。（3）利用雙線性變換法分別設計滿足下列指標的巴特沃思型、切比雪夫型數(shù)字低通濾波器，并作圖驗證設計結果：fc=1.2kHz，0.5dB，fr2kHz, At40dB,fs=8kHz。比較這三種濾波器的階數(shù)。wc=2*1000*tan(2*pi*1200/(2*8000);wt=2*1000*tan(2*pi*2000/(2*8000); %對臨界頻率進行預畸變；N,wn= buttord (wc,wt,0.5,4

52、0,'s'); %給定通帶邊界頻率wc和阻帶邊界頻率wr、通帶波動0.5dB和最小阻帶衰減40dB，求滿足指標的模擬濾波器的最低階數(shù)N和通帶邊界頻率wn；B,A= butter (N, wn, 's'); %給定模擬濾波器的最低階數(shù)N、通帶邊界頻率wn和通帶波動0.5dB，設計N階模擬巴特沃斯低通濾波器，B和A分別表示系統(tǒng)函數(shù)的分子和分母多項式的系數(shù)；num,den=bilinear(B,A,8000); %給定模擬濾波器系統(tǒng)函數(shù)H(s)=B(s)/A(s)和采樣頻率8000Hz，根據(jù)雙線性變換法求出數(shù)字濾波器的系統(tǒng)函數(shù)H(z)=B(z)/A(z)，num和de

53、n分別為數(shù)字濾波器的分子分母多項式系數(shù)；h,w=freqz(num,den); %計算以(num,den)為參數(shù)的濾波器的頻率向量w和復頻率響應向量h；f=w/pi*500;plot(f,20*log10(abs(h);axis(0,400,-200,10);grid;xlabel('頻率/Hz');ylabel('幅度/dB');title（巴特沃思濾波器）;wc=2*1000*tan(2*pi*1200/(2*8000);wt=2*1000*tan(2*pi*2000/(2*8000); %對臨界頻率進行預畸變；N,wn=cheb1ord(wc,wt,0.5

54、,40,'s'); %給定通帶邊界頻率wc和阻帶邊界頻率wr、通帶波動0.5dB和最小阻帶衰減40dB，求滿足指標的模擬濾波器的最低階數(shù)N和通帶邊界頻率wn；B,A=cheby1(N,0.5,wn, 's'); %給定模擬濾波器的最低階數(shù)N、通帶邊界頻率wn和通帶波動0.5dB，設計N階模擬切比雪夫低通濾波器，B和A分別表示系統(tǒng)函數(shù)的分子和分母多項式的系數(shù)；num,den=bilinear(B,A,8000); %給定模擬濾波器系統(tǒng)函數(shù)H(s)=B(s)/A(s)和采樣頻率8000Hz，根據(jù)雙線性變換法求出數(shù)字濾波器的系統(tǒng)函數(shù)H(z)=B(z)/A(z)，num

55、和den分別為數(shù)字濾波器的分子分母多項式系數(shù)；h,w=freqz(num,den); %計算以(num,den)為參數(shù)的濾波器的頻率向量w和復頻率響應向量h；f=w/pi*500;plot(f,20*log10(abs(h);axis(0,400,-200,10);grid;xlabel('頻率/Hz');ylabel('幅度/dB');title（切比雪夫濾波器）;分析：切比雪夫濾波器的階數(shù)為5；巴特沃思濾波器的階數(shù)為9。（4）分別用脈沖響應不變法及雙線性不變法設計一巴特沃思型數(shù)字帶通濾波器，已知fs=30kHz，其等效的模擬濾波器指標為3dB，2kHzf3k

56、Hz；At5dB，f6kHz；At20dB，f1.5kHz。w1=2*30000*tan(2*pi*2000/(2*30000);w2=2*30000*tan(2*pi*3000/(2*30000);wr=2*30000*tan(2*pi*6000/(2*30000); %對臨界頻率進行預畸變；N,wn=buttord(w1 w2,1 wr,3,5,'s'); %給定通帶邊界頻率帶w1 w2和阻帶邊界頻率帶1 wr、通帶波動3dB和最小阻帶衰減5dB，求滿足指標的模擬濾波器的最低階數(shù)N和3dB邊界頻率wn；B,A=butter(N,wn,'s'); %給定模擬濾

57、波器的最低階數(shù)N和3dB邊界頻率wn，設計N階模擬巴特沃斯低通濾波器，B和A分別表示系統(tǒng)函數(shù)的分子和分母多項式的系數(shù)；num,den=bilinear(B,A,30000); %給定模擬濾波器系統(tǒng)函數(shù)H(s)=B(s)/A(s)和采樣頻率30000Hz，根據(jù)雙線性變換法求出數(shù)字濾波器的系統(tǒng)函數(shù)H(z)=B(z)/A(z)，num和den分別為數(shù)字濾波器的分子分母多項式系數(shù)；h,w=freqz(num,den); %計算以(num,den)為參數(shù)的濾波器的頻率向量w和復頻率響應向量h；f=w/pi*200;plot(f,20*log10(abs(h);axis(-40,40,-30,10); %給定幅頻響應的頻率和幅度的范圍；grid;xlabel('頻率/kHz');ylabel('幅度/dB');4、 思考題1. 雙線性變換法中和之間的關系是非線性的，在實驗中你注意到這種非線性關系了嗎？從那幾種數(shù)字濾波器的幅頻特性曲線中可以觀察到這種非線性關系？ 答：在雙線性變化法中，模擬頻率與數(shù)字頻率不再是線性關系，所以一個線性相位模擬濾波器經(jīng)雙線性變換后，得到的數(shù)字濾波器不再保持原有的線性相位了，在每一幅使用了雙線性變換的圖中，可以看到在采樣頻率一半處，幅度為零，這