數(shù)字信號處理實驗指導書

1、實驗1 時域離散信號的產(chǎn)生一、實驗目的l 學會運用MATLAB產(chǎn)生常用離散時間信號。二、實驗涉及的matlab子函數(shù)1、square功能：產(chǎn)生矩形波調用格式：x=square(t)；類似于sin（t），產(chǎn)生周期為2*pi，幅值為+1的方波。x=square(t，duty)；產(chǎn)生制定周期的矩形波，其中duty用于指定脈沖寬度與整個周期的比例。2、rand功能：產(chǎn)生rand隨機信號。調用格式：x=rand（n，m）；用于產(chǎn)生一組具有n行m列的隨機信號。三、實驗原理在時間軸的離散點上取值的信號，稱為離散時間信號。通常，離散時間信號用x（n）表示，其幅度可以在某一范圍內連續(xù)取值。由于信號處理所用的設備

2、主要是計算機或專用的信號處理芯片，均以有限的位數(shù)來表示信號的幅度，因此，信號的幅度也必須“量化”，即取離散值。我們把時間和幅度上均取離散值的信號稱為時域離散信號或數(shù)字信號。在MATLAB中，時域離散信號可以通過編寫程序直接生成，也可以通過對連續(xù)信號的等間隔抽樣獲得。 下面介紹常用的時域離散信號及其程序。1、單位抽樣序列MATLAB源程序為1) function x,n = impuls (n0,n1,n2)% Generates x(n) = delta(n-n0); n=n0 處建立一個單位抽樣序列% x,n = impuls (n0,n1,n2)if (n0 < n1) | (n0

3、> n2) | (n1 > n2)error('arguments must satisfy n1 <= n0 <= n2')endn = n1:n2;x = zeros(1,(n0-n1), 1, zeros(1,(n2-n0);將上述文件存為：impuls.m，在命令窗口輸入n0=0,n1=-10,n2=11;x,n=impuls (n0,n1,n2); stem(n,x,filled)2) n1=-5;n2=5;n0=0;n=n1:n2;x=n=n0; stem(n,x,'filled','k');axis(n1,n

4、2,1.1*min(x),1.1*max(x);title('單位脈沖序列');xlabel('時間(n)');ylabel('幅度x(n)');3) n1=-5;n2=5;k=0; n=n1:n2; nt=length(n); %求n點的個數(shù)nk=abs(k-n1)+1; %確定k在n序列中的位置x=zeros(1,nt); %對所有樣點置0x(nk)=1; %對抽樣點置1stem(n,x,'filled','k');axis(n1,n2,0,1.1*max(x);title('單位脈沖序列');

5、xlabel('時間(n)');Ylabel('幅度x(n)');2、單位階躍序列MATLAB源程序為：1) n1=-2;n2=8;n0=0; n=n1:n2; %生成離散信號的時間序列x=n>=n0; %生成離散信號x(n)stem(n, x,'filled','k'); %繪制脈沖桿圖，且圓點處用實芯圓表示axis(n1,n2,0,1.1*max(x);title('單位階躍序列');xlabel('時間(n)');Ylabel('幅度x(n)');2) n1=-2;n2=

6、8;k=0; n=n1:n2; nt=length(n); %求n點的個數(shù)nk=abs(k-n1)+1; %確定k在n序列中的位置x=zeros(1,nk-1),ones(1,nt-nk+1); %對所有樣點置0stem(n,x,'filled','k');axis(n1,n2,0,1.1*max(x);title('單位階躍序列');xlabel('時間(n)');ylabel('幅度x(n)');3、正弦序列x（n）=Um例、已知一時域周期性正弦信號的頻率為1HZ，振幅幅度為1V，在窗口中顯示兩個周期的信號波形

7、，并對該信號的一個周期進行32點采樣獲得離散信號。顯示連續(xù)信號和采樣獲得離散信號波形。MATLAB源程序為：f=1;Um=1;nt=2; %輸入信號頻率、振幅和顯示周期數(shù)N=32; T=1/f; %N為采樣點數(shù)，T為窗口顯示時間 dt=T/N; %采樣時間間隔n=0:nt*N-1; %離散信號的時間序列tn=n*dt; %時間序列樣點在時間軸上的位置x=Um*sin(2*f*pi*tn);subplot(2,1,1);plot(tn,x); %顯示原信號axis(0,nt*T,1.1*min(x) 1.1*max(x);ylabel('x(t)');subplot(2,1,2)

8、;stem(tn,x); %顯示經(jīng)采樣的信號axis(0,nt*T,1.1*min(x),1.1*max(x);ylabel('x(n)');4、矩形序列將square表示式中的t換成n，且n取整數(shù)，則可以獲得矩形序列。例、一個連續(xù)的周期性矩形信號頻率為5kHZ，信號幅度為0-2V之間，脈沖寬度與周期的比例為1:4,且要求在窗口上顯示其2個周期的信號波形，并對信號的一個周期進行16點采樣來獲得離散信號，顯示原連續(xù)信號與采樣獲得的離散信號。MATLAB源程序為：f=5000;nt=2; N=16;T=1/f; dt=T/N; n=0:nt*N-1;tn=n*dt;x=square

9、(2*f*pi*tn,25)+1; %產(chǎn)生時域信號,且幅度在02V之間 subplot(2,1,1);stairs(tn,x,'k'); axis(0 nt*T 1.1*min(x) 1.1*max(x);ylabel('x(t)');subplot(2,1,2);stem(tn,x,'filled','k'); axis(0 nt*T 1.1*min(x) 1.1*max(x);ylabel('x(n)');注意：直接用square子函數(shù)產(chǎn)生的信號波形，其幅度為-11之間。5、rand函數(shù)在實際系統(tǒng)的研究和處理中

10、，常常需要產(chǎn)生隨機信號，MATLAB提供的rand函數(shù)可以為我們生成隨機信號。例、生成一組41點構成的連續(xù)隨機信號和與之相應的隨機序列。MATLAB源程序為：tn=0:40;N=length(tn);x=rand(1,N);subplot(1,2,1),plot(tn,x,'k');subplot(1,2,2),stem(tn,x,'filled','k');四、實驗任務1、產(chǎn)生離散序列：（1）f（n）= （顯示-3<n<4）（2）f（n）=u（n-1） （顯示-5<n<5）2、一個連續(xù)的周期性正弦信號頻率為50HZ，信號

11、幅度在02V之間，在窗口上顯示2個周期的信號波形，對信號的一個周期進行16點采樣來獲得離散信號，顯示原連續(xù)信號和采樣獲得的離散信號波形。3、一個連續(xù)的周期性方波信號頻率為200HZ,信號幅度在02V之間，在窗口上顯示兩個周期的信號波形，用4Khz的頻率對連續(xù)信號進行采樣，顯示原連續(xù)信號和采樣獲得的離散信號波形。實驗2 離散序列的基本運算一、實驗目的l 學會運用MATLAB進行離散序列的運算，并掌握程序的編寫方法。二、實驗涉及的matlab子函數(shù)1、find功能：尋找非零元素的索引號。調用格式：find（n>=min(n1)）&(n<=max(n1)）；在符合關系運算條件的范

12、圍內尋找非零元素的索引號。2、fliplr功能：對矩陣行元素進行左右翻轉。調用格式：x1= fliplr（x）；將x的行元素左右翻轉，賦給變量x1。三、實驗原理離散序列的時域運算包括信號的相加、相乘，信號的時域變換包括信號的移位、反折、尺度變換等。在MATLAB中，離散序列的相加、相乘等運算是兩個向量之間的運算，因此參加運算的兩個序列向量必須具有相同的維數(shù)，否則應進行相應的處理。下面介紹各種離散序列的時域運算和時域變換的性質。1、 序列的移位x1（n）=x（n-m）例、x0=u（n），（顯示-10<n<10）x1=u（n+6）（顯示-10<n<10）,x2= u（n-4

13、）（顯示-10<n<10）MATLAB源程序為：n1=-10;n2=10;k0=0;k1=-6;k2=4; n=n1:n2; %生成離散信號的時間序列x0=n>=k0; %生成離散信號x0(n)x1=(n-k1)>=0; %生成離散信號x1(n)x2=(n-k2)>=0; %生成離散信號x2(n)subplot(3,1,1),stem(n, x0,'filled','k'); axis(n1,n2,1.1*min(x0),1.1*max(x0);ylabel('u(n)');subplot(3,1,2),stem(n

14、, x1,'filled','k'); axis(n1,n2,1.1*min(x1),1.1*max(x1);ylabel('u(n+6)');subplot(3,1,3),stem(n, x2,'filled','k'); axis(n1,n2,1.1*min(x2),1.1*max(x2);ylabel('u(n-4)');例、已知x（n）=，求x（n-2）和x（n+2）在-210區(qū)間的波形。MATLAB源程序為：n=-2:10;n0=2;n1=-2;x=2*sin(2*pi*n/10); %建

15、立原信號x(n)x1=2*sin(2*pi*(n-n0)/10); %建立x(n-2)信號x2=2*sin(2*pi*(n-n1)/10); %建立x(n+2)信號subplot(3,1,1),stem(n,x,'filled','k'); ylabel('x(n)');subplot(3,1,2),stem(n,x1,'filled','k'); ylabel('x(n-2)');subplot(3,1,3),stem(n,x2,'filled','k'); yla

16、bel('x(n+2)');2、 序列的相加情況1 兩序列具有相同的維數(shù)例、求x(n)= (0<n<10)MATLAB源程序為：n1=0;n2=10;n01=2;n02=4;n=n1:n2;x1=(n-n01)=0;x2=(n-n02)=0;x3=x1+x2;subplot(3,1,1);stem(n,x1,'filled','k');axis(n1,n2,1.1*min(x1),1.1*max(x1);ylabel('(n-2)');subplot(3,1,2);stem(n,x2,'filled',

17、'k');axis(n1,n2,1.1*min(x2),1.1*max(x2);ylabel('(n-4)');subplot(3,1,3);stem(n,x3,'filled','k');axis(n1,n2,1.1*min(x3),1.1*max(x3);ylabel('(n-2)+(n-4)');情況2 兩序列具有不同的維數(shù)例、x1（n）=u（n+2） （-4<n<6）X2(n)=u(n-4) (-5<n<8)求x（n）= x1（n）+ X2(n)=u(n-4)MATLAB源程序為：n

18、1=-4:6;n01=-2; x1=(n1-n01)>=0; %建立x1信號n2=-5:8;n02=4;x2=(n2-n02)>=0; %建立x2信號n=min(n1,n2):max(n1,n2); %為x信號建立時間序列nN=length(n); %求時間序列n的點數(shù)Ny1=zeros(1,N); %新建一維N列的y1全0數(shù)組y2=zeros(1,N); %新建一維N列的y2全0數(shù)組y1(find(n>=min(n1)&(n<=max(n1)=x1; %為y1賦值y2(find(n>=min(n2)&(n<=max(n2)=x2; %為y2

19、賦值x=y1+y2;subplot(3,1,1),stem(n1,x1,'filled','k');ylabel('x1(n)');axis(min(n),max(n),1.1*min(x),1.1*max(x);subplot(3,1,2),stem(n2,x2,'filled','k');ylabel('x2(n)');axis(min(n),max(n),1.1*min(x),1.1*max(x);subplot(3,1,3),stem(n,x,'filled','k&

20、#39;);ylabel('x(n)');axis(min(n),max(n),1.1*min(x),1.1*max(x);3、 序列相乘例、x1（n）= （-4<n<10）X2(n)=u(n+1) (-2<n<6) 求x（n）= x1（n）xX2(n)MATLAB源程序為：n1=-4:10; x1=3*exp(-0.25*n); %建立x1信號n2=-2:6;n02=-1;x2=(n2-n02)>=0; %建立x2信號n=min(n1,n2):max(n1,n2); %為x信號建立時間序列nN=length(n); %求時間序列n的點數(shù)Ny1=z

21、eros(1,N); %新建一維N列的y1全0數(shù)組y2=zeros(1,N); %新建一維N列的y2全0數(shù)組y1(find(n>=min(n1)&(n<=max(n1)=x1; %為y1賦值y2(find(n>=min(n2)&(n<=max(n2)=x2; %為y2賦值x=y1.*y2;subplot(3,1,1),stem(n1,x1,'filled','k');ylabel('x1(n)');axis(min(n),max(n),1.1*min(x1),1.1*max(x1);subplot(3,1,

22、2),stem(n2,x2,'filled','k');ylabel('x2(n)');axis(min(n),max(n),1.1*min(x2),1.1*max(x2);subplot(3,1,3),stem(n,x,'filled','k');ylabel('x(n)');axis(min(n),max(n),1.1*min(x),1.1*max(x);4、 序列反折例、x（n）= （-4<n<4）求x（-n）。MATLAB源程序為：n=-4:4;x=exp(-0.3*n);x1=

23、fliplr(x);n1=-fliplr(n);subplot(1,2,1),stem(n,x,'filled','k');title('x(n)');subplot(1,2,2),stem(n1,x1,'filled','k');title('x(-n)');5、 序列的尺度變換例、x（n）=sin（2*pi*n），求x（2n）和x（n/2）MATLAB源程序為：n=0:40;tn=n./20; %每周期取20個點x=sin(2*pi*tn); %建立原信號x(n)x1=sin(2*pi*tn*2

24、); %建立x(2n)信號x2=sin(2*pi*tn/2); %建立x(n/2)信號subplot(3,1,1),stem(tn,x,'filled','k'); ylabel('x(n)');axis(0,2,1.1*min(x),1.1*max(x);subplot(3,1,2),stem(tn,x1,'filled','k'); ylabel('x(2n)');axis(0,2,1.1*min(x),1.1*max(x);subplot(3,1,3),stem(tn,x2,'fill

25、ed','k'); ylabel('x(n/2)');axis(0,2,1.1*min(x),1.1*max(x); 四、實驗任務1、x（n）= （-5<n<5）2、x（n）=u（n-2）+u（n+2） （-5<n<5）3、已知x（n）=3cos（2*pi*n/10）,顯示x（n-3）和x（n+3）在020區(qū)間的波形。4、已知x1=exp（- n /16），x2=5sin（2*pi*n/10）,顯示x1乘以x2在024區(qū)間的波形。5、已知x（n）=nsin（n），顯示020區(qū)間的波形 y1(n)=x(n-3),y2（n）=x（-n

26、），y3=-x(n),y4=x(-n+3),y5（n）=x（n/2）.實驗3 模擬原型濾波器的設計一、實驗目的l 學會運用MATLAB設計模擬低通濾波器原型的設計方法。二、實驗涉及的matlab子函數(shù)1、buttord功能：確定巴特沃斯濾波器的階數(shù)和3dB截止頻率。調用格式：N,Wc= buttord(wp,ws,Rp,As,s)，計算巴特沃斯濾波器的階數(shù)和3dB截止頻率。Rp為通帶最大衰減，As為阻帶最小衰減。2、che1ord功能：確定切比雪夫1型濾波器的階數(shù)和通帶截止頻率。調用格式：N,Wp= buttord(wp,ws,Rp,As,s)3、che2ord功能：確定切比雪夫2型濾波器的階

27、數(shù)和阻帶截止頻率。調用格式： N,Wn= buttord(wp,ws,Rp,As,s) 4、ellipord功能：確定橢圓濾波器的階數(shù)和通帶截止頻率。調用格式： N,Wn= buttord(wp,ws,Rp,As,s)。5、buttap功能：巴特沃斯模擬低通濾波器原型（即歸一化的濾波器）。調用格式：z，p，k= buttap(n)，設計巴特沃斯模擬低通濾波器原型，其系統(tǒng)函數(shù)為6、cheb1ap功能：切比雪夫1型模擬低通濾波器原型。調用格式：z，p，k= cheb1ap (n，Rp)，設計切比雪夫1型模擬低通濾波器原型，其通帶最大衰減為Rp，系統(tǒng)函數(shù)為7、cheb2ap功能：切比雪夫2型模擬低通

28、濾波器原型。調用格式：z，p，k= cheb2ap (n，As)，設計切比雪夫2型模擬低通濾波器原型，其阻帶最小衰減為As，系統(tǒng)函數(shù)為8、ellipap功能：橢圓模擬低通濾波器原型。調用格式：z，p，k= ellipap (n，Rp，As)，設計橢圓模擬低通濾波器原型，其系統(tǒng)函數(shù)為9、poly功能：求某向量指定根所對應的特征多項式。調用格式：P=poly（），多項式P是一個特征多項式，的元素是多項式P的根。例如：>>a=1 2 ；P=poly(a)；則P = 1 -3 210、poly2strPa= poly2str(a,s)例如：>> P= 1,-3,2; PA=po

29、ly2str(P,'s')則PA = s2 - 3 s + 211、pzmap功能：顯示連續(xù)系統(tǒng)的零極點分布圖。調用格式：pzmap（b，a）；繪制由行向量b和a構成的系統(tǒng)函數(shù)確定的零極點分布圖pzmap(p，z)；繪制由列向量z確定的零點和p確定的極點構成的零極點分布圖。三、實驗原理由于IIR數(shù)字濾波器是在已知的歸一化的低通模擬濾波器的基礎上設計的，主要包括巴特沃斯、切比雪夫、橢圓低通濾波器，因此把這些低通濾波器稱為濾波器原型。下面介紹各種濾波器的設計。例1、 進行巴特沃斯濾波器原型的設計，獲得任意階數(shù)N的歸一化的系統(tǒng)函數(shù)公式和幅頻響應。MATLAB源程序為：%巴特沃斯模擬濾

30、波器原型n=input('N= '); %輸入濾波器階數(shù)N%計算n階模擬低通原型，得到左半平面零極點z0,p0,k0=buttap(n); b0=k0*poly(z0) %求濾波器系數(shù)b0a0=poly(p0) %求濾波器系數(shù)a0h,w=freqs(b0,a0); %顯示系統(tǒng)的頻率特性plot(w,abs(h),'r'), %畫幅頻特性圖axis(0,5,0,1.1),ylable('幅度');xlabel('f（HZ）'); Pb=poly2str(b0,'p') %給出b0決定的關于p多項式Pa=poly2st

31、r(a0,'p') %給出a0決定的關于p多項式輸入階數(shù)N=2則運行結果為b0 = 1a0 = 1.0000 1.4142 1.0000Pb = 1Pa = p2 + 1.4142 p + 1即歸一化的2階巴特沃斯濾波器原型的系統(tǒng)函數(shù)為例2、 通過模擬濾波器原型設計一個巴特沃斯模擬低通濾波器，要求通帶截止頻率2khz，通帶最大衰減1db，阻帶截止頻率5khz，阻帶最小衰減20db。MATLAB源程序為：%巴特沃斯模擬濾波器fp=2000;wp=2*pi*fp %輸入濾波器的通帶截止頻率fs=5000;ws=2*pi*fs %輸入濾波器的阻帶截止頻率Rp=1;As=20; %輸入

32、濾波器的通阻帶衰減指標%計算濾波器的階數(shù)和3dB截止頻率N,wc=buttord(wp,ws,Rp,As,'s') %計算n階模擬低通原型z0,p0,k0=buttap(n) %字母后加0表明這是原型濾波器的各指標，而不是所求的濾波器的b0=k0*poly(z0) %求歸一化濾波器分子系數(shù)b0a0=poly(p0) %求歸一化濾波器分母系數(shù)a0H,w=freqs(b0,a0); %求歸一化系統(tǒng)的頻率特性dbH=20*log10(abs(H)/max(abs(H); %將歸一化系統(tǒng)的幅頻特性化為分貝值，注意此時的分貝值為負的，為了使分貝圖和幅頻響應圖一致subplot(2,2,1

33、),plot(w*wc/(2*pi),abs(H),grid %畫所求濾波器的幅頻響應圖axis(0,6000,0,1.1);ylabel('幅度');xlabel('f(Hz)');subplot(2,2,2),plot(w*wc/(2*pi),angle(H),grid %畫所求濾波器的相頻響應圖axis(0,6000,-4,4);ylabel('相位');xlabel('f(Hz)');subplot(2,2,3),plot(w*wc/(2*pi),dbH),grid %畫所求濾波器的幅頻響應分貝圖axis(0,6000,-

34、30,2);ylabel('幅度(dB)');xlabel('f(Hz)');subplot(2,2,4),plot(p0*wc,'x'); %畫所求濾波器的極點圖，用x號標出axis square,axis equal,grid on %使x、y軸等比例顯示%以下程序可不寫wx0=wp,ws/wc; %設置頻率向量Hx=freqs(b0,a0,wx0); %計算通帶阻帶該兩點的頻率特性，即H(jwp)、H(jws)dbHx=-20*log10(abs(Hx)/max(abs(H) %將H(jwp)、H(jws)化為分貝值，即通帶最大衰減、阻帶最

35、小衰減運行結果為wp = 1.2566e+004ws = 3.1416e+004N = 4wc = 1.7689e+004z0 = p0 = -0.3827 + 0.9239i （此為歸一化濾波器的極點，第四幅圖所畫的為所求濾波器的極點p0*wc） -0.3827 - 0.9239i -0.9239 + 0.3827i -0.9239 - 0.3827ik0 = 1b0 = 1 a0 = 1.0000 2.6131 3.4142 2.6131 1.0000dbHx = 0.2730 20.0000 （實際求出來的濾波器的通帶最大衰減、阻帶最小衰減，可見阻帶指標剛好，通帶指標有富余）所求濾波器的

36、系統(tǒng)函數(shù)為例3、 進行切比雪夫1型濾波器原型的設計，獲得任意階數(shù)N的歸一化的系統(tǒng)函數(shù)公式和幅頻響應。MATLAB源程序為：n=input('N= '); %由使用者輸入濾波器階數(shù)NRp=input('Rp='); %輸入濾波器通帶衰減指標%計算n階模擬低通原型，得到左半平面零極點z0,p0,k0=cheb1ap(n,Rp); b0=k0*poly(z0) %求濾波器系數(shù)b0a0=poly(p0) %求濾波器系數(shù)a0h,w=freqs(b0,a0); %顯示系統(tǒng)的頻率特性plot(w,abs(h),'r'),axis(0,5,0,1.1),ylab

37、el('幅度');xlabel('f（HZ）'); title('切比雪夫1型原型模擬濾波器');Pb=poly2str(b0,'p') %給出b0決定的關于s多項式Pa=poly2str(a0,'p') %給出a0決定的關于s多項式輸入N=5,Rp=1運行結果為b0 = 0.1228a0 = 1.0000 0.9368 1.6888 0.9744 0.5805 0.1228Pb = 0.12283Pa = p5 + 0.93682 p4 + 1.6888 p3 + 0.9744 p2 + 0.58053 p +

38、0.12283即5階通帶最大衰減1db的切比雪夫1型原形濾波器系統(tǒng)函數(shù)為例4、設計一個切比雪夫1型模擬低通濾波器，要求通帶截止頻率2khz，通帶最大衰減1db，阻帶截止頻率5khz，阻帶最小衰減20db。MATLAB源程序為：%切比雪夫型模擬濾波器fp=2000;wp=2*pi*fp %輸入濾波器的通帶截止頻率fs=5000;ws=2*pi*fs %輸入濾波器的阻帶截止頻率Rp=1;As=20; %輸入濾波器的通阻帶衰減指標%計算濾波器的階數(shù)和通帶截止頻率n,wp=cheb1ord(wp,ws,Rp,As,'s')%計算n階模擬低通原型，得到左半平面零極點z0,p0,k0=ch

39、eb1ap(n,Rp)b0=k0*poly(z0) %求濾波器系數(shù)b0a0=poly(p0) %求濾波器系數(shù)a0H,w=freqs(b0,a0); %求系統(tǒng)的頻率特性dbH=20*log10(abs(H)/max(abs(H); %化為分貝值subplot(2,2,1),plot(w*wp/(2*pi),abs(H),'k'),grid axis(0,6000,0,1.1);ylabel('幅度');xlabel('f(Hz)');subplot(2,2,2),plot(w*wp/(2*pi),angle(H),'k'),grid

40、axis(0,6000,-4,4);ylabel('相位');xlabel('f(Hz)');subplot(2,2,3),plot(w*wp/(2*pi),dbH,'k'),grid axis(0,6000,-50,2);ylabel('幅度(dB)');xlabel('f(Hz)');subplot(2,2,4),plot(p0*wp,'xk') %顯示系統(tǒng)的極點分布axis square,axis equal,grid on %使X,Y軸等比例顯示wx0=wp,ws/wp; %設置頻率向量Hx

41、=freqs(b0,a0,wx0); %計算該兩點的頻率特性dbHx=-20*log10(abs(Hx)/max(abs(H) %化為分貝值運行結果為wp = 1.2566e+004ws = 3.1416e+004n = 3wp = 1.2566e+004z0 = p0 = -0.2471 + 0.9660i -0.4942 -0.2471 - 0.9660ik0 = 0.4913b0 = 0.4913a0 = 1.0000 0.9883 1.2384 0.4913dbHx = 0.9996 28.9442 （可見實現(xiàn)同樣的指標，階數(shù)比巴特沃斯低，但性能卻更好）例5、進行橢圓型濾波器原型的設計

42、，獲得任意階數(shù)N的歸一化的系統(tǒng)函數(shù)公式和幅頻響應。MATLAB源程序為：n=input('N= '); %由使用者輸入濾波器階數(shù)NRp=input('Rp='); %輸入濾波器通帶衰減指標As=input('As=');%計算n階模擬低通原型，得到左半平面零極點z0,p0,k0=ellipap(n,Rp,As); b0=k0*poly(z0) %求濾波器系數(shù)b0a0=poly(p0) %求濾波器系數(shù)a0h,w=freqs(b0,a0); %顯示系統(tǒng)的頻率特性plot(w,abs(h),'r'),axis(0,5,0,1.1),yl

43、abel('幅度');xlabel('f（HZ）'); title('橢圓型模擬濾波器原型');Pb=poly2str(b0,'p') %給出b0決定的關于s多項式Pa=poly2str(a0,'p') %給出a0決定的關于s多項式輸入N= 4 Rp=1 As=20運行結果為b0 = 0.1000 0 0.5422 0 0.5256a0 = 1.0000 0.9038 1.6765 0.8688 0.5897Pb = 0.1 p4 + 0.5422 p2 + 0.52559Pa = p4 + 0.90377 p3

44、+ 1.6765 p2 + 0.86881 p + 0.58972例6、設計一個切比雪夫1型模擬低通濾波器，要求通帶截止頻率2khz，通帶最大衰減2db，阻帶截止頻率5khz，阻帶最小衰減50db。MATLAB源程序為：%橢圓模擬濾波器fp=2000;wp=2*pi*fp %輸入濾波器的通帶截止頻率fs=5000;ws=2*pi*fs %輸入濾波器的阻帶截止頻率Rp=2;As=50; %輸入濾波器的通阻帶衰減指標%計算濾波器的階數(shù)和通帶截止頻率n,wn=ellipord(wp,ws,Rp,As,'s') %當不知道對哪個頻率歸一化時，用wn表示，最后在相應的指標處去除歸一化即可

45、%計算n階模擬低通原型，得到左半平面零極點z0,p0,k0=ellipap(n,Rp,As) b0=k0*poly(z0) %求濾波器系數(shù)b0Pb=poly2str(b0,'s') %給出b0決定的關于s多項式a0=(poly(p0) %求濾波器系數(shù)a0Pa=poly2str(a0,'s') %給出a0決定的關于s多項式H,w=freqs(b0,a0); %求系統(tǒng)的頻率特性dbH=20*log10(abs(H)/max(abs(H); %化為分貝值subplot(2,2,1),plot(w*wn/(2*pi),abs(H),grid %去除歸一化axis(0,1

46、0000,-0.1,1.1);ylabel('幅度');xlabel('f(Hz)');subplot(2,2,2),plot(w*wn/(2*pi),angle(H),grid %去除歸一化axis(0,10000,-4,4);ylabel('相位');xlabel('f(Hz)');subplot(2,2,3),plot(w*wn/(2*pi),dbH),grid axis(0,10000,-100,2);ylabel('幅度(dB)');xlabel('f(Hz)');subplot(2,2,

47、4),pzmap(b0,a0); %此為歸一化原型濾波器的零極點，不是所求濾波器的零極點wx0=wp,ws/wn; %設置頻率向量Hx=freqs(b0,a0,wx0); %計算該兩點的頻率特性dbHx=-20*log10(abs(Hx)/max(abs(H) %化為分貝值運行結果為wp = 1.2566e+004ws = 3.1416e+004n = 4wn = 1.2566e+004 （wn=wp可見橢圓濾波器是對通帶截止頻率進行歸一化的）z0 = 0 - 4.2226i 0 + 4.2226i 0 - 1.8719i 0 + 1.8719ip0 = -0.2673 - 0.4381i -

48、0.2673 + 0.4381i -0.0879 - 0.9678i -0.0879 + 0.9678ik0 = 0.0032b0 = 0.0032 0 0.0675 0 0.1976Pb = 0.0031623 s4 + 0.067465 s2 + 0.19757a0 = 1.0000 0.7104 1.3017 0.5512 0.2487Pa = s4 + 0.71038 s3 + 1.3017 s2 + 0.5512 s + 0.24872dbHx = 1.9993 50.2300四、實驗任務設計一個模擬低通濾波器，通帶截止頻率為6khz，通帶最大衰減為1db，阻帶截止頻率為15khz，

49、阻帶最小衰減為30db。要求：分別實現(xiàn)符合以上指標的巴特沃斯、切比雪夫1型、橢圓濾波器，繪制幅頻特性、幅頻特性分布圖、相頻特性、零極點分布圖，并列寫系統(tǒng)函數(shù)表示式，注意：是所求系統(tǒng)的系統(tǒng)函數(shù)，而不是歸一化的系統(tǒng)函數(shù)。實驗4 模擬域的頻率變換法一、實驗目的l 學會運用MATLAB設計模擬低通、帶通、高通、帶阻濾波器的設計方法。二、實驗涉及的matlab子函數(shù)1、lp2lp功能：低通到低通模擬濾波器變換。調用格式：bt,at= lp2lp(b,a,w0)，將系統(tǒng)函數(shù)表示的截止頻率為1rad/s的模擬低通濾波器原型變換為截止頻率為w0的低通濾波器。2、lp2hp功能：低通到高通模擬濾波器變換。調用格

50、式：bt,at= lp2hp(b,a,w0)，將系統(tǒng)函數(shù)表示的截止頻率為1rad/s的模擬低通濾波器原型變換為截止頻率為w0的高通濾波器。3、lp2bp功能：低通到帶通模擬濾波器變換。調用格式：bt,at= lp2bp(b,a,w0,bw)，將系統(tǒng)函數(shù)表示的截止頻率為1rad/s的模擬低通濾波器原型變換為中心頻率為w0、帶寬為bw的帶通濾波器。如果被設計的濾波器低端截止頻率為w1，高端截止頻率為w2，則W0=sqrt（w1*w2），bw=w2-w14、lp2bs功能：低通到帶阻模擬濾波器變換。調用格式：bt,at= lp2bs(b,a,w0,bw)，將系統(tǒng)函數(shù)表示的截止頻率為1rad/s的模擬

51、低通濾波器原型變換為中心頻率為w0、帶寬為bw的帶阻濾波器。如果被設計的濾波器低端截止頻率為w1，高端截止頻率為w2，則W0=sqrt（w1*w2），bw=w2-w1。三、實驗原理下面舉例說明用頻率變換法設計各種模擬濾波器的過程。例1、通過頻率變換法設計一個巴特沃斯模擬低通濾波器，要求通帶截止頻率2khz，通帶最大衰減1db，阻帶截止頻率5khz，阻帶最小衰減20db。MATLAB源程序為：%將巴特沃斯模擬濾波器原型變換成實際模擬低通濾波器fp=2000;wp=2*pi*fp; %輸入實際濾波器的通帶截止頻率fs=5000;ws=2*pi*fs; %輸入實際濾波器的阻帶截止頻率Rp=1;As=

52、20; %輸入濾波器的通阻帶衰減指標%計算濾波器的階數(shù)和3dB截止頻率n,wc=buttord(wp,ws,Rp,As,'s')%計算n階模擬低通原型，得到左半平面零極點z0,p0,k0=buttap(n); b0=k0*poly(z0) %求歸一化的濾波器系數(shù)b0a0=poly(p0) %求歸一化的濾波器系數(shù)a0 H,w0=freqs(b0,a0); %求歸一化的濾波器頻率特性dbH=20*log10(abs(H)/max(abs(H); %幅度化為分貝值%變換為實際模擬低通濾波器bt,at=lp2lp(b0,a0,wc); %從歸一化低通變換到實際低通Ht,wt=freqs

53、(bt,at); %求實際系統(tǒng)的頻率特性dbHt=20*log10(abs(Ht)/max(abs(Ht); %幅度化為分貝值%實際模擬低通頻率特性作圖subplot(2,2,1),plot(wt/2/pi,abs(Ht);axis(0,6000,0,1.1);ylabel('幅度');xlabel('f(Hz)');gridsubplot(2,2,2),plot(wt/2/pi,dbHt);axis(0,6000,-30,1);ylabel('幅度dB');xlabel('f(Hz)');gridsubplot(2,2,3),p

54、lot(wt/2/pi,angle(Ht);axis(0,6000,-4,4);ylabel('相位');xlabel('f(Hz)');grid運行結果為n = 4wc = 1.7689e+004b0 = 1a0 = 1.0000 2.6131 3.4142 2.6131 1.0000例2、通過頻率變換法設計一個巴特沃斯模擬高通濾波器，要求通帶截止頻率5khz，通帶最大衰減1db，阻帶截止頻率2khz，阻帶最小衰減20db。MATLAB源程序為：%將巴特沃斯模擬濾波器原型變換成實際模擬高通濾波器fp=5000;wp=2*pi*fp; %輸入實際濾波器的通帶截止頻率fs=2000;ws=2*pi*fs; %輸入實際濾波器的阻帶截止頻率Rp=1;As=20; %輸入濾波器的通阻帶衰減指標%計算濾波器的階數(shù)和3dB截止頻率n,wc=buttord(wp,ws,Rp,As,'s')%計算n階模擬低通原型，得到左半平面零極點z0,p0,k0=buttap(n); b0=k0*poly(z0) %求歸一化的濾波器系數(shù)b0a0=pol