第一篇：數(shù)字信號處理第四章介紹

第四章 線性時不變離散時間系統(tǒng)的頻域分析

一、傳輸函數(shù)和頻率響應(yīng) 例4.1傳輸函數(shù)分析 Q4.1 clear;M = input('Enter the filter length M: ');w = 0:2*pi/1023:2*pi;num =(1/M)*ones(1,M);den = [1];h = freqz(num, den, w);subplot(2,1,1)plot(w/pi,abs(h));grid title('Magnitude Spectrum |H(e^{jomega})|')xlabel('omega /pi');ylabel('Amplitude');subplot(2,1,2)plot(w/pi,angle(h));grid title('Phase Spectrum arg[H(e^{jomega})]')xlabel('omega /pi');ylabel('Phase in radians');

M=2

M=10

M=15

幅度譜為偶對稱，相位譜為奇對稱，這是一個低通濾波器。M越大，通帶越窄且過渡帶越陡峭。

Q4.2使用修改后的程序P3.1，計算并畫出當(dāng)w=[0,pi]時傳輸函數(shù)因果線性時不變離散時間系統(tǒng)的頻率響應(yīng)。它表示哪種類型的濾波器？ w = 0:pi/511:pi;

num = [0.15 0-0.15];den = [1-0.5 0.7];如下圖1這是一個帶通濾波器。的圖1

圖2 Q4.3對下面的傳輸函數(shù)重做習(xí)題Q4.2：，式（4.36）和式（4.37）給出的兩個濾波器之間的區(qū)別是什么？你將選擇哪一個濾波器來濾波，為什么？ w = 0:pi/511:pi;num = [0.15 0-0.15];den = [0.7-0.5 1];如上圖2也是一個帶通濾波器，這兩個濾波器的幅度譜是一樣的，相位譜不太一樣，我會選擇第一個帶通濾波器，因為它的相位譜更加平滑，相位失真小。

Q4.4 使用MATLAB計算并畫出當(dāng)w=[0,pi]時因果線性時不變離散時間系統(tǒng)的群延遲。系統(tǒng)的傳輸函數(shù)為clf;w = 0:pi/511:pi;num = [1-1.2 1];den = [1-1.3 1.04-0.222];h= grpdelay(num,den,w);plot(w/pi,h);xlabel('w/pi');ylabel('群延遲')。

Q4.5 使用Q3.50中編寫的程序，分別計算并畫出式（4.36）和式（4.37）確定的兩個濾波器的沖激響應(yīng)中的前一百個樣本。討論你的結(jié)果。clf;num = [0.15 0-0.15];den = [0.7-0.5 1];L = input('輸入樣本數(shù) L: ');[g t] = impz(num,den,L);stem(t,g);title(['前 ',num2str(L),' 脈沖響應(yīng)的樣本']);xlabel('時間序號 n');ylabel('h[n]');

（4.36）式（4.37）式

由圖可知：這些情節(jié)由impz給生成的因果的脈沖響應(yīng)實現(xiàn)的H(z)。我們觀察到Q4.3因果濾波器與H(z)在(4.36)穩(wěn)定,這意味著H[n]是絕對可和,我們看到交替和指數(shù)衰減的脈沖響應(yīng)。在另一方面,因果編檔人員與H(z)在(4.37)極點以外的單位圓,是不穩(wěn)定的。不足為奇的是,相應(yīng)的h[n]上圖顯示與n指數(shù)增長。

Q4.6 傳輸函數(shù)的極零點圖同樣能分析線性時不變離散時間系統(tǒng)的性質(zhì)。使用命令zplane可以很容易地得到系統(tǒng)的極零點圖。使用zplane分別生成式（4.36）和式（4.37）確定的兩個濾波器的極零點圖。討論你的結(jié)果。clf;num = [0.15 0-0.15];den = [1-0.5 0.7];[z p k] = tf2zpk(num,den);disp('Zeros:');disp(z);disp('Poles:');disp(p);input('Hit to continue...');[sos k] = zp2sos(z,p,k)input('Hit to continue...');zplane(z,p);式（4.36）

式（4.37）

由圖可知：過濾器在(4.36)在單位圓和兩極因此它的因果實現(xiàn)穩(wěn)定；較低的圖顯示過濾器(4.37)極點在單位圓外,其因果關(guān)系的實現(xiàn)是不穩(wěn)定的。

二、傳輸函數(shù)的類型 例4.2濾波器 Q4.7 clf;fc = 0.25;n = [-6.5:1:6.5];y = 2*fc*sinc(2*fc*n);k = n+6.5;stem(k,y);title('N = 14');axis([0 13-0.2 0.6]);xlabel('Time index n');ylabel('Amplitude');grid;

圖1 圖2 如圖1低通有限沖激濾波器的長度為14，決定濾波器長度的語句為n = [-6.5:1:6.5]，而控制截止頻率的參數(shù)是fc = 0.25。Q4.8 fc = 0.45;n = [-9.5:1:9.5];y = 2*fc*sinc(2*fc*n);k = n+9.5;stem(k,y);title('N = 20');axis([0 19-0.2 0.7]);xlabel('Time index n');ylabel('Amplitude');grid;修改參數(shù)fc和n，得到如上圖2，可知低通有限沖激濾波器的長度變?yōu)?0.Q4.9 clf;fc = 0.65;n = [-7.0:1:7.0];y = 2*fc*sinc(2*fc*n);k = n+7.0;stem(k,y);title('N = 14');axis([0 14-0.4 1.4]);xlabel('Time index n');ylabel('Amplitude');grid;

Q4.10 clear;N = input('Enter the filter time shift N: ');No2 = N/2;fc = 0.25;n = [-No2:1:No2];y = 2*fc*sinc(2*fc*n);w = 0:pi/511:pi;h = freqz(y, [1], w);plot(w/pi,abs(h));grid;title(strcat('|H(e^{jomega})|, N=',num2str(N)));xlabel('omega /pi');ylabel('Amplitude');

上圖依次分別為N=5,10,30,100的四幅圖，從這四幅圖可以看出隨著階數(shù)N的增大，低通濾波器的過渡帶越來越窄，阻帶衰減越來越快，濾波器越來越接近理想低通濾波器。Q4.11 clf;M = 2;num = ones(1,M)/M;[g,w] = gain(num,1);plot(w/pi,g);grid axis([0 1-50 0.5])xlabel('omega /pi');ylabel('Gain in dB');title(['M = ',num2str(M)])

可以驗證3dB截止頻率在π/2處。Q4.12 clear;K = input('Enter the number of sections K: ');Hz = [1];for i=1:K;Hz = conv(Hz,[1 1]);end;Hz =(0.5)^K * Hz;[g,w] = gain(Hz,1);ThreedB =-3*ones(1,length(g));t1 = 2*acos((0.5)^(1/(2*K)))*ones(1,512)/pi;t2 =-50:50.5/511:0.5;plot(w/pi,g,w/pi,ThreedB,t1,t2);grid;axis([0 1-50 0.5])xlabel('omega /pi');ylabel('Gain in dB');title(['K = ',num2str(K),';Theoretical omega_{c} = ',num2str(t1(1))]);

Q4.13 clear;M = input('Enter the filter length M: ');n = 0:M-1;num =(-1).^n.* ones(1,M)/M;[g,w] = gain(num,1);plot(w/pi,g);grid;axis([0 1-50 0.5]);xlabel('omega /pi');ylabel('Gain in dB');title(['M = ', num2str(M)]);

其3dB截止頻率約為0.82pi Q4.14 設(shè)計一個在0.45pi處具有3dB截止頻率wc的一階無限沖激響應(yīng)低通濾波器和一階無限沖激響應(yīng)高通濾波器。用MATLAB計算并畫出它們的增益響應(yīng)，驗證設(shè)計的濾波器是否滿足指標(biāo)。用MATLAB證明兩個濾波器是全通互補和功率互補的。

Q4.15 級聯(lián)10個式（4.15）所示一階無限沖激響應(yīng)低通濾波器，設(shè)計一個在0.3pi處具有3dB截止頻率wc的無限沖激響應(yīng)低通濾波器。把它與一個具有相同截止頻率的一階無限沖激響應(yīng)低通濾波器的增益響應(yīng)作比較。

Q4.16 設(shè)計一個中心頻率wo在0.61pi處、3dB帶寬為0.51pi的二階無限沖激響應(yīng)帶通濾波器。由于式（4.20）是α的二次方程，為了產(chǎn)生相同的3dB帶寬，參數(shù)α將有兩個數(shù)值，得到的傳輸函數(shù)HBP（z）也會有兩個不同的表達式。使用函數(shù)zplane可產(chǎn)生兩個傳輸函數(shù)的極零點圖，從中可以選擇一個穩(wěn)定的傳輸函數(shù)。用MATLAB計算并畫出你所設(shè)計的濾波器的增益響應(yīng)，并驗證它確實滿足給定的條件。用設(shè)計的穩(wěn)定無限沖激響應(yīng)帶通濾波器的傳輸函數(shù)的參數(shù)α和β，生成一個二階無限沖激響應(yīng)帶阻濾波器的傳輸函數(shù)HBS（z）。用MATLAB證明HBP(z)和HBS(z)都是全通互補和功率互補的。

Q4.17 用MATLAB計算并畫出一個梳狀濾波器的幅度響應(yīng)，該梳狀濾波器是在L取不同值的情況下，由式（4.40）給出的原型有限沖激響應(yīng)低通濾波器得到的。證明新濾波器的幅度響應(yīng)在k=0,1,2,3......，L-1.處有L個極小值，在處有L個極大值，Q4.18 用MATLAB計算并畫出一個梳狀濾波器的幅度響應(yīng)，該梳狀濾波器是在L取不同值的情況下，由式（4.42）在M=2時給出的原型有限沖激響應(yīng)低通濾波器得到的。確定這種梳狀濾波器沖激響應(yīng)的極大值和極小值的位置。

從這些情節(jié)我們觀察,梳狀濾波器極距為1kπ/L,山峰為(2k+1)π/L.Q4.19 clf;b = [1-8.5 30.5-63];num1 = [b 81 fliplr(b)];num2 = [b 81 81 fliplr(b)];num3 = [b 0-fliplr(b)];num4 = [b 81-81-fliplr(b)];n1 = 0:length(num1)-1;n2 = 0:length(num2)-1;subplot(2,2,1);stem(n1,num1);xlabel('Time index n');ylabel('Amplitude');grid;title('Type 1 FIR Filter');subplot(2,2,2);stem(n2,num2);xlabel('Time index n');ylabel('Amplitude');grid;title('Type 2 FIR Filter');subplot(2,2,3);stem(n1,num3);xlabel('Time index n');ylabel('Amplitude');grid;title('Type 3 FIR Filter');subplot(2,2,4);stem(n2,num4);xlabel('Time index n');ylabel('Amplitude');grid;title('Type 4 FIR Filter');pause subplot(2,2,1);zplane(num1,1);title('Type 1 FIR Filter');subplot(2,2,2);zplane(num2,1);title('Type 2 FIR Filter');subplot(2,2,3);zplane(num3,1);title('Type 3 FIR Filter');subplot(2,2,4);zplane(num4,1);title('Type 4 FIR Filter');disp('Zeros of Type 1 FIR Filter are');disp(roots(num1));disp('Zeros of Type 2 FIR Filter are');disp(roots(num2));disp('Zeros of Type 3 FIR Filter are');disp(roots(num3));disp('Zeros of Type 4 FIR Filter are');disp(roots(num4));

1型有限沖激響應(yīng)濾波器的零點是 Zeros of Type 1 FIR Filter are 2.9744 2.0888 0.9790 + 1.4110i 0.97900.4784i 0.4787 0.3362 2型有限沖激響應(yīng)濾波器的零點是 Zeros of Type 2 FIR Filter are 3.7585 + 1.5147i 3.75852.6623i-1.0000 0.0893 + 0.3530i 0.08930.0922i 3型有限沖激響應(yīng)濾波器的零點是 Zeros of Type 3 FIR Filter are 4.7627 1.6279 + 3.0565i 1.62790.2549i 0.2100 4型有限沖激響應(yīng)濾波器的零點是 Zeros of Type 4 FIR Filter are 3.4139 1.6541 + 1.5813i 1.65410.9973i 1.0000 0.3159 + 0.3020i 0.31592.0140i-1.2659 + 2.0135i-1.26590.3559i 0.2457 + 0.2126i 0.24572.0392i-1.0101 + 2.1930i-1.01010.3762i 0.2397 + 0.1934i 0.23971.2263i 1.0000 0.6576 + 0.7534i 0.65760.7803i 4型有限沖激響應(yīng)濾波器的零點是 Zeros of Type 4 FIR Filter are 2.0841 + 2.0565i 2.08411.9960i 1.0000-0.2408 + 0.3197i-0.24080.2399i Q4.21 用MATLAB 確定如下傳輸函數(shù)是否是有界實函數(shù)：一個有界實函數(shù)，求一個與

有著相同幅度的有界實函數(shù)。

它若不是由下圖可知：H1(z)不是有界實函數(shù)。

故H2(z)為

Q4.22 用MATLAB 確定如下傳輸函數(shù)是否是有界實函數(shù)：有界實函數(shù)，求一個與

有著相同幅度的有界實函數(shù)。

