數(shù)字信號(hào)處理上機(jī)實(shí)驗(yàn)指導(dǎo)書

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

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

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

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

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

6、幅度相位譜加深理解系統(tǒng)的特性二、實(shí)驗(yàn)內(nèi)容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('時(shí)間(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('時(shí)間(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('復(fù)指數(shù)信號(hào)的實(shí)部');subplot(2,2,3);stem(n,real(x),'filled');title('復(fù)指數(shù)序列的實(shí)部');s

8、ubplot(2,2,2);plot(n,imag(x);title('復(fù)指數(shù)信號(hào)的虛部');subplot(2,2,4);stem(n,imag(x),'filled');title('復(fù)指數(shù)序列的虛部');box %box on 加邊框/ box off不加邊框/ box 加或不加切換，注意只對(duì)上面一個(gè)圖有效（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 %對(duì)當(dāng)前圖形右邊及上邊加邊框（5）一個(gè)連續(xù)的周期性方波信號(hào)頻率為200Hz，信號(hào)幅度在-1+1V之間，要求在圖形窗口上顯示其兩個(gè)周期的波形。以4kHz的頻率對(duì)連續(xù)信號(hào)進(jìn)行采樣，編寫程序生成連續(xù)信號(hào)和其采樣獲得的離散信號(hào)波形。f=200;nt=2; %顯示周期數(shù)

10、N=20;T=1/f;dt=T/N; %每個(gè)周期顯示20個(gè)離散值，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;%清除所有的繪圖內(nèi)容(如果本次程序執(zhí)行前已經(jīng)有繪圖窗口存在，則可能將本程序?qū)⒁L制的圖形繪制到之前的窗口中，可能導(dǎo)致疑惑)w=0:0.1:2*pi;%將連續(xù)w分成極小的間隔，%(1)取消下面一行可以做第1個(gè)幅度相位譜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);%第一個(gè)窗口subplot(211)%將第一個(gè)窗口分成2行1列的子圖，分

12、的方法是從左到右，從上到下plot(w,mag);%以w為橫坐標(biāo)，幅度譜為縱坐標(biāo)繪制幅度譜set(gca,'ytick',-0.5:0.7:3);%設(shè)定y軸顯示范圍，從-0.5到3，每一個(gè)刻度間隔是0.4set(gca,'Xtick',0:pi/8:2*pi);%同y軸含義，此處對(duì)X軸標(biāo)注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');%設(shè)定橫軸上坐標(biāo)如何顯示，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('幅度譜');%設(shè)定第1幅子圖的標(biāo)題subplot(212);%第2個(gè)子圖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('相位譜');%設(shè)定第2幅子圖的標(biāo)題grid on %可以在同一個(gè)子圖里面疊繪多張圖（7）已知系統(tǒng)函數(shù)，用MATLAB繪出8階系統(tǒng)函數(shù)的零極點(diǎn)圖、幅頻響應(yīng)和相頻響應(yīng)曲線。b=1 0 0 0 0

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

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

18、和高分辨率譜，并比較二者的不同。三、 實(shí)驗(yàn)報(bào)告要求 見各程序要求%以下為4個(gè)擴(kuò)展函數(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) 實(shí)序列的分解 % 實(shí)序列可分解為共扼對(duì)稱分

19、量 % 和共扼反對(duì)稱分量 function xec,xoc=circevod(x)N=length(x);n=0:(N-1);xec=0.5*(x+x(mod(-n,N)+1); %根據(jù)上面的公式計(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') %要求移位周期大于信號(hào)長(zhǎng)度endx=x zeros(1,N-length(x);n=0:1:N-1;n=mod(n-m

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

21、xec(n)');axis(-0.5,10.5,-1,11)subplot(2,1,2);stem(n,xoc); %畫出序列的共扼反對(duì)稱分量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); %求序列的共扼對(duì)稱分量的DFTXoc=dft(xoc,11); %求序列的共扼反對(duì)稱分量的DFTsubplot(2,2,1);stem(n,r

22、eal(X);axis(-0.5,10.5,-5,50)title('RealDFTx(n)');xlabel('k'); %畫出序列DFT的實(shí)部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 本例為計(jì)算序列的圓周卷積程序% 運(yùn)行之前應(yīng)在命令窗口輸入 x1,x2,N 的值%實(shí)驗(yàn)報(bào)告要求：自己選擇2個(gè)序列進(jìn)行計(jì)算，將實(shí)驗(yàn)結(jié)果寫出。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補(bǔ)0成為N長(zhǎng)序列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); %和我們手工計(jì)算圓周卷積得到的表是一致的endy=x1*H' %用矩陣相乘的方法得到結(jié)果% 例3 本例驗(yàn)證采樣定理%令，繪制其傅立葉變換。用不同頻率對(duì)其進(jìn)行采樣，分別

25、畫出離散時(shí)間傅立葉變換。已給出采樣頻率為時(shí)的的程序 %實(shí)驗(yàn)報(bào)告要求：（1）請(qǐng)寫出時(shí)的程序 (2)將實(shí)驗(yàn)結(jié)果畫出 (3)實(shí)驗(yàn)結(jié)果說明什么（采樣頻率不同結(jié)果有何不同）。Dt=0.00005; %步長(zhǎng)為0.00005st=-0.005:Dt:0.005; xa=exp(-1000*abs(t); %取時(shí)間從-0.005s到0.005s這段模擬信號(hào)Wmax=2*pi*2000; %信號(hào)最高頻率為2*2000K=500; %頻域正半軸取500個(gè)點(diǎn)進(jìn)行計(jì)算k=0:1:K;W=k*Wmax/K; % 求模擬角頻率Xa=xa*exp(-j*t'*W)*Dt; %計(jì)算連續(xù)時(shí)間傅立葉變換（利用矩陣運(yùn)算實(shí)現(xiàn)

26、） Xa=real(Xa); %取實(shí)部W=-fliplr(W),W(2:501); %將角頻率范圍擴(kuò)展為從-到+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); %畫出模擬信號(hào)，橫坐標(biāo)為時(shí)間（毫秒），縱坐標(biāo)為幅度xlabel('time(millisecond)');ylabel('xa(t)'); title('anolog signal');subplot(2,2,2);plot(W/(2*pi

27、*1000),Xa*1000); %畫出連續(xù)時(shí)間傅立葉變換 xlabel('frequency(kHZ)'); %橫坐標(biāo)為頻率（kHz）ylabel('xa(jw)'); %縱坐標(biāo)為幅度title('FT');%下面為采樣頻率5kHz時(shí)的程序Ts=0.0002; %采樣間隔為n=-25:1:25;x=exp(-1000*abs(n*Ts); %離散時(shí)間信號(hào)K=500;k=0:1:K;w=pi*k/K; %w為數(shù)字頻率X=x*exp(-j*n'*w); %計(jì)算離散時(shí)間傅立葉變換（序列的傅立葉變換）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); %畫出采樣信號(hào)（離散時(shí)間信號(hào)）xlabel('time(millisecond)');gtext('Ts=0.2ms'); %該語句可以將引號(hào)中的內(nèi)容放置在figure中的任何地方，只需 %將十字的中心放在想放置內(nèi)容的地方，然后按鼠標(biāo)即可。ylabel('x1(n)');title('discrete signal');subplot(2,2,4);plot(w/pi,X); %畫出離散時(shí)間傅立葉變換

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

30、散時(shí)間傅立葉變換Xe=abs(X); %取模subplot(3,2,1);stem(n,x);ylabel('x(n)'); %畫出矩形序列subplot(3,2,2);plot(w/pi,Xe);ylabel('|X(ejw)|'); %畫出離散時(shí)間傅立葉變換N=10;x=ones(1,5),zeros(1,N-5); %將原序列補(bǔ)零為10長(zhǎng)序列n=0:1:N-1;X=dft(x,N); %進(jìn)行DFTmagX=abs(X); k=(0:length(magX)'-1)*N/length(magX);subplot(3,2,3);stem(n,x);yl

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

32、DFT%（3）采樣n=0：39，計(jì)算有40個(gè)有效采樣點(diǎn)的序列的DFT%實(shí)驗(yàn)報(bào)告要求： （1）請(qǐng)編寫將有10個(gè)有效采樣點(diǎn)的序列補(bǔ)零到40長(zhǎng)后計(jì)算DFT的程序 (2) 請(qǐng)編寫計(jì)算有40個(gè)有效采樣點(diǎn)的序列的DFT的程序 (3) 將實(shí)驗(yàn)結(jié)果畫出并分析實(shí)驗(yàn)結(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('沒有足夠采樣點(diǎn)的信號(hào)');Y=dft(x,M);k1=0:M-1;w1=2*pi/M*k1;subplot(2,1,2);stem(w1/pi,abs(Y);title(

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

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

35、分別打印其幅頻特性曲線。并進(jìn)行對(duì)比、分析和討論。（3）對(duì)模擬周期信號(hào)進(jìn)行譜分析 選擇采樣頻率，變換區(qū)間N=16,32,64 三種情況進(jìn)行譜分析。分別打印其幅頻特性，并進(jìn)行分析和討論。 四、思考題（1）對(duì)于周期序列，如果周期不知道，如何用FFT進(jìn)行譜分析？（2）如何選擇FFT的變換區(qū)間？（包括非周期信號(hào)和周期信號(hào)）（3）當(dāng)N=8時(shí)，和的幅頻特性會(huì)相同嗎？為什么？N=16 呢？五、實(shí)驗(yàn)報(bào)告要求（1）完成各個(gè)實(shí)驗(yàn)任務(wù)和要求。附上程序清單和有關(guān)曲線。（2）簡(jiǎn)要回答思考題。六、實(shí)驗(yàn)程序清單 % 用FFT對(duì)信號(hào)作頻譜分析clear all;close all%實(shí)驗(yàn)內(nèi)容(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)生長(zhǎng)度為8的三角波序列x2(n)x3n=xb,xa;X1k8=fft(x1n,8); %計(jì)算x1n的8點(diǎn)DFTX1k16=fft(x1n,16); %計(jì)算x1n的16點(diǎn)DFTX2k8=fft(x2n,8); %計(jì)算x2n的8點(diǎn)DFTX2k16=fft(x2n,16); %計(jì)算x2n的16點(diǎn)DFTX3k8=fft(x3n,8); %計(jì)算x3n的8點(diǎn)DFTX3k16=fft(x3n,16); %計(jì)算x3n的16點(diǎn)DFT%以下繪制幅頻特性曲線subplot(2,2,1);mst

37、em(X1k8); %繪制8點(diǎn)DFT的幅頻特性圖title('(1a) 8點(diǎn)DFTx_1(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X1k8)subplot(2,2,3);mstem(X1k16); %繪制16點(diǎn)DFT的幅頻特性圖title('(1b)16點(diǎn)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點(diǎn)DFT的幅頻特性圖title('(2a) 8點(diǎn)DFTx_2(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X2k8)subplot(2,2,2);mstem(X2k16); %繪制16點(diǎn)DFT的幅頻特性圖title('(2b)16點(diǎn)DFTx_2(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X2k16)subplot(2,2,3);mstem(X3

39、k8); %繪制8點(diǎn)DFT的幅頻特性圖title('(3a) 8點(diǎn)DFTx_3(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X3k8)subplot(2,2,4);mstem(X3k16); %繪制16點(diǎn)DFT的幅頻特性圖title('(3b)16點(diǎn)DFTx_3(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X3k16)%實(shí)驗(yàn)內(nèi)容(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); %計(jì)算x4n的8點(diǎn)DFTX5k8=fft(x5n,8); %計(jì)算x5n的8點(diǎn)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); %計(jì)算x4n的16點(diǎn)DFTX5k16=fft(x5n); %計(jì)算x5n的16點(diǎn)DFTfigure(3)subplot(2,2,1);mstem(X4k8); %繪制8點(diǎn)DFT的幅頻特性圖title(

41、'(4a) 8點(diǎn)DFTx_4(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X4k8)subplot(2,2,3);mstem(X4k16); %繪制16點(diǎn)DFT的幅頻特性圖title('(4b)16點(diǎn)DFTx_4(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X4k16)subplot(2,2,2);mstem(X5k8); %繪制8點(diǎn)DFT的幅頻特性圖title('(5

42、a) 8點(diǎn)DFTx_5(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X5k8)subplot(2,2,4);mstem(X5k16); %繪制16點(diǎn)DFT的幅頻特性圖title('(5b)16點(diǎn)DFTx_5(n)');xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(X5k16)%實(shí)驗(yàn)內(nèi)容(3) 模擬周期信號(hào)譜分析=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); %對(duì)x6(t)16點(diǎn)采樣X6k16=fft(x6nT); %計(jì)算x6nT的16點(diǎn)DFTX6k16=fftshift(X6k16); %將零頻率移到頻譜中心 Tp=N*T;F=1/Tp; %頻率分辨率Fk=-N/2:N/2-1;fk=k*F; %產(chǎn)生16點(diǎn)DFT對(duì)應(yīng)的采樣點(diǎn)頻率（以零頻率為中心）subplot(3,1,1);stem(fk,abs(X6k16),'.');box on %繪制8點(diǎn)DFT的幅頻特性圖title('(6a) 16點(diǎn)|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); %對(duì)x6(t)32點(diǎn)采樣X6k32=fft(x6nT); %計(jì)算x6nT的32點(diǎn)DFTX6k32=fftshift(X6k32); %將零頻率移到頻譜中心 Tp=N*T;F=1/Tp; %頻率分辨率Fk=-N/2:N/2-1;fk=k*F; %產(chǎn)生16點(diǎn)

45、DFT對(duì)應(yīng)的采樣點(diǎn)頻率（以零頻率為中心）subplot(3,1,2);stem(fk,abs(X6k32),'.');box on %繪制8點(diǎn)DFT的幅頻特性圖title('(6b) 32點(diǎn)|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); %對(duì)x6(t)64點(diǎn)

46、采樣X6k64=fft(x6nT); %計(jì)算x6nT的64點(diǎn)DFTX6k64=fftshift(X6k64); %將零頻率移到頻譜中心 Tp=N*T;F=1/Tp; %頻率分辨率Fk=-N/2:N/2-1;fk=k*F; %產(chǎn)生16點(diǎn)DFT對(duì)應(yīng)的采樣點(diǎn)頻率（以零頻率為中心）subplot(3,1,3);stem(fk,abs(X6k64),'.'); box on%繪制8點(diǎn)DFT的幅頻特性圖title('(6a) 64點(diǎn)|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點(diǎn)DFT對(duì)應(yīng)的采樣點(diǎn)頻率(關(guān)于歸一化值)stem(wk,abs(Xk),'.');box on %繪制M點(diǎn)DFT的幅頻特性圖xlabel('/');ylabel('幅度');axis(0,2,0,1.2*max(abs(Xk)七、實(shí)驗(yàn)程序運(yùn)行結(jié)果運(yùn)行結(jié)果如圖3所示。圖3程序運(yùn)行結(jié)果分析討論：請(qǐng)讀者注意

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

49、N=8不是其周期的整數(shù)倍，得到的頻譜不正確，如圖（5a）所示。N=16是其一個(gè)周期，得到正確的頻譜，僅在0.25和0.125處有2根單一譜線, 如圖（5b）所示。3、實(shí)驗(yàn)內(nèi)容（3），對(duì)模擬周期信號(hào)譜分析有3個(gè)頻率成分，。所以的周期為0.5s。 采樣頻率。變換區(qū)間N=16時(shí)，觀察時(shí)間Tp=16T=0.25s，不是的整數(shù)倍周期，所以所得頻譜不正確，如圖（6a）所示。變換區(qū)間N=32,64 時(shí)，觀察時(shí)間Tp=0.5s，1s，是的整數(shù)周期，所以所得頻譜正確，如圖（6b）和（6c）所示。圖中3根譜線正好位于處。變換區(qū)間N=64 時(shí)頻譜幅度是變換區(qū)間N=32 時(shí)2倍，這種結(jié)果正好驗(yàn)證了用DFT對(duì)中期序列譜

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

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

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

53、數(shù)去逼近。一個(gè)理想的頻率響應(yīng)的傅立葉反變換所得到的理想單位脈沖響應(yīng)往往是一個(gè)無限長(zhǎng)序列。對(duì)經(jīng)過適當(dāng)?shù)募訖?quán)、截短處理才能得到一個(gè)所需要的有限長(zhǎng)脈沖響應(yīng)序列。對(duì)應(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ù)長(zhǎng)度n 產(chǎn)生一個(gè)矩形窗w。（2）三角窗（Triangular Window）調(diào)用格式：w=triang(n

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

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

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