第一篇：數字信號處理第六章介紹

第六章

數字濾波器結構

6.1：級聯的實現

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

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

級聯框圖：

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

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

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 級聯框圖：

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

6.2：級聯和并聯實現

Q6.3使用程序P6.1生成如下因果無限沖激響應傳輸函數的級聯實現： 畫出級聯實現的框圖。答：

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

級聯實現框圖：

Q6.4使用程序P6.1生成如下因果無限沖激響應傳輸函數的級聯實現：

畫出級聯實現的框圖。答：級聯實現框圖：

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

并聯II型框圖：

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

并聯II型框圖：

6.3：全通傳輸函數的實現

Q6.7使用程序P4.4生成如下全通傳輸函數的級聯格型實現:

As(z)是一個穩定的傳輸函數嗎? 答：

運行結果：

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

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

Q6.8使用程序P4.4生成如下全通傳輸函數的級聯格型實現:

A6(z)足一個穩定的傳輸函數嗎? 答：得到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}的值可以得到傳輸函數A6(z)是穩定的，因為反饋系數的平均幅值小于整體。

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

使用1型和2型全通項生成所示全通函數的典范級聯實現，實現的結構框圖如下：

整體結構中乘法器的總數是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)的典范級聯實現框圖如下：

整體結構中乘法器的總數是6。6.4：無限沖激響應傳輸函數的Gary-Markel實現

num = input('分子系數向量 = ');den = input('分母系數向量 = ');N = length(den)-1;% 分母多項式的階數 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('格型參數是');disp(k)disp('前饋乘法器是');disp(alpha)Q6.11 使用程序 P6_3我們通過IIR將Q6.3給的正向傳輸函數H1(z)的 Gray-Markel級聯格型實現參數如下： 晶格參數和前饋乘數分別如下：

對應Gray-Markel的結構框圖如下：

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

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

對應Gray-Markel的結構框圖如下：

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

Q6.13使用函數tf2latc編寫出一個MATLAB程序，以生成一個因果無限沖激響應傳輸函數的GrayMarkel實現。用該程序實現式(6.27)所示的傳輸函數。你的結果與習題6.11中得到的結果相符嗎？使用函數1atc2tf由向量k和alpha確定傳輸函數。所得到的傳輸函數和式(6.27)給出的傳輸函數相同嗎？ 答：程序如下： 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)運行結果如下： k = 0.62459686089013 0.68373782742919 0.48111942348398 0.21960784313725 0.06250000000000 alpha =-0.01982100623522-0.09085169508677 0.***849 0.16053921568627 0.31250000000000-0.12500000000000 結果與習題6.11中得到的結果相符。Q6.14使用在習題6.13中生成的程序，實現式(6.28)給出的傳輸函數。你的結果與習題6.12中得到的結果相符嗎？使用函數latc2tf由向量k和alpha確定傳輸函數。所得到的傳輸函數和式(6.28)給出的傳輸函數相同嗎? 答：運行結果： 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中得到的結果相符。

6.5：無限沖激響應傳輸函數的并聯全通實現

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

G(z)全通和的分解：

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

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

G(z)全通和的分解：

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

兩個全通傳輸函數的階數是3和2.

第二篇：數字信號處理第四章介紹

第四章 線性時不變離散時間系統的頻域分析

一、傳輸函數和頻率響應 例4.1傳輸函數分析 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，計算并畫出當w=[0,pi]時傳輸函數因果線性時不變離散時間系統的頻率響應。它表示哪種類型的濾波器？ w = 0:pi/511:pi;

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

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

Q4.4 使用MATLAB計算并畫出當w=[0,pi]時因果線性時不變離散時間系統的群延遲。系統的傳輸函數為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）確定的兩個濾波器的沖激響應中的前一百個樣本。討論你的結果。clf;num = [0.15 0-0.15];den = [0.7-0.5 1];L = input('輸入樣本數 L: ');[g t] = impz(num,den,L);stem(t,g);title(['前 ',num2str(L),' 脈沖響應的樣本']);xlabel('時間序號 n');ylabel('h[n]');

（4.36）式（4.37）式

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

Q4.6 傳輸函數的極零點圖同樣能分析線性時不變離散時間系統的性質。使用命令zplane可以很容易地得到系統的極零點圖。使用zplane分別生成式（4.36）和式（4.37）確定的兩個濾波器的極零點圖。討論你的結果。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)在單位圓和兩極因此它的因果實現穩定；較低的圖顯示過濾器(4.37)極點在單位圓外,其因果關系的實現是不穩定的。

