數(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í)驗(yàn)一離散信號(hào)產(chǎn)生及頻譜的繪制4實(shí)驗(yàn)二時(shí)域采樣與頻域采樣9實(shí)驗(yàn)三離散傅立葉變換及譜分析16實(shí)驗(yàn)四用FFT對(duì)信號(hào)作頻譜分析21實(shí)驗(yàn)五雙線性變換法設(shè)計(jì)IIR數(shù)字濾波器32實(shí)驗(yàn)六 FIR數(shù)字濾波器設(shè)計(jì)與軟件實(shí)現(xiàn)38附錄：實(shí)驗(yàn)用MATLAB語言工具箱函數(shù)簡介45引言“數(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)容，抽樣定

2、理、時(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á)到處理目的” 這一重要概念的理解。這樣便可以使學(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ù)字濾

3、波器的基本設(shè)計(jì)方法。IIR濾波器的單位沖激響應(yīng)是無限長的，設(shè)計(jì)方法是先設(shè)計(jì)模擬濾波器，然后再通過SZ平面轉(zhuǎn)換，求出相應(yīng)的數(shù)字濾波器的系統(tǒng)函數(shù)。這里的平面轉(zhuǎn)換有兩種方法，即沖激響應(yīng)不變法和雙線性變換法，后者沒有頻率混疊的缺點(diǎn)，且轉(zhuǎn)換簡單，是一種普遍應(yīng)用的方法。FIR濾波器的單位沖激響應(yīng)是有限長的，設(shè)計(jì)濾波器的目的即是求出符合要求的單位沖激響應(yīng)。窗函數(shù)法是一種基本的，也是一種重要的設(shè)計(jì)方法。學(xué)習(xí)完第七章后可以進(jìn)行實(shí)驗(yàn)四。由于數(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語言高好幾

4、倍的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í)，不要求明確的維數(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)

5、，實(shí)驗(yàn)內(nèi)容可增加幾倍。例如，C語言FFT子程序有70多行，而用MATLAB只調(diào)用一個(gè)fft函數(shù)即可實(shí)現(xiàn)對(duì)序列進(jìn)行FFT計(jì)算。另外，MATLAB的工具箱及圖形顯示（打?。┕δ埽蓾M足各層次人員直觀、方便的進(jìn)行分析、計(jì)算和設(shè)計(jì)工作，從而可大大節(jié)省時(shí)間。例如，序列的卷積、濾波，系統(tǒng)函數(shù)H(z)的幅頻特性和相頻特性等計(jì)算，均有現(xiàn)成的工具箱函數(shù)。而用其它算法語言完成這些計(jì)算的編程比較麻煩，且程序較長。由于上述特點(diǎn)，在美國一些大學(xué)里，MATLAB已成為輔助教學(xué)的有益工具。MATLAB已成功地用于數(shù)字信號(hào)處理課程中的問題分析、實(shí)驗(yàn)、濾波器設(shè)計(jì)及計(jì)算機(jī)模擬。附錄中所介紹的信號(hào)處理工具箱函數(shù)及繪圖函數(shù)基本可滿足

6、本教材所要求的上機(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）通過編程繪制的幅度相位譜加深理解系統(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)&#

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

8、(n,real(x);title('復(fù)指數(shù)信號(hào)的實(shí)部');subplot(2,2,3);stem(n,real(x),'filled');title('復(fù)指數(shù)序列的實(shí)部');subplot(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)

9、=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);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 on %對(duì)當(dāng)前圖形右邊及上邊加邊框（5）一個(gè)

10、連續(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ù)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)

11、;stem(tn,x);axis(0,nt*T,1.1*min(x),1.1*max(x);ylabel('x(n)');box（6）繪制的幅度譜和相位譜。clc;%本語句的作用是清除命令執(zhí)行界面中所有的輸出信息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*

12、w)+exp(-5*j*w)+exp(-6*j*w);mag=abs(hjw);%取模/幅度值ang=angle(hjw);%相位figure(1);%第一個(gè)窗口subplot(211)%將第一個(gè)窗口分成2行1列的子圖，分的方法是從左到右，從上到下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'

