信號(hào)處理實(shí)驗(yàn)四離散傅里葉變換

1、哈爾濱工程大學(xué)實(shí) 驗(yàn) 報(bào) 告實(shí) 驗(yàn) 名 稱： 實(shí)驗(yàn)四：離散傅里葉變換 班 級(jí)： 電子信息工程4班 學(xué) 號(hào)： 姓 名： 實(shí) 驗(yàn) 時(shí) 間： 2016年10月19日 成 績(jī)：_指 導(dǎo) 教 師： 欒曉明 實(shí)驗(yàn)室名稱： 數(shù)字信號(hào)處理實(shí)驗(yàn)室 哈爾濱工程大學(xué)實(shí)驗(yàn)室與資產(chǎn)管理處 制實(shí)驗(yàn)四 離散傅里葉變換一、 實(shí)驗(yàn)原理1. 由DFT定義式： k=0, 1 , , N-1,將其寫(xiě)成矩陣方程表示為利用MATLAB的矩陣運(yùn)算功能，可編寫(xiě)出計(jì)算DFT的函數(shù)文件。function Xk = dft(xn,N)%計(jì)算離散傅里葉變換%Xk = 在0<=k<=N-1間的DFT系數(shù)數(shù)組%xn = N點(diǎn)有限長(zhǎng)序列% N

2、 = DFT的長(zhǎng)度n = 0:1:N-1;%n的行向量k = 0:1:N-1;%k的行向量WN = exp(-j*2*pi/N);%Wn因子nk = n'*k;%產(chǎn)生一個(gè)含bk值的N乘N維矩陣WNnk = WN.nk;%DFT矩陣Xk = xn*WNnk;%DFT系數(shù)的行向量由IDFT定義式： ，n = 0, 1, 2, ,N-1，利用MATLAB的矩陣運(yùn)算功能，可編寫(xiě)出計(jì)算傅里葉反變換的函數(shù)文件。function xn = idft(Xk,N)%計(jì)算離散傅里葉反變換%-%xn = 在0<=n<=N-1%Xk = N點(diǎn)有限長(zhǎng)序列% N = IDFT的長(zhǎng)度k = 0:1:N-

3、1;%k的行向量n = 0:1:N-1;%n的行向量WN = exp(-j*2*pi/N);%Wn因子nk = n'*k;%產(chǎn)生一個(gè)含bk值的N乘N維矩陣WNnk = WN.nk;%DFT矩陣xn = Xk*WNnk;%傅里葉反變換計(jì)算序列值 DFT的快速算法FFT利用了的三個(gè)固有特性：(1)對(duì)稱性，(2)周期性，(3)可約性，和。FFT算法基本上可以分為兩大類，即按時(shí)間抽選法(DIT，Decimation-In-Time)和按頻率抽選法(DIF，Decimation-In-Frequency )。MATLAB中提供了進(jìn)行快速傅里葉變換的fft函數(shù)：X=fft(x)，基2時(shí)間抽取FFT

