DFT頻譜分析作業(yè)

1、DFT頻譜分析作業(yè)通信工程 2013040150作業(yè)一設周期序列xn=cos0.48n+cos0.52n截取 N 點長得到 x(n)=xnRN(n)（1）N=10，做10點DFT，得到X1(k);（2）N=10，做100點補零DFT,得到X2(k);（3）N=100，做100點DFT,得到X3(k)。要求：針對以上三種情況，分別輸出| X1(k)|、| X2(k)|、| X3(k)|的圖形，并進行比較、分析和討論。解：代碼如下%雙頻周期信號的DFT分析n=0:1000;x=cos(0.48*pi*n)+cos(0.52*pi*n);%時域的離散時間信號%-N=10 10點DFT-%N=10;%

2、截取的序列長度N0=10;%DFT點數(shù) N0Nxn=x(1:N);%取一個周期的長度Xk=fft(xn,N0);subplot(3,1,1)stem(0:1:N0-1,abs(Xk),fill,MarkerSize,3.75)%取模畫圖title(N=10 10點DFT)xlabel(k);ylabel(|X1(k)|)grid%-N=10 100點補零DFT-%N=10;%截取的序列長度N0=100;%DFT點數(shù) N0Nxn=x(1:N);%取一個周期的長度Xk=fft(xn,N0);subplot(3,1,2)stem(0:1:N0-1,abs(Xk),fill,MarkerSize,3.

3、75)%取模畫圖gridtitle(N=10 100點補零DFT)xlabel(k);ylabel(|X2(k)|)%-N=100 100點DFT-%N=100;%截取的序列長度N0=100;%DFT點數(shù) N0Nxn=x(1:N);%取一個周期的長度Xk=fft(xn,N0);subplot(3,1,3)stem(0:1:N0-1,abs(Xk),fill,MarkerSize,3.75)%取模畫圖gridtitle(N=100 100點DFT)xlabel(k);ylabel(|X3(k)|)運行,得結(jié)果如圖1所示：圖 1 DFT頻譜觀察對比：1) 10點DFT的| X1(k)|出現(xiàn)明顯的柵

4、欄效應，相鄰譜線間隔很大，而100點DFT的間隔較小，說明通過增大DFT點數(shù)N0即頻域抽樣點數(shù)，可以減小柵欄效應；2) | X2(k)|與| X3(k)|相比出現(xiàn)了明顯的頻譜泄漏，即截斷效應。由于DFT需要對時域序列進行截斷，對x(n)乘一個N點長的矩形窗函數(shù)RN(n) (主瓣寬度為4/N)， N越大主瓣越小引起的泄漏就越小。作業(yè)二設xat =cos200t+sin100t+cos(50t)。用FFT分析其頻譜結(jié)構(gòu)，選擇不同的截取長度Tp，觀察截斷效應，試用加窗的方法降低旁瓣。解：加矩形窗的FFT代碼如下%單頻周期信號的DFT分析 cos(pi * t) 0.5Hzfs=400;% 采樣頻率H