13、,'0','pi/8','pi/4','3*pi/8','pi/2','5*pi/8','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');%設(shè)定橫軸上坐標(biāo)如何顯示，w1=0:2*pi

14、/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 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(gc

15、a,'XtickLabel','0','pi/8','pi/4','3*pi/8','pi/2','5*pi/8','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');

16、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 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

17、/pi');ylabel('|H(ejomega)|'); subplot(1,3,3); plot(w/pi,angle(H); %繪制相頻響應(yīng)曲線xlabel('omega/pi');ylabel('phi(omega)');實(shí)驗(yàn)二時(shí)域采樣與頻域采樣一、實(shí)驗(yàn)?zāi)康模?）掌握模擬信號(hào)采樣前后頻譜的變化（2）要求掌握頻率域采樣會(huì)引起時(shí)域周期化的概念（3）頻率域采樣定理及其對(duì)頻域采樣點(diǎn)數(shù)選擇的指導(dǎo)作用。二、實(shí)驗(yàn)原理與方法時(shí)域采樣定理的要點(diǎn)是：a) 對(duì)模擬信號(hào)以間隔T進(jìn)行時(shí)域等間隔理想采樣，形成的采樣信號(hào)的頻譜是原模擬信號(hào)頻譜以采樣角頻率（）為

18、周期進(jìn)行周期延拓。公式為：b) 采樣頻率必須大于等于模擬信號(hào)最高頻率的兩倍以上，才能使采樣信號(hào)的頻譜不產(chǎn)生頻譜混疊。利用計(jì)算機(jī)計(jì)算上式并不方便，下面我們導(dǎo)出另外一個(gè)公式，以便用計(jì)算機(jī)上進(jìn)行實(shí)驗(yàn)。理想采樣信號(hào)和模擬信號(hào)之間的關(guān)系為：對(duì)上式進(jìn)行傅立葉變換，得到：在上式的積分號(hào)內(nèi)只有當(dāng)時(shí)，才有非零值，因此：上式中，在數(shù)值上，再將代入，得到：上式的右邊就是序列的傅立葉變換，即上式說明理想采樣信號(hào)的傅立葉變換可用相應(yīng)的采樣序列的傅立葉變換得到，只要將自變量用代替即可。頻域采樣定理的要點(diǎn)是：a) 對(duì)信號(hào)x(n)的頻譜函數(shù)X(ej)在0，2上等間隔采樣N點(diǎn)，得到則N點(diǎn)IDFT得到的序列就是原序列x(n)以N

19、為周期進(jìn)行周期延拓后的主值區(qū)序列，公式為：b) 由上式可知，頻域采樣點(diǎn)數(shù)N必須大于等于時(shí)域離散信號(hào)的長度M(即NM)，才能使時(shí)域不產(chǎn)生混疊，則N點(diǎn)IDFT得到的序列就是原序列x(n),即=x(n)。如果N>M，比原序列尾部多N-M個(gè)零點(diǎn)；如果N<M，z則=IDFT發(fā)生了時(shí)域混疊失真，而且的長度N也比x(n)的長度M短，因此。與x(n)不相同。在數(shù)字信號(hào)處理的應(yīng)用中，只要涉及時(shí)域或者頻域采樣，都必須服從這兩個(gè)采樣理論的要點(diǎn)。 對(duì)比上面敘述的時(shí)域采樣原理和頻域采樣原理，得到一個(gè)有用的結(jié)論，這兩個(gè)采樣理論具有對(duì)偶性：“時(shí)域采樣頻譜周期延拓，頻域采樣時(shí)域信號(hào)周期延拓”。因此放在一起進(jìn)行實(shí)驗(yàn)

