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

第四章 線性時(shí)不變離散時(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

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

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

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

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

Q4.4 使用MATLAB計(jì)算并畫出當(dāng)w=[0,pi]時(shí)因果線性時(shí)不變離散時(shí)間系統(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中編寫的程序，分別計(jì)算并畫出式（4.36）和式（4.37）確定的兩個(gè)濾波器的沖激響應(yīng)中的前一百個(gè)樣本。討論你的結(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('時(shí)間序號(hào) n');ylabel('h[n]');

（4.36）式（4.37）式

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

Q4.6 傳輸函數(shù)的極零點(diǎn)圖同樣能分析線性時(shí)不變離散時(shí)間系統(tǒng)的性質(zhì)。使用命令zplane可以很容易地得到系統(tǒng)的極零點(diǎn)圖。使用zplane分別生成式（4.36）和式（4.37）確定的兩個(gè)濾波器的極零點(diǎn)圖。討論你的結(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)在單位圓和兩極因此它的因果實(shí)現(xiàn)穩(wěn)定；較低的圖顯示過濾器(4.37)極點(diǎn)在單位圓外,其因果關(guān)系的實(shí)現(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低通有限沖激濾波器的長(zhǎng)度為14，決定濾波器長(zhǎng)度的語句為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，可知低通有限沖激濾波器的長(zhǎng)度變?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)])

可以驗(yàn)證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è)計(jì)一個(gè)在0.45pi處具有3dB截止頻率wc的一階無限沖激響應(yīng)低通濾波器和一階無限沖激響應(yīng)高通濾波器。用MATLAB計(jì)算并畫出它們的增益響應(yīng)，驗(yàn)證設(shè)計(jì)的濾波器是否滿足指標(biāo)。用MATLAB證明兩個(gè)濾波器是全通互補(bǔ)和功率互補(bǔ)的。

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

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

Q4.17 用MATLAB計(jì)算并畫出一個(gè)梳狀濾波器的幅度響應(yīng)，該梳狀濾波器是在L取不同值的情況下，由式（4.40）給出的原型有限沖激響應(yīng)低通濾波器得到的。證明新濾波器的幅度響應(yīng)在k=0,1,2,3......，L-1.處有L個(gè)極小值，在處有L個(gè)極大值，Q4.18 用MATLAB計(jì)算并畫出一個(gè)梳狀濾波器的幅度響應(yīng)，該梳狀濾波器是在L取不同值的情況下，由式（4.42）在M=2時(shí)給出的原型有限沖激響應(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)濾波器的零點(diǎn)是 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)濾波器的零點(diǎn)是 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)濾波器的零點(diǎn)是 Zeros of Type 3 FIR Filter are 4.7627 1.6279 + 3.0565i 1.62790.2549i 0.2100 4型有限沖激響應(yīng)濾波器的零點(diǎn)是 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)濾波器的零點(diǎn)是 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ù)：一個(gè)有界實(shí)函數(shù)，求一個(gè)與

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

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

