實(shí)驗(yàn)1 離散系統(tǒng)的時域分析

1、實(shí)驗(yàn)1 離散系統(tǒng)的時域分析一、實(shí)驗(yàn)?zāi)康模杭由顚﹄x散系統(tǒng)的差分方程、單位抽樣響應(yīng)和卷積分析方法的理解。 二、實(shí)驗(yàn)原理：離散系統(tǒng)Discrete-time systmex(n)y(n)其輸入、輸出關(guān)系可用以下差分方程描述：輸入信號分解為沖激信號，系統(tǒng)單位抽樣序列h（n），則系統(tǒng)響應(yīng)為如下的卷積計算式：當(dāng)時，h(n)是有限長度的（n：0，M），稱系統(tǒng)為FIR系統(tǒng)；反之，稱系統(tǒng)為IIR系統(tǒng)。三 、實(shí)驗(yàn)內(nèi)容編制程序求解下列兩個系統(tǒng)的單位抽樣響應(yīng)，并繪出其圖形。（1）（2）（1）源程序：N=21;b=1 -1;a=1 0.75 0.125;x=1 zeros(1,N-1);n=0:1:N-1;y=filt

2、er(b,a,x);stem(n,y);xlabel('n');ylabel('幅度');title('單位抽樣響應(yīng)');圖形：（2）源程序：N=20;d=0 0.25 0.25 0.25 0.25;c=1;x=1 zeros(1,N-1);n=0:1:N-1;y=filter(d,c,x);stem(n,y);xlabel('n');ylabel('幅度');title('單位抽樣響應(yīng)');圖形如下： 實(shí)驗(yàn)2 離散系統(tǒng)的頻率響應(yīng)分析和零、極點(diǎn)分布一、實(shí)驗(yàn)?zāi)康模杭由顚﹄x散系統(tǒng)的頻率響應(yīng)分析和零、極點(diǎn)

3、分布的概念理解。二、實(shí)驗(yàn)原理：離散系統(tǒng)的時域方程為其變換域分析方法如下：系統(tǒng)函數(shù)為 分解因式 ，其中 和 稱為零、極點(diǎn)。 在MATLAB中，可以用函數(shù)z，p，K=tf2zp（num，den）求得有理分式形式的系統(tǒng)函數(shù)的零、極點(diǎn)，用函數(shù)zplane（z，p）繪出零、極點(diǎn)分布圖；也可以用函數(shù)zplane（num，den）直接繪出有理分式形式的系統(tǒng)函數(shù)的零、極點(diǎn)分布圖。使h=freqz(num,den,w)函數(shù)可求系統(tǒng)的頻率響應(yīng)，w是頻率的計算點(diǎn)，如w=0:pi/255:pi, h是復(fù)數(shù)，abs(h)為幅度響應(yīng)，angle(h)為相位響應(yīng)。另外，在MATLAB中，可以用函數(shù) r，p，k=residu