20、。二、實(shí)驗(yàn)內(nèi)容及步驟（1）時(shí)域采樣理論的驗(yàn)證。給定模擬信號(hào)，式中A=444.128，=50，=50rad/s，它的幅頻特性曲線如圖3.1圖2.1 的幅頻特性曲線現(xiàn)用DFT(FFT)求該模擬信號(hào)的幅頻特性，以驗(yàn)證時(shí)域采樣理論。按照的幅頻特性曲線，選取三種采樣頻率，即=1kHz，300Hz，200Hz。觀測時(shí)間選。為使用DFT，首先用下面公式產(chǎn)生時(shí)域離散信號(hào)，對(duì)三種采樣頻率，采樣序列按順序用，表示。因?yàn)椴蓸宇l率不同，得到的，的長度不同， 長度（點(diǎn)數(shù)）用公式計(jì)算。選FFT的變換點(diǎn)數(shù)為M=64，序列長度不夠64的尾部加零。X(k)=FFTx(n) ， k=0,1,2,3,-,M-1式中k代表的頻率為

21、。要求： 編寫實(shí)驗(yàn)程序，計(jì)算、和的幅度特性，并繪圖顯示。觀察分析頻譜混疊失真。（2）頻域采樣理論的驗(yàn)證。給定信號(hào)如下：編寫程序分別對(duì)頻譜函數(shù)在區(qū)間上等間隔采樣32和16點(diǎn)，得到：再分別對(duì)進(jìn)行32點(diǎn)和16點(diǎn)IFFT，得到：分別畫出、的幅度譜，并繪圖顯示x(n)、的波形，進(jìn)行對(duì)比和分析，驗(yàn)證總結(jié)頻域采樣理論。提示：頻域采樣用以下方法容易變程序?qū)崿F(xiàn)。 直接調(diào)用MATLAB函數(shù)fft計(jì)算就得到在的32點(diǎn)頻率域采樣 抽取的偶數(shù)點(diǎn)即可得到在的16點(diǎn)頻率域采樣，即。 當(dāng)然也可以按照頻域采樣理論，先將信號(hào)x(n)以16為周期進(jìn)行周期延拓，取其主值區(qū)（16點(diǎn)），再對(duì)其進(jìn)行16點(diǎn)DFT(FFT),得到的就是在的1

22、6點(diǎn)頻率域采樣。四、思考題： 如果序列x(n)的長度為M，希望得到其頻譜在上的N點(diǎn)等間隔采樣，當(dāng)N<M時(shí)， 如何用一次最少點(diǎn)數(shù)的DFT得到該頻譜采樣？五、實(shí)驗(yàn)報(bào)告及要求a）運(yùn)行程序打印要求顯示的圖形，。b) 分析比較實(shí)驗(yàn)結(jié)果，簡述由實(shí)驗(yàn)得到的主要結(jié)論c) 簡要回答思考題d) 附上程序清單和有關(guān)曲線。六、實(shí)驗(yàn)程序清單1 時(shí)域采樣理論的驗(yàn)證程序清單% 時(shí)域采樣理論驗(yàn)證程序exp2a.mTp=64/1000;%觀察時(shí)間Tp=64微秒%產(chǎn)生M長采樣序列x(n)% Fs=1000;T=1/Fs;Fs=1000;T=1/Fs;M=Tp*Fs;n=0:M-1;A=444.128;alph=pi*50*