它若不是一個使用zplane我們觀察到的G1(z)在單位圓,因此傳遞函數(shù)是穩(wěn)定的。

Q4.23 用MATLAB產(chǎn)生如下兩個因果系統(tǒng)傳輸函數(shù)的極零點圖：，研究生成的極零點圖，你可以推斷它們的穩(wěn)定性么？

用Q4.6的程序做H1(z)

，Q4.24 用程序P4.4檢測Q4,23中兩個傳輸函數(shù)的穩(wěn)定性。這兩個傳輸函數(shù)哪一個是穩(wěn)定的？ % Program P4_4 clf;den = input('分母系數(shù) = ');ki = poly2rc(den);disp('穩(wěn)定性測試參數(shù)是');disp(ki);

由此我們可以總結(jié)出H1(z)穩(wěn)定，H2(z)不穩(wěn)定

Q4.25 用程序P4.4確定下面這個多項式的所有根是否都在單位圓內(nèi)：

由此看出，都在單位圓內(nèi)。

Q4.26 用程序P4.4確定下面這個多項式的所有根是否都在單位圓內(nèi)：

由此看出，都在單位圓內(nèi)。

第二篇：數(shù)字信號處理第七章介紹

第七章 數(shù)字濾波器設(shè)計

7.1：無限沖激響應(yīng)濾波器的階數(shù)的估計

Q7.1用MATTAB確定一個數(shù)字無限沖激響應(yīng)低通濾波器所有四種類型的最低階數(shù)。指標(biāo)如下:40 kHz的抽樣率，,4 kHz的通帶邊界頻率，8 kHz的阻帶邊界頻率，0.5 dB的通帶波紋，40 dB的最小阻帶衰減。評論你的結(jié)果。答：標(biāo)準(zhǔn)通帶邊緣角頻率Wp是：

標(biāo)準(zhǔn)阻帶邊緣角頻率Ws是：

理想通帶波紋Rp是0.5dB 理想阻帶波紋Rs是40dB 1.使用這些值得到巴特沃斯低通濾波器最低階數(shù)N=8，相應(yīng)的標(biāo)準(zhǔn)通帶邊緣頻率Wn是0.2469.2.使用這些值得到切比雪夫1型低通濾波器最低階數(shù)N=5，相應(yīng)的標(biāo)準(zhǔn)通帶邊緣頻率Wn是0.2000.3/使用這些值得到切比雪夫2型低通濾波器最低階數(shù)N=5，相應(yīng)的標(biāo)準(zhǔn)通帶邊緣頻率Wn是0.4000.4.使用這些值得到橢圓低通濾波器最低階數(shù)N=8，相應(yīng)的標(biāo)準(zhǔn)通帶邊緣頻率Wn是0.2000.從以上結(jié)果中觀察到橢圓濾波器的階數(shù)最低，并且符合要求。

Q7.2用MATLAB確定一個數(shù)字無限沖激響應(yīng)高通濾波器所有四種類型的最低階數(shù)。指標(biāo)如下:3500Hz的抽樣率，1050 Hz的通帶邊界頻率,600 Hz的阻帶邊界頻率，1 dB的通帶波紋，50 dB的最小阻帶衰減。評論你的結(jié)果 答：標(biāo)準(zhǔn)通帶邊緣角頻率Wp是：

標(biāo)準(zhǔn)阻帶邊緣角頻率Ws是：

理想通帶波紋Rp是1dB 理想阻帶波紋Rs是50dB 1.使用這些值得到巴特沃斯高通濾波器最低階數(shù)N=8，相應(yīng)的標(biāo)準(zhǔn)通帶邊緣頻率Wn是0.5646.2.使用這些值得到切比雪夫1型高通濾波器最低階數(shù)N=5，相應(yīng)的標(biāo)準(zhǔn)通帶邊緣頻率Wn是0.6000.3.使用這些值得到切比雪夫2型高通濾波器最低階數(shù)N=5，相應(yīng)的標(biāo)準(zhǔn)通帶邊緣頻率Wn是0.3429.4.使用這些值得到橢圓低通濾波器最低階數(shù)N=4，相應(yīng)的標(biāo)準(zhǔn)通帶邊緣頻率Wn是0.6000.從以上結(jié)果中觀察到橢圓濾波器的階數(shù)最低，并且符合要求。

Q7.3用MATLAB確定一個數(shù)字無限沖激響應(yīng)帶通濾波器所有四種類型的最低階數(shù)。指標(biāo)如下:7 kHz的抽樣率,1.4 kHz和2.1 kHz的通帶邊界頻率，1.05 kHz和2.45 kHz的阻帶邊界頻率，,0.4 dB的通帶波紋,50 dB的最小阻帶衰減。評論你的結(jié)果。答：標(biāo)準(zhǔn)通帶邊緣角頻率Wp是：

標(biāo)準(zhǔn)阻帶邊緣角頻率Ws是：

理想通帶波紋Rp是0.4dB 理想阻帶波紋Rs是50dB 1.使用這些值得到巴特沃斯帶通濾波器最低階數(shù)2N=18，相應(yīng)的標(biāo)準(zhǔn)通帶邊緣頻率Wn是 [0.3835 0.6165].2.使用這些值得到切比雪夫1型帶通濾波器最低階數(shù)2N=12，相應(yīng)的標(biāo)準(zhǔn)通帶邊緣頻率Wn是[0.4000 0.6000].3.使用這些值得到切比雪夫2型帶通濾波器最低階數(shù)2N=12，相應(yīng)的標(biāo)準(zhǔn)通帶邊緣頻率Wn是[0.3000 0.7000].4.使用這些值得到橢圓帶通濾波器最低階數(shù)2N=8，相應(yīng)的標(biāo)準(zhǔn)通帶邊緣頻率Wn是[0.4000 0.6000].從以上結(jié)果中觀察到橢圓濾波器的階數(shù)最低，并且符合要求。

Q7.4用MATLAB確定一個數(shù)字無限沖激響應(yīng)帶阻濾波器所有四種類型的最低階數(shù)。指標(biāo)如下:12 kHz的抽樣率，2.1 kHz和4.5 kHz的通帶邊界頻率,2.7 kHz和3.9 kHz的阻帶邊界頻率，0.6 dB的通帶波紋,45 dB的最小阻帶衰減。評論你的結(jié)果。

答：標(biāo)準(zhǔn)通帶邊緣角頻率Wp是：

標(biāo)準(zhǔn)阻帶邊緣角頻率Ws是：

理想通帶波紋Rp是0.6dB 理想阻帶波紋Rs是45dB 1.使用這些值得到巴特沃斯帶阻濾波器最低階數(shù)2N=18，相應(yīng)的標(biāo)準(zhǔn)通帶邊緣頻率Wn是[0.3873 0.7123].2.使用這些值得到切比雪夫1型帶阻濾波器最低階數(shù)2N=10，相應(yīng)的標(biāo)準(zhǔn)通帶邊緣頻率Wn是[0.3500 0.7500].3.使用這些值得到切比雪夫2型帶阻濾波器最低階數(shù)2N=10，相應(yīng)的標(biāo)準(zhǔn)通帶邊緣頻率Wn是[0.4500 0.6500].4.使用這些值得到橢圓帶阻濾波器最低階數(shù)2N=8，相應(yīng)的標(biāo)準(zhǔn)通帶邊緣頻率Wn是[0.3500 0.7500].從以上結(jié)果中觀察到橢圓濾波器的階數(shù)最低，并且符合要求。

7.2：無限沖激響應(yīng)濾波器設(shè)計

程序P7.1說明巴特沃斯帶阻濾波器的設(shè)計。% 巴特沃斯帶阻濾波器的設(shè)計

Ws = [0.4 0.6];Wp = [0.2 0.8];Rp = 0.4;Rs = 50;% 估計濾波器的階數(shù)

[N1, Wn1] = buttord(Wp, Ws, Rp, Rs);% 設(shè)計濾波器

[num,den] = butter(N1,Wn1,'stop');% 顯示傳輸函數(shù)

disp('分子系數(shù)是 ');disp(num);disp('分母系數(shù)是 ');disp(den);% 計算增益響應(yīng)

[g, w] = gain(num,den);% 繪制增益響應(yīng) plot(w/pi,g);grid axis([0 1-60 5]);xlabel('omega /pi');ylabel('增益，dB');title('巴特沃斯帶阻濾波器的增益響應(yīng)');Q7.5通過運行程序P7.1來設(shè)計巴特沃茲帶阻濾波器。寫出所產(chǎn)生的傳輸函數(shù)的準(zhǔn)確表達式。濾波器的指標(biāo)是什么，你的設(shè)計符合指標(biāo)嗎，使用MATLAB，計算并繪制濾波器的未畸變的相位響應(yīng)及群延遲響應(yīng)。答：表達式是：

濾波器參數(shù)是：

Wp1=0.2π，Ws1=0.4π，Ws2=0.6π，Wp2=0.8π，Rp=0.4dB,Rs=50dB.設(shè)計的濾波器增益響應(yīng)如下：

從圖中可以總結(jié)出設(shè)計符合指標(biāo)。

濾波器的未畸變的相位響應(yīng)及群延遲響應(yīng)如下：

Q7.6修改程序P7.1來設(shè)計符合習(xí)題Q7.1所給指標(biāo)的切比雪夫1型低通濾波器。寫出所產(chǎn)生的傳輸函數(shù)的準(zhǔn)確表達式。你的設(shè)計符合指標(biāo)嗎?使用MATLAB,計算并繪制濾波器的未畸變的相位響應(yīng)及群延遲響應(yīng)。答：表達式如下：

設(shè)計的濾波器增益響應(yīng)如下：

從圖中可以總結(jié)出設(shè)計符合指標(biāo)。

濾波器的未畸變的相位響應(yīng)及群延遲響應(yīng)如下：

Q7.7修改程序P7.1來設(shè)計符合習(xí)題Q7.2所給指標(biāo)的切比雪夫2型高通濾波器。寫出所產(chǎn)生的傳輸函數(shù)的準(zhǔn)確表達式。你的設(shè)計符合指標(biāo)嗎?使用MATLAB,計算并繪制濾波器的未畸變的相位響應(yīng)及群延遲響應(yīng)。答：表達式如下：

設(shè)計的濾波器增益響應(yīng)如下：

從圖中可以總結(jié)出設(shè)計符合指標(biāo)。濾波器的未畸變的相位響應(yīng)及群延遲響應(yīng)如下：

Q7.8修改程序P7.1來設(shè)計符合習(xí)題Q7.3所給指標(biāo)的橢圓帶通濾波器。寫出所產(chǎn)生的傳輸函數(shù)的準(zhǔn)確表達式。你的設(shè)計符合指標(biāo)嗎，使用MATLAB，計算井繪制濾波器的未畸變的相位響應(yīng)及群延遲響應(yīng)。答：表達式如下：

設(shè)計的濾波器增益響應(yīng)如下：

從圖中可以總結(jié)出設(shè)計符合指標(biāo)。

濾波器的未畸變的相位響應(yīng)及群延遲響應(yīng)如下：

7.3：吉布斯現(xiàn)象

Q7.9使用函數(shù)sinc編寫一個MATLAB程序，以產(chǎn)生截止頻率在Wc= 0.4π處、長度分別為81,61,41和21的四個零相位低通濾波器的沖激響應(yīng)系數(shù)，然后計算并畫出它們的幅度響應(yīng)。使用冒號“:”運算符從長度為81的濾波器的沖激響應(yīng)系數(shù)中抽出較短長度濾波器的沖激響應(yīng)系數(shù)。在每一個濾波器的截止頻率兩邊研究頻率響應(yīng)的擺動行為。波紋的數(shù)量與濾波器的長度之間有什么關(guān)系?最大波紋的高度與濾波器的長度之間有什么關(guān)系?你將怎樣修改上述程序以產(chǎn)生一個偶數(shù)長度的零相位低通濾波器的沖激響應(yīng)系數(shù)? 答：長度為81時幅度響應(yīng)如下：

長度分別為61,41和21的幅度響應(yīng)如下：

從中可以觀察到由于吉布斯現(xiàn)象產(chǎn)生的幅度響應(yīng)的擺動行為。

波紋的數(shù)量與濾波器的長度之間的關(guān)系——波紋的數(shù)量減少與長度成正比。最大波紋的高度與濾波器的長度之間的關(guān)系——最大波紋的高度與長度無關(guān)。

Q7.10使用函數(shù)sinc編寫一個MATLAB程序，以產(chǎn)生一個截止頻率在Wc= 0.4π處、長度為45的零相位高通濾波器的沖激響應(yīng)系數(shù)，計算并畫出其幅度響應(yīng)。在每一個濾波器的截止頻率兩邊研究頻率響應(yīng)的擺動行為。你將怎樣修改上述程序以產(chǎn)生一個偶數(shù)長度的零相位高通濾波器的沖激響應(yīng)系數(shù)? 答：長度為45時幅度響應(yīng)如下：

從中可以觀察到由于吉布斯現(xiàn)象產(chǎn)生的幅度響應(yīng)擺動行為。

在這種情況下你不能改變長度。原因：這是一個零相位濾波器，這意味著它也是一個線性相位濾波器，因為零相是一種特殊的線性相位的子集?，F(xiàn)在，理想的有限脈沖響應(yīng)長度甚至有對稱的中點h[n]。使其成了一個線性相位FIR濾波器。二型濾波器不可能是高通濾波器，因為必須在z=-1處有零點，意味著w=+-π。