4、算法，x是表示離散信號(hào)的向量；X是系數(shù)向量；X=fft(x,N)，補(bǔ)零或截?cái)嗟腘點(diǎn)DFT，當(dāng)x的長(zhǎng)度小于N時(shí)，對(duì)x補(bǔ)零使其長(zhǎng)度為N，當(dāng)x的長(zhǎng)度大于N時(shí)，對(duì)x截?cái)嗍蛊溟L(zhǎng)度為N。Ifft函數(shù)計(jì)算IDFT，其調(diào)用格式與fft函數(shù)相同，參考help文件。2利用DFT做連續(xù)信號(hào)的頻譜分析DFT(實(shí)際中用FFT計(jì)算)可用來(lái)對(duì)連續(xù)信號(hào)和數(shù)字信號(hào)進(jìn)行譜分析。在實(shí)際分析過(guò)程中，要對(duì)連續(xù)信號(hào)采樣和截?cái)啵纱丝赡芤鸱治稣`差。(1) 混疊效應(yīng)對(duì)連續(xù)信號(hào)進(jìn)行頻譜分析時(shí)，首先要對(duì)其采樣，變成時(shí)域離散信號(hào)后才能用DFT（FFT）進(jìn)行譜分析。采樣速率fs必須滿足采樣定理，否則會(huì)在w=（對(duì)應(yīng)模擬頻率f=fs/2）附近發(fā)生頻譜

5、混疊現(xiàn)象。(2) 截?cái)嘈?yīng)處理實(shí)際信號(hào)序列x(n)時(shí)，一般總要將它截?cái)酁橐挥邢揲L(zhǎng)序列，長(zhǎng)為N點(diǎn)，相當(dāng)于乘以一個(gè)矩形窗形成有限長(zhǎng)序列y(n)=x(n)w(n)。矩形窗函數(shù)其頻譜有主瓣，也許許多副瓣，窗口越大，主瓣越窄，當(dāng)窗口趨于無(wú)窮大時(shí)，就是一個(gè)沖擊函數(shù)。時(shí)域的乘積對(duì)應(yīng)于頻域的卷積，所以，加窗后的頻域?qū)嶋H是原信號(hào)頻譜與矩形窗函數(shù)頻譜的卷積，卷積的結(jié)果使頻譜延伸到了主瓣以外，且一直延時(shí)到無(wú)窮。當(dāng)窗口無(wú)窮大時(shí)，與沖擊函數(shù)的卷積才是其本身，這是無(wú)需畸變。由此可見(jiàn)，階段后頻譜Y(ejw)與原序列頻譜X(ejw)必然有差別，這種差別表現(xiàn)在：a. 頻譜泄露。原來(lái)序列x(n)的頻譜是離散譜線，經(jīng)截?cái)嗪?，原?lái)離

6、散譜線向附近展寬，成為泄露。顯然，泄露使頻譜變模糊，使譜分辨率降低。b. 譜間干擾。主譜線兩邊又很多旁譜，引起不同頻率分量間干擾，這使譜分析產(chǎn)生較大偏差。程度與窗函數(shù)幅度譜主瓣寬度直接相關(guān)。(3) 柵欄效應(yīng)N點(diǎn)DFT是頻率區(qū)間0,2上對(duì)時(shí)域離散信號(hào)的頻譜進(jìn)行N點(diǎn)等間隔采樣，而采樣點(diǎn)之間的頻譜函數(shù)是看不見(jiàn)的。這就好像從N個(gè)柵欄縫隙中觀看信號(hào)的頻譜情況，僅得到N個(gè)縫隙中看到的函數(shù)值。由于柵欄效應(yīng)，有可能漏掉大的頻譜分量。例3.1 對(duì)連續(xù)的單一頻率周期信號(hào)按采樣頻率fs=8fa采樣，截取長(zhǎng)度N分別選N=20和N=16，觀察其DFT結(jié)果的幅度譜。解：此時(shí)離散序列，用MATLAB計(jì)算并作圖，函數(shù)fft(

7、)用于計(jì)算離散傅里葉變換DFT，程序如下：k=8;n1=0:19; %序列點(diǎn)數(shù)20xa1=sin(2*pi*n1/k);subplot(221)stem(n1,xa1)xlabel('t/T');ylabel('x(n)');%圖一：序列圖像xk1=fft(xa1);%調(diào)用fft函數(shù)xk1=abs(xk1);%取xk1絕對(duì)值subplot(222)stem(n1,xk1)xlabel('k');ylabel('X(k)');圖二：序列DFT圖像n2=0:1:15;%序列點(diǎn)數(shù)16xa2=sin(2*pi*n2/k);subplot(

8、223)stem(n2,xa2)xlabel('t/T');ylabel('x(n)');xk2=fft(xa2);xk2=abs(xk2);subplot(224)stem(n2,xk2)xlabel('k');ylabel('X(k)'); 二、 實(shí)驗(yàn)內(nèi)容1已知連續(xù)周期信號(hào) (1) 確定信號(hào)的基頻和基本周期.(2) 當(dāng)信號(hào)長(zhǎng)度取0.5Tp,1.5Tp,2Tp,時(shí),對(duì)x(t)采樣，利用FFT計(jì)算其幅度譜；對(duì)所得結(jié)果進(jìn)行比較，總結(jié)應(yīng)如何選取分析長(zhǎng)度。解：(1) 基頻=2，基本周期T=1s(2) 采樣頻率fs=10fa采樣，0.5T

9、p程序：f=10; %設(shè)置fs=10fan=0:4; %設(shè)置分析長(zhǎng)度為0.5Tpxa=cos(10*pi*n/f)+2*sin(18*pi*n/f); %設(shè)置采樣函數(shù)subplot(211)stem(n,xa) xlabel('t/T');ylabel('x(n)'); %圖一：采樣的結(jié)果title('采樣結(jié)果')xk=fft(xa); %進(jìn)行FFT變換xk=abs(xk);%取X(k)的絕對(duì)值subplot(212)stem(n,xk)xlabel('k');ylabel('X(k)'); %繪制出DFT結(jié)果title('幅度譜') 采樣頻率fs=10fa采樣，1.5Tp時(shí)將上述程序中n=0,4改為n=0,14即可采樣頻率fs=10fa采樣，2Tp時(shí)將上述程序中n=0,4改為n=0,19即可結(jié)果分析：有三圖比較可知，當(dāng)長(zhǎng)度為2Tp時(shí)，采樣信號(hào)頻譜不發(fā)生泄露，故應(yīng)取2Tp2對(duì)模擬信號(hào) 以采樣間隔T=0.01s采樣，分別取N=40，N=50，N=60得到x(n)，用N點(diǎn)DFT得到xa(t)幅度譜的估計(jì)并比較結(jié)果解：T=0.01s,N=40時(shí)程序：T=0.01;%采樣間隔N=40;%采樣點(diǎn)數(shù)n=0:N-1;xa=2*sin(4*pi*n*T)+5*cos(8*pi*n*T);%對(duì)信