23、20.5;omega=pi*50*20.5;xnt=A*exp(-alph*n*T).*sin(omega*n*T);Xk=T*fft(xnt,M); %M點(diǎn)FFTxnt)yn='xa(nT)'subplot(3,2,1);tstem(xnt,yn);%調(diào)用自編繪圖函數(shù)tstem繪制序列圖box on;title('(a) Fs=1000Hz');k=0:M-1;fk=k/Tp;subplot(3,2,2);plot(fk,abs(Xk);title('(a) T*FTxa(nT),Fs=1000Hz');xlabel('f(Hz)

24、9;);ylabel('幅度');axis(0,Fs,0,1.2*max(abs(Xk)%=% Fs=300Hz和 Fs=200Hz的程序與上面Fs=1000Hz完全相同。2 頻域采樣理論的驗(yàn)證程序清單%頻域采樣理論驗(yàn)證程序exp2b.mM=27;N=32;n=0:M;%產(chǎn)生M長三角波序列x(n)xa=0:floor(M/2); xb= ceil(M/2)-1:-1:0; xn=xa,xb;Xk=fft(xn,1024);%1024點(diǎn)FFTx(n), 用于近似序列x(n)的TFX32k=fft(xn,32);%32點(diǎn)FFTx(n)x32n=ifft(X32k);%32點(diǎn)IFFT

25、X32(k)得到x32(n)X16k=X32k(1:2:N);%隔點(diǎn)抽取X32k得到X16(K)x16n=ifft(X16k,N/2);%16點(diǎn)IFFTX16(k)得到x16(n)subplot(3,2,2);stem(n,xn,'.');box ontitle('(b) 三角波序列x(n)');xlabel('n');ylabel('x(n)');axis(0,32,0,20)k=0:1023;wk=2*k/1024;%subplot(3,2,1);plot(wk,abs(Xk);title('(a)FTx(n)'

26、;);xlabel('omega/pi');ylabel('|X(ejomega)|');axis(0,1,0,200)k=0:N/2-1;subplot(3,2,3);stem(k,abs(X16k),'.');box ontitle('(c) 16點(diǎn)頻域采樣');xlabel('k');ylabel('|X_1_6(k)|');axis(0,8,0,200)n1=0:N/2-1;subplot(3,2,4);stem(n1,x16n,'.');box ontitle('(

27、d) 16點(diǎn)IDFTX_1_6(k)');xlabel('n');ylabel('x_1_6(n)');axis(0,32,0,20)k=0:N-1;subplot(3,2,5);stem(k,abs(X32k),'.');box ontitle('(e) 32點(diǎn)頻域采樣');xlabel('k');ylabel('|X_3_2(k)|');axis(0,16,0,200)n1=0:N-1;subplot(3,2,6);stem(n1,x32n,'.');box ontitl

28、e('(f) 32點(diǎn)IDFTX_3_2(k)');xlabel('n');ylabel('x_3_2(n)');axis(0,32,0,20)%=tstem 程序清單=function tstem(xn,yn)%時(shí)域序列繪圖函數(shù)% xn:信號(hào)數(shù)據(jù)序列，yn:繪圖信號(hào)的縱坐標(biāo)名稱（字符串）n=0:length(xn)-1;stem(n,xn,'.');box onxlabel('n');ylabel(yn);axis(0,n(end),min(xn),1.2*max(xn)七、實(shí)驗(yàn)程序運(yùn)行結(jié)果1 時(shí)域采樣理論的驗(yàn)證程

29、序運(yùn)行結(jié)果exp2a.m如圖3.2所示。由圖可見，采樣序列的頻譜的確是以采樣頻率為周期對(duì)模擬信號(hào)頻譜的周期延拓。當(dāng)采樣頻率為1000Hz時(shí)頻譜混疊很??；當(dāng)采樣頻率為300Hz時(shí)，在折疊頻率150Hz附近頻譜混疊很嚴(yán)重；當(dāng)采樣頻率為200Hz時(shí)，在折疊頻率110Hz附近頻譜混疊更很嚴(yán)重。圖2.22 時(shí)域采樣理論的驗(yàn)證程序exp2b.m運(yùn)行結(jié)果如圖3.3所示。圖2.3該圖驗(yàn)證了頻域采樣理論和頻域采樣定理。對(duì)信號(hào)x(n)的頻譜函數(shù)X(ej)在0，2上等間隔采樣N=16時(shí)， N點(diǎn)IDFT得到的序列正是原序列x(n)以16為周期進(jìn)行周期延拓后的主值區(qū)序列：由于N<M，所以發(fā)生了時(shí)域混疊失真，因此。

30、與x(n)不相同，如圖圖3.3(c)和(d)所示。當(dāng)N=32時(shí)，如圖圖3.3(c)和(d)所示，由于N>M，頻域采樣定理，所以不存在時(shí)域混疊失真，因此與x(n)相同。實(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ù)信

31、號(hào)傅里葉變換以及由不同采樣頻率采樣得到的離散信號(hào)的傅里葉變換。5實(shí)現(xiàn)補(bǔ)零序列的傅里葉變換。6實(shí)現(xiàn)高密度譜和高分辨率譜，并比較二者的不同。三、 實(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=W

32、N.(-nk);xn=(Xk*WNnk)/N;% (3) 實(shí)序列的分解 % 實(shí)序列可分解為共扼對(duì)稱分量% 和共扼反對(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)

33、長度endx=x zeros(1,N-length(x);n=0:1:N-1;n=mod(n-m,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 com

34、ponent')xlabel('n');ylabel('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,

35、11); %求序列的共扼反對(duì)稱分量的DFTsubplot(2,2,1);stem(n,real(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(

36、'DFTxec(n)');xlabel('k');subplot(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)>Nerror('N mustbe >= the length of x1')endif length(x2)>

37、;Nerror('N mustbe >= the length of x2')endx1=x1 zeros(1,N-length(x1); %將x1,x2補(bǔ)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); %和我們手工計(jì)算圓周卷積得到的表是一致的endy=x1*H' %用矩陣相乘的方法得到結(jié)果

38、% 例3 本例驗(yàn)證采樣定理%令，繪制其傅立葉變換。用不同頻率對(duì)其進(jìn)行采樣，分別畫出離散時(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; %步長為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(

39、-j*t'*W)*Dt; %計(jì)算連續(xù)時(shí)間傅立葉變換（利用矩陣運(yùn)算實(shí)現(xiàn)） 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 signa

40、l');subplot(2,2,2);plot(W/(2*pi*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í)間傅立葉

41、變換（序列的傅立葉變換）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); %畫出采樣信號(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(

42、2,2,4);plot(w/pi,X); %畫出離散時(shí)間傅立葉變換xlabel('frequency(radian)'); %橫坐標(biāo)為弧度ylabel('x1(jw)');title('DTFT');%例4 本例說明補(bǔ)零序列的離散傅立葉變換%序列，已給出序列的傅立葉變換程序和將原序列補(bǔ)零到10長序列的DFT%實(shí)驗(yàn)報(bào)告要求： （1）編寫將序列補(bǔ)零到20長后計(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

43、*(exp(-j*pi/500).(n'*k); %計(jì)算離散時(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長序列n=0:1:N-1;X=dft(x,N); %進(jìn)行DFTmagX=abs(X);k=(0:length(magX)'-1)*N/length(m

44、agX);subplot(3,2,3);stem(n,x);ylabel('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)有效長度不變，所以其分辨率也不變，因此還是很難看出信號(hào)的頻譜成分。高分辨率譜是將信號(hào)有效長度加長，因此分辨率提高，可以看出信號(hào)的成分。%有一個(gè)序列為 （該序列周期計(jì)算可得40）%（1）下面給出有10個(gè)有效采樣點(diǎn)序列的DF

45、T程序%（2）請(qǐng)寫出將第一問中的10長序列補(bǔ)零到40長，計(jì)算其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長后計(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(

46、2,1,2);stem(w1/pi,abs(Y);title('信號(hào)的頻譜');M=10;n=0:M-1;b=40;s=0:b-1;x=2*cos(0.35*pi*n)+cos(0.5*pi*n), zeros(1,30);subplot(2,1,1);stem(s,x);title('有足夠采樣點(diǎn)的信號(hào)');Y=dft(x,b);k1=0:b-1;w1=2*pi/b*k1;subplot(2,1,2);stem(w1/pi,abs(Y);title('信號(hào)的頻譜');實(shí)驗(yàn)四用FFT對(duì)信號(hào)作頻譜分析一、實(shí)驗(yàn)?zāi)康?、理解用FFT對(duì)周期序列進(jìn)行頻譜分

47、析時(shí)所面臨的問題并掌握其解決方法。2、掌握用時(shí)域窗函數(shù)加權(quán)處理的技術(shù)。3、理解用FFT對(duì)非周期信號(hào)進(jìn)行頻譜分析所面臨的問題并掌握其解決方法。二、實(shí)驗(yàn)原理與計(jì)算方法、對(duì)周期序列進(jìn)行頻譜分析應(yīng)注意的問題k k(a)時(shí)域周期整數(shù)倍截?cái)?(b)時(shí)域非周期整數(shù)倍截?cái)鄨D 4-1 周期函數(shù)的幅頻曲線k圖 4-2 矩形窗的頻譜對(duì)時(shí)間序列作FFT時(shí)，實(shí)際上要作周期延拓（如果取長序列的一段進(jìn)行計(jì)算還要先作截?cái)啵?。周期序列是無限長時(shí)間序列，如果截?cái)鄥^(qū)間剛好就是該序列周期的整數(shù)倍，那么在進(jìn)行周期延拓后，將還原出原來的周期序列，由此可以較精確地計(jì)算出的該周期序列的頻譜。反之，如果截?cái)鄥^(qū)間并不是該序列周期的整數(shù)倍，那么在

48、進(jìn)行周期延拓后，就不可能還原出原來的周期序列，由此計(jì)算出的頻譜與該周期序列的頻譜存在誤差，而且誤差的大小與截?cái)鄥^(qū)間的選取直接相關(guān)，如圖4-1所示，其中幅度頻譜的量值表示為，以dB(分貝)為單位。這種誤差是由于周期序列與矩形截?cái)嘈蛄邢喑嗽陬l域產(chǎn)生二者的頻譜卷積形成的。矩形窗的頻譜是抽樣函數(shù)序列，如圖4-2所示。除了k = 0處主瓣內(nèi)集中了大部分的能量外，兩旁的較小峰值處的旁瓣也分散了一部分能量，它與周期序列頻譜卷積的結(jié)果使原來集中的頻譜展寬，稱為頻率泄漏。如果對(duì)已知周期的信號(hào)作頻譜分析，在進(jìn)行時(shí)域截?cái)鄷r(shí)，完全可以選取其周期的整數(shù)倍裁取，從而可以避免這種頻率泄漏的發(fā)生。不過，通常需要進(jìn)行頻譜分析的

49、信號(hào)是周期未知的信號(hào)，或隨機(jī)信號(hào)，無法判斷它的周期值，為了盡量避免頻率泄漏對(duì)結(jié)果的影響，在作時(shí)間截?cái)鄷r(shí)，就應(yīng)選取其頻譜的旁瓣較小的截?cái)嗪瘮?shù)，以減輕泄漏問題。2、時(shí)域窗函數(shù)的應(yīng)用作為截?cái)嗪瘮?shù)，矩形窗在作時(shí)間截?cái)鄷r(shí)，對(duì)所截取區(qū)間內(nèi)的信號(hào)不加以任何影響，而其它的窗函數(shù)都將對(duì)所截取區(qū)間內(nèi)的信號(hào)作加權(quán)處理。除三角窗、Hanning窗和Hamming窗外，常用的窗函數(shù)還有很多，例如Parzen窗、Kaiser窗、Chebyshev窗、Tukey窗、Poisson窗、Caushy窗、Gaussian窗和Blackman窗等等。本次實(shí)驗(yàn)利用窗函數(shù)作時(shí)域加權(quán)截?cái)唷?0 t 0 k(a) 正弦函數(shù)的加權(quán)的非周期時(shí)

50、域截?cái)?b)減小了泄漏的頻譜圖 4-3 采用Hanning窗加權(quán)后的時(shí)域截?cái)嗪皖l譜圖 4-3 中給出了采用Hanning窗對(duì)正弦函數(shù)作非整周期的時(shí)域加權(quán)截?cái)嗪蟮牟ㄐ魏皖l譜，與圖4-1(b)比較，泄漏已明顯減少。應(yīng)該指出，前面所給出的窗函數(shù)都是定義為以0點(diǎn)為中心、寬度為N +1的加權(quán)函數(shù)，在這里應(yīng)用時(shí)，需要將其右移，成為區(qū)間內(nèi)的加權(quán)函數(shù)。3、對(duì)非周期序列進(jìn)行頻譜分析應(yīng)注意的問題（）混疊一般非周期信號(hào)作FFT之前要進(jìn)行時(shí)域采樣和周期延拓(無限長時(shí)間信號(hào)還應(yīng)先截?cái)嘣傺油?。根據(jù)Fourier變換理論，經(jīng)等周期的沖激采樣后，離散序列的頻譜是原信號(hào)頻譜以為周期的周期延拓。按照Nyquist采樣定理，由采

