DFT和FFT實(shí)驗(yàn)(上傳)

1、DFT和FFT實(shí)驗(yàn)一、實(shí)驗(yàn)?zāi)康暮鸵?、掌握DFT變換2、掌握DFT性質(zhì)3、掌握快速傅立葉變換(FFT)二、實(shí)驗(yàn)容和原理1、實(shí)驗(yàn)容,1)求有限長(zhǎng)離散時(shí)間信號(hào)的離散時(shí)間傅立葉變換X(ej )并繪圖。裝n訂？已知 x(n) 0.9e 30 n 10線？已知 x(n) 2n10 n 10J 2)已知序列x(n) cos(0.82 n) 2sin( n), 0 n 50,繪制x(n)及其離散傅立葉變換 X(k)的幅度、：相位圖。r:3)設(shè) x(n) sin(0.2 n) randn(n) , 0 n N 1,其中，randn(n)為高斯白噪聲。求出 N4m , m=2,3,4的matlab采用不同算法

2、的執(zhí)行時(shí)間。94)研究高密度頻譜和高分辨率頻譜。：設(shè)有連續(xù)信號(hào)x(t) cos(2 6.5 103t) cos(2 7 103t) cos(2 9 103t)?以采樣頻率fs 32kHz對(duì)信號(hào)x(t)采樣，分析下列三種情況的幅頻特性。9?采集數(shù)據(jù)長(zhǎng)度N=16點(diǎn)，做N=16點(diǎn)的FFT,并畫出幅頻特性。?采集數(shù)據(jù)長(zhǎng)度N=16點(diǎn)，補(bǔ)零到256點(diǎn)，做N=256點(diǎn)的FFT,并畫出幅頻特性。:?采集數(shù)據(jù)長(zhǎng)度N=256點(diǎn)，做N=256點(diǎn)的FFT,并畫出幅頻特性。:觀察三種不同頻率特性圖，分析和比較它們的特點(diǎn)以及形成的原因。2、實(shí)驗(yàn)原理1) DFT序列x(n)的離散時(shí)間傅里葉變換(DTFT)表示為X(ej )

3、x( n)e jnN 1如果x(n)為因果有限長(zhǎng)序列，n=0,1,.,N-1 ,則x(n)的DTFT表示為X(ej ) x(n)e n 0N 1j 2 nkx(n)的離散傅里葉變換(DFT)表達(dá)式為X(k) x(n)e N (k 0,1，N 1)n 0序列的N點(diǎn)DFT是序列DTFT在頻率區(qū)間0,2可上的N點(diǎn)燈間隔采樣，采樣間隔為2”N。通過DFT ,可以完成由一組有限個(gè)信號(hào)采樣值x(n)直接計(jì)算得到一組有限個(gè)頻譜采樣值X(k)o X(k)的幅度譜為X(k) VXR(k)X2(k),其中下標(biāo)R和I分別表示取實(shí)部、虛部的運(yùn)算。X(k)的相位譜為(k) arctanXR(k)一、4，1 N離散傅里葉

4、反變換(IDFT)定義為x(n) 一N n1j_2_nkX(k)e N (n 0,1,., N 1)。02) FFTj2-n快速傅里葉變換(FFT)是DFT的快速算法，并不是一個(gè)新的映射。FFT利用了 e N函數(shù)的周期性和對(duì)稱性以及一些特殊值來減少DFT的運(yùn)算量，可使 DFT的運(yùn)算量下降幾個(gè)數(shù)量級(jí)，從而使數(shù)字信號(hào)處理的速度大大提高。若信號(hào)是連續(xù)信號(hào),用FFT進(jìn)行譜分析時(shí)，首先必須對(duì)信號(hào)進(jìn)行采樣，使之變成離散信號(hào)， 然后就可以用FFT來對(duì)連續(xù)信號(hào)進(jìn)行譜分析。為了滿足采樣定理，一般在采樣之前要設(shè)置一個(gè)抗混疊低通濾波器, 且抗混疊濾波器的截止頻率不得高于與采樣頻率的一半。比較DFT和IDFT的定義，

