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

1、數(shù)字信號處理實驗指導書數(shù)字信號處理上機實驗指導書陳純鍇電子與信息工程學院一、引言“數(shù)字信號處理”是一門理論和實驗密切結合的課程，為了深入地掌握課程內容，應當在學習理論的同時，做習題和上機實驗。上機實驗不僅可以幫助學生深入地理解和消化基本理論，而且能鍛煉初學者的獨立解決問題的能力。所以，根據(jù)本課程的重點要求編寫了四個實驗。第一章是全書的基礎內容，抽樣定理、時域離散系統(tǒng)的時域和頻域分析以及系統(tǒng)對輸入信號的響應是重要的基本內容。由于第一章大部分內容已經(jīng)在前期信號與系統(tǒng)課程中學習完，所以可通過實驗一幫助學生溫習以上重要內容，加深學生對“數(shù)字信號處理是通過對輸入信號的一種運算達到處理目的” 這一重要概念

2、的理解。這樣便可以使學生從信號與系統(tǒng)課程順利的過渡到本課程的學習上來。第二章、三章DFT、FFT是數(shù)字信號處理的重要數(shù)學工具，它有廣泛的使用內容。采用實驗二、實驗三加深理解DFT的基本概念、基本性質。FFT是它的快速算法，必須學會使用。數(shù)字濾波器的基本理論和設計方法是數(shù)字信號處理技術的重要內容。學習這一部分時，應重點掌握IIR和FIR兩種不同的數(shù)字濾波器的基本設計方法。IIR濾波器的單位沖激響應是無限長的，設計方法是先設計模擬濾波器，然后再通過SZ平面轉換，求出相應的數(shù)字濾波器的系統(tǒng)函數(shù)。這里的平面轉換有兩種方法，即沖激響應不變法和雙線性變換法，后者沒有頻率混疊的缺點，且轉換簡單，是一種普遍應

3、用的方法。FIR濾波器的單位沖激響應是有限長的，設計濾波器的目的即是求出符合要求的單位沖激響應。窗函數(shù)法是一種基本的，也是一種重要的設計方法。學習完第七章后可以進行實驗四。二、關于使用計算機語言由于數(shù)字信號處理實驗的主要目的是驗證數(shù)字信號處理的有關理論，進一步理解鞏固所學理論知識。所以，實現(xiàn)實驗用的算法語言可以有許多種，但為了提高實驗效率，要求學生用編程效率比C語言高好幾倍的MATLAB語言。下面介紹MATLAB的主要特點。（有關MATLAB的啟動、程序運行和有關信號處理工具箱函數(shù)等內容將放到最后附錄中介紹。）MATLAB是一種交互式的以矩陣為基本數(shù)據(jù)結構的系統(tǒng)。在生成矩陣對象時，不要求明確的

4、維數(shù)說明。所謂交互式，是指MATLAB的草稿紙編程環(huán)境。即用戶每輸入一條命令并按回車鍵，MATLAB系統(tǒng)便解釋執(zhí)行之，并顯示執(zhí)行結果。根據(jù)該結果，用戶立即知道剛輸入的命令的正確性，或利用中間結果進行其他處理等。與C語言或FORTRON語言做科學數(shù)值計算的程序設計相比較，利用MATLAB可節(jié)省大量的編程時間。將其用于數(shù)字信號處理實驗，則可大大提高實驗效率，在有限的上機時間內，實驗內容可增加幾倍。例如，C語言FFT子程序有70多行，而用MATLAB只調用一個fft函數(shù)即可實現(xiàn)對序列進行FFT計算。另外，MATLAB的工具箱及圖形顯示（打?。┕δ?，可滿足各層次人員直觀、方便的進行分析、計算和設計工作