51、樣引起周期延拓后頻譜之間不發(fā)生混疊的條件是：(1)原信號(hào)應(yīng)該是有限帶寬信號(hào)，設(shè)其頻帶寬度為fm；(2)頻譜的周期，即采樣周期應(yīng)滿足Nyquist條件。0 n 0 N/2 N(a)時(shí)域按周期Ts采樣(b)頻域一個(gè)周期內(nèi)在N/2附近出現(xiàn)混疊圖 4-4非周期函數(shù)采樣后的幅頻曲線由于實(shí)際上有限長時(shí)間信號(hào)的FT是頻域的無限函數(shù)，因此采樣所得的離散序列的頻譜必定產(chǎn)生混疊，減小采樣周期只能減小而不能消除混疊。對(duì)于時(shí)間有限函數(shù)，當(dāng)采樣周期較大時(shí)，也會(huì)在FFT得到的頻域出現(xiàn)混疊，形成頻譜失真，造成頻譜分析結(jié)果與原信號(hào)的實(shí)際頻譜的差異，也無法恢復(fù)出原信號(hào)。當(dāng)然，實(shí)際工作中只要采樣和截?cái)喈a(chǎn)生的誤差在許可的范圍內(nèi)就行

52、了，但應(yīng)該認(rèn)識(shí)到混疊是引起頻譜分析誤差的一個(gè)主要原因。還應(yīng)該注意的是，離散Fourier變換的頻域也是周期化的，區(qū)間內(nèi)的樣點(diǎn)實(shí)際上是負(fù)頻率區(qū)的量值，因此如果出現(xiàn)混疊，就將在一個(gè)周期內(nèi)出現(xiàn)，并發(fā)生在附近的區(qū)域，如圖4-4所示。要減少混疊，就要盡量減小采樣周期。（）泄漏周期函數(shù)截?cái)嘁鸬念l率泄漏問題，在非周期函數(shù)截?cái)嗵幚砗笸瑯哟嬖?，這種誤差是由于采樣后的離散序列與矩形截?cái)嘈蛄邢喑嗽陬l域造成二者的頻譜卷積形成的。矩形窗的頻譜是抽樣函數(shù)序列，它與離散序列頻譜卷積的結(jié)果使原來集中于每一個(gè)樣點(diǎn)處的頻譜展寬，其影響在高頻區(qū)(接近N/2的樣點(diǎn))特別明顯，如圖4-5所示。同樣，為了盡量避免頻率泄漏對(duì)結(jié)果的影響，