Q7.11編寫一個MATLAB程序，以產(chǎn)生長度分別為81,61,41和21的四個零相位微分器的沖激響應(yīng)系數(shù)，計算并畫出它們的幅度響應(yīng)。下面的代碼段顯示了怎樣產(chǎn)生一個長度為2M+1的微分器。n=1:M;b=cos(pi*n)./n;num=[-fliplr(b)0 b];對于每種情況,研究微分器的頻率響應(yīng)的擺動行為。波紋的數(shù)量與微分器的長度之間有什么關(guān)系，最大波紋的高度與濾波器的長度之間有什么關(guān)系? 答：幅度響應(yīng)分別如下：

從中可以觀察到由于吉布斯現(xiàn)象產(chǎn)生的幅度響應(yīng)的擺動行為。波紋的數(shù)量與微分器的長度之間的關(guān)系——兩者成正比。

最大波紋的高度與濾波器的長度之間的關(guān)系——兩者間沒有關(guān)系。

Q7.12編寫一個MA11AB程序，以產(chǎn)生長度分別為81,61.41和21的四個離散時間希爾伯特變換器的沖激響應(yīng)系數(shù)，計算并畫出它們的幅度響應(yīng)。下面的代碼段顯示了怎樣產(chǎn)生一個長度為2M十1的希爾伯特變換器。n=1:M;c=sin(pi*n)./2;b=2*(c.*c)./(pi*n);num=[-fliplr(b)0 b];對于每種情況，研究希爾伯特變換器的頻率響應(yīng)的擺動行為。波紋的數(shù)量與希爾伯特變換器的長度之間有什么關(guān)系?最大波紋的高度與濾波器的長度之間有什么關(guān)系? 答：幅度響應(yīng)如下：

從中可以觀察到由于吉布斯現(xiàn)像產(chǎn)生的幅度響應(yīng)的擺動行為。波紋的數(shù)量與希爾伯特變換器的長度之間的關(guān)系——兩者成正比。最大波紋的高度與濾波器的長度之間的關(guān)系——兩者無關(guān)系。7.4：有限沖激響應(yīng)濾波器的階數(shù)估計

Q7.13 線性相位低通FIR濾波器的階數(shù)估算，參數(shù)如下:

p =2 kHz, s =2.5 kHz，p = 0.005，s = 0.005, FT = 10kHz 使用 kaiord 的結(jié)果為N = 46 使用 ceil 命令的目的是朝正方向最接近整數(shù)方向取整。使用nargin命令的目的是表明函數(shù)M文件體內(nèi)變量的數(shù)目。

Q7.14(a)線性相位FIR濾波器的階數(shù)估算，其中采樣頻率改為FT = 20 kHz，則結(jié)果為 N=91。(b)線性相位FIR(c)線性相位FIR從上述結(jié)果和7.13的對比我們可以觀察到：

濾波器階數(shù)和采樣頻率的關(guān)系為–對于一個給定的模擬過渡帶寬，采樣頻率的增加導(dǎo)致估算階數(shù)也相應(yīng)增加，朝下一個整數(shù)取整。

其中模擬過渡帶寬|Fp-Fs|和Δω的關(guān)系：Δω=2pi*|Fp-Fs|/FT。因此增加FT會減小Δω。

濾波器階數(shù)和通帶波紋寬度的關(guān)系為估計的階數(shù)大致和log（底數(shù)為10）成比例的擴散。濾波器階數(shù)和過渡帶寬度的關(guān)系為在舍入的時候，階數(shù)隨著過渡帶寬成比例的改變。有兩個因素增加過渡帶寬來分割順序。

Q7.15 線性相位FIR低通濾波器階數(shù)的估算，其中濾波器滿足7.13給的規(guī)格，使用kaiserord的結(jié)果為N=54 正確結(jié)果：kaiserord([2000 2500],[1 0],[0.005 0.005],10000)將上述結(jié)果和7.13比較我們觀察到：用凱瑟來估算階數(shù)是較小的。因為凱瑟使用了一個不同的近似估計。這種估計經(jīng)常和FIR設(shè)計的凱瑟窗一起用。

Q7.16 線性相位FIR低通濾波器的階數(shù)估算滿足的規(guī)格和7.13中的一樣，使用remezord函數(shù)的結(jié)果為N=47.正確結(jié)果：firpmord([2000 2500],[1 0],[0.005 0.005],10000)通過和7.13和7.15比較我們可以觀察到：在這里，firpmord給了一個比凱爾更大比凱瑟更小一點的結(jié)果。使用凱爾則更接近與一般情況。而使用凱瑟和firpmord則有專門的用途。Q7.17 線性相位帶通FIR濾波器的階數(shù)估算滿足如下規(guī)格：通帶邊界為1.8和3.6kHz，阻帶邊界為1.2kHz到4.2kHz

p = 0.01,阻帶波紋

s = 0.02,FT = 12 kHz。= 0.002p

= 0.002 結(jié)果為 N=57。s

s = 2.3 kHz，結(jié)果為N=76.使用kaiord 函數(shù)求得的結(jié)果為：通帶波紋δp= 0.1，得到的結(jié)果為：kaiord([1800 3600],[1200 4200],0.1,0.02,12000)，N=20。但是當(dāng)δp= 0.01時結(jié)果為：kaiord([1800 3600],[1200 4200],0.01,0.02,12000)，得到的N=33。所以答案不唯一，可以選擇其中一個。Q7.18 線性相位帶通FIR濾波器的階數(shù)估算，其中FIR濾波器的規(guī)格和7.17一樣，則使用kaiserord的結(jié)果為同樣，它也有矛盾。當(dāng)使用δp= 0.1時，得到的結(jié)果為：kaiserord([1200 1800 3600 4200],[0 1 0],[0.02 0.1 0.02],12000)，則N=37.當(dāng)用δp= 0.01時，結(jié)果為：kaiserord([1200 1800 3600 4200],[0 1 0],[0.02 0.01 0.02],12000)，此時N=45.和7.17的結(jié)果比較我們觀察到通過kaiserord函數(shù)估計的階數(shù)要更高，但如果你要設(shè)計Kaiser窗的話則結(jié)果更精確。

Q7.19 線性相位帶通FIR濾波器的階數(shù)估算，其中FIR濾波器的規(guī)格和7.17一樣，使用函數(shù)remezord。

當(dāng)取δp= 0.01時，結(jié)果為firpmord([1200 1800 3600 4200],[0 1 0],[0.02 0.1 0.02],12000)，此時N=22.而如果δp= 0.01，則結(jié)果為：firpmord([1200 1800 3600 4200],[0 1 0],[0.020.01 0.02],12000)，此時N=35.可以從中任意選擇。

和7.17和7.18的結(jié)果比較我們可以觀察到通過firpmord來估算的階數(shù)在另外兩個的中間，在設(shè)計Parks-McClellan時更準(zhǔn)確。

7.5有限沖激響應(yīng)濾波器設(shè)計

Q7.20 使用matlab程序設(shè)計并畫出線性相位FIR濾波器增益和相位反應(yīng)，使用fir1如下。通過使用函數(shù)kaiserord.來估計濾波器階數(shù)，輸出結(jié)果為濾波器的系數(shù)。低通濾波器滿足7.20所要求的規(guī)格的系數(shù)如下：

增益和相位響應(yīng)如下:

從增益圖像我們可以知道這個設(shè)計不能滿足規(guī)格.這個濾波器滿足規(guī)格的階數(shù)為N=66.為了滿足規(guī)格，圖如下：

Q7.21 漢寧窗：

布萊克曼窗：

切比雪夫窗：

Q7.22 程序如下：

% Program Q7_22 % Use Parks-McClellan to design a linear phase Lowpass % FIR Digital Filter meeting the design specification given % in Q7.13.%Compute and plot the gain function.%Compute and plot the unwrapped phase response.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear;% Design spec as given in Q7.13.Fp = 2*10^3;Fs = 2.5*10^3;FT = 10*10^3;Rp = 0.005;Rs = 0.005;% Estimate the filter order and print to console N = kaiord(Fp,Fs,Rp,Rs,FT)% Design the filter using Parks-McClellan Wp = 2*Fp/FT;% These freqs are normalized: they go Ws = 2*Fs/FT;% zero to one, not zero to pi.F = [0 Wp Ws 1];A = [1 1 0 0];h = firpm(N,F,A);% Show the Numerator Coefficients disp('Numerator Coefficients are ');disp(h);% Compute and plot the gain response [g, w] = gain(h,[1]);figure(1);plot(w/pi,g);grid;xlabel('omega /pi');ylabel('Gain in dB');title('Gain Response');% Compute the frequency response w2 = 0:pi/511:pi;Hz = freqz(h,[1],w2);% TEST: did we meet the spec? MagH = abs(Hz);T1 = 1.005*ones(1,length(w2));T2 = 0.995*ones(1,length(w2));T3 = 0.005*ones(1,length(w2));figure(4);plot(w2/pi,MagH,w2/pi,T1,w2/pi,T2,w2/pi,T3);grid;% Find and plot the phase figure(2);Phase = angle(Hz);plot(w2/pi,Phase);grid;xlabel('omega /pi');ylabel('Phase(rad)');title('Phase Response');figure(3);UPhase = unwrap(Phase);plot(w2/pi,UPhase);grid;xlabel('omega /pi');ylabel('Unwrapped Phase(rad)');title('Unwrapped Phase Response');低通濾波器系數(shù)：

0.0028-0.0022-0.0046-0.0006 0.0053 0.0019-0.0073-0.0058 0.0079 0.0106-0.0069-0.0170 0.0032 0.0243 0.0045-0.0319-0.0182 0.0390 0.0422-0.0448-0.0924 0.0486 0.3136 0.4501 0.3136 0.0486-0.0924-0.0448 0.0422 0.0390-0.0182-0.0319 0.0045 0.0243 0.0032-0.0170-0.0069 0.0106 0.0079-0.0058-0.0073 0.0019 0.0053-0.0006-0.0046-0.0022 0.0028 增益和相位響

從圖中可以看出此時的濾波器不滿足指標(biāo)。欲滿足指標(biāo)，應(yīng)調(diào)節(jié)N=47.Q7.23 用凱澤窗設(shè)計一個有限沖激響應(yīng)低通濾波器。程序：

% Program Q7_23 % Use Kaiser window to design a linear phase Lowpass % FIR Digital Filter meeting the design specification given % in Q7.23.% % It is clear from the statement of the question that Mitra % wants us to use(7.36)and(7.37)for this problem.That % isn't the greatest thing to try because kaiserord already does % exactly what we need....but that's Q7_24!So here goes!%Compute and plot the gain function.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear;% Design spec as given in Q7.23.Wp = 0.31;Ws = 0.41;Wn = Wp +(Ws-Wp)/2;As = 50;Ds = 10^(-As/20);Dp = Ds;%Kaiser window design has equal ripple in % passband and stopband.% estimate order using(7.37)if As > 21 N = ceil((As-7.95)*2/(14.36*(abs(Wp-Ws)))+1)else N = ceil(0.9222*2/abs(Wp-Ws)+1)end % Use(7.36)to get Beta if As > 50 BTA = 0.1102*(As-8.7);elseif As >= 21 BTA = 0.5842*(As-21)^0.4+0.07886*(As-21);else BTA = 0;end Win = kaiser(N+1,BTA);h = fir1(N,Wn,Win);% Show the Numerator Coefficients disp('Numerator Coefficients are ');disp(h);% Compute and plot the gain response [g, w] = gain(h,[1]);figure(1);plot(w/pi,g);grid;axis([0 1-80 5]);xlabel('omega /pi');ylabel('Gain in dB');title('Gain Response');% Compute the frequency response w2 = 0:pi/511:pi;Hz = freqz(h,[1],w2);% Find and plot the phase figure(2);Phase = angle(Hz);plot(w2/pi,Phase);grid;xlabel('omega /pi');ylabel('Phase(rad)');title('Phase Response');figure(3);UPhase = unwrap(Phase);plot(w2/pi,UPhase);grid;xlabel('omega /pi');ylabel('Unwrapped Phase(rad)');title('Unwrapped Phase Response');低通濾波器系數(shù)：

0.0003 0.0008 0.0003-0.0011-0.0017 0.0000 0.0026 0.0027-0.0010-0.0049-0.0035 0.0033 0.0080 0.0034-0.0074-0.0119-0.0018 0.0140 0.0161-0.0027-0.0241-0.0201 0.0127 0.0406 0.0236-0.0354-0.0754-0.0258 0.1214 0.2871 0.3597 0.2871 0.1214-0.0258-0.0754-0.0354 0.0236 0.0406 0.0127-0.0201-0.0241-0.0027 0.0161 0.0140-0.0018-0.0119-0.0074 0.0034 0.0080 0.0033-0.0035-0.0049-0.0010 0.0027 0.0026 0.0000-0.0017-0.0011 0.0003 0.0008 0.0003 增益和相位響應(yīng)如下：

從圖中可以看出設(shè)計的濾波器滿足要求。N=60.Q7.24 用函數(shù)kaiserord和firl重做習(xí)題Q7.23 程序：