二、傳輸函數的類型 例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]，而控制截止頻率的參數是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;修改參數fc和n，得到如上圖2，可知低通有限沖激濾波器的長度變為20.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的四幅圖，從這四幅圖可以看出隨著階數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 設計一個在0.45pi處具有3dB截止頻率wc的一階無限沖激響應低通濾波器和一階無限沖激響應高通濾波器。用MATLAB計算并畫出它們的增益響應，驗證設計的濾波器是否滿足指標。用MATLAB證明兩個濾波器是全通互補和功率互補的。

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

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

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

從這些情節我們觀察,梳狀濾波器極距為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型有限沖激響應濾波器的零點是 Zeros of Type 1 FIR Filter are 2.9744 2.0888 0.9790 + 1.4110i 0.97900.4784i 0.4787 0.3362 2型有限沖激響應濾波器的零點是 Zeros of Type 2 FIR Filter are 3.7585 + 1.5147i 3.75852.6623i-1.0000 0.0893 + 0.3530i 0.08930.0922i 3型有限沖激響應濾波器的零點是 Zeros of Type 3 FIR Filter are 4.7627 1.6279 + 3.0565i 1.62790.2549i 0.2100 4型有限沖激響應濾波器的零點是 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型有限沖激響應濾波器的零點是 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 確定如下傳輸函數是否是有界實函數：一個有界實函數，求一個與

有著相同幅度的有界實函數。

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

故H2(z)為

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

有著相同幅度的有界實函數。

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

Q4.23 用MATLAB產生如下兩個因果系統傳輸函數的極零點圖：，研究生成的極零點圖，你可以推斷它們的穩定性么？

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

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

由此我們可以總結出H1(z)穩定，H2(z)不穩定

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

由此看出，都在單位圓內。

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

由此看出，都在單位圓內。

第三篇：數字信號處理第七章介紹

第七章 數字濾波器設計

7.1：無限沖激響應濾波器的階數的估計

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

標準阻帶邊緣角頻率Ws是：

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

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

標準阻帶邊緣角頻率Ws是：

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

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

標準阻帶邊緣角頻率Ws是：

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

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

答：標準通帶邊緣角頻率Wp是：

標準阻帶邊緣角頻率Ws是：

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

7.2：無限沖激響應濾波器設計

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

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

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

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

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

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

濾波器參數是：

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

從圖中可以總結出設計符合指標。

濾波器的未畸變的相位響應及群延遲響應如下：

Q7.6修改程序P7.1來設計符合習題Q7.1所給指標的切比雪夫1型低通濾波器。寫出所產生的傳輸函數的準確表達式。你的設計符合指標嗎?使用MATLAB,計算并繪制濾波器的未畸變的相位響應及群延遲響應。答：表達式如下：

設計的濾波器增益響應如下：

從圖中可以總結出設計符合指標。

濾波器的未畸變的相位響應及群延遲響應如下：

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

設計的濾波器增益響應如下：

從圖中可以總結出設計符合指標。濾波器的未畸變的相位響應及群延遲響應如下：

Q7.8修改程序P7.1來設計符合習題Q7.3所給指標的橢圓帶通濾波器。寫出所產生的傳輸函數的準確表達式。你的設計符合指標嗎，使用MATLAB，計算井繪制濾波器的未畸變的相位響應及群延遲響應。答：表達式如下：

設計的濾波器增益響應如下：

從圖中可以總結出設計符合指標。

濾波器的未畸變的相位響應及群延遲響應如下：

7.3：吉布斯現象

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

長度分別為61,41和21的幅度響應如下：

從中可以觀察到由于吉布斯現象產生的幅度響應的擺動行為。

波紋的數量與濾波器的長度之間的關系——波紋的數量減少與長度成正比。最大波紋的高度與濾波器的長度之間的關系——最大波紋的高度與長度無關。

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

從中可以觀察到由于吉布斯現象產生的幅度響應擺動行為。

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

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

從中可以觀察到由于吉布斯現象產生的幅度響應的擺動行為。波紋的數量與微分器的長度之間的關系——兩者成正比。

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

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

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

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

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

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

濾波器階數和采樣頻率的關系為–對于一個給定的模擬過渡帶寬，采樣頻率的增加導致估算階數也相應增加，朝下一個整數取整。

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

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

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

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

p = 0.01,阻帶波紋

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

= 0.002 結果為 N=57。s

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

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

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

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

7.5有限沖激響應濾波器設計

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

增益和相位響應如下:

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

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');低通濾波器系數：

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 增益和相位響

從圖中可以看出此時的濾波器不滿足指標。欲滿足指標，應調節N=47.Q7.23 用凱澤窗設計一個有限沖激響應低通濾波器。程序：