4、ez（num，den）完成部分分式展開計算；可以用函數(shù)sos=zp2sos（z，p，K）完成將高階系統(tǒng)分解為2階系統(tǒng)的串聯(lián)。三、實(shí)驗(yàn)內(nèi)容練習(xí)1 求下列直接型系統(tǒng)函數(shù)的零、極點(diǎn)，并將它轉(zhuǎn)換成二階節(jié)形式 解 用MATLAB計算程序如下： num=1 -0.1 -0.3 -0.3 -0.2;den=1 0.1 0.2 0.2 0.5;z,p,k=tf2zp(num,den);disp('零點(diǎn)');disp(z);disp('極點(diǎn)');disp(p);disp('增益系數(shù)');disp(k);sos=zp2sos(z,p,k);disp('二階節(jié)

5、');disp(real(sos);zplane(num,den)輸入到“num”和“den”的分別為分子和分母多項(xiàng)式的系數(shù)。計算求得零、極點(diǎn)增益系數(shù)和二階節(jié)的系數(shù)：零點(diǎn) 0.9615 -0.5730 -0.1443 + 0.5850i -0.1443 - 0.5850i  極點(diǎn) 0.5276 + 0.6997i 0.5276 - 0.6997i -0.5776 + 0.5635i -0.5776 - 0.5635i  增益系數(shù) 1  二階節(jié) 1.0000 -0.3885 -0.5509 1.0000 1.1552 0.6511 1.0000 0.2885

6、 0.3630 1.0000 -1.0552 0.7679  系統(tǒng)函數(shù)的二階節(jié)形式為： 極點(diǎn)圖如右圖。    2 差分方程 所對應(yīng)的系統(tǒng)的頻率響應(yīng)。 解：差分方程所對應(yīng)的系統(tǒng)函數(shù)為用MATLAB計算的程序如下：k=256;num=0.8 -0.44 0.36 0.02;den=1 0.7 -0.45 -0.6;w=0:pi/k:pi;h=freqz(num,den,w);subplot(2,2,1);plot(w/pi,real(h);gridtitle('實(shí)部')xlabel('omega/pi');ylabel('幅度&

7、#39;)subplot(2,2,2);plot(w/pi,imag(h);gridtitle('虛部')xlabel('omega/pi');ylabel('Amplitude')subplot(2,2,3);plot(w/pi,abs(h);gridtitle('幅度譜')xlabel('omega/pi');ylabel('幅值')subplot(2,2,4);plot(w/pi,angle(h);gridtitle('相位譜')xlabel('omega/pi'

8、;);ylabel('弧度')                   3求系統(tǒng) 的零、極點(diǎn)和幅度頻率響應(yīng)和相位響應(yīng)。要求：繪出零、極點(diǎn)分布圖，幅度頻率響應(yīng)和相位響應(yīng)曲線。 源程序：num=0.0528 0.0797 0.1295 0.1295 0.0797 0.0528;den=1 -1.8107 2.4947 -1.8801 0.9537 -0.2336;subplot(3,1,1); z,p,k=tf2zp(

9、num,den);disp('零點(diǎn)');disp(z);disp('極點(diǎn)');disp(p);disp('增益系數(shù)');disp(k);sos=zp2sos(z,p,k);disp('二階節(jié)');disp(real(sos);zplane(num,den)title('零極點(diǎn)分布')k=256;m=0.0528 0.0797 0.1295 0.1295 0.0797 0.0528;n=1 -1.8107 2.4947 -1.8801 0.9537 -0.2336;w=0:pi/k:pi;h=freqz(m,n,w)

10、;subplot(3,2,3);plot(w/pi,real(h);gridtitle('實(shí)部')xlabel('omega/pi');ylabel('幅度')subplot(3,2,4);plot(w/pi,imag(h);gridtitle('虛部')xlabel('omega/pi');ylabel('Amplitude')subplot(3,2,5);plot(w/pi,abs(h);gridtitle('幅度譜')xlabel('omega/pi');yla

11、bel('幅值')subplot(3,2,6);plot(w/pi,angle(h);gridtitle('相位譜')xlabel('omega/pi');ylabel('弧度') 實(shí)驗(yàn)結(jié)果：零點(diǎn) -1.0000 -0.3018 + 0.9534i -0.3018 - 0.9534i 0.0471 + 0.9989i 0.0471 - 0.9989i極點(diǎn) 0.2788 + 0.8973i 0.2788 - 0.8973i 0.3811 + 0.6274i 0.3811 - 0.6274i 0.4910 增益系數(shù) 0.0528二階節(jié)

12、 0.0528 0.0528 0 1.0000 -0.4910 0 1.0000 0.6036 1.0000 1.0000 -0.7622 0.53891.0000 -0.0941 1.0000 1.0000 -0.5575 0.8829圖形：實(shí)驗(yàn)3 FFT算法的應(yīng)用一 、實(shí)驗(yàn)?zāi)康模杭由顚﹄x散信號的DFT的理解及其FFT算法的運(yùn)用。二 、實(shí)驗(yàn)原理：N點(diǎn)序列的DFT和IDFT變換定義式如下：， 利用旋轉(zhuǎn)因子具有周期性，可以得到快速算法（FFT）。 在MATLAB中，可以用函數(shù)X=fft（x，N）和x=ifft（X，N）計算N點(diǎn)序列的DFT正、反變換。三、實(shí)驗(yàn)內(nèi)容1. 對連續(xù)的單一頻率周期信號 按

13、采樣頻率采樣，截取長度N分別選N =20和N =16，觀察其DFT結(jié)果的幅度譜。 解 此時離散序列： 即N=8。用MATLAB計算并作圖，函數(shù)fft用于計算離散傅里葉變換DFT，程序如下： q=8;k1=0:1:19;xa1=sin(2*pi*k1/q);subplot(2,2,1)plot(k1,xa1)xlabel('k');ylabel('x1(k)');xm1=fft(xa1);xm1=abs(xm1);subplot(2,2,2)stem(k1,xm1)xlabel('m');ylabel('X1(m)');k2=0:1

14、:15;xa2=sin(2*pi*k2/q);subplot(2,2,3)plot(k2,xa2)xlabel('k');ylabel('x2(k)');xm2=fft(xa2);xm2=abs(xm2);subplot(2,2,4)stem(k2,xm2)xlabel('m');ylabel('X2(m)');(d)(c)(b)(a)         圖2.1計算結(jié)果示于圖2.1，(c) 和(d) 分別是N=16時的截取信號和DFT結(jié)果，由于截取了兩

15、個整周期，得到單一譜線的頻譜；而(a) 和(b)分別是N=20時的截取信號和DFT結(jié)果，由于截取了兩個半周期，不再是單一譜線的頻譜，而是含有豐富的諧波分量，與N=16時的差別很大；上述頻譜的誤差主要是由于有限長序列的離散傅里葉變換是把有限長序列作為周期序列的一個周期來表示的，隱含有周期性意義。 2. 分別計算16點(diǎn)序列 的16點(diǎn)和32點(diǎn)DFT,繪出幅度譜圖形。源程序：q=16;n1=0:1:15;xn1=sin(5*pi*n1/q);subplot(2,2,1)plot(n1,xn1)xlabel('n');ylabel('x1(n)');xm1=fft(xn1

16、);xm1=abs(xm1);subplot(2,2,2)stem(n1,xm1)xlabel('m');ylabel('X1(m)');n2=0:1:31;xn2=sin(5*pi*n2/q);subplot(2,2,3)plot(n2,xn2)xlabel('n');ylabel('x2(n)');xm2=fft(xn2);xm2=abs(xm2);subplot(2,2,4)stem(n2,xm2)xlabel('m');ylabel('X2(m)');實(shí)驗(yàn)結(jié)果：實(shí)驗(yàn)4 無限沖激響應(yīng)數(shù)字濾波器

17、設(shè)計一、實(shí)驗(yàn)?zāi)康恼莆针p線性變換法及沖激響應(yīng)不變法設(shè)計IIR數(shù)字濾波器的具體設(shè)計方法及其原理；熟悉用雙線性變換法及沖激響應(yīng)不變法設(shè)計低通、高通和帶通IIR數(shù)字濾波器的計算機(jī)編程。二、實(shí)驗(yàn)編程函數(shù) 在MATLAB中，可以用下列函數(shù)輔助設(shè)計IIR數(shù)字濾波器：1) 利用buttord和cheblord可以確定低通原型巴特沃斯和切比雪夫?yàn)V波器的階數(shù)和截止頻率；2) num,den=butter（N,Wn）(巴特沃斯)；3) num,den=cheby1（N,Wn）,num,den=cheby2（N,Wn）(切比雪夫1型和2型)可以進(jìn)行濾波器的設(shè)計；4) lp2hp，lp2bp，lp2bs可以完成低通濾波

