南郵2013DSP實驗報告(英文mitra版)

1、南京郵電大學實 驗 報 告實驗名稱：離散時間信號與系統(tǒng)的時、頻域表示離散傅立葉變換和z變換 數(shù)字濾波器的頻域分析和實現(xiàn)數(shù)字濾波器的設計課程名稱 數(shù)字信號處理A(雙語) 班級學號_ 姓 名_開課時間 2013 /2014 學年， 第一學期實驗名稱：離散時間信號與系統(tǒng)的時、頻域表示實驗目的和任務：熟悉Matlab基本命令，理解和掌握離散時間信號與系統(tǒng)的時、頻域表示及簡單應用。在Matlab環(huán)境中，按照要求產(chǎn)生序列，對序列進行基本運算；對簡單離散時間系統(tǒng)進行仿真，計算線性時不變（LTI）系統(tǒng)的沖激響應和卷積輸出；計算和觀察序列的離散時間傅立葉變換（DTFT）幅度譜和相位譜。實驗內(nèi)容：基本序列產(chǎn)生和運

2、算： Q1.11.3，Q1.23，Q1.301.33離散時間系統(tǒng)仿真： Q2.12.3LTI系統(tǒng)：Q2.19，Q2.21，Q2.28DTFT：Q3.1，Q3.2，Q3.4實驗過程與結果分析：Q1.1 運行P1_1產(chǎn)生單位樣本序列 un 的程序與顯示的波形如下:clf;n=-10:20;u=zeros(1,10) 1 zeros(1,20);stem(n,u);xlabel('時間序號n');ylabel('振幅');title('單位樣本序列');axis(-10 20 0 1.2);Q1.2 clf 命令的作用是- 清除圖形窗口上的圖形axis

3、命令的作用是- 設置坐標軸的范圍和顯示方式title命令的作用是- 給圖形加名字xlabel命令的作用是- 添加x坐標標注ylabel命令的作用是- 添加y坐標標注Q1.3產(chǎn)生有延時11個樣本udn的程序及其運行結果如下： clf;n=-10:20;u=zeros(1,21) 1 zeros(1,9);stem(n,u);xlabel('時間序號n');ylabel('振幅');title('單位樣本序列');axis(-10 20 0 1.2);Q1.23 修改上述程序，以長生長度為50、頻率為0.08、振幅為2.5、相移為90度的一個正弦序列

4、并顯示它。該序列的周期是多少？n=0:50;f=0.1;f=0.08;phase=pi;A=2.5;arg=2*pi*f*n-phase;x=A*cos(arg);clf;stem(n,x);axis(0 50 -3 3);grid;周期為：2/ = 1/f = 1/0.08 = 1/(8/100) = 100/8 = 25/2.Q1.30 未污染的信號sn 是什么樣的形式？加性噪聲dn 是什么樣的形式？加性噪聲dn 是均勻分布在 -0.4和+0.4之間的隨機序列Q1.31 使用語句s=s+d能產(chǎn)生被噪聲污染的信號嗎？若不能，為什么？不能。因為d是列向量，s是行向量Q1.32 信號x1、x2、

5、x3與x之間的關系是什么？x1是x的延時，x2和x相等，x3超前于xQ1.33 legend的作用是什么legend用于產(chǎn)生圖例說明Q2.1 對于M = 2 和輸入 xn = s1n+s2n，程序P2.1的輸出為： 輸入 xn 被該離散時間系統(tǒng)抑制的分量為- Signal #2的高頻分量 Q2.2 程序P2.1 中 LTI system 被修改為 yn = 0.5(xnxn1)后， 輸入 xn = s1n+s2n 導致的輸出為:對于輸入的影響是- 該系統(tǒng)現(xiàn)在是一個高通濾波器。它通過高頻率的輸入分量S2，而不是低頻分量輸入S1。Q2.3程序 P2_1 對于不同M（M=4,6）取值和不同正弦分量(