% Use Kaiser window to design a linear phase Lowpass % FIR Digital Filter meeting the design specification given % in Q7.23.Use kaiserord and fir1.%Compute and plot the gain function.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear;% Design spec as given in Q7.23.Wp = 0.31;Ws = 0.41;As = 50;Ds = 10^(-As/20);% Design the Filter F = [Wp Ws];A = [1 0];DEV = [Ds Ds];[N,Wn,BTA,Ftype] = kaiserord(F,A,DEV);Win = kaiser(N+1,BTA);h = fir1(N,Wn,Ftype,Win);% Show the Numerator Coefficients disp('Numerator Coefficients are ');disp(h);% Compute and plot the gain response [g, w] = gain(h,[1]);figure(1);plot(w/pi,g);grid;axis([0 1-80 5]);xlabel('omega /pi');ylabel('Gain in dB');title('Gain Response');% Compute the frequency response w2 = 0:pi/511:pi;Hz = freqz(h,[1],w2);% Find and plot the phase figure(2);Phase = angle(Hz);plot(w2/pi,Phase);grid;xlabel('omega /pi');ylabel('Phase(rad)');title('Phase Response');figure(3);UPhase = unwrap(Phase);plot(w2/pi,UPhase);grid;xlabel('omega /pi');ylabel('Unwrapped Phase(rad)');title('Unwrapped Phase Response');參數(shù)如下：

Wp ? 0.31;Ws ? 0.41;As ? 50 dB

增益和相位響應(yīng)如下：

從圖中可以看出設(shè)計的濾波器滿足要求。N=59.Q7.25 用fir2設(shè)計一個95階有限沖激響應(yīng)濾波器。

程序：

% Program Q7_25 % Use fir2 to design a linear phase Lowpass % FIR Digital Filter meeting the design specification given % in Q7.23.%Print out the numerator coefficients % for the transfer function.%-Compute and plot the gain function.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear;% Design spec as given in Q7.17.F = [1200 1800 3600 4200];A = [0 1 0];DEV = [0.02 0.1 0.02];Fs = 12000;Dp = 0.1;Ds = 0.02;[N,Wn,BTA,FILTYPE] = kaiserord(F,A,DEV,Fs);N % firpm setup F2 = 2*[0 1200 1800 3600 4200 6000]/Fs;A2 = [0 0 1 1 0 0];wgts = max(Dp,Ds)*[1/Ds 1/Dp 1/Ds];h = firpm(N,F2,A2,wgts);% Show the Numerator Coefficients disp('Numerator Coefficients are ');disp(h);% Compute and plot the gain response [g, w] = gain(h,[1]);figure(1);plot(w/pi,g);grid;axis([0 1-80 5]);xlabel('omega /pi');ylabel('Gain in dB');title('Gain Response');% Compute the frequency response w2 = 0:pi/511:pi;Hz = freqz(h,[1],w2);% Find and plot the phase figure(2);Phase = angle(Hz);plot(w2/pi,Phase);grid;xlabel('omega /pi');ylabel('Phase(rad)');title('Phase Response');figure(3);UPhase = unwrap(Phase);plot(w2/pi,UPhase);grid;xlabel('omega /pi');ylabel('Unwrapped Phase(rad)');title('Unwrapped Phase Response');

增益響應(yīng)： 相位響應(yīng)：

從增益響應(yīng)的圖像中可以看出，此濾波器滿足指標(biāo)。N=37.Q7.27 用remez設(shè)計具有如下指標(biāo)的有限沖激響應(yīng)帶通濾波器。

程序：

% Program Q7_27 % Use kaiserord and firpm to design the linear phase bandpass % FIR Digital Filter specified in Q7.27.% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear;% Design spec as given in Q7.27.Fs1 = 1500;Fp1 = 1800;Fp2 = 3000;Fs2 = 4200;Fs = 12000;Dp = 0.1;Ds = 0.02;F = [Fs1 Fp1 Fp2 Fs2];A = [0 1 0];DEV = [Ds Dp Ds];[N,Wn,BTA,FILTYPE] = kaiserord(F,A,DEV,Fs);% firpm setup ws1 = 2*Fs1/Fs;wp1 = 2*Fp1/Fs;wp2 = 2*Fp2/Fs;ws2 = 2*Fs2/Fs;F2 = [0 ws1 wp1 wp2 ws2 1];A2 = [0 0 1 1 0 0];wgts = max(Dp,Ds)*[1/Ds 1/Dp 1/Ds];h = firpm(N,F2,A2,wgts);% Show the Numerator Coefficients disp('Numerator Coefficients are ');disp(h);% Compute and plot the gain response [g, w] = gain(h,[1]);figure(1);plot(w/pi,g);grid;axis([0 1-80 5]);xlabel('omega /pi');ylabel('Gain in dB');title('Gain Response');% Compute the frequency response w2 = 0:pi/511:pi;Hz = freqz(h,[1],w2);% Find and plot the phase figure(2);Phase = angle(Hz);plot(w2/pi,Phase);grid;xlabel('omega /pi');ylabel('Phase(rad)');title('Phase Response');figure(3);UPhase = unwrap(Phase);plot(w2/pi,UPhase);grid;xlabel('omega /pi');ylabel('Unwrapped Phase(rad)');title('Unwrapped Phase Response');figure(4);% Add lines to the plot to help determine if the spec was met.hold on;tmpY = 0:1.4/4:1.4;tmpX = ones(1,length(tmpY))*wp1;plot(tmpX,tmpY, 'r-');% vertical line at passband edge freq tmpX = ones(1,length(tmpY))*wp2;plot(tmpX,tmpY, 'r-');% vertical line at passband edge freq tmpX = ones(1,length(tmpY))*ws1;plot(tmpX,tmpY, 'g-');% vertical line at stopband edge freq tmpX = ones(1,length(tmpY))*ws2;plot(tmpX,tmpY, 'g-');% vertical line at stopband edge freq tmpY = ones(1,length(w))*(Dp);plot(w/pi,tmpY, 'r-');% horizontal line at Dp tmpY = ones(1,length(w))*(Ds);plot(w/pi,tmpY, 'g-');% horizontal line at Ds % now plot the Frequency Response plot(w2/pi,abs(Hz));grid;hold off;增益響應(yīng)：

從增益響應(yīng)的圖像可以看出，該濾波器不滿足設(shè)計指標(biāo)。用kaiserord估計濾波器的階數(shù)N=73。

第三篇：數(shù)字信號處理第六章介紹

第六章

數(shù)字濾波器結(jié)構(gòu)

6.1：級聯(lián)的實現(xiàn)

num = input('分子系數(shù)向量 = ');den = input('分母系數(shù)向量 = ');[z,p,k] = tf2zp(num,den);sos = zp2sos(z,p,k)Q6.1使用程序P6.1，生成如下有限沖激響應(yīng)傳輸函數(shù)的一個級聯(lián)實現(xiàn): H1(z)=2+10z^(-1)+23z^(-2)+34z^(-3)+31z^(-4)+16 z^(-5)+4z^(-6)畫出級聯(lián)實現(xiàn)的框圖。H1(z)是一個線性相位傳輸函數(shù)嗎? 答：運行結(jié)果：

sos = zp2sos(z,p,k)Numerator coefficient vector = [2,10,23,34,31,16,4] Denominator coefficient vector = [1] sos = 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)框圖：

H1（z）不是一個線性相位傳輸函數(shù)，因為系數(shù)不對稱。

Q6.2使用程序P6.1，生成如下有限沖激響應(yīng)傳輸函數(shù)的一個級聯(lián)實現(xiàn): H2(z)=6+31z^(-1)+74z^(-2)+102z^(-3)+74z^(-4)+31 z^(-5)+6z^(-6)畫出級聯(lián)實現(xiàn)的框圖。H2(z)是一個線性相位傳輸函數(shù)嗎?只用4個乘法器生成H2(z)的一級聯(lián)實現(xiàn)。顯示新的級聯(lián)結(jié)構(gòu)的框圖。

Numerator coefficient vector = [6,31,74,102,74,31,6] Denominator coefficient vector = [1] sos = 6.0000 15.0000 6.0000 1.0000 0 0 1.0000 2.0000 3.0000 1.0000 0 0 1.0000 0.6667 0.3333 1.0000 0 0 級聯(lián)框圖：

H2（z）是一個線性相位傳輸函數(shù)。只用四個乘法器生成級聯(lián)框圖：

6.2：級聯(lián)和并聯(lián)實現(xiàn)

Q6.3使用程序P6.1生成如下因果無限沖激響應(yīng)傳輸函數(shù)的級聯(lián)實現(xiàn)： 畫出級聯(lián)實現(xiàn)的框圖。答：

Numerator coefficient vector = [3,8,12,7,2,-2] Denominator coefficient vector = [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.2500 1.0000 1.0000 1.0000 1.0000 0.5000 0.5000

級聯(lián)實現(xiàn)框圖：

Q6.4使用程序P6.1生成如下因果無限沖激響應(yīng)傳輸函數(shù)的級聯(lián)實現(xiàn)：

畫出級聯(lián)實現(xiàn)的框圖。答：級聯(lián)實現(xiàn)框圖：

程序P6.2生成兩種類型的并聯(lián)實現(xiàn) num = input('分子系數(shù)向量 = ');den = input('分母系數(shù)分量 = ');[r1,p1,k1] = residuez(num,den);[r2,p2,k2] = residue(num,den);disp('并聯(lián)I型')disp('留數(shù)是');disp(r1);disp('極點在');disp(p1);disp('常數(shù)');disp(k1);disp('并聯(lián)II型')disp('留數(shù)是');disp(r2);disp('極點在');disp(p2);disp('常數(shù)');disp(k2);Q6.5使用程序P6.2生成式(6.27)所示因果無限沖激響應(yīng)傳輸函數(shù)的兩種不同并聯(lián)形式實現(xiàn)。畫出兩種實現(xiàn)的框圖。答：并聯(lián)I型框圖：

并聯(lián)II型框圖：

Q6.6使用程序P6.2生成式(6.28)所示因果無限沖激響應(yīng)傳輸函數(shù)的兩種不同并聯(lián)形式實現(xiàn)。畫出兩種實現(xiàn)的框圖。答：并聯(lián)I型框圖：

并聯(lián)II型框圖：

6.3：全通傳輸函數(shù)的實現(xiàn)

Q6.7使用程序P4.4生成如下全通傳輸函數(shù)的級聯(lián)格型實現(xiàn):

As(z)是一個穩(wěn)定的傳輸函數(shù)嗎? 答：

運行結(jié)果：

k(5)= 0.0625 k(4)= 0.2196 k(3)= 0.4811 k(2)= 0.6837 k(1)= 0.6246

2從{ki}的值我們可以得到傳輸函數(shù)A5(z)是穩(wěn)定的，因為對所有的1

Q6.8使用程序P4.4生成如下全通傳輸函數(shù)的級聯(lián)格型實現(xiàn):

A6(z)足一個穩(wěn)定的傳輸函數(shù)嗎? 答：得到A6(z)的{ki}值如下：

k(6)= 0.0278 k(5)= 0.1344 k(4)= 0.3717 k(3)= 0.5922 k(2)= 0.7711 k(1)= 0.8109

從{ki}的值可以得到傳輸函數(shù)A6(z)是穩(wěn)定的，因為反饋系數(shù)的平均幅值小于整體。

Q6.9 使用l型和2型全通項生成式(6.29)所示全通傳輸函數(shù)的典范級聯(lián)實現(xiàn)。顯示實現(xiàn)的框圖。在最終的結(jié)構(gòu)中，乘法器的總數(shù)是多少? 答：全通因子如下所示：

使用1型和2型全通項生成所示全通函數(shù)的典范級聯(lián)實現(xiàn)，實現(xiàn)的結(jié)構(gòu)框圖如下：

整體結(jié)構(gòu)中乘法器的總數(shù)是5.Q6.10 用zp2sos 我們可以得到 A6(z)的因子如下：

sos = 0.0278 0.0556 0.1111 1.0000 0.5000 0.2500 1.0000 2.0000 3.0000 1.0000 0.6667 0.3333 1.0000 3.0000 3.0000 1.0000 1.0000 0.3333 從上面因子可以分解 A6(z)為低階的全通因子:

使用2型的全通項生成A6(z)的典范級聯(lián)實現(xiàn)框圖如下：

整體結(jié)構(gòu)中乘法器的總數(shù)是6。6.4：無限沖激響應(yīng)傳輸函數(shù)的Gary-Markel實現(xiàn)

num = input('分子系數(shù)向量 = ');den = input('分母系數(shù)向量 = ');N = length(den)-1;% 分母多項式的階數(shù) k = ones(1,N);a1 = den/den(1);alpha = num(N+1:-1:1)/den(1);for ii = N:-1:1, alpha(N+2-ii:N+1)= alpha(N+2-ii:N+1)-alpha(N-ii+1)*a1(2:ii+1);k(ii)= a1(ii+1);a1(1:ii+1)=(a1(1:ii+1)-k(ii)*a1(ii+1:-1:1))/(1-k(ii)*k(ii));end disp('格型參數(shù)是');disp(k)disp('前饋乘法器是');disp(alpha)Q6.11 使用程序 P6_3我們通過IIR將Q6.3給的正向傳輸函數(shù)H1(z)的 Gray-Markel級聯(lián)格型實現(xiàn)參數(shù)如下： 晶格參數(shù)和前饋乘數(shù)分別如下：

對應(yīng)Gray-Markel的結(jié)構(gòu)框圖如下：

使用程序P6_3，從這些格型參數(shù)可以得到傳輸函數(shù)H1(z)是穩(wěn)定的，因為所有格型參數(shù)的平方值比整體的小。

Q6.12 使用程序 P6_3我們通過IIR將Q6.4給的正向傳輸函數(shù)H2(z)的 Gray-Markel級聯(lián)格型實現(xiàn)參數(shù)如下：

對應(yīng)Gray-Markel的結(jié)構(gòu)框圖如下：

使用程序P6_3，從這些格型參數(shù)可以得到傳輸函數(shù)H2(z)是穩(wěn)定的，因為所有格型參數(shù)的平方值比整體的小。

Q6.13使用函數(shù)tf2latc編寫出一個MATLAB程序，以生成一個因果無限沖激響應(yīng)傳輸函數(shù)的GrayMarkel實現(xiàn)。用該程序?qū)崿F(xiàn)式(6.27)所示的傳輸函數(shù)。你的結(jié)果與習(xí)題6.11中得到的結(jié)果相符嗎？使用函數(shù)1atc2tf由向量k和alpha確定傳輸函數(shù)。所得到的傳輸函數(shù)和式(6.27)給出的傳輸函數(shù)相同嗎？ 答：程序如下： format long num = input('Numerator coefficient vector = ');den = input('Denominator coefficient vector = ');num = num/den(1);% normalize upstairs and down by d0.den = den/den(1);% here is the lattice/ladder realization from the transfer fcn: [k,alpha] = tf2latc(num,den)% now check inversion disp('Check of Lattice/Ladder Inversion:');[num2,den2] = latc2tf(k,alpha)運行結(jié)果如下： k = 0.62459686089013 0.68373782742919 0.48111942348398 0.21960784313725 0.06250000000000 alpha =-0.01982100623522-0.09085169508677 0.***849 0.16053921568627 0.31250000000000-0.12500000000000 結(jié)果與習(xí)題6.11中得到的結(jié)果相符。Q6.14使用在習(xí)題6.13中生成的程序，實現(xiàn)式(6.28)給出的傳輸函數(shù)。你的結(jié)果與習(xí)題6.12中得到的結(jié)果相符嗎？使用函數(shù)latc2tf由向量k和alpha確定傳輸函數(shù)。所得到的傳輸函數(shù)和式(6.28)給出的傳輸函數(shù)相同嗎? 答：運行結(jié)果： k = 0.81093584641352 0.77112772506402 0.592*** 0.37169052478550 0.***293 0.02777777777778 alpha =-0.01112037033486 0.02345313662512-0.0***79-0.04739265773254 0.***485 0.20370370370370 0.11111111111111 與題6.12中得到的結(jié)果相符。