故H2(z)為

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Q7.1用MATTAB確定一個(gè)數(shù)字無限沖激響應(yīng)低通濾波器所有四種類型的最低階數(shù)。指標(biāo)如下:40 kHz的抽樣率，,4 kHz的通帶邊界頻率，8 kHz的阻帶邊界頻率，0.5 dB的通帶波紋，40 dB的最小阻帶衰減。評(píng)論你的結(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確定一個(gè)數(shù)字無限沖激響應(yīng)高通濾波器所有四種類型的最低階數(shù)。指標(biāo)如下:3500Hz的抽樣率，1050 Hz的通帶邊界頻率,600 Hz的阻帶邊界頻率，1 dB的通帶波紋，50 dB的最小阻帶衰減。評(píng)論你的結(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確定一個(gè)數(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的最小阻帶衰減。評(píng)論你的結(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確定一個(gè)數(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的最小阻帶衰減。評(píng)論你的結(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è)計(jì)

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

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

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

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

disp('分子系數(shù)是 ');disp(num);disp('分母系數(shù)是 ');disp(den);% 計(jì)算增益響應(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通過運(yùn)行程序P7.1來設(shè)計(jì)巴特沃茲帶阻濾波器。寫出所產(chǎn)生的傳輸函數(shù)的準(zhǔn)確表達(dá)式。濾波器的指標(biāo)是什么，你的設(shè)計(jì)符合指標(biāo)嗎，使用MATLAB，計(jì)算并繪制濾波器的未畸變的相位響應(yīng)及群延遲響應(yīng)。答：表達(dá)式是：

濾波器參數(shù)是：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

在這種情況下你不能改變長(zhǎng)度。原因：這是一個(gè)零相位濾波器，這意味著它也是一個(gè)線性相位濾波器，因?yàn)榱阆嗍且环N特殊的線性相位的子集?，F(xiàn)在，理想的有限脈沖響應(yīng)長(zhǎng)度甚至有對(duì)稱的中點(diǎn)h[n]。使其成了一個(gè)線性相位FIR濾波器。二型濾波器不可能是高通濾波器，因?yàn)楸仨氃趜=-1處有零點(diǎn)，意味著w=+-π。

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

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

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

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

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

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的對(duì)比我們可以觀察到：

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

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

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

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ù)是較小的。因?yàn)閯P瑟使用了一個(gè)不同的近似估計(jì)。這種估計(jì)經(jīng)常和FIR設(shè)計(jì)的凱瑟窗一起用。

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給了一個(gè)比凱爾更大比凱瑟更小一點(diǎn)的結(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時(shí)結(jié)果為：kaiord([1800 3600],[1200 4200],0.01,0.02,12000)，得到的N=33。所以答案不唯一，可以選擇其中一個(gè)。Q7.18 線性相位帶通FIR濾波器的階數(shù)估算，其中FIR濾波器的規(guī)格和7.17一樣，則使用kaiserord的結(jié)果為同樣，它也有矛盾。當(dāng)使用δp= 0.1時(shí)，得到的結(jié)果為：kaiserord([1200 1800 3600 4200],[0 1 0],[0.02 0.1 0.02],12000)，則N=37.當(dāng)用δp= 0.01時(shí)，結(jié)果為：kaiserord([1200 1800 3600 4200],[0 1 0],[0.02 0.01 0.02],12000)，此時(shí)N=45.和7.17的結(jié)果比較我們觀察到通過kaiserord函數(shù)估計(jì)的階數(shù)要更高，但如果你要設(shè)計(jì)Kaiser窗的話則結(jié)果更精確。

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

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

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

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

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

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

從增益圖像我們可以知道這個(gè)設(shè)計(jì)不能滿足規(guī)格.這個(gè)濾波器滿足規(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 增益和相位響

從圖中可以看出此時(shí)的濾波器不滿足指標(biāo)。欲滿足指標(biāo)，應(yīng)調(diào)節(jié)N=47.Q7.23 用凱澤窗設(shè)計(jì)一個(gè)有限沖激響應(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è)計(jì)的濾波器滿足要求。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è)計(jì)的濾波器滿足要求。N=59.Q7.25 用fir2設(shè)計(jì)一個(gè)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è)計(jì)具有如下指標(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è)計(jì)指標(biāo)。用kaiserord估計(jì)濾波器的階數(shù)N=73。

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

第六章

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

6.1：級(jí)聯(lián)的實(shí)現(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ù)的一個(gè)級(jí)聯(lián)實(shí)現(xiàn): H1(z)=2+10z^(-1)+23z^(-2)+34z^(-3)+31z^(-4)+16 z^(-5)+4z^(-6)畫出級(jí)聯(lián)實(shí)現(xiàn)的框圖。H1(z)是一個(gè)線性相位傳輸函數(shù)嗎? 答：運(yùn)行結(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

級(jí)聯(lián)框圖：

H1（z）不是一個(gè)線性相位傳輸函數(shù)，因?yàn)橄禂?shù)不對(duì)稱。

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

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

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

Q6.3使用程序P6.1生成如下因果無限沖激響應(yīng)傳輸函數(shù)的級(jí)聯(lián)實(shí)現(xiàn)： 畫出級(jí)聯(lián)實(shí)現(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

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

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

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

程序P6.2生成兩種類型的并聯(lián)實(shí)現(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('極點(diǎn)在');disp(p1);disp('常數(shù)');disp(k1);disp('并聯(lián)II型')disp('留數(shù)是');disp(r2);disp('極點(diǎn)在');disp(p2);disp('常數(shù)');disp(k2);Q6.5使用程序P6.2生成式(6.27)所示因果無限沖激響應(yīng)傳輸函數(shù)的兩種不同并聯(lián)形式實(shí)現(xiàn)。畫出兩種實(shí)現(xiàn)的框圖。答：并聯(lián)I型框圖：

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

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

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

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

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

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

運(yùn)行結(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)定的，因?yàn)閷?duì)所有的1

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

A6(z)足一個(gè)穩(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)定的，因?yàn)榉答佅禂?shù)的平均幅值小于整體。

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