% 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');低通濾波器系數：

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 增益和相位響應如下：

從圖中可以看出設計的濾波器滿足要求。N=60.Q7.24 用函數kaiserord和firl重做習題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');參數如下：

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

增益和相位響應如下：

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

程序：

% 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');

增益響應： 相位響應：

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

程序：

% 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;增益響應：

從增益響應的圖像可以看出，該濾波器不滿足設計指標。用kaiserord估計濾波器的階數N=73。

第四篇：數字信號處理課程設計..

課程設計報告

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

姓 名： 學 號： 院 系： 專業班級： 指導教師： 完成日期： 2013年7月2日

目錄

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

第1部分 課程設計報告

一．設計目的

綜合運用本課程的理論知識進行頻譜分析以及濾波器設計，通過理論推導得出相應結論，并利用MATLAB作為工具進行實現，從而復習鞏固課堂所學的理論知識，提高對所學知識的綜合應用能力，并從實踐上初步實現對數字信號的處理。

二．設計內容

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

三．設計原理

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

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

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

4.在Matlab中，FIR濾波器利用函數fftfilt對信號進行濾波。

5.抽樣定理

連續信號經理想抽樣后時域、頻域發生的變化（理想抽樣信號與連續信號頻譜之間的關系）

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

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

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

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?)的關系：理想抽樣信號的頻譜是連續信號頻譜的Xaa

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

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

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

四．具體實現

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('錄音結束');%文字提示 fprintf('按Enter鍵回放錄音');pause;%暫停命令

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

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

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

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

grid;%加網格

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

1.2濾除無效點

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

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

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

圖1.2 去除無效點

2.巴特沃斯濾波器的設計

2.1設計巴特沃思低通濾波器

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設計巴特沃思高通濾波器

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設計巴特沃思帶通濾波器

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數字信號sound.wav');%播放原始信號

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

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

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

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('低通濾波后信號的時域波形');%加標題 xlabel('頻率w/pi');ylabel('幅值/A');%顯示標表示的意義 grid;%網格

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

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

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');%顯示標表示的意義 title('高通濾波后信號的時域波形');%標題 ylabel('幅值/A');%顯示縱坐標的表示意義 grid;%網格

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

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

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('帶通濾波后信號的時域波形');%加標題 xlabel('頻率w/pi');ylabel('幅值/A');%顯示縱坐標表示的意義 grid;%網格

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

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

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

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

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

4.語音信號的回放

sound(xlow,Fs,bits);%在Matlab中，函數sound可以對聲音進行回放，其調用格式: 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('錄音結束');%文字提示 fprintf('按Enter鍵回放錄音');pause;%暫停命令 wavplay(x,Fs);%按任意鍵播放語音信號

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

5.2 分析男女聲音的頻譜

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

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

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

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

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

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

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

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

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

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

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('錄音結束');%文字提示 fprintf('按Enter鍵回放錄音');pause;%暫停命令 wavplay(x,Fs);%按任意鍵播放語音信號 wavwrite(x,Fs,'C:UsersacerDesktop數字信號噪音.wav');%存儲語音信號

6.2 對噪聲進行頻譜的分析

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

wavplay(x1,fs);%播放噪聲信號 N=length(x1);%返回采樣點數 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數字信號sound.wav ');%對錄入信號進行采樣

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

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

6.4疊加信號的頻譜分析：

wavplay(yy,fs);%播放疊加后信號 N=length(yy);%返回采樣點數 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 設計一個合適的濾波器將噪聲濾除 fs=18000;%采樣頻率 Wp=2*1000/fs;%通帶截至頻率 Ws=2*2000/fs;%阻帶截至頻率 Rp=1;%最大衰減 Rs=100;%最小衰減

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

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

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

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

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

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.結果分析

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

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

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

頻成分較多。

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

第2部分 課程設計總結

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

八．參考文獻

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

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

[5] 數字信號處理及其MATLAB實驗 趙紅怡 張常 化學工業出版社

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

錢同惠 機械工業出版社 29

第五篇：數字信號處理實驗報告

JIANGSU

UNIVERSITY OF TECHNOLOGY

數字信號處理實驗報告

學院名稱： 電氣信息工程學院

專 業：

班 級： 姓 名： 學 號： 指導老師： 張維璽（教授）

2013年12月20日

實驗一 離散時間信號的產生

一、實驗目的

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

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

二、實驗原理

離散時間信號是指在離散時刻才有定義的信號，簡稱離散信號，或者序列。離散序列通常用x（n）來表示，自變量必須是整數。常見的離散信號如下：（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）矩形序列

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

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

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

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

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