5、，從而可大大節(jié)省時間。例如，序列的卷積、濾波，系統(tǒng)函數(shù)H(z)的幅頻特性和相頻特性等計算，均有現(xiàn)成的工具箱函數(shù)。而用其它算法語言完成這些計算的編程比較麻煩，且程序較長。由于上述特點，在美國一些大學里，MATLAB已成為輔助教學的有益工具。MATLAB已成功地用于數(shù)字信號處理課程中的問題分析、實驗、濾波器設計及計算機模擬。附錄中所介紹的信號處理工具箱函數(shù)及繪圖函數(shù)基本可滿足本教材所要求的上機實驗需要。對序列進行譜分析的MATLAB程序及運行結果見附錄。實驗一 離散信號產(chǎn)生及頻譜的繪制一、實驗目的: （1）熟悉Matlab環(huán)境。（2）掌握 Matlab 中一些基本函數(shù)的建立方法（2）通過編程繪制的

6、幅度相位譜加深理解系統(tǒng)的特性二、實驗內容1、編寫程序，產(chǎn)生以下離散序列：（1）f(n)=(n) (-3<n<4)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) (-5<n<5)n1=-5;n2=5;n0=0;n=n1:n2;x=n>=n0;stem(n,x,'filled&

7、#39;);axis(n1,n2,0,1.1*max(x);xlabel('時間(n)');ylabel('幅度x(n)');title('單位階躍序列');box（3） (0<n<16)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('復指數(shù)信號的實部');subplot(2,2,3);stem(n,real(x),'filled');title('復指數(shù)序列的實部');s

8、ubplot(2,2,2);plot(n,imag(x);title('復指數(shù)信號的虛部');subplot(2,2,4);stem(n,imag(x),'filled');title('復指數(shù)序列的虛部');box %box on 加邊框/ box off不加邊框/ box 加或不加切換，注意只對上面一個圖有效（4）f(n)=3sin(n/4) (0<n<20)f= 1/8;Um=3;nt=2; %顯示的周期數(shù)目N=32; T=1/f;dt=T/N;n=0:nt*N-1;tn=n*dt;x=Um*sin(2*f*pi*tn);sub

9、plot(2,1,1);plot(tn,x);axis(0,nt*T,1.1*min(x),1.1*max(x);ylabel('x(t)');subplot(2,1,2);stem(tn,x);axis(0,nt*T,1.1*min(x),1.1*max(x);ylabel('x(n)');box on %對當前圖形右邊及上邊加邊框（5）一個連續(xù)的周期性方波信號頻率為200Hz，信號幅度在-1+1V之間，要求在圖形窗口上顯示其兩個周期的波形。以4kHz的頻率對連續(xù)信號進行采樣，編寫程序生成連續(xù)信號和其采樣獲得的離散信號波形。f=200;nt=2; %顯示周期數(shù)

10、N=20;T=1/f;dt=T/N; %每個周期顯示20個離散值，4kHz的頻率n=0:nt*N-1;tn=n*dt;x=square(2*f*pi*tn,25); %其中25為占空比subplot(2,1,1);plot(tn,x);axis(0,nt*T,1.1*min(x),1.1*max(x);ylabel('x(t)');subplot(2,1,2);stem(tn,x);axis(0,nt*T,1.1*min(x),1.1*max(x);ylabel('x(n)');box（6）繪制的幅度譜和相位譜。clc;%本語句的作用是清除命令執(zhí)行界面中所有的輸

11、出信息clear ;%清除workspace中所有的變量clf;%清除所有的繪圖內容(如果本次程序執(zhí)行前已經(jīng)有繪圖窗口存在，則可能將本程序將要繪制的圖形繪制到之前的窗口中，可能導致疑惑)w=0:0.1:2*pi;%將連續(xù)w分成極小的間隔，%(1)取消下面一行可以做第1個幅度相位譜hjw=1+exp(-j*w)+exp(-2*j*w)+exp(-3*j*w)+exp(-4*j*w)+exp(-5*j*w)+exp(-6*j*w);mag=abs(hjw);%取模/幅度值ang=angle(hjw);%相位figure(1);%第一個窗口subplot(211)%將第一個窗口分成2行1列的子圖，分

12、的方法是從左到右，從上到下plot(w,mag);%以w為橫坐標，幅度譜為縱坐標繪制幅度譜set(gca,'ytick',-0.5:0.7:3);%設定y軸顯示范圍，從-0.5到3，每一個刻度間隔是0.4set(gca,'Xtick',0:pi/8:2*pi);%同y軸含義，此處對X軸標注set(gca,'XtickLabel','0','pi/8','pi/4','3*pi/8','pi/2','5*pi/8','3*pi/4',&#

13、39;7*pi/8','pi','9*pi/8','10*pi/8','11*pi/8','3*pi/2','13*pi/8','7*pi/4','15*pi/8','2*pi');%設定橫軸上坐標如何顯示，w1=0:2*pi/7:2*pi;hjw1=1+exp(-j*w1)+exp(-2*j*w1)+exp(-3*j*w1)+exp(-4*j*w1)+exp(-5*j*w1)+exp(-6*j*w1);mag1=abs(hjw1);hold

14、on;stem(w1,mag1,'r');title('幅度譜');%設定第1幅子圖的標題subplot(212);%第2個子圖plot(w,ang);%顯示相位%stem(w,mag)set(gca,'ytick',-pi:pi/4:pi);set(gca,'Xtick',0:pi/8:2*pi);set(gca,'XtickLabel','0','pi/8','pi/4','3*pi/8','pi/2','5*pi/8

15、9;,'3*pi/4','7*pi/8','pi','9*pi/8','10*pi/8','11*pi/8','3*pi/2','13*pi/8','7*pi/4','15*pi/8','2*pi');title('相位譜');%設定第2幅子圖的標題grid on %可以在同一個子圖里面疊繪多張圖（7）已知系統(tǒng)函數(shù)，用MATLAB繪出8階系統(tǒng)函數(shù)的零極點圖、幅頻響應和相頻響應曲線。b=1 0 0 0 0

16、0 0 0 -1; %H(z)的分子多項式系數(shù)矢量a=1; %H(z)的分母多項式系數(shù)矢量subplot(1,3,1);, zplane(b,a); %繪制H(z)的零極點圖H,w=freqz(b,a); %計算系統(tǒng)的頻率響應 subplot(1,3,2); plot(w/pi,abs(H); %繪制幅頻響應曲線axis(0,1,0,2.5);xlabel('omega/pi');ylabel('|H(ejomega)|'); subplot(1,3,3); plot(w/pi,angle(H); %繪制相頻響應曲線xlabel('omega/pi

17、9;);ylabel('phi(omega)');實驗二 離散傅立葉變換及譜分析一、 實驗目的1掌握離散傅里葉變換的計算機實現(xiàn)方法。2檢驗實序列傅里葉變換的性質。3掌握計算序列的圓周卷積的方法。4熟悉連續(xù)信號經(jīng)理想采樣前后的頻譜變化關系，加深對時域采樣定理的理解。5學習用DFT對連續(xù)信號和時域離散信號進行譜分析的方法，了解可能出現(xiàn)的分析誤差，以便在實際中正確應用DFT。二、 實驗內容1實現(xiàn)離散傅里葉變換。2計算序列圓周卷積。3計算實序列傅里葉變換并檢驗DFT性質。4實現(xiàn)連續(xù)信號傅里葉變換以及由不同采樣頻率采樣得到的離散信號的傅里葉變換。5實現(xiàn)補零序列的傅里葉變換。6實現(xiàn)高密度譜

18、和高分辨率譜，并比較二者的不同。三、 實驗報告要求 見各程序要求%以下為4個擴展函數(shù)% （1）離散傅立葉變換 采用矩陣相乘的方法function Xk=dft(xn,N)n=0:1:N-1;k=0:1:N-1;WN=exp(-j*2*pi/N);nk=n'*k;WNnk=WN.(nk);Xk=xn*WNnk;%（2）逆離散傅立葉變換 function xn=idft(Xk,N)n=0:1:N-1;k=0:1:N-1;WN=exp(-j*2*pi/N);nk=n'*k;WNnk=WN.(-nk);xn=(Xk*WNnk)/N;% (3) 實序列的分解 % 實序列可分解為共扼對稱分

19、量 % 和共扼反對稱分量 function xec,xoc=circevod(x)N=length(x);n=0:(N-1);xec=0.5*(x+x(mod(-n,N)+1); %根據(jù)上面的公式計算，mod函數(shù)為取余xoc=0.5*(x-x(mod(-n,N)+1);% (4) 序列的循環(huán)移位 function y=cirshftt(x,m,N)if length(x)>N error('N mustbe >= the length of x') %要求移位周期大于信號長度endx=x zeros(1,N-length(x);n=0:1:N-1;n=mod(n-m

20、,N);y=x(n+1);%例1 本例檢驗實序列的性質DFTxec(n)=ReX(k) DFTxoc(n)=ImX(k)% 設 x(n)=10*(0.8).n 0<=n<=10 將x(n)分解為共扼對稱及共扼反對稱部分 %實驗報告要求：（1）將實驗結果畫出 （2）實驗結果說明什么n=0:10;x=10*(0.8).n;xec,xoc=circevod(x);subplot(2,1,1);stem(n,xec); %畫出序列的共扼對稱分量title('Circular -even component') xlabel('n');ylabel('

21、xec(n)');axis(-0.5,10.5,-1,11)subplot(2,1,2);stem(n,xoc); %畫出序列的共扼反對稱分量title('Circular -odd component') xlabel('n');ylabel('xoc(n)');axis(-0.5,10.5,-4,4)figure(2)X=dft(x,11); %求出序列的DFTXec=dft(xec,11); %求序列的共扼對稱分量的DFTXoc=dft(xoc,11); %求序列的共扼反對稱分量的DFTsubplot(2,2,1);stem(n,r

22、eal(X);axis(-0.5,10.5,-5,50)title('RealDFTx(n)');xlabel('k'); %畫出序列DFT的實部subplot(2,2,2);stem(n,imag(X);axis(-0.5,10.5,-20,20)title('ImagDFTx(n)');xlabel('k'); %畫出序列DFT的虛部subplot(2,2,3);stem(n,Xec);axis(-0.5,10.5,-5,50) title('DFTxec(n)');xlabel('k');su

23、bplot(2,2,4);stem(n,imag(Xoc);axis(-0.5,10.5,-20,20)title('DFTxoc(n)');xlabel('k');% 例2 本例為計算序列的圓周卷積程序% 運行之前應在命令窗口輸入 x1,x2,N 的值%實驗報告要求：自己選擇2個序列進行計算，將實驗結果寫出。if length(x1)>N error('N must be >= the length of x1')endif length(x2)>N error('N must be >= the length

24、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' %用矩陣相乘的方法得到結果% 例3 本例驗證采樣定理%令，繪制其傅立葉變換。用不同頻率對其進行采樣，分別

25、畫出離散時間傅立葉變換。已給出采樣頻率為時的的程序 %實驗報告要求：（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; %計算連續(xù)時間傅立葉變換（利用矩陣運算實現(xiàn)

26、） Xa=real(Xa); %取實部W=-fliplr(W),W(2:501); %將角頻率范圍擴展為從-到+Xa=fliplr(Xa),Xa(2:501); %A = 1 3 5 7 9 則 fliplr(A)= 9 7 5 3 1 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

27、*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

28、),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); %畫出離散時間傅立葉變換

29、xlabel('frequency(radian)'); %橫坐標為弧度ylabel('x1(jw)');title('DTFT');%例4 本例說明補零序列的離散傅立葉變換%序列，已給出序列的傅立葉變換程序和將原序列補零到10長序列的DFT%實驗報告要求： （1）編寫將序列補零到20長后計算DFT的程序(2)給出實驗結果(3)實驗結果說明什么（即序列補零后進行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); %計算離

30、散時間傅立葉變換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);yl

31、abel('x(n)'); %畫出補零序列subplot(3,2,4);stem(k,magX); %畫出DFT結果axis(0,10,0,5);ylabel('|X(k)|'); %例5 本題說明高密度譜和高分辨率譜之間的區(qū)別，高密度譜是信號補零后得到的，雖然譜線相當密，但是因為信號有效長度不變，所以其分辨率也不變，因此還是很難看出信號的頻譜成分。高分辨率譜是將信號有效長度加長，因此分辨率提高，可以看出信號的成分。%有一個序列為 （該序列周期計算可得40）%（1）下面給出有10個有效采樣點序列的DFT程序%（2）請寫出將第一問中的10長序列補零到40長，計算其

32、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(

33、'信號的頻譜');實驗三 用FFT對信號作頻譜分析一、實驗目的（1）學習用FFT對連續(xù)信號和時域離散信號進行譜分析的方法（2）了解可能出現(xiàn)的分析誤差及其原因，以便正確應用FFT。二、實驗原理用FFT對信號作頻譜分析是學習數(shù)字信號處理的重要內容。經(jīng)常需要進行譜分析的信號是模擬信號和時域離散信號。對信號進行譜分析的重要問題是頻譜分辨率D和分析誤差。頻譜分辨率直接和FFT的變換區(qū)間N有關，因為FFT能夠實現(xiàn)的頻率分辨率是，因此要求?？梢愿鶕?jù)此式選擇FFT的變換區(qū)間N。誤差主要來自于用FFT作頻譜分析時，得到的是離散譜，而信號（周期信號除外）是連續(xù)譜，只有當N較大時離散譜的包絡才能逼近

34、于連續(xù)譜，因此N要適當選擇大一些。周期信號的頻譜是離散譜，只有用整數(shù)倍周期的長度作FFT，得到的離散譜才能代表周期信號的頻譜。如果不知道信號周期，可以盡量選擇信號的觀察時間長一些。對模擬信號進行譜分析時，首先要按照采樣定理將其變成時域離散信號。如果是模擬周期信號，也應該選取整數(shù)倍周期的長度，經(jīng)過采樣后形成周期序列，按照周期序列的譜分析進行。三、實驗步驟及內容（1）對以下序列進行譜分析。 選擇FFT的變換區(qū)間N為8和16 兩種情況進行頻譜分析。分別打印其幅頻特性曲線。 并進行對比、分析和討論。（2）對以下周期序列進行譜分析。 選擇FFT的變換區(qū)間N為8和16 兩種情況分別對以上序列進行頻譜分析。

35、分別打印其幅頻特性曲線。并進行對比、分析和討論。（3）對模擬周期信號進行譜分析 選擇采樣頻率，變換區(qū)間N=16,32,64 三種情況進行譜分析。分別打印其幅頻特性，并進行分析和討論。 四、思考題（1）對于周期序列，如果周期不知道，如何用FFT進行譜分析？（2）如何選擇FFT的變換區(qū)間？（包括非周期信號和周期信號）（3）當N=8時，和的幅頻特性會相同嗎？為什么？N=16 呢？五、實驗報告要求（1）完成各個實驗任務和要求。附上程序清單和有關曲線。（2）簡要回答思考題。六、實驗程序清單 % 用FFT對信號作頻譜分析clear all;close all%實驗內容(1)=x1n=ones(1,4);

36、%產(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);mst

37、em(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);m

38、stem(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(X3

39、k8); %繪制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;

40、%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(

41、'(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('(5

42、a) 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的變換

43、區(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對應的采樣點頻率（以零頻率為中心）subplot(3,1,1);stem(fk,abs(X6k16),'.');box on %繪制8點DFT的幅頻特性圖title('(6a) 16點|DFTx_6(nT

44、)|');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點

45、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的變換區(qū)間N=16x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %對x6(t)64點

46、采樣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對應的采樣點頻率（以零頻率為中心）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-

47、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對應的采樣點頻率(關于歸一化值)stem(wk,abs(Xk),'.');box on %繪制M點DFT的幅頻特性圖xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(Xk)七、實驗程序運行結果運行結果如圖3所示。圖3程序運行結果分析討論：請讀者注意

48、，用DFT（或FFT）分析頻譜，繪制頻譜圖時，最好將X(k)的自變量k換算成對應的頻率，作為橫坐標便于觀察頻譜。為了便于讀取頻率值，最好關于歸一化，即以作為橫坐標。1、實驗內容（1）圖（1a）和（1b）說明的8點DFT和16點DFT分別是的頻譜函數(shù)的8點和16點采樣；因為，所以，與的8點DFT的模相等，如圖（2a）和（3a）。但是，當N=16時，與不滿足循環(huán)移位關系，所以圖（2b）和（3b）的模不同。2、實驗內容（2），對周期序列譜分析的周期為8，所以N=8和N=16均是其周期的整數(shù)倍，得到正確的單一頻率正弦波的頻譜，僅在0.25處有1根單一譜線。如圖（4b）和（4b）所示。的周期為16，所以

49、N=8不是其周期的整數(shù)倍，得到的頻譜不正確，如圖（5a）所示。N=16是其一個周期，得到正確的頻譜，僅在0.25和0.125處有2根單一譜線, 如圖（5b）所示。3、實驗內容（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倍，這種結果正好驗證了用DFT對中期序列譜

50、分析的理論。注意：（1）用DFT（或FFT）對模擬信號分析頻譜時，最好將X(k)的自變量k換算成對應的模擬頻率fk，作為橫坐標繪圖，便于觀察頻譜。這樣，不管變換區(qū)間N取信號周期的幾倍，畫出的頻譜圖中有效離散諧波譜線所在的頻率值不變，如圖（6b）和（6c）所示。（2）本程序直接畫出采樣序列N點DFT的模值，實際上分析頻譜時最好畫出歸一化幅度譜，這樣就避免了幅度值隨變換區(qū)間N變化的缺點。本實驗程序這樣繪圖只要是為了驗證了用DFT對中期序列譜分析的理論。實驗四 FIR數(shù)字濾波器設計與軟件實現(xiàn)一、實驗目的（1）掌握用窗函數(shù)法設計FIR數(shù)字濾波器的原理和方法。（2）掌握用等波紋最佳逼近法設計FIR數(shù)字濾

51、波器的原理和方法。（3）掌握FIR濾波器的快速卷積實現(xiàn)原理。（4）學會調用MATLAB函數(shù)設計與實現(xiàn)FIR濾波器。二、實驗原理與方法（1）FIR濾波器的設計在前面的實驗中，我們介紹了IIR濾波器的設計方法并實踐了其中的雙線性變換法，IIR具有許多誘人的特性；但如此同時，也具有一些缺點。例如：若想利用快速傅立葉變換技術進行快速卷積實現(xiàn)濾波器，則要求單位脈沖響應是有限長的。此外，IIR濾波器的優(yōu)異幅度響應，一般是以相位的非線性為代價的，非線性相位會引起頻率色散。FIR濾波器具有嚴格的相位特性，這對于語音信號處理和數(shù)據(jù)傳輸是很重要的。目前FIR濾波器的設計方法主要有三種：窗函數(shù)法、頻率取樣法和切比雪

52、夫等波紋逼近的最優(yōu)化設計方法。常用的是窗函數(shù)法和切比雪夫等波紋逼近的最優(yōu)化設計方法。本實驗中的窗函數(shù)法比較簡單，可應用現(xiàn)成的窗函數(shù)公式，在技術指標要求不高的時候是比較靈活方便的。它是從時域出發(fā)，用一個窗函數(shù)截取理想的得到，以有限長序列近似理想的；如果從頻域出發(fā)，用理想的在單位圓上等角度取樣得到，根據(jù)得到將逼近理想的，這就是頻率取樣法。（2）窗函數(shù)設計法同其它的數(shù)字濾波器的設計方法一樣，用窗函數(shù)設計濾波器也是首先要對濾波器提出性能指標。一般是給定一個理想的頻率響應，使所設計的FIR濾波器的頻率響應去逼近所要求的理想的濾波器的響應。窗函數(shù)法設計的任務在于尋找一個可實現(xiàn)（有限長單位脈沖響應）的傳遞函

53、數(shù)去逼近。一個理想的頻率響應的傅立葉反變換所得到的理想單位脈沖響應往往是一個無限長序列。對經(jīng)過適當?shù)募訖唷⒔囟烫幚聿拍艿玫揭粋€所需要的有限長脈沖響應序列。對應不同的加權、截短，就有不同的窗函數(shù)。所要尋找的濾波器脈沖響應就等于理想脈沖響應和窗函數(shù)的乘積，即 由此可見，窗函數(shù)的形狀就決定了濾波器的性質。例如：窗函數(shù)的主瓣寬度決定了濾波器的過渡帶寬；窗函數(shù)的旁瓣大小決定了濾波器的阻帶衰減。1、幾種常見的窗函數(shù)：（1）矩形窗（Rectangle Window）調用格式：w=boxcar(n)，根據(jù)長度n 產(chǎn)生一個矩形窗w。（2）三角窗（Triangular Window）調用格式：w=triang(n

54、) ，根據(jù)長度n 產(chǎn)生一個三角窗w。（3）漢寧窗（Hanning Window）調用格式：w=hanning(n) ，根據(jù)長度n 產(chǎn)生一個漢寧窗w。（4）海明窗（Hamming Window）調用格式：w=hamming(n) ，根據(jù)長度n 產(chǎn)生一個海明窗w。（5）布拉克曼窗（Blackman Window）調用格式：w=blackman(n) ，根據(jù)長度n 產(chǎn)生一個布拉克曼窗w。（6）愷撒窗（Kaiser Window）調用格式：w=kaiser(n,beta) ，根據(jù)長度n 和影響窗函數(shù)旁瓣的參數(shù)產(chǎn)生一個愷撒窗w。三、實驗內容及步驟（1）認真復習第六章中用窗函數(shù)法和等波紋最佳逼近法設計FI

55、R數(shù)字濾波器的原理；（2）調用信號產(chǎn)生函數(shù)xtg產(chǎn)生具有加性噪聲的信號xt，并自動顯示xt及其頻譜，如圖4.1所示；圖4.1 具有加性噪聲的信號x(t)及其頻譜如圖（3）請設計低通濾波器，從高頻噪聲中提取xt中的單頻調幅信號，要求信號幅頻失真小于0.1dB，將噪聲頻譜衰減60dB。先觀察xt的頻譜，確定濾波器指標參數(shù)。（4）根據(jù)濾波器指標選擇合適的窗函數(shù)，計算窗函數(shù)的長度N，調用MATLAB函數(shù)fir1設計一個FIR低通濾波器。并編寫程序，調用MATLAB快速卷積函數(shù)fftfilt實現(xiàn)對xt的濾波。繪圖顯示濾波器的頻響特性曲線、濾波器輸出信號的幅頻特性圖和時域波形圖。（5）重復（3），濾波器指

56、標不變，但改用等波紋最佳逼近法，調用MATLAB函數(shù)remezord和remez設計FIR數(shù)字濾波器。并比較兩種設計方法設計的濾波器階數(shù)。提示：1）、MATLAB函數(shù)fir1和fftfilt的功能及其調用格式請查閱本書第6章；2）、采樣頻率Fs=1000Hz，采樣周期T=1/Fs；3）、根據(jù)圖4.1(b)和實驗要求，可選擇濾波器指標參數(shù)：通帶截止頻率fp=120Hz，阻帶截至頻率fs=150Hz，換算成數(shù)字頻率，通帶截止頻率，通帶最大衰為0.1dB，阻帶截至頻率，阻帶最小衰為60dB。4）、實驗程序框圖如圖4.2所示，供讀者參考。Fs=1000，T=1/Fsxt=xtg產(chǎn)生信號xt, 并顯示xt及其頻譜用窗函數(shù)法或等波紋最佳逼近法設計FIR濾波器hn對信號xt濾波：yt=fftfilt(hn,xt)1、計算并繪圖顯示濾波器損耗函數(shù)2、繪圖顯示濾波器輸出信號ytEnd圖6.2 實驗程序框圖四、思考題(1)如果給定通帶截止頻率和阻帶截止頻率以及阻帶最小衰減,如何用窗函數(shù)法設計線性相位低通濾波器?請寫出設計步驟。(2)如果要求用窗函數(shù)法設計帶通濾波器,且給定通帶上、下截止頻率為和，阻帶上、下截止頻率為和，試求理想帶通濾波器的截止頻率。（3）解釋為什么對同樣的技術指