使用1型和2型全通項(xiàng)生成所示全通函數(shù)的典范級(jí)聯(lián)實(shí)現(xiàn)，實(shí)現(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型的全通項(xiàng)生成A6(z)的典范級(jí)聯(lián)實(shí)現(xiàn)框圖如下：

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

num = input('分子系數(shù)向量 = ');den = input('分母系數(shù)向量 = ');N = length(den)-1;% 分母多項(xiàng)式的階數(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級(jí)聯(lián)格型實(shí)現(xiàn)參數(shù)如下： 晶格參數(shù)和前饋乘數(shù)分別如下：

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

使用程序P6_3，從這些格型參數(shù)可以得到傳輸函數(shù)H1(z)是穩(wěn)定的，因?yàn)樗懈裥蛥?shù)的平方值比整體的小。

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

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

使用程序P6_3，從這些格型參數(shù)可以得到傳輸函數(shù)H2(z)是穩(wěn)定的，因?yàn)樗懈裥蛥?shù)的平方值比整體的小。

Q6.13使用函數(shù)tf2latc編寫出一個(gè)MATLAB程序，以生成一個(gè)因果無限沖激響應(yīng)傳輸函數(shù)的GrayMarkel實(shí)現(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)運(yùn)行結(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中生成的程序，實(shí)現(xiàn)式(6.28)給出的傳輸函數(shù)。你的結(jié)果與習(xí)題6.12中得到的結(jié)果相符嗎？使用函數(shù)latc2tf由向量k和alpha確定傳輸函數(shù)。所得到的傳輸函數(shù)和式(6.28)給出的傳輸函數(shù)相同嗎? 答：運(yùn)行結(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)全通實(shí)現(xiàn)

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

G(z)全通和的分解：

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

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

G(z)全通和的分解：

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

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

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

課程設(shè)計(jì)報(bào)告

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

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

目錄

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

第1部分 課程設(shè)計(jì)報(bào)告

一．設(shè)計(jì)目的

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

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

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

三．設(shè)計(jì)原理

1.在Matlab軟件平臺(tái)下，利用函數(shù)wavrecord()，wavwrite(),wavread(),wavplay()對(duì)語音信號(hào)進(jìn)行錄制，存儲(chǔ)，讀取，回放。

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

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

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

5.抽樣定理

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

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

理想采樣過程描述： 時(shí)域描述：

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

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)系：理想抽樣信號(hào)的頻譜是連續(xù)信號(hào)頻譜的Xaa

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

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

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

四．具體實(shí)現(xiàn)

1.錄制一段聲音

1.1錄制并分析

在MATLAB中用wavrecord、wavread、wavplay、wavwrite對(duì)聲音進(jìn)行錄制、讀取、回放、存儲(chǔ)。

程序如下：

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

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

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

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

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

grid;%加網(wǎng)格