5、Zf1=100;f2=50;f3=25;%模擬頻率HZn=0:1000;x=cos(2*pi*f1*n/fs)+sin(2*pi*f2*n/fs)+cos(2*pi*f3*n/fs);%時域的離散時間信號%-Tp=0.04Tp=0.04;N=fs*Tp;%截取的序列長度N0=2048;%DFT點數(shù) N0Nxn=x(1:N);%取一個周期的長度Xk=fft(xn,N0);Xk=abs(Xk)./max(abs(Xk);subplot(3,1,1)stem(0:1:N0-1,Xk,fill,MarkerSize,3.75)%取模畫圖gridtitle(Tp=0.04s);xlabel(k);yla

6、bel(|X(k)|)axis(0 2200 0 1.2)%-Tp=4*0.04Tp=4*0.04;N=fs*Tp;%截取的序列長度N0=2048;%DFT點數(shù) N0Nxn=x(1:N);%取一個周期的長度Xk=fft(xn,N0);Xk=abs(Xk)./max(abs(Xk);subplot(3,1,2)stem(0:1:N0-1,Xk,fill,MarkerSize,3.75)%取模畫圖gridtitle(Tp=4*0.04s);xlabel(k);ylabel(|X(k)|)axis(0 2200 0 1.2)%-Tp=16*0.04Tp=16*0.04;N=fs*Tp;%截取的序列長

7、度N0=2048;%DFT點數(shù) N0Nxn=x(1:N);%取一個周期的長度Xk=fft(xn,N0);Xk=abs(Xk)./max(abs(Xk);subplot(3,1,3)stem(0:1:N0-1,Xk,fill,MarkerSize,3.75)%取模畫圖gridtitle(Tp=16*0.04s);xlabel(k);ylabel(|X(k)|)axis(0 2200 0 1.2)運行,得圖 2 DFT頻譜圖 3 DFT頻譜（模擬頻率）由模擬頻率圖可以看出，在做到Tp=160.04S 也即N取256時，出現(xiàn)了三個頻點，就是原信號的三個頻點100、50、25。當Tp=0.04S、40

8、.04S也即N取16、64時頻率分辨率較低。加Hamming窗進行分析：%單頻周期信號的DFT分析 cos(pi * t) 0.5Hzfs=400;% 采樣頻率HZf1=100;f2=50;f3=25;%模擬頻率HZn=0:1000;x=cos(2*pi*f1*n/fs)+sin(2*pi*f2*n/fs)+cos(2*pi*f3*n/fs);%時域的離散時間信號%-Tp=0.04Tp=0.04;N=fs*Tp;%截取的序列長度N0=2048;%DFT點數(shù) N0Nxn=x(1:N).*hamming(N);%取一個周期的長度Xk=fft(xn,N0);Xk=abs(Xk)./max(abs(X

9、k);subplot(3,1,1)stem(0:1:N0-1,Xk,fill,MarkerSize,3.75)%取模畫圖gridtitle(Tp=0.04s);xlabel(k);ylabel(|X(k)|)axis(0 2200 0 1.2)%-Tp=4*0.04Tp=4*0.04;N=fs*Tp;%截取的序列長度N0=2048;%DFT點數(shù) N0Nxn=x(1:N).*hamming(N);%取一個周期的長度Xk=fft(xn,N0);Xk=abs(Xk)./max(abs(Xk);subplot(3,1,2)stem(0:1:N0-1,Xk,fill,MarkerSize,3.75)%取

10、模畫圖gridtitle(Tp=4*0.04s);xlabel(k);ylabel(|X(k)|)axis(0 2200 0 1.2)%-Tp=16*0.04Tp=16*0.04;N=fs*Tp;%截取的序列長度N0=2048;%DFT點數(shù) N0Nxn=x(1:N).*hamming(N);%取一個周期的長度Xk=fft(xn,N0);Xk=abs(Xk)./max(abs(Xk);subplot(3,1,3)stem(0:1:N0-1,Xk,fill,MarkerSize,3.75)%取模畫圖gridtitle(Tp=16*0.04s);xlabel(k);ylabel(|X(k)|)axi

11、s(0 2200 0 1.2)圖 4 加Hamming窗的DFT頻譜對比可以明顯看出，加窗后得到的頻譜圖上非信號頻點處幅度很低。這說明旁瓣確實低了。下圖對比更明顯： 右側(cè)的為加矩形窗，左側(cè)為加海明窗。也可以看出雖然海明窗降低了旁瓣幅度、但是主瓣寬度變寬，主瓣泄漏厲害，本來能夠分開的兩個頻率分量幅值很相近。這顯示了加窗的矛盾。作業(yè)三讀取一個音樂文件并播放，取其中一段數(shù)據(jù)顯示其時域波形，利用DFT進行頻譜分析并顯示其幅度譜。解：clear all%聲音文件的DFT分析%data,fs,bits=wavread(warning);data,fs,bits=wavread(clock);sound(d

12、ata,fs)subplot(2,1,1)plot(data);axis(0 length(data) -1.2 1.2)title(聲音信號波形圖)N=65536;%截取的序列長度N0=65536*2;%DFT點數(shù) N0Nxn=data(1:N);%取一個周期的長度Xk=fft(xn,N0);subplot(2,1,2)plot(0:1:N0-1,abs(Xk)%fill,MarkerSize,3.75)%取模畫圖title(N=,num2str(N), ,num2str(N0), 點DFT)xlabel(k);ylabel(|X1(k)|)grid% subplot(3,1,3)% Xa=abs(Xk);% Xkdb=20*log10(Xa/max(Xa)% stem(0:1:N0-1,Xkdb,fill,MarkerSize,3.75)%取模畫圖% title(N=,num2str(N), ,num2str(N0), 點DFT)% axis(0 N0+5 -50 0)% xlabel(k);ylabe