5、兩者的區(qū)別僅在于指數(shù)因子的指數(shù)部分的符號(hào)差異和幅度尺度變換， 可用FFT算法來計(jì)算IDFT 。因此三、主要儀器設(shè)備Matlab四、操作方法和實(shí)驗(yàn)步驟1、認(rèn)真分析原函數(shù)，取點(diǎn)2、用matlab編寫程序，運(yùn)行程序得出結(jié)果五、實(shí)驗(yàn)數(shù)據(jù)記錄、處理和分析1、求有限長(zhǎng)離散時(shí)間信號(hào)的離散時(shí)間傅立葉變換X(ej )并繪圖。n? 已知 x(n) 0.9e 30 n 10?已知 x(n) 2n10 n 10【解答】思路：這是一道 DFT的題，按照題目要求只需要取11個(gè)點(diǎn)即可。第(1)小題M文件源代碼N=11;% 取點(diǎn)個(gè)數(shù)為11個(gè)j=sqrt(-1);%定義j為復(fù)數(shù)單位f=inline('(0.9*exp(

6、j*pi/3)An','n');W=0:2*pi/1000:2*pi;%Xw=zeros(size(W);%for n=0:N-1Xw=Xw+f(n)*exp(-j*W*n);end%定義一個(gè)函數(shù)f(n)定義離散域的基本頻率W為數(shù)組，間距為2*pi/1000定義一個(gè)與W位數(shù)相等的數(shù)組對(duì)f(n) 函數(shù)做 DFT變換xn=;for n=0:N-1xn(n+1)=f(n);endmagXw=abs(Xw);angleXw= angle(Xw);的值放進(jìn)數(shù)組xn里面，便于最后畫出 xn的圖像定義一個(gè)數(shù)組magXw,將abw(Xw)的值賦給它定義數(shù)組angleXw, 將angle

7、(Xw)的值賦給它figure。)；plot(xn,'.-');xlabel('n'); ylabel('x(n)'); %figure(2);k=0:1:N-1;plot(W,magXw,'-') xlabel('W); ylabel('|X(W)|'); figure(3);plot(W,angleXw,'-');xlabel(W); ylabel('angle(X(W)'); %運(yùn)行結(jié)果xn圖像：0.8-0.6畫出xn的圖畫出magXw的圖像畫出angleXw 的圖像0

8、.60.40.20-0.2-0.4-0.8-0.6-0.4-0.200.20.40.60.81nX(W)的幅度圖X(W)的相位圖)a【分析】11-1=10可見整的幅度頻譜有ii-i=i0個(gè)極大，ii-i=i0個(gè)極小。而的相位則有ii-i=i0個(gè)極大, 個(gè)極小，并且相位在-區(qū)和區(qū)之間擺動(dòng)。第(2)小題M文件源代碼N=10;%j=sqrt(-1);%f=inline('2An','n'); %W=0:2*pi/1000:2*pi; %Xw=zeros(size(W); % for n=-N:N取點(diǎn)個(gè)數(shù)為11個(gè)定義j為復(fù)數(shù)單位定義一個(gè)函數(shù)f(n)定義離散域的基本頻率，

9、將其設(shè)置為間距為2*pi/1000 的數(shù)組定義一個(gè)數(shù)組Xw,位數(shù)與W相等Xw=Xw+f(n)*exp(-j*W*n);end %對(duì)f(n) 函數(shù)做 DTFT變換xn=;%for n=-N:Nxn(n+1+N)=f(n);end%magXw=abs(Xw);angleXw=angle(Xw);figure。);plot(xn,'.-');xlabel('n'); ylabel('x(n)');figure(2);plot(W,magXw,'-')定義一個(gè)數(shù)組xn將f(n)的值放進(jìn)數(shù)組xn里面，便于最后畫出 xn的圖像%將X(W)的模

10、值放進(jìn)數(shù)組magXw%將X(W)的相位放進(jìn)數(shù)組 magXw%畫出xn的圖xlabel('W); ylabel('|X(w)|'); % figure(3);畫出magXw的圖像plot(W,angleXw,'-');xlabel('W); ylabel('angle(X(w)'); %畫出angleXw 的圖像xn圖像1200100080060040020010152025nX(W)的幅度圖2200200018001600140012001000800600X(W)的相位圖4-401234567【分析】曲卬） 1110x（n）=

11、y 2“ 川= n -10W(n)的相位在一凡到江之間來回振動(dòng)，并且中間出現(xiàn)突變的可見*(。的幅度有一個(gè)極大值，一個(gè)極小值。情況。2)已知序列x(n)cos(0. 82n) 2sin( n) , 0 n 50 ,繪制x(n)及其離散傅立葉變換 X(k)的幅度、相位圖。【解答】思路：這是一道 DFT的題，按照題目要求只需要取51個(gè)點(diǎn)即可。M文件源代碼N=51;%取點(diǎn)個(gè)數(shù)為50 個(gè)定義一個(gè)函數(shù)f(n)j=sqrt(-1);%定義 j 為復(fù)數(shù)單位f=inline('cos(0.82*pi*n)+2*sin(pi*n)','n');%Xk=;%定義一個(gè)數(shù)組XkW=2*p

12、i/N;%定義離散域的基本頻率for k=1:NXk(k)=0;for n=0:N-1Xk(k)=Xk(k)+f(n)*exp(-j*(k-1)*W*n);endend%xn=;%for n=0:N-1xn(n+1)=f(n);對(duì) f(n) 函數(shù)做DFT 變換定義一個(gè)數(shù)組xnend%magXk=;%for k=1:NmagXk(k)=abs(Xk(k);end%angleXk=;%for k=1:NangleXk(k)=angle(Xk(k);end%figure(1);plot(xn,'.-');xlabel('n'); ylabel('x(n)

13、9;); % figure(2);k=0:1:N-1;plot(k,magXk,'+-') xlabel('k'); ylabel('|X(k)|'); % figure(3);plot(k,angleXk,'x-');將 f(n) 的值放進(jìn)數(shù)組xn 里面，便于最后畫出定義一個(gè)數(shù)組magXk將 X(kW) 的模值放進(jìn)數(shù)組magXk定義數(shù)組angleXk將 X(kW) 的相位放進(jìn)數(shù)組magXk畫出xn 的圖畫出magXk 的圖像xnxlabel('k'); ylabel('angle(X(k)');

14、命令窗口中的運(yùn)行及其結(jié)果：%畫出angleXk 的圖像xn 圖像的圖像10.80.60.40.20-0.2-0.4-0.6-0.8-1102030405060Xk的幅度圖30252015100T11 l lin l【l I1_l I U L U L 05101520253035404550Xk的相位圖【分析】ejt).B2Jm + e-p.82mi 一四色/c 51部)=y(情等cosO.827m + 2simm)e = J小勃+ 卦1 1 (J,必/3)>1:小如十碟5171。必+或)可見的幅度頻譜擁有兩個(gè)峰值，尺(的的相位頻譜在一汽到冗之間來回振動(dòng)，且中間存在3個(gè)臺(tái)階式的向 下跳變，

15、一個(gè)臺(tái)階式的向上跳變。3 .設(shè) x(n) sin(0.2 n) randn(n) , 0 n N 1,其中，randn(n)為高斯白噪聲。求出 N4m , m=2,3,4的matlab采用不同算法的執(zhí)行時(shí)間?！窘獯稹克悸罚河?jì)算DFT算法和FFT算法的運(yùn)行時(shí)間可以使用的etime函數(shù)M文件源代碼m=input('m=:');%輸入m值N=4Am;%求出所取得x(n)的點(diǎn)數(shù)j=sqrt(-1);%定義j為復(fù)數(shù)單位arr=;%定義一個(gè)數(shù)組arrW=2*pi/N;%定義離散域的基本頻率dft_time=0;%定義dft_time為0t1=clock;%此處為dft計(jì)算的時(shí)間起點(diǎn)for