圖1.1 原始信號(hào)的時(shí)域與頻譜圖

1.2濾除無效點(diǎn)

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

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

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

圖1.2 去除無效點(diǎn)

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

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

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è)計(jì)巴特沃思高通濾波器

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è)計(jì)巴特沃思帶通濾波器

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．將聲音信號(hào)送入濾波器濾波

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

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

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

y0=fft(x);%快速傅立葉變換 subplot(4,2,3);plot(f,abs(y0(n1)));%原始信號(hào)的頻譜圖 title('原始信號(hào)的頻譜圖');%加標(biāo)題 xlabel('頻率w/pi');%顯示橫坐標(biāo)表示的意義 ylabel('幅值 ');%顯示縱坐標(biāo)表示的意義 title('原始信號(hào)的頻譜圖');%加標(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);%信號(hào)送入濾波器濾波，ys為輸出 fftwave=fft(ys);%將濾波后的語音信號(hào)進(jìn)行快速傅立葉變換 figure(4);subplot(4,2,2);%在四行兩列的第二個(gè)窗口顯示圖形 plot(ys);%信號(hào)的時(shí)域波形

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

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

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

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

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);%將信號(hào)送入高通濾波器濾波 subplot(4,2,5);%在四行兩列的第五個(gè)窗口顯示圖形 plot(ys);%信號(hào)的時(shí)域波形 xlabel('頻率w/pi');ylabel('幅值/A');%顯示標(biāo)表示的意義 title('高通濾波后信號(hào)的時(shí)域波形');%標(biāo)題 ylabel('幅值/A');%顯示縱坐標(biāo)的表示意義 grid;%網(wǎng)格

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

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

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

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);%將信號(hào)送入高通濾波器濾波 figure(4);subplot(4,2,6);%在四行兩列的第六個(gè)窗口顯示圖形 plot(ys1);%繪制后信號(hào)的時(shí)域的圖形 title('帶通濾波后信號(hào)的時(shí)域波形');%加標(biāo)題 xlabel('頻率w/pi');ylabel('幅值/A');%顯示縱坐標(biāo)表示的意義 grid;%網(wǎng)格

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

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

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

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

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

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

4.語音信號(hào)的回放

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

5.1 錄制一段異性的聲音進(jìn)行頻譜分析

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

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

5.2 分析男女聲音的頻譜

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

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

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

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

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

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

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

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

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

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

5.3男女聲音的頻譜圖

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

圖5.3男女聲音信號(hào)波形與頻譜對(duì)比

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

6.噪聲的疊加和濾除

6.1錄制一段背景噪聲

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

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

6.2 對(duì)噪聲進(jìn)行頻譜的分析

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

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

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

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

6.3原始信號(hào)與噪音的疊加

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

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

yy=x+x1;%將兩個(gè)聲音疊加

6.4疊加信號(hào)的頻譜分析：

wavplay(yy,fs);%播放疊加后信號(hào) N=length(yy);%返回采樣點(diǎn)數(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);%信號(hào)的時(shí)域波形

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

6.5 設(shè)計(jì)一個(gè)合適的濾波器將噪聲濾除 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);%信號(hào)送入濾波器濾波，ys為輸出 fftwave=fft(ys);%將濾波后的語音信號(hào)進(jìn)行快速傅立葉變換 figure(5);subplot(3,2,5);plot(ys);%信號(hào)的時(shí)域波形

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

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

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

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

7.結(jié)果分析

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

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

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

頻成分較多。

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

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

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

八．參考文獻(xiàn)

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

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

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

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

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

第五篇：數(shù)字信號(hào)處理實(shí)驗(yàn)報(bào)告

JIANGSU

UNIVERSITY OF TECHNOLOGY

數(shù)字信號(hào)處理實(shí)驗(yàn)報(bào)告

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

專 業(yè)：

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

2013年12月20日

實(shí)驗(yàn)一 離散時(shí)間信號(hào)的產(chǎn)生

一、實(shí)驗(yàn)?zāi)康?/p>

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

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