53、在對(duì)非周期函數(shù)作時(shí)間截?cái)鄷r(shí)，除盡量增加截?cái)嘈蛄械膶挾韧猓矐?yīng)選取其頻譜的旁瓣較小的截?cái)嗪瘮?shù)，以減輕泄漏問題。0 n 0 N/2 N(a)時(shí)域截?cái)?b)頻域一個(gè)周期內(nèi)在N/2附近出現(xiàn)泄漏圖4-5函數(shù)采樣后作截?cái)嗟姆l曲線在選取了適當(dāng)?shù)拇昂瘮?shù)后，應(yīng)當(dāng)使窗函數(shù)的寬度與被處理的序列長度相同，如果作變換前還需要補(bǔ)零（例如為了作卷積運(yùn)算或避免柵欄效應(yīng)），則應(yīng)將原序列與窗函數(shù)相乘后再補(bǔ)零，即補(bǔ)零的樣點(diǎn)不用窗函數(shù)加權(quán)處理。（）柵欄效應(yīng)非周期信號(hào)應(yīng)具有連續(xù)的頻譜，在對(duì)作抽樣后進(jìn)行DFT，得到的是離散的頻譜。如果排除混疊和泄漏等誤差的影響，所得的結(jié)果也只是的連續(xù)頻譜上的個(gè)樣值。這就象通過柵欄上的等間距縫隙觀看到的