6.5：無限沖激響應(yīng)傳輸函數(shù)的并聯(lián)全通實現(xiàn)

Q6.15 生成下式給出的只階因果有界實低通1型切比雪夫傳輸函數(shù)G(z)的全通和的分解。使用zplane獲得G(z)的零極點分布圖：

G(z)全通和的分解：

G(z)的功率補充傳輸函數(shù)H(z)的表達式如下：

兩個全通傳輸函數(shù)的階數(shù)是1和2.Q6.15 生成一個五階因果有界實低通橢圓傳輸函數(shù)G(z)的全通和的分解。使用zplane獲得G(z)的零極點分布圖：

G(z)全通和的分解：

G(z)的功率補充傳輸函數(shù)H(z)的表達式如下：

兩個全通傳輸函數(shù)的階數(shù)是3和2.

第四篇：數(shù)字信號處理課程設(shè)計..

課程設(shè)計報告

課程名稱： 數(shù)字信號處理 課題名稱： 語音信號的處理與濾波

姓 名： 學(xué) 號： 院 系： 專業(yè)班級： 指導(dǎo)教師： 完成日期： 2013年7月2日

目錄

第1部分 課程設(shè)計報告………………………………………3 一．設(shè)計目的……………………………………………3 二．設(shè)計內(nèi)容……………………………………………3 三．設(shè)計原理……………………………………………3 四．具體實現(xiàn)……………………………………………5 1.錄制一段聲音…………………………………5 2.巴特沃斯濾波器的設(shè)計………………………8 3.將聲音信號送入濾波器濾波…………………13 4.語音信號的回放………………………………19 5.男女語音信號的頻譜分析……………………19 6.噪聲的疊加和濾除……………………………22 五． 結(jié)果分析……………………………………………27 第2部分 課程設(shè)計總結(jié)………………………………28 一． 參考文獻……………………………………………28

第1部分 課程設(shè)計報告

一．設(shè)計目的

綜合運用本課程的理論知識進行頻譜分析以及濾波器設(shè)計，通過理論推導(dǎo)得出相應(yīng)結(jié)論，并利用MATLAB作為工具進行實現(xiàn)，從而復(fù)習(xí)鞏固課堂所學(xué)的理論知識，提高對所學(xué)知識的綜合應(yīng)用能力，并從實踐上初步實現(xiàn)對數(shù)字信號的處理。

二．設(shè)計內(nèi)容

錄制一段個人自己的語音信號，并對錄制的信號進行采樣；畫出采樣后語音信號的時域波形和頻譜圖；給定濾波器的性能指標(biāo)，采用窗函數(shù)法和雙線性變換法設(shè)計濾波器，并畫出濾波器的頻率響應(yīng)；然后用自己設(shè)計的濾波器對采集的信號進行濾波，畫出濾波后信號的時域波形和頻譜，并對濾波前后的信號進行對比，分析信號的變化；回放語音信號；換一個與你性別相異的人錄制同樣一段語音內(nèi)容，分析兩段內(nèi)容相同的語音信號頻譜之間有什么特點；再錄制一段同樣長時間的背景噪聲疊加到你的語音信號中，分析疊加前后信號頻譜的變化，設(shè)計一個合適的濾波器，能夠把該噪聲濾除；

三．設(shè)計原理

1.在Matlab軟件平臺下，利用函數(shù)wavrecord()，wavwrite(),wavread(),wavplay()對語音信號進行錄制，存儲，讀取，回放。

2.用y=fft(x)對采集的信號做快速傅立葉變換，并用[h1,w]=freqz(h)進行DTFT變換。

3.掌握FIR DF線性相位的概念，即線性相位對h(n)、H(?)及零點的約束，了解四種FIR DF的頻響特點。

4.在Matlab中，F(xiàn)IR濾波器利用函數(shù)fftfilt對信號進行濾波。

5.抽樣定理

連續(xù)信號經(jīng)理想抽樣后時域、頻域發(fā)生的變化（理想抽樣信號與連續(xù)信號頻譜之間的關(guān)系）

理想抽樣信號能否代表原始信號、如何不失真地還原信號即由離散信號恢復(fù)連續(xù)信號的條件（抽樣定理）

理想采樣過程描述： 時域描述：

?a(t)?xa(t)?T(t)??xa(t)?(t?nT)??xa(nT)?(t?nT)xn???n??????T(t)?頻域描述：利用傅氏變換的性質(zhì)，時域相乘頻域卷積，若

n?????(t?nT)??a(t)Xa(j?)?xXa(j?)?xa(t)?T(j?)??T(t)

則有

?(j?)?1X(j?)??(j?)XaaT2?1?2?1??Xa(j?)??Xa(j??jk)??Xa(j??jk?s)Tk???TTk????(j?)與X(j?)的關(guān)系：理想抽樣信號的頻譜是連續(xù)信號頻譜的Xaa

周期延拓,重復(fù)周期為?s(采樣角頻率)。如果：

?X(j?)?Xa(j?)??a??0???s/2???s/2即連續(xù)信號是帶限的，且信號最高頻率不超過抽樣頻率的二分之一，則可不失真恢復(fù)。

奈奎斯特采樣定理：要使實信號采樣后能夠不失真還原，采樣頻率必須大于信號最高頻率的兩倍：?s?2?h 或 fs?2fh

四．具體實現(xiàn)

1.錄制一段聲音

1.1錄制并分析

在MATLAB中用wavrecord、wavread、wavplay、wavwrite對聲音進行錄制、讀取、回放、存儲。

程序如下：

Fs=8000;%抽樣頻率 time=3;%錄音時間 fprintf('按Enter鍵錄音%ds',time);%文字提示 pause;%暫停命令 fprintf('錄音中......');x=wavrecord(time*Fs,Fs,'double');%錄制語音信號 fprintf('錄音結(jié)束');%文字提示 fprintf('按Enter鍵回放錄音');pause;%暫停命令

wavplay(x,Fs);%按任意鍵播放語音信號

wavwrite(x,Fs,'C:UsersacerDesktop數(shù)字信號sound.wav');%存儲語音信號

N=length(x);%返回采樣點數(shù) df=fs/N;%采樣間隔 n1=1:N/2;f=[(n1-1)*(2*pi/N)]/pi;%頻帶寬度 figure(2);subplot(2,1,1);plot(x);%錄制信號的時域波形 title('原始信號的時域波形');%加標(biāo)題 ylabel('幅值/A');%顯示縱坐標(biāo)的表示意義 grid;%加網(wǎng)格

y0=fft(x);%快速傅立葉變換 figure(2);subplot(2,1,2);plot(f,abs(y0(n1)));%原始信號的頻譜圖 title('原始信號的頻譜圖');%加標(biāo)題 xlabel('頻率w/pi');%顯示橫坐標(biāo)表示的意義 ylabel('幅值 ');%顯示縱坐標(biāo)表示的意義 title('原始信號的頻譜圖');%加標(biāo)題

grid;%加網(wǎng)格

圖1.1 原始信號的時域與頻譜圖

1.2濾除無效點

針對實際發(fā)出聲音落后錄制動作半拍的現(xiàn)象，如何拔除對無效點的采樣的問題: 出現(xiàn)這種現(xiàn)象的原因主要是錄音開始時，人的反應(yīng)慢了半拍，導(dǎo)致出現(xiàn)了一些無效點，而后而出現(xiàn)的無效的點，主要是已經(jīng)沒有聲音的動作，先讀取聲音出來，將原始語音信號時域波形圖畫出來，根據(jù)己得到的信號，可以在第二次讀取聲音的后面設(shè)定采樣點，取好有效點，畫出濾除無效點后的語音信號時域波形圖，對比可以看出。這樣就可以解決這個問題。

x=wavread('C:UsersacerDesktop數(shù)字信號sound.wav', 7

[4000,24000]);%從4000點截取到24000結(jié)束 plot(x);%畫出截取后的時域圖形 title('截取后的聲音時域圖形');%標(biāo)題 xlabel('頻率');ylabel('振幅');grid;%畫網(wǎng)格

圖1.2 去除無效點

2.巴特沃斯濾波器的設(shè)計

2.1設(shè)計巴特沃思低通濾波器

MATLAB程序如下。濾波器圖如圖3.3所示。

%低通濾波

fp=1000;fs=1200;Fs=22050;rp=1;rs=100;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;Fs1=1;wap=2*tan(wp/2);was=2*tan(ws/2);[N,wc]=buttord(wap,was,rp,rs,'s');[B,A]=butter(N,wc,'s');[Bz,Az]=bilinear(B,A,Fs1);figure(1);[h,w]=freqz(Bz,Az,512,Fs1*22050);plot(w,abs(h));title('巴特沃斯低通濾波器');xlabel('頻率（HZ）');ylabel('耗損（dB）');gridon;9

圖2.1 巴特沃思低通濾波器

2.2設(shè)計巴特沃思高通濾波器

MATLAB程序如下。濾波器圖如圖3.5所示。%高通濾波

fp=4800;fs=5000;Fs=22050;rp=1;rs=100;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;T=1;Fs1=1;wap=2*tan(wp/2);was=2*tan(ws/2);10

[N,wc]=buttord(wap,was,rp,rs,'s');[B,A]=butter(N,wc,'high','s');[Bz,Az]=bilinear(B,A,Fs1);figure(1);[h,w]=freqz(Bz,Az,512,Fs1*22050);plot(w,abs(h));title('巴特沃斯高通濾波器');xlabel('頻率（HZ）');ylabel('耗損（dB）');grid on;

圖2.2巴特沃思高通濾波器

2.3設(shè)計巴特沃思帶通濾波器

MATLAB程序如下。濾波器圖如圖3.7所示。%帶通濾波

fp=[1200,3000];fs=[1000,3200];Fs=8000;rp=1;rs=100;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;T=1;Fs1=1;wap=2*tan(wp/2);was=2*tan(ws/2);[N,wc]=buttord(wap,was,rp,rs,'s');[B,A]=butter(N,wc,'s');[Bz,Az]=bilinear(B,A,Fs1);figure(4);[h,w]=freqz(Bz,Az,512,Fs1*1000);plot(w,abs(h));title('巴特沃斯帶通濾波器');xlabel('頻率（HZ）');ylabel('耗損（dB）');grid on;12

圖2.3巴特沃思帶通濾波器

3．將聲音信號送入濾波器濾波

x=wavread('C:UsersacerDesktop數(shù)字信號sound.wav');%播放原始信號

wavplay(x,fs);%播放原始信號 N=length(x);%返回采樣點數(shù) df=fs/N;%采樣間隔 n1=1:N/2;f=[(n1-1)*(2*pi/N)]/pi;%頻帶寬度 figure(4);subplot(4,2,1);plot(x);%錄制信號的時域波形