6、任取2個)取值的運行結果如下： M=4 f1=0.05 f2=0.10M=6 f1=0.30 f2=0.30Q2.19 運行 P2_5 生成的結果如下：:Q2.21 生成的MATLAB代碼如下:clf;N = 40;num = 0.9 -0.45 0.35 0.002;den = 1.0 0.71 -0.46 -0.62;x = 1 zeros(1,N-1);y = filter(num,den,x);stem(y);xlabel('時間序號n'); ylabel('振幅');title('沖擊響應'); grid;程序產(chǎn)生的40個樣本如下所示:

7、Q2.28 程序P2_7產(chǎn)生的序列 yn and y1n 如下所示： yn 和 y1n 的差別為- 它們無差別。將xn補零后得到 x1n作為輸入,產(chǎn)生y1n的原因是-對于長度N1和N2的兩個序列，轉(zhuǎn)化率返回得到的序列長度N1 + N2-1。與此相反，過濾器接受一個輸入信號和一個系統(tǒng)規(guī)范。返回的結果是相同的長度作為輸入信號。因此，為了從轉(zhuǎn)化率和濾波器得到直接比較的結果，有必要供應濾波器的輸入已經(jīng)零填充為長度L(x)+L(h)-1。Q3.1 程序P3_1計算離散時間傅里葉變換的原始序列為-H(e)=pause 命令的作用為- 不加參數(shù),直接用pause的話,就是程序暫停,直至用戶按任意一個按鍵。如