54、另一邊的景象，故此稱柵欄效應(yīng)。被柵欄遮住的景象中有可能存在與顯現(xiàn)出的頻譜差異較大的變化，即顯示信號(hào)特征的頻譜分量。為了使被柵欄遮住的部分能盡可能地顯現(xiàn)出來，可以采用增加頻域樣點(diǎn)密度的方法，即在不增加信號(hào)采樣點(diǎn)的情況下，用時(shí)域補(bǔ)零加寬變換尺度N來實(shí)現(xiàn)，稱為補(bǔ)零重構(gòu)。例如原來信號(hào)采樣得到12個(gè)樣點(diǎn)，在其后面再加上4個(gè)零，使序列的總長度為16個(gè)樣點(diǎn)。這樣處理的結(jié)果原來信號(hào)的采樣間隔和頻率都沒有改變，設(shè)采樣頻率為，經(jīng)補(bǔ)零重構(gòu)之后，采樣頻率仍然為，但是原來頻域樣點(diǎn)間寬度為/12，經(jīng)補(bǔ)零重構(gòu)之后頻域樣點(diǎn)間寬度為/16。這就使補(bǔ)零重構(gòu)之后頻域樣點(diǎn)密度增加，而且顯示出原來沒有顯露的一些頻率位置的頻譜。三、實(shí)驗(yàn)