title('原始信號的時域波形');%加標(biāo)題 ylabel('幅值/A');%顯示縱坐標(biāo)的表示意義 grid;%加網(wǎng)格

y0=fft(x);%快速傅立葉變換 subplot(4,2,3);plot(f,abs(y0(n1)));%原始信號的頻譜圖 title('原始信號的頻譜圖');%加標(biāo)題 xlabel('頻率w/pi');%顯示橫坐標(biāo)表示的意義 ylabel('幅值 ');%顯示縱坐標(biāo)表示的意義 title('原始信號的頻譜圖');%加標(biāo)題 grid;%加網(wǎng)格

3．1低通濾波器濾波 fs=8000;beta=10.056;wc=2*pi*1000/fs;ws=2*pi*1200/fs;width=ws-wc;wn=(ws+wc)/2;n=ceil(12.8*pi /width);h=fir1(n,wn/pi,'band',kaiser(n+1,beta));[h1,w]=freqz(h);

ys=fftfilt(h,x);%信號送入濾波器濾波，ys為輸出 fftwave=fft(ys);%將濾波后的語音信號進行快速傅立葉變換 figure(4);subplot(4,2,2);%在四行兩列的第二個窗口顯示圖形 plot(ys);%信號的時域波形

title('低通濾波后信號的時域波形');%加標(biāo)題 xlabel('頻率w/pi');ylabel('幅值/A');%顯示標(biāo)表示的意義 grid;%網(wǎng)格

subplot(4,2,4);%在四行兩列的第四個窗口顯示圖形 plot(f, abs(fftwave(n1)));%繪制模值 xlabel('頻率w/pi');ylabel('幅值/A');%顯示標(biāo)表示的意義

title('低通濾波器濾波后信號的頻譜圖');%標(biāo)題 grid;%加網(wǎng)格

wavplay(ys,8000);%播放濾波后信號

3．2高通濾波器濾波 fs=8000;beta=10.056;ws=2*5000/fs;wc=2*4800/fs;

width=ws-wc;wn=(ws+wc)/2;n=ceil(12.8*pi/width);h=fir1(n,wn/pi, 'high',kaiser(n+2,beta));[h1,w]=freqz(h);ys=fftfilt(h,x);%將信號送入高通濾波器濾波 subplot(4,2,5);%在四行兩列的第五個窗口顯示圖形 plot(ys);%信號的時域波形 xlabel('頻率w/pi');ylabel('幅值/A');%顯示標(biāo)表示的意義 title('高通濾波后信號的時域波形');%標(biāo)題 ylabel('幅值/A');%顯示縱坐標(biāo)的表示意義 grid;%網(wǎng)格

fftwave=fft(ys);%將濾波后的語音信號進行快速傅立葉變換 subplot(4,2,7);%在四行兩列的第七個窗口顯示圖形 plot(f,abs(fftwave(n1)));%繪制模值 axis([0 1 0 50]);xlabel('頻率w/pi');ylabel('幅值/A');%顯示標(biāo)表示的意義

title('高通濾波器濾波后信號的頻譜圖');%標(biāo)題 grid;%加網(wǎng)格

wavplay(ys,8000);%播放濾波后信號

3．3帶通濾波器 fs=8000;beta=10.056;wc1=2*pi*1000/fs;wc2=2*pi*3200/fs;ws1=2*pi*1200/fs;ws2=2*pi*3000/fs;width=ws1-wc1;wn1=(ws1+wc1)/2;wn2=(ws2+wc2)/2;wn=[wn1 wn2];n=ceil(12.8/width*pi);h=fir1(n,wn/pi,'band',kaiser(n+1,beta));[h1,w]=freqz(h);ys1= fftfilt(h,x);%將信號送入高通濾波器濾波 figure(4);subplot(4,2,6);%在四行兩列的第六個窗口顯示圖形 plot(ys1);%繪制后信號的時域的圖形 title('帶通濾波后信號的時域波形');%加標(biāo)題 xlabel('頻率w/pi');ylabel('幅值/A');%顯示縱坐標(biāo)表示的意義 grid;%網(wǎng)格

fftwave=fft(ys1);%對濾波后的信號進行快速傅立葉變換 subplot(4,2,8);%在四行兩列的第八個窗口顯示圖形

plot(f, abs(fftwave(n1)));%繪制模值 axis([0 1 0 50]);xlabel('頻率w/pi');ylabel('幅值/A');%顯示標(biāo)表示的意義 title('帶通濾波器濾波后信號的頻譜圖');%加標(biāo)題 grid;%網(wǎng)格

wavplay(ys1,8000);%播放濾波后信號 圖形如下：

原始信號的時域波形幅值/A0-1012x 10原始信號的頻譜圖34幅值/A1低通濾波后信號的時域波形0.50-0.5012頻率w/pi3400.51頻率w/pi高通濾波后信號的時域波形幅值/A0幅值/A0幅值/Ax 10高通濾波器濾波后信號的頻譜圖5012頻率w/pi34幅值/A0.20-0.2幅值/A2001000x 10低通濾波器濾波后信號的頻譜圖200100000.51頻率w/pi帶通濾波后信號的時域波形0.50-0.501234頻率w/pix 10帶通濾波器濾波后信號的頻譜圖50幅值 00.5頻率w/pi1000.5頻率w/pi1

分析：三個濾波器濾波后的聲音與原來的聲音都發(fā)生了變化。其中低

通的濾波后與原來聲音沒有很大的變化，其它兩個都又明顯的變化

4.語音信號的回放

sound(xlow,Fs,bits);%在Matlab中，函數(shù)sound可以對聲音進行回放，其調(diào)用格式: sound(xhigh, Fs,bits);%sound(x, Fs, bits);sound(xdaitong, Fs,bits);5.男女語音信號的頻譜分析

5.1 錄制一段異性的聲音進行頻譜分析

Fs=8000;%抽樣頻率 time=3;%錄音時間 fprintf('按Enter鍵錄音%ds',time);%文字提示 pause;%暫停命令 fprintf('錄音中......');x=wavrecord(time*Fs,Fs,'double');%錄制語音信號 fprintf('錄音結(jié)束');%文字提示 fprintf('按Enter鍵回放錄音');pause;%暫停命令 wavplay(x,Fs);%按任意鍵播放語音信號

wavwrite(x,Fs,'C:UsersacerDesktop數(shù)字信號sound2.wav');%存儲語音信號

5.2 分析男女聲音的頻譜

x=wavread(' C:UsersacerDesktop數(shù)字信號sound2.wav ');%播放原始信號,解決落后半拍

wavplay(x,fs);%播放原始信號 N=length(x);%返回采樣點數(shù) df=fs/N;%采樣間隔 n1=1:N/2;

f=[(n1-1)*(2*pi/N)]/pi;%頻帶寬度 figure(1);subplot(2,2,1);plot(x);%錄制信號的時域波形

title('原始女生信號的時域波形');%加標(biāo)題 ylabel('幅值/A');%顯示縱坐標(biāo)的表示意義 grid;%加網(wǎng)格

y0=fft(x);%快速傅立葉變換 subplot(2,2,2);plot(f,abs(y0(n1)));%原始信號的頻譜圖 title('原始女生信號的頻譜圖');%加標(biāo)題 xlabel('頻率w/pi');%顯示橫坐標(biāo)表示的意義 ylabel('幅值 ');%顯示縱坐標(biāo)表示的意義 grid;%加網(wǎng)格

[y,fs,bits]=wavread(' C:UsersacerDesktop數(shù)字信號sound.wav ');% 對語音信號進行采樣

wavplay(y,fs);%播放原始信號 N=length(y);%返回采樣點數(shù) df=fs/N;%采樣間隔 n1=1:N/2;f=[(n1-1)*(2*pi/N)]/pi;%頻帶寬度 subplot(2,2,3);plot(y);%錄制信號的時域波形

title('原始男生信號的時域波形');%加標(biāo)題 ylabel('幅值/A');%顯示縱坐標(biāo)的表示意義 grid;%加網(wǎng)格

y0=fft(y);%快速傅立葉變換

subplot(2,2,4);%在四行兩列的第三個窗口顯示圖形 plot(f,abs(y0(n1)));%原始信號的頻譜圖 title('原始男生信號的頻譜圖');%加標(biāo)題 xlabel('頻率w/pi');%顯示橫坐標(biāo)表示的意義 ylabel('幅值 ');%顯示縱坐標(biāo)表示的意義 grid;%加網(wǎng)格

5.3男女聲音的頻譜圖

原始女生信號的時域波形0.50-0.5-1150100原始女生信號的頻譜圖幅值/A幅值 012345000x 10原始男生信號的時域波形0.50.5頻率w/pi原始男生信號的頻譜圖1300200幅值/A0幅值 012x 1034100-0.5000.5頻率w/pi1

圖5.3男女聲音信號波形與頻譜對比

分析：就時域圖看，男生的時域圖中振幅比女生的高，對于頻譜圖女生的高頻成分比較多

6.噪聲的疊加和濾除

6.1錄制一段背景噪聲

Fs=8000;%抽樣頻率 time=3;%錄音時間 fprintf('按Enter鍵錄音%ds',time);%文字提示 pause;%暫停命令 fprintf('錄音中......');x=wavrecord(time*Fs,Fs,'double');%錄制語音信號

fprintf('錄音結(jié)束');%文字提示 fprintf('按Enter鍵回放錄音');pause;%暫停命令 wavplay(x,Fs);%按任意鍵播放語音信號 wavwrite(x,Fs,'C:UsersacerDesktop數(shù)字信號噪音.wav');%存儲語音信號

6.2 對噪聲進行頻譜的分析

[x1,fs,bits]=wavread(' C:UsersacerDesktop數(shù)字信號噪音.wav ');%對語音信號進行采樣

wavplay(x1,fs);%播放噪聲信號 N=length(x1);%返回采樣點數(shù) df=fs/N;%采樣間隔

n1=1:N/2;f=[(n1-1)*(2*pi/N)]/pi;%頻帶寬度 figure(5);subplot(3,2,1);plot(x1);%信號的時域波形 title('噪聲信號的時域波形');grid;ylabel('幅值/A');y0=fft(x1);%快速傅立葉變換

subplot(3,2,2);plot(f,abs(y0(n1)));%噪聲信號的頻譜圖 ylabel('幅值');title('噪聲信號的頻譜圖');

6.3原始信號與噪音的疊加

fs=8000;[x,fs,bits]=wavread(' C:UsersacerDesktop數(shù)字信號sound.wav ');%對錄入信號進行采樣

[x1,fs,bits]=wavread(' C:UsersacerDesktop數(shù)字信號噪音.wav ');%對噪聲信號進行采樣

yy=x+x1;%將兩個聲音疊加

6.4疊加信號的頻譜分析：

wavplay(yy,fs);%播放疊加后信號 N=length(yy);%返回采樣點數(shù) df=fs/N;%采樣間隔 n1=1:N/2;f=[(n1-1)*(2*pi/N)]/pi;%頻帶寬度 figure(5);subplot(3,2,3);plot(yy,'LineWidth',2);%信號的時域波形

title('疊加信號的時域波形');xlabel('時間/t');ylabel('幅值/A');grid;y0=fft(yy);%快速傅立葉變換 subplot(3,2,4);plot(f,abs(y0(n1)));%疊加信號的頻譜圖 title('疊加信號的頻譜圖');xlabel('頻率w/pi');ylabel('幅值/db');grid;

6.5 設(shè)計一個合適的濾波器將噪聲濾除 fs=18000;%采樣頻率 Wp=2*1000/fs;%通帶截至頻率 Ws=2*2000/fs;%阻帶截至頻率 Rp=1;%最大衰減 Rs=100;%最小衰減

[N,Wn]=buttord(Wp,Ws,Rp,Rs);%buttord函數(shù)（n為階數(shù)，Wn為截至頻率）

[num,den]=butter(N,Wn);%butter函數(shù)（num為分子系數(shù)den為分母系數(shù)）

[h,w]=freqz(num,den);%DTFT變換

ys=filter(num,den,yy);%信號送入濾波器濾波，ys為輸出 fftwave=fft(ys);%將濾波后的語音信號進行快速傅立葉變換 figure(5);subplot(3,2,5);plot(ys);%信號的時域波形

title('低通濾波后信號的時域波形');%加標(biāo)題 ylabel('幅值/A');%顯示標(biāo)表示的意義 grid;%網(wǎng)格 subplot(3,2,6);plot(f, abs(fftwave(n1)));%繪制模值 title('低通濾波器濾波后信號的頻譜圖');%標(biāo)題 xlabel('頻率w/pi');ylabel('幅值/A');%顯示標(biāo)表示的意義 grid;%加網(wǎng)格

wavplay(ys,8000);%播放濾波后信號 grid;圖形如下：

噪聲信號的時域波形1100噪聲信號的頻譜圖幅值/A0-1幅值0123450000.5疊加信號的頻譜圖1x 10疊加信號的時域波形10-101時間/t2200幅值/db34幅值/A100000.5頻率w/pi1x 10低通濾波后信號的時域波形0.5低通濾波器濾波后信號的頻譜圖200幅值/A0-0.5幅值/A012x 1034100000.5頻率w/pi1

圖6.1噪音的疊加與濾除前后頻譜對比