8、果加參數(shù)，比如pause(1.5)就是程序暫停1.5秒。Q3.2 程序 P3_1 運行結果為:DTFT 是關于 w的周期函數(shù)么？答：DTFT是關于的周期函數(shù) 周期是 2四個圖形的對稱性為：實部是2周期和偶對稱； 虛部是2周期和奇對稱；幅度 是2周期和偶對稱；相位是2周期和奇對稱性。Q3.4 修改程序 P3_1 重做Q3.2的程序如下: clf;w = -4*pi:8*pi/511:4*pi;num = 1 3 5 7 9 11 13 15 17;den = 1;h = freqz(num, den, w);subplot(2,1,1)plot(w/pi,real(h);gridtitle(&#

9、39;H(ejomega)的實部')xlabel('omega /pi');ylabel('振幅');subplot(2,1,2)plot(w/pi,imag(h);gridtitle('H(ejomega)的虛部')xlabel('omega /pi');ylabel('振幅');pausesubplot(2,1,1)plot(w/pi,abs(h);gridtitle('|H(ejomega)|幅度譜')xlabel('omega /pi');ylabel('振幅

10、');subplot(2,1,2)plot(w/pi,angle(h);gridtitle('相位譜argH(ejomega)')xlabel('omega /pi');ylabel('以弧度為單位的相位');修改程序后的運行結果為: DTFT 是關于 w的周期函數(shù)么？ 答：DTFT 是關于 w的周期函數(shù)。周期是 - 2 相位譜中跳變的原因是 - 角度返回到arctan的本值實驗名稱：離散傅立葉變換和z變換實驗目的和任務：掌握離散傅立葉變換（DFT）及逆變換（IDFT）、z變換及逆變換的計算和分析。利用Matlab語言，完成DFT和IDF

11、T的計算及常用性質(zhì)的驗證，用DFT實現(xiàn)線性卷積，實現(xiàn)z變換的零極點分析，求有理逆z變換。實驗內(nèi)容：DFT和IDFT計算： Q3.233.24DFT的性質(zhì)： Q3.263.29，Q3.36，Q3.38，Q3.40z變換分析：Q3.463.48逆z變換：Q3.50實驗過程與結果分析：Q3.23 編寫一個MATLAB程序，計算并畫出長度為N的L點離散傅里葉變換Xk的值，其中LN，然后計算并畫出L點離散傅里葉變換Xk。對不同長度N和不同的離散傅里葉變換長度L，運行程序。討論你的結果。編寫的MATLAB程序：clf;N=200; % length of signalL=256; % length of

12、DFTnn = 0:N-1;kk = 0:L-1;xR = 0.1*(1:100) zeros(1,N-100); xI = zeros(1,N); x = xR + i*xI;XF = fft(x,L);subplot(3,2,1);grid;plot(nn,xR);grid;title('Rexn');xlabel('Time index n');ylabel('Amplitude');subplot(3,2,2);plot(nn,xI);grid;title('Imxn');xlabel('Time index n&

13、#39;);ylabel('Amplitude');subplot(3,2,3);plot(kk,real(XF);grid;title('ReXk');xlabel('Frequency index k');ylabel('Amplitude');subplot(3,2,4);plot(kk,imag(XF);grid;title('ImXk');xlabel('Frequency index k');ylabel('Amplitude');% IDFTxx = ifft(XF,L

14、);subplot(3,2,5);plot(kk,real(xx);grid;title('Real part of IDFTXk');xlabel('Time index n');ylabel('Amplitude');subplot(3,2,6);plot(kk,imag(xx);grid;title('Imag part of IDFTXk');xlabel('Time index n');ylabel('Amplitude');Q3.24 寫一個MATLAB程序，用一個N點復數(shù)離散傅里葉計算

15、兩個長度為N的實數(shù)序列的N點離散傅里葉變換，并將結果同直接使用兩個N點離散傅里葉變換得到的結果進行比較。 編寫的MATLAB程序：clf;N=256; % length of signalnn = 0:N-1;ntime = -N/2:N/2-1;g = (0.75).abs(ntime); % signal gh = (-0.9).ntime; % signal hGF = fft(g);HF = fft(h);x = g + i*h;XF = fft(x);XFstar = conj(XF);XFstarmod = XFstar(1) fliplr(XFstar(2:N);GF2 = 0.

16、5*(XF + XFstarmod);HF2 = -i*0.5*(XF - XFstarmod);abs(max(GF-GF2)abs(max(HF-HF2)figure(1);clf;subplot(2,2,1);grid;plot(nn,real(GF);grid;title('Two N-point DFT''s');xlabel('Frequency index k');ylabel('ReGk');subplot(2,2,2);plot(nn,imag(GF);grid;title('Two N-point DF

17、T''s');xlabel('Frequency index k');ylabel('ImGk');subplot(2,2,3);grid;plot(nn,real(GF2);grid;title('Single N-point DFT');xlabel('Frequency index k');ylabel('ReGk');subplot(2,2,4);plot(nn,imag(GF2);grid;title('Single N-point DFT');xlabel(

18、9;Frequency index k');ylabel('ImGk');figure(2);clf;subplot(2,2,1);grid;plot(nn,real(HF);grid;title('Two N-point DFT''s');xlabel('Freq index k');ylabel('ReHk');subplot(2,2,2);plot(nn,imag(HF);grid;title('Two N-point DFT''s');xlabel('Freq

19、 index k');ylabel('ImHk');subplot(2,2,3);grid;plot(nn,real(HF2);grid;title('Single N-point DFT');xlabel('Freq index k');ylabel('ReHk');subplot(2,2,4);plot(nn,imag(HF2);grid;title('Single N-point DFT');xlabel('Freq index k');ylabel('ImHk');

20、Q3.26 在函數(shù)circshift中，命令rem的作用是什么？答：rem(x,y)是用y對x求余數(shù)函數(shù)。Q3.27 解釋函數(shù)circshift怎樣實現(xiàn)圓周移位運算。答：在輸入序列x由M的位置開始被循環(huán)移位。如果M> 0，則circshift刪除從矢量x最左邊開始的M個元素和它們附加在右側(cè)的剩余元素，以獲得循環(huán)移位序列。如果如果M<0，則circshift首先通過x的長度來彌補M，即序列x最右邊的長度的M樣品從x中刪除和所附在其余的M個樣本的右側(cè)，以獲得循環(huán)移位序列。Q3.28 在函數(shù)circshift中，運算符=的作用是什么？ 答：=是不等于的意思。Q3.29 解釋函數(shù)circo

21、nv怎樣實現(xiàn)圓周卷積運算。 答：輸入是兩個長度都為L的向量x1和x2，它是非常有用的定期延長X2的函數(shù)。讓x2p成為x2延長無限長的周期的序列。從概念上講，在定點時間上通過時序交換后的x2p的長度L交換x2p序列和x2tr等于1的元素。然后元素1至L的輸出向量y是通過取x1和獲得的長度為L的sh矢量之間的內(nèi)積得到通過循環(huán)右移的時間反轉(zhuǎn)向量x2tr。對于輸出樣本Yn的1NL時，右循環(huán)移位的量為n-1個位置上。Q3.36 運行程序P3.9并驗證離散傅里葉變換的圓周卷積性質(zhì)。g1 = 1 2 3 4 5 6; g2 = 1 -2 3 3 -2 1;ycir = circonv(g1,g2);disp

22、('Result of circular convolution = ');disp(ycir)G1 = fft(g1); G2 = fft(g2);yc = real(ifft(G1.*G2);disp('Result of IDFT of the DFT products = ');disp(yc)Result of circular convolution =12 28 14 0 16 14Result of IDFT of the DFT products =12 28 14 0 16 14一個圓周卷積的DTF是DTF的逐點產(chǎn)物。Q3.38 運行程序P3

23、.10并驗證線性卷積可通過圓周卷積得到。g1 = 1 2 3 4 5;g2 = 2 2 0 1 1;g1e = g1 zeros(1,length(g2)-1);g2e = g2 zeros(1,length(g1)-1);ylin = circonv (g1e,g2e);disp('Linear convolution via circular convolution = ');disp(ylin);y = conv(g1, g2);disp('Direct linear convolution = ');disp(y)Linear convolution v

24、ia circular convolution =2 6 10 15 21 15 7 9 5Direct linear convolution =2 6 10 15 21 15 7 9 5使用圓周卷積確實有可能得到線性卷積Q3.40 編寫一個MATLAB程序，對兩個序列做離散傅里葉變換，已生成他們的線性卷積。用此程序驗證Q3.38和Q3.39的結果編寫的MATLAB程序：% Program Q3.40g1 = 1 2 3 4 5;g2 = 2 2 0 1 1;g1e = g1 zeros(1,length(g2)-1);g2e = g2 zeros(1,length(g1)-1);G1EF =

25、 fft(g1e);G2EF = fft(g2e);ylin = real(ifft(G1EF.*G2EF);disp('Linear convolution via DFT = ');disp(ylin);Linear convolution via DFT =2.0000 6.0000 10.0000 15.0000 21.0000 15.0000 7.0000 9.0000 5.0000Q3.46 使用程序P3.1在單位圓上求下面的z變換： G(z)=Q3.47 編寫一個MATLAB程序，計算并顯示零點和極點，計算并顯示其因式形式，并產(chǎn)生z的兩個多項式之比的形式表示的z變

26、換的極零點圖。使用該程序，分析式(3.32)的z變換G(z)。 編寫的MATLAB程序：% Program Q3_47clf;num = 2 5 9 5 3;den = 5 45 2 1 1;z p k = tf2zpk(num,den);disp('Zeros:');disp(z);disp('Poles:');disp(p);input('Hit <return> to continue.');sos k = zp2sos(z,p,k)input('Hit <return> to continue.')

27、;zplane(z,p);運行結果：Zeros: -1.0000 + 1.4142i -1.0000 - 1.4142i -0.2500 + 0.6614i -0.2500 - 0.6614iPoles: -8.9576 -0.2718 0.1147 + 0.2627i 0.1147 - 0.2627isos = 1.0000 2.0000 3.0000 1.0000 9.2293 2.4344 1.0000 0.5000 0.5000 1.0000 -0.2293 0.0822k =0.4000Q3.48 通過習題Q3.47產(chǎn)生的極零點圖，求出G(z)的收斂域的數(shù)目。清楚地顯示所有的收斂域。

28、由極零點圖說明離散時間傅里葉變換是否存在。R1 : | z | < 0.2718 (left-sided, not stable)R2 : 0.2718 < | z | < 0.2866 (two-sided, not stable)R 3: 0.2866 < | z | < 8.9576 (two-sided, stable)R4 : | z | > 8.9576 (right-sided, not stable) 不能從極零點圖肯定的說DTFT是否存在，因為其收斂域一定要指定。當收斂域在上述R 3內(nèi)所獲得的序列卻是證明了DTFT的存在，它是一個具有雙面沖

29、激響應的穩(wěn)定系統(tǒng)。Q3.50 編寫一個MATLAB程序，計算一個有理逆z變換的前L個樣本，其中L的值由用戶通過命令input提供。用該程序計算并畫出式(3.32)中G(z)的逆變換的前50個樣本。使用命令stem畫出由逆變換產(chǎn)生的序列。編寫的MATLAB程序：% Program Q3.50clf;% initializenum = 2 5 9 5 3;den = 5 45 2 1 1;% Query user for parameter LL = input('Enter the number of samples L: ');% find impulse responseg

30、t = impz(num,den,L);%plot the impulse responsestem(t,g);title('First ',num2str(L),' samples of impulse response');xlabel('Time Index n');ylabel('hn');Enter the number of samples L: 50實驗名稱：數(shù)字濾波器的頻域分析和實現(xiàn)實驗目的和任務：（1）求濾波器的幅度響應和相位響應，觀察對稱性，判斷濾波器類型。（2）用Matlab函數(shù)實現(xiàn)系統(tǒng)的級聯(lián)型和并聯(lián)型結構，

31、并對濾波器進行結構仿真。實驗內(nèi)容：系統(tǒng)傳遞函數(shù)的級聯(lián)和并聯(lián)實現(xiàn)：Q6.1, Q6.3, Q6.5, Q8.3,Q8.5,濾波器的幅頻特性分析：Q8.1, Q8.9,Q8.10,Q8.14實驗過程與結果分析： Q6.1 使用程序P6.1，生成如下有限沖激響應傳輸函數(shù)的一個級聯(lián)實現(xiàn)：H(z)=2+10z+23z+34z+31z+16z編寫的MATLAB程序：% Program P6_1 num=input('分子系數(shù)向量=');den=input('分母系數(shù)向量=');z,p,k=tf2zp(num,den);sos=zp2sos(z,p,k);分子系數(shù)向量=2 1

32、0 23 34 31 4分母系數(shù)向量=1 1 1 1 1 1sos = 2.0000 6.0000 4.0000 1.0000 0 0 1.0000 1.0000 2.0000 1.0000 0 0 1.0000 1.0000 0.5000 1.0000 0 0畫出級聯(lián)實現(xiàn)的框圖：H(z)不是一個線性相位傳輸函數(shù)。Q6.3 使用程序P6.1，H(z)=分子系數(shù)向量 = 3 8 12 7 2 2分母系數(shù)向量= 16 24 24 14 5 1 sos = 0.1875 -0.0625 0 1.0000 0.5000 0 1.0000 2.0000 2.0000 1.0000 0.5000 0.25

33、00 1.0000 1.0000 1.0000 1.0000 0.5000 0.5000畫出級聯(lián)實現(xiàn)的框圖：Q6.5 使用程序P6.2，式H(z)=所示因果無限沖激響應傳輸函數(shù)的兩種不同并聯(lián)形式實現(xiàn)。畫出兩種實現(xiàn)的框圖：并聯(lián) I 形 留數(shù)是 -0.4219 + 0.6201i -0.4219 - 0.6201i 2.3437 0.3438 - 2.5079i 0.3438 + 2.5079i極點在 -0.2500 + 0.6614i -0.2500 - 0.6614i -0.5000 -0.2500 + 0.4330i -0.2500 - 0.4330i常數(shù) -2 并聯(lián)II型留數(shù)是 -0.30

34、47 - 0.4341i -0.3047 + 0.4341i -1.1719 1.0000 + 0.7758i 1.0000 - 0.7758i極點在 -0.2500 + 0.6614i -0.2500 - 0.6614i -0.5000 -0.2500 + 0.4330i -0.2500 - 0.4330i常數(shù) 0.1875Q8.1 程序P8.1設計的濾波器類型為 IIR帶阻濾波器 ，其指標為 =0.4 0.5 =0.1 0.8 , R=1dB，R=30dB，濾波器的階數(shù)為 4 。為了驗證仿真需計算的沖激響應樣本為 5個。仿真是正確的嗎？Q8.3 生成Q8.1中產(chǎn)生的傳輸函數(shù)的一個級聯(lián)實現(xiàn)，

35、并編寫一個程序來仿真它。其中每個單獨的部分用一個直接型實現(xiàn)。驗證仿真。format shortnum = input('Numerator coefficients = ');den = input('Denominator coefficients = ');Numfactors = factorize(num);Denfactors = factorize(den);disp('Numerator Factors'),disp(Numfactors)disp('Denominator Factors'),disp(Denfac

36、tors)Numerator coefficients = 0.0571 0 -0.1143 0 0.0571Denominator coefficients = 1.0000 -0.5099 1.2862 -0.3350 0.4479Numerator Factors Columns 1 through 3 1.000000000000000 1.021143238995107 0 1.000000000000000 0.979294541463240 0 1.000000000000000 -1.021143238995070 0 1.000000000000000 -0.97929454

37、1463277 0Denominator Factors Columns 1 through 3 1.000000000000000 -0.597633166080704 0.678485370310279 1.000000000000000 0.087733166080704 0.660146879504817Q8.5 生成Q8.1中傳輸函數(shù)的一個并聯(lián)型實現(xiàn)，并編寫一個程序來仿真它。每一個單獨的部分用直接型實現(xiàn)。驗證仿真。num = input('Numerator coefficient vector = ');den = input('Denominator co

38、efficient vector = ');r1,p1,k1 = residuez(num,den);r2,p2,k2 = residue(num,den);disp('Parallel Form I')disp('Residues are');disp(r1);disp('Poles are at');disp(p1);disp('Constant value');disp(k1);disp('Parallel Form II')disp('Residues are');disp(r2);

39、disp('Poles are at');disp(p2);disp('Constant value');disp(k2); b1,a1=residuez(R1,P1,0); R1=r1(1) r1(2) r1(3) r1(4); P1=p1(1) p1(2) p1(3) p1(4);disp('b1=');disp(b1); disp('a1=');disp(a1); Numerator coefficient vector = 0.0571 0 -0.1143 0 0.0571Denominator coefficient v

40、ector = 1.0000 -0.5099 1.2862 -0.3350 0.4479Parallel Form IResidues are -0.023475623417519 + 0.197736845648730i -0.023475623417519 - 0.197736845648730i -0.011716283258078 - 0.213037612139943i -0.011716283258078 + 0.213037612139943iPoles are at 0.298816583040352 + 0.767589747202480i 0.298816583040352

41、 - 0.767589747202480i -0.043866583040352 + 0.811309190381311i -0.043866583040352 - 0.811309190381311iConstant value 0.127483813351194Parallel Form IIResidues are -0.158795680938490 + 0.041067400713457i -0.158795680938490 - 0.041067400713457i 0.173353325938490 - 0.000160296180734i 0.173353325938490 +

42、 0.000160296180734iPoles are at 0.298816583040352 + 0.767589747202480i 0.298816583040352 - 0.767589747202480i -0.043866583040352 + 0.811309190381311i -0.043866583040352 - 0.811309190381311iConstant value 0.057100000000000b1= -0.046951246835039 -0.289531739579519 0a1= 1.000000000000000 -0.59763316608

43、0704 0.678485370310279Q8.9 程序P8.3設計的濾波器類型為 IIR低通橢圓濾波器 ，指標為=0.25 =0.55，R=0.5dB，R=50dB，階數(shù)為 8，輸入的正弦序列的頻率為f=0.7kHz。Q8.10 運行P8.3并產(chǎn)生兩個圖形。哪種輸入成分會在濾波器輸出出現(xiàn)？為什么輸出序列的開始部分不是一種理想的正弦曲線？修改P8.3，以便紙過濾序列X2n。產(chǎn)生的輸出序列和預料的是一樣的嗎？答：產(chǎn)生的輸出序列和預料的是不一樣的Q8.14 程序P8.4設計的濾波器類型為 FIR低通濾波器，指標為=0 0.3 =0.5 ，階數(shù)為 9 ，為了驗證仿真需計算的沖激響應樣本為 10個。

44、仿真是正確的嗎？實驗名稱：數(shù)字濾波器的設計實驗目的和任務：（1）用窗口法設計滿足指標的FIR數(shù)字濾波器。（2）以模擬低通濾波器為原型設計IIR數(shù)字濾波器（3）選定一個信號濾波問題，設計數(shù)字濾波器，驗證濾波效果。實驗內(nèi)容：閱讀Page 91-93相應的函數(shù)和程序P7.1,完成Q7.1, Q7.5,Q7.6 閱讀Page 94-96相應的函數(shù), 完成Q7.9,Q7.13,Q7.14,Q7.20 sinc函數(shù)的功能與使用，可通過help查詢： Matlab-Help-Search-Function Name-輸入sinc 實驗指導書Page 49有sinc函數(shù)使用實例 幅度響應的分析： 通過DTFT

45、定義計算幅度響應 通過freqz函數(shù)分析實驗過程與結果分析： Q7.1 用MATLAB確定一個數(shù)字無線沖激響應低通濾波器所有四種類型的最低階數(shù)。指標如下：40kHz的抽樣率，4kHZ的通帶邊界頻率，8kHz的阻帶邊界頻率，0.5dB的通帶波紋，40dB的最小阻帶衰減。評論你的結果。根據(jù)題意：F=40kHz，F(xiàn)=4 kHz，F(xiàn)=8kHz，通帶增益R=0.5dB，阻帶增益R=40dB 可以得出= W = W (1) 根據(jù)上述數(shù)據(jù)和buttord函數(shù)N, Wn = buttord(0.2,0.4,0.5,40)得到巴特沃茲濾波器的最低階數(shù)N=8。Wn=0.2469 =0.2469(2) 根據(jù)上述數(shù)據(jù)

46、和cheb1ord函數(shù)N, Wn = cheb1ord(0.2,0.4,0.5,40)得到切比雪夫1型濾波器的最低階數(shù)N=5。Wn=0.2000 =0.2000(3) 根據(jù)上述數(shù)據(jù)和cheb2ord函數(shù)N, Wn = cheb1ord(0.2,0.4,0.5,40)得到切比雪夫2型濾波器的最低階數(shù)N=5。Wn=0.4000 =0.4000(4) 根據(jù)上述數(shù)據(jù)和ellipord函數(shù)N, Wn = ellipord (0.2,0.4,0.5,40)得到橢圓濾波器的最低階數(shù)N=4。Wn=0.2000 =0.2000Q7.5 通過運行P7.1設計巴特沃茲帶阻濾波器，寫出所產(chǎn)生的傳輸函數(shù)的準確表達式。濾

47、波器的指標是什么？使用MATLAB計算并繪制濾波器未畸變的相位相應及群延遲相應。MATLAB程序為：% Program P7_1Ws = 0.4 0.6; Wp = 0.2 0.8; Rp = 0.4; Rs = 50;N1, Wn1 = buttord(0.2, 0.4, 0.5, 40); num,den = butter(N1,Wn1,'stop'); disp('Numerator Coefficients are ');disp(num);disp('Denominator Coefficients are ');disp(den);%

48、 Compute the gain responseg, w = gain(num,den); % Plot the gain responseplot(w/pi,g);grid axis(0 1 -60 5);xlabel('omega /pi'); ylabel('Gain in dB');title('Gain Response of a Butterworth Bandstop Filter');得到：Numerator Coefficients areColumns 1 through 90.0493 0.0000 0.2465 0.0

49、000 0.4930 0.0000 0.4930 0.0000 0.2465Columns 10 through 110.0 0.0493Denominator Coefficients areColumns 1 through 91.0000 0.0000 -0.0850 0.0000 0.6360 0.0000 -0.0288 0.0000 0.0561Columns 10 through 110.0000 -0.0008因此表達式為：H (z)=指標為=0.2，=0.4，=0.6，=0.8，R=0.4dB，R=50dBQ7.6 修改P7.1后來設計符合習題Q7.1所給指標的切比雪夫1型低

50、通濾波器。寫出產(chǎn)生的傳輸函數(shù)的準確表達式。使用MATLAB計算并繪制濾波器未畸變的相位相應及群延遲響應。修改P7.1后的程序（切比雪夫1型低通濾波器）：% Program Q7_6% Design spec as given in Q7.1.FT = 40*103; % sampling freqFp = 4*103; % analog passband edge freqFs = 8*103; % analog stopband edge freqRp = 0.5; % max passband ripple, dBRs = 40; % min stopband attenuation, d

51、B% Convert spec to normalized digital frequenciesomega_p = 2*pi*Fp/FT;Wp = 2*Fp/FT; % omega_p/piomega_s = 2*pi*Fs/FT;Ws = 2*Fs/FT; % omega_s/pi% Estimate the Filter OrderN, Wn = cheb1ord(Wp, Ws, Rp, Rs);% Design the Filternum,den = cheby1(N,Rp,Wn);% Display the transfer functiondisp('Numerator C

52、oefficients are ');disp(num);disp('Denominator Coefficients are ');disp(den);% Compute the gain responseg, w = gain(num,den);% Plot the gain responsefigure(1);plot(w/pi,g);grid;axis(0 1 -60 5);xlabel('omega /pi'); ylabel('Gain in dB');title('Gain Response of a Type 1

53、Chebyshev Lowpass Filter');% Find and plot the phasefigure(2);w2 = 0:pi/511:pi;Hz = freqz(num,den,w2);Phase = unwrap(angle(Hz);plot(w2/pi,Phase);grid;xlabel('omega /pi'); ylabel('Unwrapped Phase (rad)');title('Unwrapped Phase Response of a Type 1 Chebyshev Lowpass Filter'

54、);% Find and plot the group delayfigure(3);GR = grpdelay(num,den,w2);plot(w2/pi,GR);grid;xlabel('omega /pi'); ylabel('Group Delay (sec)');title('Group Delay of a Type 1 Chebyshev Lowpass Filter');編寫的MATLAB程序如下（計算未畸變的相位響應及群延遲響應）：% Program Q7_6Ws = 0.4 0.6; Wp = 0.2 0.8; Rp = 0

55、.4; Rs = 50;% Estimate the Filter OrderN1, Wn1 = buttord(Wp, Ws, Rp, Rs);% Design the Filternum,den = butter(N1,Wn1,'stop');% Find the frequency response; find and plot unwrapped phasewp = 0:pi/1023:pi;wg = 0:pi/511:pi;Hz = freqz(num,den,wp);Phase = unwrap(angle(Hz);figure(1);plot(wp/pi,Phase);grid;% axis(0 1 a b);xlabel('omega /pi'); ylabel('Unwrapped Phase (rad)');t