55、內(nèi)容（1）將余弦函數(shù)cos(2pt)以Ts=1/53 s抽樣，對(duì)余弦序列做樣點(diǎn)數(shù)為N=128的FFT，畫出頻譜曲線，觀察并記錄頻率泄漏現(xiàn)象，然后用Hamming窗和三角窗作加權(quán)截?cái)?，觀察并記錄泄漏的衰減。實(shí)驗(yàn)代碼：n=0:127;N=128;Ts=1./53;t=n.*Ts;xn=cos(2.*pi.*t);w1=hamming(N);w2=bartlett(N);H=xn.*w1'B=xn.*w2'figure(1);subplot(2,1,1);stem(n,xn,'.');title('xn函數(shù)')subplot(2,1,2);stem(a

56、bs(fft(xn,N),'.');title('xn頻譜曲線')figure(2);subplot(2,1,1);stem(n,H,'.');title('Hamming窗加權(quán)')subplot(2,1,2);stem(abs(fft(H,N),'.');title('Hamming窗加權(quán)頻譜曲線')figure(3);subplot(2,1,1);stem(n,B,'.');title('三角窗加權(quán)')subplot(2,1,2);stem(abs(fft(B,N),'.');title('三角窗加權(quán)頻譜曲線')實(shí)驗(yàn)截圖：（2）將幅度為1，周期為2的方波信號(hào)，按Ts=1/37 s的間距抽樣，做樣點(diǎn)數(shù)N=128的FFT，畫出頻譜曲線，然后用Hamming窗和三角窗作加權(quán)截?cái)?，觀察并記錄作不同的加權(quán)截?cái)嘁鸬念l譜差異。實(shí)驗(yàn)代碼：n=0:1:127;N=12