18、器到高通、帶通、帶阻濾波器的轉(zhuǎn)換；5) 使用bilinear可以對模擬濾波器進(jìn)行雙線性變換，求得數(shù)字濾波器的傳輸函數(shù)系數(shù)；6) 利用impinvar可以完成沖激響應(yīng)不變法的模擬濾波器到數(shù)字濾波器的轉(zhuǎn)換。注： 雙線性變換法通過將數(shù)字頻率的取值范圍從0到對應(yīng)到模擬頻率，也就對應(yīng)于模擬域中所有可能的頻率值。雙線性變換法不會出現(xiàn)頻率混疊，但非線性關(guān)系卻導(dǎo)致數(shù)字濾波器的頻率響應(yīng)不能逼真地模仿模擬濾波器的頻率響應(yīng)。沖激響應(yīng)不變法通過選擇滿足設(shè)計要求的模擬濾波器單位沖激響應(yīng)h（t）的采樣值的h（n），得到的被采樣的沖激響應(yīng)將給出與原模擬濾波器非常相近的濾波器形狀。由于該方法不可避免的要發(fā)生頻率混疊現(xiàn)象，所以

19、只適合設(shè)計低通和帶通濾波器。 三、實(shí)驗(yàn)內(nèi)容 1.設(shè)采樣周期T=250s，用沖激響應(yīng)不變法和雙線性變換法設(shè)計一個三階巴特沃斯濾波器，其3dB邊界頻率為fc =1kHz。 B,A=butter(3,2*pi*1000,'s'); num1,den1=impinvar(B,A,4000); h1,w=freqz(num1,den1); B,A=butter(3,2/0.00025,'s'); num2,den2=bilinear(B,A,4000); h2,w=freqz(num2,den2); f=w/pi*2000; plot(f,abs(h1),'-.&

20、#39;,f,abs(h2),'-'); grid; xlabel('頻率/Hz ') ylabel('幅值/dB') 程序中第一個butter的邊界頻率2×1000，為沖激響應(yīng)不變法原型低通濾波器的邊界頻率；第二個butter的邊界頻率2/T=2/0.00025，為雙線性變換法原型低通濾波器的邊界頻率.圖1給出了這兩種設(shè)計方法所得到的頻響，虛線為沖激響應(yīng)不變法的結(jié)果；實(shí)線為雙線性變換法的結(jié)果。沖激響應(yīng)不變法由于混疊效應(yīng)，使得過渡帶和阻帶的衰減特性變差，并且不存在傳輸零點(diǎn)。 圖12. 設(shè)計一巴特沃斯帶通濾波器，其通帶邊界頻率分別為f2=

21、110kHz和f1=90kHz處的最大衰減小于3dB，在阻帶邊界頻率分別為f3 =80kHz 和f4=120kHz處的最小衰減大于10dB，采樣頻率fs=400Hz。 w1=2*400*tan(2*pi*90/(2*400); w2=2*400*tan(2*pi*110/(2*400); wr1=2*400*tan(2*pi*80/(2*400); wr2=2*400*tan(2*pi*120/(2*400); N,wn=buttord(w1 w2, wr1 wr2,3,10,'s'); B,A=butter(N,wn,'s'); num,den=bilinea

22、r(B,A,400); h,w=freqz(num,den); f=w/pi*200; plot(f,20*log10(abs(h); axis(40,160,-30,10); grid; xlabel('頻率/kHz') ; ylabel('幅度/dB') 圖23. 利用MATLAB編程設(shè)計一個數(shù)字帶通濾波器，指標(biāo)要求如下：通帶邊緣頻率：，通帶衰減。阻帶邊緣頻率：，最小阻帶衰減 。給出IIR數(shù)字濾波器參數(shù)，繪出幅度和相位頻響曲線，討論實(shí)現(xiàn)形式和特點(diǎn)。源程序：wp=0.45 0.65; ws=0.3 0.75; Ap=1;As=40;n1,wn1=buttord

23、(wp,ws,Ap,As);num,den=butter(n1,wn1); w=0:pi/255:pi; h=freqz(num,den,w); g=20*log10(abs(h); pha=angle(h); subplot(1,2,1)plot(w/pi,g); grid xlabel('頻率/kHz');ylabel('幅度');title('數(shù)字帶通濾波器幅頻曲線')subplot(1,2,2);plot(w/pi,pha);gridxlabel('頻率/kHz');ylabel('相位');title(&

24、#39;數(shù)字帶通濾波器相位曲線'); 實(shí)驗(yàn)結(jié)果：實(shí)現(xiàn)形式和特點(diǎn)分析：實(shí)驗(yàn)5 有限沖激響應(yīng)數(shù)字濾波器設(shè)計一 、實(shí)驗(yàn)?zāi)康模海?） 掌握用窗函數(shù)法設(shè)計FIR數(shù)字濾波器的原理和方法。（2） 了解各種窗函數(shù)對濾波特性的影響。二、實(shí)驗(yàn)原理：如果所希望的濾波器的理想頻率響應(yīng)函數(shù)為 ,則其對應(yīng)的單位抽樣響應(yīng)為窗函數(shù)設(shè)計法的基本原理是用有限長單位抽樣響應(yīng)序列h(n)逼近hd(n)，由于hd(n)往往是無限長序列，且是非因果的，所以用窗函數(shù)w(n)將hd(n)截斷，并進(jìn)行加權(quán)處理，得到:h(n)就作為實(shí)際設(shè)計的FIR數(shù)字濾波器的單位抽樣響應(yīng)序列，其頻率響應(yīng)函數(shù)為：式中，N為所選窗函數(shù)的長度 在MATLAB

25、中，可以用b=fir1(M,Wc,ftype,taper) 等函數(shù)輔助設(shè)計FIR數(shù)字濾波器。M代表濾波器階數(shù)(單位抽樣響應(yīng)h（n）的長度N=M+1)；Wc代表濾波器的截止頻率(對歸一化頻率)，當(dāng)設(shè)計帶通和帶阻濾波器時，Wc為雙元素相量；ftype代表濾波器類型，如high高通，stop帶阻等；taper為窗函數(shù)類型，默認(rèn)為海明窗；窗函數(shù)用blackman, hamming, hanning ， kaiser產(chǎn)生。三、實(shí)驗(yàn)內(nèi)容1.設(shè)計一FIR低通濾波器，通帶截止頻率 ，通帶衰減不大于0.25dB，阻帶截止頻率 ，阻帶衰減不小于50dB。選擇一個合適的窗函數(shù)，繪出所設(shè)計的濾波器的頻率響應(yīng)圖。wp=

26、0.2*pi;ws=0.3*pi;tr_width=ws-wp;M=ceil(6.6*pi/tr_width);n=0:1:M-1;wc=(ws+wp)/2;hd=ideal_lp(wc,M);w_ham=(hamming(M)'h=hd.*w_ham;H,w=freqz(h,1,1000,'whole');H=(H(1:1:501)'w=(w(1:1:501)'mag=abs(H);db=20*log10(mag+eps)/max(mag);delta_w=2*pi/1000;Rp=-(min(db(1:1:wp/delta_w+1);disp(

27、9;實(shí)際通帶波動為',num2str(Rp);As=-round(max(db(ws/delta_w+1:1:501);disp('最小阻帶衰減為',num2str(As); subplot(2,2,1);stem(n,hd);title('理想抽樣響應(yīng)');axis(0 M-1 -0.1 0.3);xlabel('n');ylabel('hd(n)');subplot(2,2,2);stem(n,w_ham);title('海明窗');axis(0 M-1 0 1.1);xlabel('n');ylabel('w(n)');subplot(2,2,3);stem(n,h);title('實(shí)際脈沖響應(yīng)');axis(0 M-1 -0.1 0.3);xlabel('n');ylabel('h(n)');subplot(2,2,4);pl