二、實(shí)驗(yàn)原理

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

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

n=-（N-1）：N-1 x=[zeros(1,N-1)1 zeros(1,N-1)]； stem(n,x)延遲K個(gè)采樣點(diǎn)的長(zhǎng)度為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）在時(shí)間軸上延遲了k個(gè)單位，得到u（n-k），即長(zhǎng)度為N的單位階躍序列u（n）可以通過下面的MATLAB命令獲得。

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

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

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

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

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

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

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

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

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

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

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

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

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

stem(n,imag(x))

三、實(shí)驗(yàn)內(nèi)容及分析

?1n?01、編制程序產(chǎn)生單位沖激序列??n???“?并繪出其圖及??n?”學(xué)號(hào)后兩位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é)號(hào)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個(gè)單位。

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

0n?0?u?n??u?n?“學(xué)號(hào)后兩位”?，并繪出其圖形。程序: 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é)號(hào)01 n=-(N-1):N-1;

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

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

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

M=1;%學(xué)號(hào)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é)號(hào)后兩位?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('余弦信號(hào)');grid on;

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

A=1;

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

;

分析:該序列具有周期性,且輸出為余弦信號(hào).(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('正弦信號(hào)');grid on;

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

n=0:0.2:N-1;

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

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

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

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

(1)N=10;

n=0:N-1;

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

(2)N=10;

n=0:N-1;

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

實(shí)驗(yàn)三 離散時(shí)間信號(hào)的頻域分析

一、實(shí)驗(yàn)?zāi)康?/p>

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

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

二、實(shí)驗(yàn)原理

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

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

三、實(shí)驗(yàn)內(nèi)容及分析

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

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

（２）計(jì)算32點(diǎn)序列

5?n16，0≦n≦31的32點(diǎn)和64點(diǎn)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('實(shí)部 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點(diǎn)');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點(diǎn)');

（１）

（２）

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

一、實(shí)驗(yàn)?zāi)康?/p>

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

二、實(shí)驗(yàn)原理

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

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

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

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

Zplane（B，A）

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

離散系統(tǒng)中z變化建立了時(shí)域函數(shù)h（n）與z域函數(shù)H（z）之間的對(duì)應(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)對(duì)各種頻率成分的響應(yīng)特性，并推出系統(tǒng)的特性（高通、低通、帶通、帶阻等）。

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

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

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

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

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

已知LTI離散時(shí)間系統(tǒng)，要求由鍵盤實(shí)現(xiàn)系統(tǒng)參數(shù)輸入，并繪出幅頻和相頻響應(yīng)曲線和零極點(diǎn)分布圖，進(jìn)而分析系統(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('零極點(diǎn)圖');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)曲線、零極點(diǎn)分布圖

實(shí)驗(yàn)六 IIR數(shù)字濾波器的設(shè)計(jì)

一、實(shí)驗(yàn)?zāi)康?/p>

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

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

二、實(shí)驗(yàn)原理

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

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

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

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

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

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

（1）butter。

（2）cheby1、cheby2。

三、實(shí)驗(yàn)內(nèi)容及分析

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

1、基于chebyshev I型模擬濾波器原型使用沖激不變轉(zhuǎn)換方法設(shè)計(jì)數(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è)計(jì)數(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)');

實(shí)驗(yàn)七 FIR數(shù)字濾波器的設(shè)計(jì)

一、實(shí)驗(yàn)?zāi)康?/p>

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

二、實(shí)驗(yàn)原理

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

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

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

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

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

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

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

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

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

三、實(shí)驗(yàn)內(nèi)容及分析

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

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è)計(jì)

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è)計(jì)一個(gè)線性相位高通濾波器。其中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è)計(jì)的系統(tǒng)脈沖響應(yīng) w=(0:1:500)*2*pi/500;

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

axis([0,8,-0.4,0.4]);title('矩形窗設(shè)計(jì)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è)計(jì)振幅譜');

實(shí)驗(yàn)心得體會(huì)：

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

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

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