7.結(jié)果分析

1.錄制剛開始時，常會出現(xiàn)實際發(fā)出聲音落后錄制動作半拍，可在[x,fs,bits]=wavread('d:matlavworkwomamaaiwo.wav')加 窗[x,fs,bits]=wavread('d:matlavworkwomamaaiwo.wav'，[100 10000])，窗的長度可根據(jù)需要定義。

2.語音信號通過低通濾波器后，把高頻濾除，聲音變得比較低沉。當(dāng)通過高通濾波器后，把低頻濾除，聲音變得比較就尖銳。通過帶通濾波器后，聲音比較適中。

3.通過觀察男生和女生圖像知：時域圖的振幅大小與性別無關(guān)，只與說話人音量大小有關(guān)，音量越大，振幅越大。頻率圖中，女生高 27

頻成分較多。

4.疊加噪聲后，噪聲與原信號明顯區(qū)分,但通過低通濾波器后，噪聲沒有濾除，信號產(chǎn)生失真。原因可能為噪聲與信號頻率相近無法濾除。

第2部分 課程設(shè)計總結(jié)

通過本次課程設(shè)計，使我們對數(shù)字信號處理相關(guān)知識有了更深刻的理解，尤其是對各種濾波器的設(shè)計。在設(shè)計的過程中遇到了很多問題，剛剛開始時曾天真的認為只要把以前的程序改了參數(shù)就可以用了，可是問題沒有我想象中的那么簡單，單純的搬程序是不能解決問題的。通過查閱資料和請教同學(xué)收獲了很多以前不懂的理論知識。再利用所學(xué)的操作，發(fā)現(xiàn)所寫的程序還是沒有能夠運行，通過不斷地調(diào)試，運行，最終得出了需要的結(jié)果。整個過程中學(xué)到了很多新的知識，特別是對Matlab的使用終于有些了解。在以后的學(xué)習(xí)中還需要深入了解這方面的內(nèi)容。在這次的課程設(shè)計中讓我體會最深的是：知識來不得半點的馬虎。也認識到自己的不足，以后要進一步學(xué)習(xí)。

八．參考文獻

[1]數(shù)字信號處理教程(第三版)程佩青 清華大學(xué)出版社 [2]MATLAB信號處理 劉波 文忠 電子工業(yè)出版社 [3]MATLAB7.1及其在信號處理中的應(yīng)用 王宏 清華大學(xué)出版社

[4]MATLAB基礎(chǔ)與編程入門 張威 西安電子科技大學(xué)出版社

[5] 數(shù)字信號處理及其MATLAB實驗 趙紅怡 張常 化學(xué)工業(yè)出版社

[6]MATLAB信號處理詳解 陳亞勇等 人民郵電出版社 [7] 數(shù)字信號處理

錢同惠 機械工業(yè)出版社 29

第五篇：數(shù)字信號處理實驗報告

JIANGSU

UNIVERSITY OF TECHNOLOGY

數(shù)字信號處理實驗報告

學(xué)院名稱： 電氣信息工程學(xué)院

專 業(yè)：

班 級： 姓 名： 學(xué) 號： 指導(dǎo)老師： 張維璽（教授）

2013年12月20日

實驗一 離散時間信號的產(chǎn)生

一、實驗?zāi)康?/p>

數(shù)字信號處理系統(tǒng)中的信號都是以離散時間形態(tài)存在的，所以對離散時間信號的研究是數(shù)字信號的基本所在。而要研究離散時間信號，首先需要產(chǎn)生出各種離散時間信號。使用MATLAB軟件可以很方便地產(chǎn)生各種常見的離散時間信號，而且它還具有強大繪圖功能，便于用戶直觀地處理輸出結(jié)果。

通過本實驗，學(xué)生將學(xué)習(xí)如何用MATLAB產(chǎn)生一些常見的離散時間信號，實現(xiàn)信號的卷積運算，并通過MATLAB中的繪圖工具對產(chǎn)生的信號進行觀察，加深對常用離散信號和信號卷積和運算的理解。

二、實驗原理

離散時間信號是指在離散時刻才有定義的信號，簡稱離散信號，或者序列。離散序列通常用x（n）來表示，自變量必須是整數(shù)。常見的離散信號如下：（1）單位沖激序列δ（n）

如果δ（n）在時間軸上延遲了k個單位，得到δ（n-k），即長度為N的單位沖激序列δ（n）可以通過下面的MATLAB命令獲得。

n=-（N-1）：N-1 x=[zeros(1,N-1)1 zeros(1,N-1)]； stem(n,x)延遲K個采樣點的長度為N的單位沖激序列δ（n-k）（k

n=0：N-1 y=[zeros(1,M)1 zeros(1,N-M-1)]； stem（n，y）

（2）單位階躍序列u（n）

如果u（n）在時間軸上延遲了k個單位，得到u（n-k），即長度為N的單位階躍序列u（n）可以通過下面的MATLAB命令獲得。

n=-（N-1）：N-1 x=[zeros(1,N-1)ones(1,N)]； stem(n,x)延遲的單位階躍序列可以使用類似于單位沖激序列的方法獲得。（3）矩形序列

矩形序列有一個重要的參數(shù)，就是序列的寬度N。矩形序列與u（n）之間的關(guān)系為矩形序列等= u（n）— u（n-N）。

因此，用MATLAB表示矩形序列可利用上面的單位階躍序列組合而成。（4）正弦序列x（n）

這里，正弦序列的參數(shù)都是實數(shù)。與連續(xù)的正弦信號不同，正弦序列的自變量n必須為整數(shù)。可以證明，只有當(dāng)2π/w為有理數(shù)時，正弦序列具有周期性。

長度為N的正弦序列x（n）可以通過下面的MATLAB命令獲得。n=0：N-1 x=A*cos（2*pi*f*n/Fs+phase）（5）單邊實指數(shù)序列x（n）

長度為N的實指數(shù)序列x（n）可以通過下面的MATLAB命令實現(xiàn)。n=0：N-1 x=a.^n stem(n,x)單邊指數(shù)序列n的取值范圍為n>=0。當(dāng)|a|>1時，單邊指數(shù)序列發(fā)散；當(dāng)|a|<1時，單邊指數(shù)序列收斂。當(dāng)a>0時，該序列均取正值；當(dāng)a<0時，序列在正負擺動。

（6）負指數(shù)序列x（n）

當(dāng)a=0時，得到虛指數(shù)序列x（n）。

與連續(xù)負指數(shù)信號一樣，我們將負指數(shù)序列實部和虛部的波形分開討論，得到如下結(jié)論：

1）當(dāng)a>0時，負指數(shù)序列x（n）的實部和虛部分別是按指數(shù)規(guī)律增長的正弦振蕩序列；

2)當(dāng)a<0時，負指數(shù)序列x（n）的實部和虛部分別是按指數(shù)規(guī)律衰減的正弦振蕩序列；

3)當(dāng)a=0時，負指數(shù)序列x（n）即為虛指數(shù)序列，其實部和虛部分別是等幅的正弦振蕩序列；

長度為N的實指數(shù)序列x（n）可以通過下面的MATLAB命令實現(xiàn)。n=0：N-1 x=exp（（a.+j*w）*n）stem(n,real(x))或

stem(n,imag(x))

三、實驗內(nèi)容及分析

?1n?01、編制程序產(chǎn)生單位沖激序列??n???“?并繪出其圖及??n?”學(xué)號后兩位0n?0?形。程序:(1)N=4;

n=-(N-1):N-1;

x=[zeros(1,N-1)1 zeros(1,N-1)];stem(n,x);

title('單位沖激序列');

grid on;

(2)N=6;

M=1;%學(xué)號01 n=-(N-1):N-1;

y=[zeros(1,N-M+1)1 zeros(1,N-M-1)];stem(n,y);

title('單位沖激序列');grid on;

分析:在上圖的基礎(chǔ)上向右平移了1個單位。

?1n?02、編制程序產(chǎn)生單位階躍序列u?n???、u?n?“學(xué)號后兩位”?及

0n?0?u?n??u?n?“學(xué)號后兩位”?，并繪出其圖形。程序: 4

(1)N=5;

n=-(N-1):N-1;

x=[zeros(1,N-1)ones(1,N)];stem(n,x);

title('單位階躍序列');grid on;

(2)N=6;

M=1;%學(xué)號01 n=-(N-1):N-1;

x=[zeros(1,N-M+1)ones(1,N-M)];stem(n,x);

title('單位階躍序列');grid on;

分析:在上圖的基礎(chǔ)上平移了1個單位.(3)N=6;

M=1;%學(xué)號01 n=-(N-1):N-1;

x=[zeros(1,N-1)ones(1,N)];y=[zeros(1,N-M+1)ones(1,N-M)];z=x-y;stem(n,z);

title('單位階躍序列');grid on;

2??

3、編制程序產(chǎn)生正弦序列x?n??cos?2?n?、x?n??cos??n?及

?學(xué)號后兩位?x?n??sin?2n?并繪出其圖形。

程序:(1)N=5;

A=1;

w=2*pi;phi=0;n=0:0.05:N-1;x=A*cos(w*n+phi);stem(n,x);title('余弦信號');grid on;

分析:該序列具有周期性,且輸出為余弦信號.(2)N=5;

A=1;

w=2*pi/1;%學(xué)號01 phi=0;n=0:0.05:N-1;x=A*cos(w*n+phi);stem(n,x);title('余弦信號');grid on;

;

分析:該序列具有周期性,且輸出為余弦信號.(3)N=5;

A=1;

w=2*pi;phi=0;

n=0:0.05:N-1;x=A*sin(w*n+phi);stem(n,x);title('正弦信號');grid on;

分析:該序列具有周期性,且輸出為正弦信號.4、編制程序產(chǎn)生復(fù)正弦序列x?n??e(2?j學(xué)號后兩位）n，并繪出其圖形。N=3;

n=0:0.2:N-1;

w=1;%學(xué)號01 x=exp((2+j*w)*n);subplot(2,1,1)

stem(n,real(x)),title('實部');grid on;subplot(2,1,2)

stem(n,imag(x)),title('虛部');grid on;

5、編制程序產(chǎn)生指數(shù)序列x?n??an，并繪出其圖形。其中a=學(xué)號后兩位、a=1/“學(xué)號后兩位”。

(1)N=10;

n=0:N-1;

a=1;%學(xué)號01 x=a.^n;stem(n,x);title('指數(shù)序列');grid on;

(2)N=10;

n=0:N-1;

a=1;%學(xué)號01 x=a.^(-n);stem(n,x);title('指數(shù)序列');grid on;

實驗三 離散時間信號的頻域分析

一、實驗?zāi)康?/p>

信號的頻域分析是信號處理中一種有效的工具。在離散信號的頻域分析中，通常將信號表示成單位采樣序列的線性組合，而在頻域中，將信號表示成復(fù)變量或的線性組合。通過這樣的表示，可以將時域的離散序列映射到頻域以便于進一步的處理。

在本實驗中，將學(xué)習(xí)利用MATLAB計算離散時間信號的DTFT和DFT，并加深對其相互關(guān)系的理解。

二、實驗原理

（１）DTFT和DFT的定義及其相互關(guān)系。

（２）使用到的MATLAB命令有基于DTFT離散時間信號分析函數(shù)以及求解序列的DFT函數(shù)。

三、實驗內(nèi)容及分析

（１）編程計算并畫出下面DTFT的實部、虛部、幅度和相位譜。

X(e)?jw0.0518?0.1553e1?1.2828ex(n)?cos?jw?jw?0.1553e?j2w?1.0388e?j2w?0.0518e?j3w?0.3418e?j3w

（２）計算32點序列

5?n16，0≦n≦31的32點和64點DFT，分別繪出幅度譜圖形，并繪出該序列的DTFT圖形。

３－１

clear;

x=[0.0518,-0.1553,0.1553,0.0518];y=[1,1.2828,1.0388,0.3418];w=[0:500]*pi/500 H=freqz(x,y,w);

magX=abs(H);angX=angle(H);realX=real(H);imagX=imag(H);subplot(221);plot(w/pi,magX);grid;

xlabel('frequency in pi unit');ylabel('magnitude');title('幅度 part');axis([0 0.9 0 1.1]);

subplot(223);plot(w/pi,angX);grid;

xlabel('frequency in pi unit');ylabel('radians');title('相位 part');axis([0 1-3.2 3.2]);

subplot(222);plot(w/pi,realX);grid;

xlabel('frequency in pi unit');ylabel('real part');title('實部 part');axis([0 1-1 1]);

subplot(224);plot(w/pi,imagX);grid;

xlabel('frequency in pi unit');ylabel('imaginary');title('虛部 part');axis([0 1-1 1.1]);

３－２

N=32;n=0:N-1;

xn=cos(5*pi*n/16);k=0:1:N-1;Xk=fft(xn,N);subplot(2,1,1);stem(n,xn);subplot(2,1,2);stem(k,abs(Xk));title('32點');figure N=64;n=0:N-1;

xn=cos(5*pi*n/16);k=0:1:N-1;Xk=fft(xn,N);subplot(2,1,1);stem(n,xn);subplot(2,1,2);stem(k,abs(Xk));title('64點');

（１）

（２）

實驗四 離散時間LTI系統(tǒng)的Z域分析

一、實驗?zāi)康?/p>

