數字信號處理實驗指導書數字信號處理上機實驗指導書陳純鍇電子與信息工程學院一、引言“數字信號處理”是一門理論和實驗密切結合的課程，為了深入地掌握課程內容，應當在學習理論的同時，做習題和上機實驗。上機實驗不僅可以幫助學生深入地理解和消化基本理論，而且能鍛煉初學者的獨立解決問題的能力。所以，根據本課程的重點要求編寫了四個實驗。第一章是全書的基礎內容，抽樣定理、時域離散系統的時域和頻域分析以及系統對輸入信號的響應是重要的基本內容。由于第一章大部分內容已經在前期信號與系統課程中學習完，所以可通過實驗一幫助學生溫習以上重要內容，加深學生對“數字信號處理是通過對輸入信號的一種運算達到處理目的” 這一重要概念的理解。這樣便可以使學生從信號與系統課程順利的過渡到本課程的學習上來。第二章、三章DFT、FFT是數字信號處理的重要數學工具，它有廣泛的使用內容。采用實驗二、實驗三加深理解DFT的基本概念、基本性質。FFT是它的快速算法，必須學會使用。數字濾波器的基本理論和設計方法是數字信號處理技術的重要內容。學習這一部分時，應重點掌握IIR和FIR兩種不同的數字濾波器的基本設計方法。IIR濾波器的單位沖激響應是無限長的，設計方法是先設計模擬濾波器，然后再通過SZ平面轉換，求出相應的數字濾波器的系統函數。這里的平面轉換有兩種方法，即沖激響應不變法和雙線性變換法，后者沒有頻率混疊的缺點，且轉換簡單，是一種普遍應用的方法。FIR濾波器的單位沖激響應是有限長的，設計濾波器的目的即是求出符合要求的單位沖激響應。窗函數法是一種基本的，也是一種重要的設計方法。學習完第七章后可以進行實驗四。二、關于使用計算機語言由于數字信號處理實驗的主要目的是驗證數字信號處理的有關理論，進一步理解鞏固所學理論知識。所以，實現實驗用的算法語言可以有許多種，但為了提高實驗效率，要求學生用編程效率比C語言高好幾倍的MATLAB語言。下面介紹MATLAB的主要特點。（有關MATLAB的啟動、程序運行和有關信號處理工具箱函數等內容將放到最后附錄中介紹。）MATLAB是一種交互式的以矩陣為基本數據結構的系統。在生成矩陣對象時，不要求明確的維數說明。所謂交互式，是指MATLAB的草稿紙編程環境。即用戶每輸入一條命令并按回車鍵，MATLAB系統便解釋執行之，并顯示執行結果。根據該結果，用戶立即知道剛輸入的命令的正確性，或利用中間結果進行其他處理等。與C語言或FORTRON語言做科學數值計算的程序設計相比較，利用MATLAB可節省大量的編程時間。將其用于數字信號處理實驗，則可大大提高實驗效率，在有限的上機時間內，實驗內容可增加幾倍。例如，C語言FFT子程序有70多行，而用MATLAB只調用一個fft函數即可實現對序列進行FFT計算。另外，MATLAB的工具箱及圖形顯示（打印）功能，可滿足各層次人員直觀、方便的進行分析、計算和設計工作，從而可大大節省時間。例如，序列的卷積、濾波，系統函數H(z)的幅頻特性和相頻特性等計算，均有現成的工具箱函數。而用其它算法語言完成這些計算的編程比較麻煩，且程序較長。由于上述特點，在美國一些大學里，MATLAB已成為輔助教學的有益工具。MATLAB已成功地用于數字信號處理課程中的問題分析、實驗、濾波器設計及計算機模擬。附錄中所介紹的信號處理工具箱函數及繪圖函數基本可滿足本教材所要求的上機實驗需要。對序列進行譜分析的MATLAB程序及運行結果見附錄。實驗一 離散信號產生及頻譜的繪制一、實驗目的: （1）熟悉Matlab環境。（2）掌握 Matlab 中一些基本函數的建立方法（2）通過編程繪制的幅度相位譜加深理解系統的特性二、實驗內容1、編寫程序，產生以下離散序列：（1）f(n)=(n) (-3n4)n1=-3;n2=4;n0=0;n=n1:n2;x=n=n0;stem(n,x,filled);axis(n1,n2,0,1.1*max(x);xlabel(時間(n);ylabel(幅度x(n);title(單位脈沖序列);（2）f(n)=u(n) (-5n=n0;stem(n,x,filled);axis(n1,n2,0,1.1*max(x);xlabel(時間(n);ylabel(幅度x(n);title(單位階躍序列);box（3） (0n16)n1=16;a=-0.1;w=1.6*pi;n=0:n1;x=exp(a+j*w)*n);subplot(2,2,1);plot(n,real(x);title(復指數信號的實部);subplot(2,2,3);stem(n,real(x),filled);title(復指數序列的實部);subplot(2,2,2);plot(n,imag(x);title(復指數信號的虛部);subplot(2,2,4);stem(n,imag(x),filled);title(復指數序列的虛部);box %box on 加邊框/ box off不加邊框/ box 加或不加切換，注意只對上面一個圖有效（4）f(n)=3sin(n/4 (0nN error(N mustbe = the length of x) %要求移位周期大于信號長度endx=x zeros(1,N-length(x);n=0:1:N-1;n=mod(n-m,N);y=x(n+1);%例1 本例檢驗實序列的性質DFTxec(n)=ReX(k) DFTxoc(n)=ImX(k)% 設 x(n)=10*(0.8).n 0=nN error(N must be = the length of x1)endif length(x2)N error(N must be = the length of x2)endx1=x1 zeros(1,N-length(x1); %將x1,x2補0成為N長序列x2=x2 zeros(1,N-length(x2);m=0:1:N-1;x2=x2(mod(-m,N)+1); %該語句的功能是將序列翻褶，延拓，取主值序列H=zeros(N,N);for n=1:1:N %該for循環的功能是得到x2序列的循環移位矩陣 H(n,:)=cirshftt(x2,n-1,N); %和我們手工計算圓周卷積得到的表是一致的endy=x1*H %用矩陣相乘的方法得到結果% 例3 本例驗證采樣定理%令，繪制其傅立葉變換。用不同頻率對其進行采樣，分別畫出離%散時間傅立葉變換。已給出采樣頻率為時的的程序 %實驗報告要求：（1）請寫出時的程序 (2)將實驗結果畫出 (3)實驗結果說明什么（采樣頻率不同結果有何不同）。Dt=0.00005; %步長為0.00005st=-0.005:Dt:0.005; xa=exp(-1000*abs(t); %取時間從-0.005s到0.005s這段模擬信號Wmax=2*pi*2000; %信號最高頻率為2*2000K=500; %頻域正半軸取500個點進行計算k=0:1:K;W=k*Wmax/K; % 求模擬角頻率Xa=xa*exp(-j*t*W)*Dt; %計算連續時間傅立葉變換（利用矩陣運算實現） Xa=real(Xa); %取實部W=-fliplr(W),W(2:501); %將角頻率范圍擴展為從-到+Xa=fliplr(Xa),Xa(2:501); subplot(2,2,1);plot(t*1000,xa); %畫出模擬信號，橫坐標為時間（毫秒），縱坐標為幅度xlabel(time(millisecond);ylabel(xa(t); title(anolog signal);subplot(2,2,2);plot(W/(2*pi*1000),Xa*1000); %畫出連續時間傅立葉變換 xlabel(frequency(kHZ); %橫坐標為頻率（kHz）ylabel(xa(jw); %縱坐標為幅度title(FT);%下面為采樣頻率5kHz時的程序Ts=0.0002; %采樣間隔為n=-25:1:25;x=exp(-1000*abs(n*Ts); %離散時間信號K=500;k=0:1:K;w=pi*k/K; %w為數字頻率X=x*exp(-j*n*w); %計算離散時間傅立葉變換（序列的傅立葉變換）X=real(X); w=-fliplr(w),w(2:K+1);X=fliplr(X),X(2:K+1);subplot(2,2,3);stem(n*Ts*1000,x); %畫出采樣信號（離散時間信號）xlabel(time(millisecond);gtext(Ts=0.2ms); %該語句可以將引號中的內容放置在figure中的任何地方，只需 %將十字的中心放在想放置內容的地方，然后按鼠標即可。ylabel(x1(n);title(discrete signal);subplot(2,2,4);plot(w/pi,X); %畫出離散時間傅立葉變換xlabel(frequency(radian); %橫坐標為弧度ylabel(x1(jw);title(DTFT);%例4 本例說明補零序列的離散傅立葉變換%序列，已給出序列的傅立葉變換程序和將原序列補零到10長序列的DFT%實驗報告要求： （1）編寫將序列補零到20長后計算DFT的程序(2)將實驗結果畫出(3)實驗結果說明什么（即序列補零后進行DFT,頻譜有何變化）n=0:4;x=ones(1,5); %產生矩形序列k=0:999;w=(pi/500)*k;X=x*(exp(-j*pi/500).(n*k); %計算離散時間傅立葉變換Xe=abs(X); %取模subplot(3,2,1);stem(n,x);ylabel(x(n); %畫出矩形序列subplot(3,2,2);plot(w/pi,Xe);ylabel(|X(ejw)|); %畫出離散時間傅立葉變換N=10;x=ones(1,5),zeros(1,N-5); %將原序列補零為10長序列n=0:1:N-1;X=dft(x,N); %進行DFTmagX=abs(X); k=(0:length(magX)-1)*N/length(magX);subplot(3,2,3);stem(n,x);ylabel(x(n); %畫出補零序列subplot(3,2,4);stem(k,magX); %畫出DFT結果axis(0,10,0,5);ylabel(|X(k)|); %例5 本題說明高密度譜和高分辨率譜之間的區別，高密度譜是信號補零后得到的，雖然譜線相當%密但是因為信號有效長度不變，所以其分辨率也不變，因此還是很難看出信號的頻譜成分。高分辨%率譜是將信號有效長度加長，因此分辨率提高，可以看出信號的成分。%有一個序列為 （該序列周期計算可得40）%（1）下面給出有10個有效采樣點序列的DFT程序%（2）請寫出將第一問中的10長序列補零到40長，計算其DFT%（3）采樣n=0：39，計算有40個有效采樣點的序列的DFT%實驗報告要求： （1）請編寫將有10個有效采樣點的序列補零到40長后計算DFT的程序 (2) 請編寫計算有40個有效采樣點的序列的DFT的程序 (3) 將實驗結果畫出并分析實驗結果說明什么M=10;n=0:M-1;x=2*cos(0.35*pi*n)+cos(0.5*pi*n);subplot(2,1,1);stem(n,x);title(沒有足夠采樣點的信號);Y=dft(x,M);k1=0:M-1;w1=2*pi/M*k1;subplot(2,1,2);stem(w1/pi,abs(Y);title(信號的頻譜);實驗四 用FFT對信號作頻譜分析一、實驗目的（1）學習用FFT對連續信號和時域離散信號進行譜分析的方法（2）了解可能出現的分析誤差及其原因，以便正確應用FFT。二、實驗原理用FFT對信號作頻譜分析是學習數字信號處理的重要內容。經常需要進行譜分析的信號是模擬信號和時域離散信號。對信號進行譜分析的重要問題是頻譜分辨率D和分析誤差。頻譜分辨率直接和FFT的變換區間N有關，因為FFT能夠實現的頻率分辨率是，因此要求。可以根據此式選擇FFT的變換區間N。誤差主要來自于用FFT作頻譜分析時，得到的是離散譜，而信號（周期信號除外）是連續譜，只有當N較大時離散譜的包絡才能逼近于連續譜，因此N要適當選擇大一些。周期信號的頻譜是離散譜，只有用整數倍周期的長度作FFT，得到的離散譜才能代表周期信號的頻譜。如果不知道信號周期，可以盡量選擇信號的觀察時間長一些。對模擬信號進行譜分析時，首先要按照采樣定理將其變成時域離散信號。如果是模擬周期信號，也應該選取整數倍周期的長度，經過采樣后形成周期序列，按照周期序列的譜分析進行。三、實驗步驟及內容（1）對以下序列進行譜分析。 選擇FFT的變換區間N為8和16 兩種情況進行頻譜分析。分別打印其幅頻特性曲線。 并進行對比、分析和討論。（2）對以下周期序列進行譜分析。 選擇FFT的變換區間N為8和16 兩種情況分別對以上序列進行頻譜分析。分別打印其幅頻特性曲線。并進行對比、分析和討論。（3）對模擬周期信號進行譜分析 選擇采樣頻率，變換區間N=16,32,64 三種情況進行譜分析。分別打印其幅頻特性，并進行分析和討論。 四、思考題（1）對于周期序列，如果周期不知道，如何用FFT進行譜分析？（2）如何選擇FFT的變換區間？（包括非周期信號和周期信號）（3）當N=8時，和的幅頻特性會相同嗎？為什么？N=16 呢？五、實驗報告要求（1）完成各個實驗任務和要求。附上程序清單和有關曲線。（2）簡要回答思考題。六、實驗程序清單 %實驗4程序exp3.m% 用FFT對信號作頻譜分析clear all;close all%實驗內容(1)=x1n=ones(1,4); %產生序列向量x1(n)=R4(n)M=8;xa=1:(M/2); xb=(M/2):-1:1; x2n=xa,xb; %產生長度為8的三角波序列x2(n)x3n=xb,xa;X1k8=fft(x1n,8); %計算x1n的8點DFTX1k16=fft(x1n,16); %計算x1n的16點DFTX2k8=fft(x2n,8); %計算x2n的8點DFTX2k16=fft(x2n,16); %計算x2n的16點DFTX3k8=fft(x3n,8); %計算x3n的8點DFTX3k16=fft(x3n,16); %計算x3n的16點DFT%以下繪制幅頻特性曲線subplot(2,2,1);mstem(X1k8); %繪制8點DFT的幅頻特性圖title(1a) 8點DFTx_1(n);xlabel(/);ylabel(幅度);axis(0,2,0,1.2*max(abs(X1k8)subplot(2,2,3);mstem(X1k16); %繪制16點DFT的幅頻特性圖title(1b)16點DFTx_1(n);xlabel(/);ylabel(幅度);axis(0,2,0,1.2*max(abs(X1k16)figure(2)subplot(2,2,1);mstem(X2k8); %繪制8點DFT的幅頻特性圖title(2a) 8點DFTx_2(n);xlabel(/);ylabel(幅度);axis(0,2,0,1.2*max(abs(X2k8)subplot(2,2,2);mstem(X2k16); %繪制16點DFT的幅頻特性圖title(2b)16點DFTx_2(n);xlabel(/);ylabel(幅度);axis(0,2,0,1.2*max(abs(X2k16)subplot(2,2,3);mstem(X3k8); %繪制8點DFT的幅頻特性圖title(3a) 8點DFTx_3(n);xlabel(/);ylabel(幅度);axis(0,2,0,1.2*max(abs(X3k8)subplot(2,2,4);mstem(X3k16); %繪制16點DFT的幅頻特性圖title(3b)16點DFTx_3(n);xlabel(/);ylabel(幅度);axis(0,2,0,1.2*max(abs(X3k16)%實驗內容(2) 周期序列譜分析=N=8;n=0:N-1; %FFT的變換區間N=8x4n=cos(pi*n/4);x5n=cos(pi*n/4)+cos(pi*n/8);X4k8=fft(x4n,8); %計算x4n的8點DFTX5k8=fft(x5n,8); %計算x5n的8點DFTN=16;n=0:N-1; %FFT的變換區間N=16x4n=cos(pi*n/4);x5n=cos(pi*n/4)+cos(pi*n/8);X4k16=fft(x4n); %計算x4n的16點DFTX5k16=fft(x5n); %計算x5n的16點DFTfigure(3)subplot(2,2,1);mstem(X4k8); %繪制8點DFT的幅頻特性圖title(4a) 8點DFTx_4(n);xlabel(/);ylabel(幅度);axis(0,2,0,1.2*max(abs(X4k8)subplot(2,2,3);mstem(X4k16); %繪制16點DFT的幅頻特性圖title(4b)16點DFTx_4(n);xlabel(/);ylabel(幅度);axis(0,2,0,1.2*max(abs(X4k16)subplot(2,2,2);mstem(X5k8); %繪制8點DFT的幅頻特性圖title(5a) 8點DFTx_5(n);xlabel(/);ylabel(幅度);axis(0,2,0,1.2*max(abs(X5k8)subplot(2,2,4);mstem(X5k16); %繪制16點DFT的幅頻特性圖title(5b)16點DFTx_5(n);xlabel(/);ylabel(幅度);axis(0,2,0,1.2*max(abs(X5k16)%實驗內容(3) 模擬周期信號譜分析=figure(4)Fs=64;T=1/Fs;N=16;n=0:N-1; %FFT的變換區間N=16x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %對x6(t)16點采樣X6k16=fft(x6nT); %計算x6nT的16點DFTX6k16=fftshift(X6k16); %將零頻率移到頻譜中心 Tp=N*T;F=1/Tp; %頻率分辨率Fk=-N/2:N/2-1;fk=k*F; %產生16點DFT對應的采樣點頻率（以零頻率為中心）subplot(3,1,1);stem(fk,abs(X6k16),.);box on %繪制8點DFT的幅頻特性圖title(6a) 16點|DFTx_6(nT)|);xlabel(f(Hz);ylabel(幅度);axis(-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k16)N=32;n=0:N-1; %FFT的變換區間N=16x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %對x6(t)32點采樣X6k32=fft(x6nT); %計算x6nT的32點DFTX6k32=fftshift(X6k32); %將零頻率移到頻譜中心 Tp=N*T;F=1/Tp; %頻率分辨率Fk=-N/2:N/2-1;fk=k*F; %產生16點DFT對應的采樣點頻率（以零頻率為中心）subplot(3,1,2);stem(fk,abs(X6k32),.);box on %繪制8點DFT的幅頻特性圖title(6b) 32點|DFTx_6(nT)|);xlabel(f(Hz);ylabel(幅度);axis(-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k32)N=64;n=0:N-1; %FFT的變換區間N=16x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %對x6(t)64點采樣X6k64=fft(x6nT); %計算x6nT的64點DFTX6k64=fftshift(X6k64); %將零頻率移到頻譜中心 Tp=N*T;F=1/Tp; %頻率分辨率Fk=-N/2:N/2-1;fk=k*F; %產生16點DFT對應的采樣點頻率（以零頻率為中心）subplot(3,1,3);stem(fk,abs(X6k64),.); box on%繪制8點DFT的幅頻特性圖title(6a) 64點|DFTx_6(nT)|);xlabel(f(Hz);ylabel(幅度);axis(-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k64)%= mstem程序清單=function mstem(Xk)% mstem(Xk)繪制頻域采樣序列向量Xk的幅頻特性圖M=length(Xk);k=0:M-1;wk=2*k/M; %產生M點DFT對應的采樣點頻率(關于歸一化值)stem(wk,abs(Xk),.);box on %繪制M點DFT的幅頻特性圖xlabel(/);ylabel(幅度);axis(0,2,0,1.2*max(abs(Xk)七、實驗程序運行結果實驗4程序exp3.m運行結果如圖4.1所示。圖4.1程序運行結果分析討論：請讀者注意，用DFT（或FFT）分析頻譜，繪制頻譜圖時，最好將X(k)的自變量k換算成對應的頻率，作為橫坐標便于觀察頻譜。為了便于讀取頻率值，最好關于歸一化，即以作為橫坐標。1、實驗內容（1）圖（1a）和（1b）說明的8點DFT和16點DFT分別是的頻譜函數的8點和16點采樣；因為，所以，與的8點DFT的模相等，如圖（2a）和（3a）。但是，當N=16時，與不滿足循環移位關系，所以圖（2b）和（3b）的模不同。2、實驗內容（2），對周期序列譜分析的周期為8，所以N=8和N=16均是其周期的整數倍，得到正確的單一頻率正弦波的頻譜，僅在0.25處有1根單一譜線。如圖（4b）和（4b）所示。的周期為16，所以N=8不是其周期的整數倍，得到的頻譜不正確，如圖（5a）所示。N=16是其一個周期，得到正確的頻譜，僅在0.25和0.125處有2根單一譜線, 如圖（5b）所示。 3、實驗內容（3），對模擬周期信號譜分析 有3個頻率成分，。所以的周期為0.5s。 采樣頻率。變換區間N=16時，觀察時間Tp=16T=0.25s，不是的整數倍周期，所以所得頻譜不正確，如圖（6a）所示。變換區間N=32,64 時，觀察時間Tp=0.5s，1s，是的整數周期，所以所得頻譜正確，如圖（6b）和（6c）所示。圖中3根譜線正好位于處。變換區間N=64 時頻譜幅度是變換區間N=32 時2倍，這種結果正好驗證了用DFT對中期序列譜分析的理論。注意：（1）用DFT（或FFT）對模擬信號分析頻譜時，最好將X(k)的自變量k換算成對應的模擬頻率fk，作為橫坐標繪圖，便于觀察頻譜。這樣，不管變換區間N取信號周期的幾倍，畫出的頻譜圖中有效離散諧波譜線所在的頻率值不變，如圖（6b）和（6c）所示。（2）本程序直接畫出采樣序列N點DFT的模值，實際上分析頻譜時最好畫出歸一化幅度譜，這樣就避免了幅度值隨變換區間N變化的缺點。本實驗程序這樣繪圖只要是為了驗證了用DFT對中期序列譜分析的理論。實驗四 FIR數字濾波器設計與軟件實現一、實驗目的（1）掌握用窗函數法設計FIR數字濾波器的原理和方法。（2）掌握用等波紋最佳逼近法設計FIR數字濾波器的原理和方法。（3）掌握FIR濾波器的快速卷積實現原理。（4）學會調用MATLAB函數設計與實現FIR濾波器。二、實驗原理與方法（1）FIR濾波器的設計在前面的實驗中，我們介紹了IIR濾波器的設計方法并實踐了其中的雙線性變換法，IIR具有許多誘人的特性；但如此同時，也具有一些缺點。例如：若想利用快速傅立葉變換技術進行快速卷積實現濾波器，則要求單位脈沖響應是有限長的。此外，IIR濾波器的優異幅度響應，一般是以相位的非線性為代價的，非線性相位會引起頻率色散。FIR濾波器具有嚴格的相位特性，這對于語音信號處理和數據傳輸是很重要的。目前FIR濾波器的設計方法主要有三種：窗函數法、頻率取樣法和切比雪夫等波紋逼近的最優化設計方法。常用的是窗函數法和切比雪夫等波紋逼近的最優化設計方法。本實驗中的窗函數法比較簡單，可應用現成的窗函數公式，在技術指標要求不高的時候是比較靈活方便的。它是從時域出發，用一個窗函數截取理想的得到，以有限長序列近似理想的；如果從頻域出發，用理想的在單位圓上等角度取樣得到，根據得到將逼近理想的，這就是頻率取樣法。（2）窗函數設計法同其它的數字濾波器的設計方法一樣，用窗函數設計濾波器也是首先要對濾波器提出性能指標。一般是給定一個理想的頻率響應，使所設計的FIR濾波器的頻率響應去逼近所要求的理想的濾波器的響應。窗函數法設計的任務在于尋找一個可實現（有限長單位脈沖響應）的傳遞函數去逼近。一個理想的頻率響應的傅立葉反變換所得到的理想單位脈沖響應往往是一個無限長序列。對經過適當的加權、截短處理才能得到一個所需要的有限長脈沖響應序列。對應不同的加權、截短，就有不同的窗函數。所要尋找的濾波器脈沖響應就等于理想脈沖響應和窗函數的乘積，即 由此可見，窗函數的形狀就決定了濾波器的性質。例如：窗函數的主瓣寬度決定了濾波器的過渡帶寬；窗函數的旁瓣大小決定了濾波器的阻帶衰減。1、幾種常見的窗函數：（1）矩形窗（Rectangle Window）調用格式：w=boxcar(n)，根據長度n 產生一個矩形窗w。（2）三角窗（Triangular Window）調用格式：w=triang(n) ，根據長度n 產生一個三角窗w。（3）漢寧窗（Hanning Window）調用格式：w=hanning(n) ，根據長度n 產生一個漢寧窗w。（4）海明窗（Hamming Window）調用格式：w=hamming(n) ，根據長度n 產生一個海明窗w。（5）布拉克曼窗（Blackman Window）調用格式：w=blackman(n) ，根據長度n 產生一個布拉克曼窗w。（6）愷撒窗（Kaiser Window）調用格式：w=kaiser(n,beta) ，根據長度n 和影響窗函數旁瓣的參數產生一個愷撒窗w。三、實驗內容及步驟（1）認真復習第七章中用窗函數法和等波紋最佳逼近法設計FIR數字濾波器的原理；（2）調用信號產生函數xtg產生具有加性噪聲的信號xt，并自動顯示xt及其頻譜，如圖6.1所示；圖6.1 具有加性噪聲的信號x(t)及其頻譜如圖（3）請設計低通濾波器，從高頻噪聲中提取xt中的單頻調幅信號，要求信號幅頻失真小于0.1dB，將噪聲頻譜衰減60dB。先觀察xt的頻譜，確定濾波器指標參數。（4）根據濾波器指標選擇合適的窗函數，計算窗函數的長度N，調用MATLAB函數fir1設計一個FIR低通濾波器。并編寫程序，調用MATLAB快速卷積函數fftfilt實現對xt的濾波。繪圖顯示濾波器的頻響特性曲線、濾波器輸出信號的幅頻特性圖和時域波形圖。（5）重復（3），濾波器指標不變，但改用等波紋最佳逼近法，調用MATLAB函數remezord和remez設計FIR數字濾波器。并比較兩種設計方法設計的濾波器階數。提示：1）、MATLAB函數fir1和fftfilt的功能及其調用格式請查閱本書第7章；2）、采樣頻率Fs=1000Hz，采樣周期T=1/Fs；3）、根據圖6.1(b)和實驗要求，可選擇濾波器指標參數：通帶截止頻率fp=120Hz，阻帶截至頻率fs=150Hz，換算成數字頻率，通帶截止頻率，通帶最大衰為0.1dB，阻帶截至頻率，阻帶最小衰為60dB。4）、實驗程序框圖如圖6.2所示，供讀者參考。Fs=1000，T=1/Fsxt=xtg產生信號xt, 并顯示xt及其頻譜用窗函數法或等波紋最佳逼近法設計FIR濾波器hn對信號xt濾波：yt=fftfilt(hn,xt)1、計算并繪圖顯示濾波器損耗函數2、繪圖顯示濾波器輸出信號ytEnd圖6.2 實驗程序框圖四、思考題(1)如果給定通帶截止頻率和阻帶截止頻率以及阻帶最小衰減,如何用窗函數法設計線性相位低通濾波器?請寫出設計步驟。(2)如果要求用窗函數法設計帶通濾波器,且給定通帶上、下截止頻率為和，阻帶上、下截止頻率為和，試求理想帶通濾波器的截止頻率。（3）解釋為什么對同樣的技術指標，用等波紋最佳逼近法設計的濾波器階數低？五、實驗報告要求（1）對兩種設計FIR濾波器的方法（窗函數法和等波紋最佳逼近法）進行分析比較，簡述其優缺點。（2）附程序清單、打印實驗內容要求繪圖顯示的曲線圖。（3）分析總結實驗結果。（4）簡要回答思考題。六、信號產生函數xtg程序清單function xt=xtg(N)%實驗五信號x(t)產生,并顯示信號的幅頻特性曲線%xt=xtg(N) 產生一個長度為N,有加性高頻噪聲的單頻調幅信號xt,采樣頻率Fs=1000Hz%載波頻率fc=Fs/10=100Hz,調制正弦波頻率f0=fc/10=10Hz.Fs=1000;T=1/Fs;Tp=N*T;t=0:T:(N-1)*T;fc=Fs/10;f0=fc/10; %載波頻率fc=Fs/10，單頻調制信號頻率為f0=Fc/10;mt=cos(2*pi*f0*t); %產生單頻正弦波調制信號mt，頻率為f0ct=cos(2*pi*fc*t); %產生載波正弦波信號ct，頻率為fcxt=mt.*ct; %相乘產生單頻調制信號xtnt=2*rand(1,N)-1; %產生隨機噪聲nt%=設計高通濾波器hn,用于濾除噪聲nt中的低頻成分,生成高通噪聲=fp=150; fs=200;Rp=0.1;As=70;% 濾波器指標fb=fp,fs;m=0,1; % 計算remezord函數所需參數f,m,devdev=10(-As/20),(10(Rp/20)-1)/(10(Rp/20)+1);n,fo,mo,W=remezord(fb,m,dev,Fs);% 確定remez函數所需參數hn=remez(n,fo,mo,W); % 調用remez函數進行設計,用于濾除噪聲nt中的低頻成分yt=filter(hn,1,10*nt); %濾除隨機噪聲中低頻成分，生成高通噪聲yt%=xt=xt+yt; %噪聲加信號fst=fft(xt,N);k=0:N-1;f=k/Tp;su