數(shù)字信號處理上機實驗指導書.docx

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