（6）負指數序列x（n）

當a=0時，得到虛指數序列x（n）。

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

1）當a>0時，負指數序列x（n）的實部和虛部分別是按指數規律增長的正弦振蕩序列；

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

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

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

stem(n,imag(x))

三、實驗內容及分析

?1n?01、編制程序產生單位沖激序列??n???“?并繪出其圖及??n?”學號后兩位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;%學號01 n=-(N-1):N-1;

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

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

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

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

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

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

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

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

M=1;%學號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、編制程序產生正弦序列x?n??cos?2?n?、x?n??cos??n?及

?學號后兩位?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;%學號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、編制程序產生復正弦序列x?n??e(2?j學號后兩位）n，并繪出其圖形。N=3;

n=0:0.2:N-1;

w=1;%學號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、編制程序產生指數序列x?n??an，并繪出其圖形。其中a=學號后兩位、a=1/“學號后兩位”。

(1)N=10;

n=0:N-1;

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

(2)N=10;

n=0:N-1;

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

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

一、實驗目的

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

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

二、實驗原理

（１）DTFT和DFT的定義及其相互關系。

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

三、實驗內容及分析

（１）編程計算并畫出下面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系統的Z域分析

一、實驗目的

本實驗通過使用MATLAB函數對離散時間系統的一些特性進行仿真分析，以加深對離散時間系統的零極點、穩定性，頻率響應等概念的理解。學會運用MATLAB分析離散時間系統的系統函數的零極點；學會運用MATLAB分析系統函數的零極點分布與其時域特性的關系；學會運用MATLAB進行離散時間系統的頻率特性分析。

二、實驗原理

離散時間系統的系統函數定義為系統零狀態響應的Z變化與激勵的Z變化之比。

在MATLAB中系統函數的零極點可通過函數roots得到，也可借助函數tf2zp得到，tf2zp的語句格式為

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

若要獲得系統函數H（z）的零極點分布圖，可直接應用zplane函數，其語句格式為

Zplane（B，A）

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

離散系統中z變化建立了時域函數h（n）與z域函數H（z）之間的對應關系。因此，z變化的函數H（z）從形式可以反映h（n）的部分內在性質。可根據系統的傳遞函數H（z）求單位沖激響應h（n）的函數impz、filter等。

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

MATLAB提供了求離散時間系統頻響特性的函數freqz，調用freqz的格式主要有兩種。一種形式為

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

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

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

三、實驗內容與結果分析

已知LTI離散時間系統，要求由鍵盤實現系統參數輸入，并繪出幅頻和相頻響應曲線和零極點分布圖，進而分析系統的濾波特性和穩定性。

（一）程序

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 幅頻、相頻響應曲線、零極點分布圖

實驗六 IIR數字濾波器的設計

一、實驗目的

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

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

二、實驗原理

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

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

IIR數字濾波器的設計方法如下：（1）沖激響應不變法。（2）雙線性變化法。

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

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

MATLAB提供了一組標準的數字濾波器設計函數，大大簡化了濾波器的設計工程。

（1）butter。

（2）cheby1、cheby2。

三、實驗內容及分析

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

1、基于chebyshev I型模擬濾波器原型使用沖激不變轉換方法設計數字濾波器，要求參數為通帶截止頻率?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型模擬濾波器原型使用雙線性變換方法設計數字濾波器的，要求參數為截止頻率?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數字濾波器的設計

一、實驗目的

掌握用窗函數設計FIR數字濾波的原理及其設計步驟；熟悉線性相位數字濾波器的特性。學習編寫數字濾波器的設計程序的方法，并能進行正確編程；根據給定的濾波器指標，給出設計步驟。

二、實驗原理

如果系統的沖激響應h（n）為已知，則系統的輸入輸出關系為

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

對于低通濾波器，只要設計出低通濾波器的沖激響應函數，就可以由式得到系統的輸出了。

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

設計時，要根據阻帶的最小衰減和過渡帶寬度來選擇恰當的窗函數類型和窗口長度N。常用的窗函數有矩形窗、海明窗和布萊克曼窗等。

窗函數設計FIR濾波器步驟如下：

（1）給定理想頻率響應的幅頻特性和相頻特性；

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

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

三、實驗內容及分析

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

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:萊克曼窗設計

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、用矩形窗設計一個線性相位高通濾波器。其中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);%矩形窗函數設計的系統脈沖響應 w=(0:1:500)*2*pi/500;

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

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

實驗心得體會：

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

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

最后，感謝張老師對我的諄諄教導！