本實驗通過使用MATLAB函數(shù)對離散時間系統(tǒng)的一些特性進行仿真分析，以加深對離散時間系統(tǒng)的零極點、穩(wěn)定性，頻率響應(yīng)等概念的理解。學(xué)會運用MATLAB分析離散時間系統(tǒng)的系統(tǒng)函數(shù)的零極點；學(xué)會運用MATLAB分析系統(tǒng)函數(shù)的零極點分布與其時域特性的關(guān)系；學(xué)會運用MATLAB進行離散時間系統(tǒng)的頻率特性分析。

二、實驗原理

離散時間系統(tǒng)的系統(tǒng)函數(shù)定義為系統(tǒng)零狀態(tài)響應(yīng)的Z變化與激勵的Z變化之比。

在MATLAB中系統(tǒng)函數(shù)的零極點可通過函數(shù)roots得到，也可借助函數(shù)tf2zp得到，tf2zp的語句格式為

[Z,P,K]=tf2zp(B,A)其中，B與A分別表示H（z）的分子與分母多項式的系數(shù)向量。它的作用是將H（z）的有理分式表示式轉(zhuǎn)換為零極點增益形式。

若要獲得系統(tǒng)函數(shù)H（z）的零極點分布圖，可直接應(yīng)用zplane函數(shù)，其語句格式為

Zplane（B，A）

其中，B與A分別表示H（z）的分子和分母多項式的系數(shù)向量。它的作用是在z平面上畫出單位圓、零點與極點。

離散系統(tǒng)中z變化建立了時域函數(shù)h（n）與z域函數(shù)H（z）之間的對應(yīng)關(guān)系。因此，z變化的函數(shù)H（z）從形式可以反映h（n）的部分內(nèi)在性質(zhì)。可根據(jù)系統(tǒng)的傳遞函數(shù)H（z）求單位沖激響應(yīng)h（n）的函數(shù)impz、filter等。

利用系統(tǒng)的頻率響應(yīng)，可以分析系統(tǒng)對各種頻率成分的響應(yīng)特性，并推出系統(tǒng)的特性（高通、低通、帶通、帶阻等）。

MATLAB提供了求離散時間系統(tǒng)頻響特性的函數(shù)freqz，調(diào)用freqz的格式主要有兩種。一種形式為

[H,w]= reqz(B,A,N)其中，B與A分別表示H（z）分子和分母多項式的系數(shù)向量；N為正整數(shù)，默認值為512；返回值w包含[0，π]范圍內(nèi)的N個頻率等分點；返回值H則是離散時間系統(tǒng)頻率響應(yīng)在0~π范圍內(nèi)N個頻率處的值。另一種形式為

[H，w]= freqz（B，A，N，‘whole’）

與第一種方式不同之處在于角頻率的范圍由[0，π]擴展到[0，2π]。

三、實驗內(nèi)容與結(jié)果分析

已知LTI離散時間系統(tǒng)，要求由鍵盤實現(xiàn)系統(tǒng)參數(shù)輸入，并繪出幅頻和相頻響應(yīng)曲線和零極點分布圖，進而分析系統(tǒng)的濾波特性和穩(wěn)定性。

（一）程序

b=[0.0528,0.797,0.1295,0.1295,0.797,0.0528];

a=[1,-1.8107,2.4947,-1.8801,0.9537,-0.2336];w=[0:20:500]*pi/500;

x1=0.0528+0.797*exp(-1*j*w)+0.1295*exp(-2*j*w)+0.1295*exp(-3*j*w)+0.797*exp(-4*j*w)+0.0528*exp(-5*j*w);

x2=1-1.8107*exp(-1*j*w)+2.4947*exp(-2*j*w)+1.8801*exp(-3*j*w)+0.9537*exp(-4*j*w)+0.2336*exp(-5*j*w);x22=x2+(x2==0)*eps;x=x1./x22;magx=abs(x);

angx=angle(x).*180/pi;

subplot(2,2,3);zplane(b,a);title('零極點圖');subplot(2,2,2);stem(w/pi,magx);title('幅度部分');ylabel('振幅');subplot(2,2,4);stem(w/pi,angx);

xlabel('以pi為單位的頻率');title('相位部分');ylabel('相位');

（二）波形圖

圖4-1 幅頻、相頻響應(yīng)曲線、零極點分布圖

實驗六 IIR數(shù)字濾波器的設(shè)計

一、實驗?zāi)康?/p>

從理論上講，任何的線性是不變（LTI）離散時間系統(tǒng)都可以看做一個數(shù)字濾波器，因此設(shè)計數(shù)字濾波器實際就是設(shè)計離散時間系統(tǒng)。數(shù)字濾波器你包括IIR（無限沖激響應(yīng)）和FIR（有限沖激響應(yīng)）型，在設(shè)計時通常采用不同的方法。

本實驗通過使用MATLAB函數(shù)對數(shù)字濾波器進行設(shè)計和和實現(xiàn)，要求掌握IIR數(shù)字巴特沃斯濾波器、數(shù)字切比雪夫濾波器的設(shè)計原理、設(shè)計方法和設(shè)計步驟；能根據(jù)給定的濾波器指標(biāo)進行濾波器設(shè)計；同時也加深學(xué)生對數(shù)字濾波器的常用指標(biāo)和設(shè)計過程的理解。

二、實驗原理

在IIR濾波器的設(shè)計中，常用的方法是：先根據(jù)設(shè)計要求尋找一個合適的模擬原型濾波器，然后根據(jù)一定的準(zhǔn)則將此模擬原型濾波器轉(zhuǎn)換為數(shù)字濾波器。

IIR濾波器的階數(shù)就等于所選的模擬原型濾波器的階數(shù)，所以其階數(shù)確定主要是在模擬原型濾波器中進行的。

IIR數(shù)字濾波器的設(shè)計方法如下：（1）沖激響應(yīng)不變法。（2）雙線性變化法。

一般來說，在要求時域沖激響應(yīng)能模仿模擬濾波器的場合，一般使用沖激響應(yīng)不變法。沖激響應(yīng)不變法一個重要特點是頻率坐標(biāo)的變化是線性的，因此如果模擬濾波器的頻率響應(yīng)帶限于折疊頻率的話，則通過變換后濾波器的頻率響應(yīng)可不失真地反映原響應(yīng)與頻率的關(guān)系。

與沖激響應(yīng)不變法比較，雙線性變化的主要優(yōu)點是靠頻率的非線性關(guān)系得到s平面與z平面的單值一一對應(yīng)關(guān)系，整個值對應(yīng)于單位圓一周。所以從模擬傳遞函數(shù)可直接通過代數(shù)置換得到數(shù)字濾波器的傳遞函數(shù)。

MATLAB提供了一組標(biāo)準(zhǔn)的數(shù)字濾波器設(shè)計函數(shù)，大大簡化了濾波器的設(shè)計工程。

（1）butter。

（2）cheby1、cheby2。

三、實驗內(nèi)容及分析

利用MATLAB編程方法或利用MATLAB中fdatool工具設(shè)計不同功能的IIR數(shù)字濾波器。

1、基于chebyshev I型模擬濾波器原型使用沖激不變轉(zhuǎn)換方法設(shè)計數(shù)字濾波器，要求參數(shù)為通帶截止頻率?p?0.4?；通帶最大衰減Ap?1dB；阻帶截止頻率?s?0.4?；阻帶最小衰減As?35dB。

程序:

wp=0.2*pi;

%通帶邊界頻率

ws=0.4*pi;

%阻帶截止頻率 rp=1;

%通帶最大衰減 rs=35;

%阻帶最小衰減

Fs=1000;

%?ùéè3é?ù??3?1000hz

[N,Wn]=cheb1ord(wp,ws,rp,rs,'s');

[Z,P,K]=cheby1(N,rp,Wn,'s');[H,W]=zp2tf(Z,P,K);

figure(1);freqs(H,W);[P,Q]=freqs(H,W);figure(2);plot(Q*Fs/(2*pi),abs(P));grid on;

xlabel('頻率/Hz');ylabel('幅度');

2、基于Butterworth型模擬濾波器原型使用雙線性變換方法設(shè)計數(shù)字濾波器的，要求參數(shù)為截止頻率?p?0.4?；通帶最大衰減Ap?1dB；阻帶截止頻率?s?0.25?；阻帶最小衰減AS?40dB。程序: wp=0.4*pi;ws=0.25*pi;rp=1;rs=40;fs=500;ts=1/fs;wp1=wp*ts;ws1=ws*ts;

wp2=2*fs*tan(wp1/2);ws2=2*fs*tan(ws1/2);

[N,Wn]=buttord(wp2,ws2,rp,rs,'s');[Z,P,K]=buttap(N);[Bap,Aap]=zp2tf(Z,P,K);[b,a]=lp2lp(Bap,Aap,Wn);[bz,az]=bilinear(b,a,fs);[H,W]=freqz(bz,az);subplot(2,1,1);plot(W/pi,abs(H));grid on;xlabel('頻率')ylabel('幅度')subplot(2,1,2);

plot(W/pi,20*log10(abs(H)));grid on;xlabel('頻率');ylabel('幅度(dB)');

實驗七 FIR數(shù)字濾波器的設(shè)計

一、實驗?zāi)康?/p>

掌握用窗函數(shù)設(shè)計FIR數(shù)字濾波的原理及其設(shè)計步驟；熟悉線性相位數(shù)字濾波器的特性。學(xué)習(xí)編寫數(shù)字濾波器的設(shè)計程序的方法，并能進行正確編程；根據(jù)給定的濾波器指標(biāo)，給出設(shè)計步驟。

二、實驗原理

如果系統(tǒng)的沖激響應(yīng)h（n）為已知，則系統(tǒng)的輸入輸出關(guān)系為

y（n）=x（n）*h（n）

對于低通濾波器，只要設(shè)計出低通濾波器的沖激響應(yīng)函數(shù)，就可以由式得到系統(tǒng)的輸出了。

但是將h（n）作為濾波器的脈沖響應(yīng)有兩個問題：一是它是無限長的；二是它是非因果的。對此，采取兩項措施：一是將h（n）截短；二是將其右移。

設(shè)計時，要根據(jù)阻帶的最小衰減和過渡帶寬度來選擇恰當(dāng)?shù)拇昂瘮?shù)類型和窗口長度N。常用的窗函數(shù)有矩形窗、海明窗和布萊克曼窗等。

窗函數(shù)設(shè)計FIR濾波器步驟如下：

（1）給定理想頻率響應(yīng)的幅頻特性和相頻特性；

（2）求理想單位脈沖響應(yīng)，在實際計算中，可對理想頻率響應(yīng)采樣。（3）根據(jù)過渡帶寬度和阻帶最小衰減，確定窗函數(shù)類型和窗口長度N；（4）求FIR濾波器單位脈沖響應(yīng)；

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

三、實驗內(nèi)容及分析

1、分別用海明窗和布萊克曼窗設(shè)計一個48階的FIR帶通濾波器，通帶為Wn??0.450.55?。程序1:海明窗設(shè)計

N=48;

Window=hamming(N+1);w1=0.45;w2=0.55;ws=[w1,w2];

b=fir1(N,ws/pi,Window);freqz(b,1,512);title('海明窗');grid on;

程序2:萊克曼窗設(shè)計

N=48;

Window=blackman(N+1);w1=0.45;w2=0.55;ws=[w1,w2];

b=fir1(N,ws/pi,Window);freqz(b,1,512);title('布萊克曼窗');grid on;

2、用矩形窗設(shè)計一個線性相位高通濾波器。其中He??jw?e?j????????00.3?????

0???0.3?程序: N=9;

alpha=(N-1)/2;Wc=0.7*pi;n=(0:8);i=n-alpha;i=i+(i==0)*eps;

h=(-1).^n.*sin((i).*Wc)./((i).*pi);%矩形窗函數(shù)設(shè)計的系統(tǒng)脈沖響應(yīng) w=(0:1:500)*2*pi/500;

H=h*exp(-j*n'*w);%矩形窗函數(shù)設(shè)計的頻響 magH=abs(H);% 矩形窗函數(shù)設(shè)計的振幅 subplot(211);stem(n,h);

axis([0,8,-0.4,0.4]);title('矩形窗設(shè)計h(n)');line([0,10],[0,0]);xlabel('n');ylabel('h');subplot(212);plot(w/pi,magH);

xlabel('以pi為單位的頻率');ylabel('H振幅');axis([0,2,0,1.7]);title('矩形窗設(shè)計振幅譜');

實驗心得體會：

這次實驗使我進一步加深了對MATLAB軟件的使用。從上次的信號系統(tǒng)實驗的初步使用到這一次的深入了解，有了更深刻的認識。對這種語言環(huán)境也有了新的了解。

在實驗的過程中，我對數(shù)字濾波器的整個過程有了很好的理解和掌握。IIR數(shù)字濾波器的設(shè)計讓我知道了巴特沃思濾波器和切比雪夫濾波器的頻率特性，還有雙線性變換及脈沖響應(yīng)不變法設(shè)計的濾波器的頻率特性。做這兩個實驗的時候程序有點困難，但經(jīng)過細心的改寫圖形最終出來了。FIR數(shù)字濾波器的設(shè)計出來的是兩種窗的圖形，通過兩種窗的比較，我了解了他們各自的特點，幅頻和相頻特性。

最后，感謝張老師對我的諄諄教導(dǎo)！