16、k=1:Narr(k)=0;for n=0:N-1arr(k)=arr(k)+(sin(0.2*pi*n)+rand(1)*exp(-j*(k-1)*W*n); endend%dft_time=etime(clock,t1) %Wn=exp(-j*2*pi/N); %xn=zeros(1,N); %fft_time=0;%t2=clock;%for n=0:N-1對(duì) x(n) 做 dft 變換得出 dft 變換所花的時(shí)間求出旋轉(zhuǎn)因子定義一個(gè)N 位數(shù)組定義 fft_time為 0此處為 fft 計(jì)算的時(shí)間起點(diǎn)xn(n+1)=sin(0.2*pi*n)+randn(1);%end將 x(n) 的值

17、放入數(shù)組xn 中n1=fliplr(dec2bin(0:N-1);%n2=bin2dec(n1);%x=zeros(2*m+1,N);%for i=1:Nx(1,i)=xn(n2(i)+1);end%for i=1:2*m%Number=2A(2*m-i);Interval_of_Unit=2A(i-1); % Interval_of_Group=2Ai; % Wnr=;%for r=1:2A(i-1)Wnr(r)=WnA(r-1)*N/2Ai);end%碼位倒置步驟1 ：將碼位轉(zhuǎn)換為二進(jìn)制，再進(jìn)行倒序碼位倒置步驟2：將碼位轉(zhuǎn)換為十進(jìn)制后翻轉(zhuǎn)定義一個(gè)(2m+1)XN的矩陣將碼位倒序后的值重新賦

18、值進(jìn)入矩陣第一列進(jìn)行第 v 級(jí)計(jì)算每組中每個(gè)計(jì)算單元的間距每組之間的間距定義一個(gè)新數(shù)組Wnr將每一級(jí)運(yùn)算的指數(shù)因子賦值給Wnrfor k=0:Number-1for l=1:2A(i-1)x(i+1,l+k*2Ai)=x(i,l+k*2Ai)+Wnr(l)*x(i,l+k*2Ai+2A(i-1);x(i+1,l+k*2Ai+2A(i-1)=x(i,l+k*2Ai)-Wnr(l)*x(i,l+k*2Ai+2A(i-1);end對(duì) x(n) 做 fft 變換得出 dft 變換所需要的時(shí)間endend%fft_time=etime(clock,t2) %命令窗口中的運(yùn)行及其結(jié)果m=:2dft_tim

19、e =0fft_time =0 m=:3dft_time = 0.000ffttime =0.000 m=:4dft_time =0.000fft_time =0.000【分析】從實(shí)驗(yàn)結(jié)果可以看出，對(duì)于同樣的x(n) ， FFT 變換的計(jì)算時(shí)間小于DFT 變換的計(jì)算時(shí)間，并且隨著取樣點(diǎn)數(shù)的增多，這種差距越來越明顯。這是由于FFT 算法利用了旋轉(zhuǎn)因子的對(duì)稱性和周期性，將長(zhǎng)序列分解為短序列，分級(jí)進(jìn)行蝶形運(yùn)算使得復(fù)數(shù)乘法的次數(shù)減少到分解前的一半。通過對(duì)短序列的計(jì)算進(jìn)行適當(dāng)?shù)慕M合，從而達(dá)到了刪除重復(fù)運(yùn)算，提供運(yùn)算速度的目的。4 ) 研究高密度頻譜和高分辨率頻譜。設(shè)有連續(xù)信號(hào)x(t) cos(2 6.5

20、 103t) cos(2 7 103t) cos(2 9 103t)?以采樣頻率fs 32kHz對(duì)信號(hào)x(t)采樣，分析下列三種情況的幅頻特性。?采集數(shù)據(jù)長(zhǎng)度N=16點(diǎn)，做N=16點(diǎn)的FFT,并畫出幅頻特性。?采集數(shù)據(jù)長(zhǎng)度N=16點(diǎn)，補(bǔ)零到256點(diǎn)，做N=256點(diǎn)的FFT,并畫出幅頻特性。?采集數(shù)據(jù)長(zhǎng)度N=256點(diǎn)，做N=256點(diǎn)的FFT,并畫出幅頻特性。觀察三種不同頻率特性圖，分析和比較它們的特點(diǎn)以及形成的原因?！窘獯稹?思路：此處為FFT變換，要注意采樣數(shù)據(jù)長(zhǎng)度 N和FFT變換點(diǎn)數(shù)M不同，因此當(dāng)M>N時(shí)，應(yīng)該將x(n)中非 0-N-1 圍的點(diǎn)補(bǔ)零。M文件源代碼N=input('

21、;N=:');%M=input('M=:');%fs=32000;v=log2(M); %輸入采樣數(shù)據(jù)長(zhǎng)度N輸入 FFT 變換點(diǎn)數(shù)M給采樣頻率fs 賦值為 32000 ， v 為 FFT 變換的級(jí)數(shù)Wn=exp(-j*2*pi/M); %j=sqrt(-1);%xn=zeros(1,M); %定義旋轉(zhuǎn)因子Wn定義復(fù)數(shù)單位j定義一個(gè)M 位數(shù)組 xnfor n=0:N-1xn(n+1)=cos(2*pi*6.5*1000*n/fs)+cos(2*pi*7*1000*n/fs)+cos(2*pi*9*1000*n/fs);end將 x(n) 的值賦值進(jìn)入xn 中n1=flip

22、lr(dec2bin(0:M-1);%碼位倒置步驟1 ：將碼位轉(zhuǎn)換為二進(jìn)制，再進(jìn)行倒序n2=bin2dec(n1);%x=zeros(v+1,M); %for i=1:M碼位倒置步驟2 ：將碼位轉(zhuǎn)換為十進(jìn)制后翻轉(zhuǎn)定義一個(gè)(v+1) XM的矩陣xx(1,i)=xn(n2(i)+1);end %將碼位倒序后的for i=1:v % 進(jìn)行第 v 級(jí)計(jì)算x(n) 值賦值進(jìn)入矩陣x 的第一行Number=2A(v-i);%Interval_of_Unit=2A(i-1); %Interval_of_Group=2Ai; %Wnr=;%每一級(jí)計(jì)算中的“ 群 ” 數(shù)每組中每個(gè)計(jì)算單元的間距每組之間的間距定義一個(gè)數(shù)組Wnrfor r=1:2A(i-1)Wnr(r)=WnA(r-1)*M/2Ai);將每一級(jí)運(yùn)算的指數(shù)因子賦值給Wnrend%for k=0:Number-1for l=1:2A(i-1)x(i+1,l+k*2Ai)=x(i,l+k*2Ai)+Wnr(l)*x(i,l+k*2Ai+2A(i-1)；x(i+1,l+k*2Ai+2A(i-1)=x(i,l+k*2Ai)-Wnr(l)*x(i,l+k*2Ai+2A(i-1)； endendend%Xk=;%for k=1:MXk(k)=x(v+1,k);end%figure。);k=