DFT頻譜分析作業(yè)

1、DFT頻譜分析作業(yè)通信工程50作業(yè)一設(shè)周期序列,?截取N點長得到k(n)=x(n)RN(n)(1)N=10,做10點DFT,得到X1(k);(2)N=10,做100點補零DFT導(dǎo)到X2(k);(3)N=100,做100點DFT得到X3(k)o要求：針對以上三種情況?分別輸出|X1(k)|、|X2(k)卜|X3(k)|的圖形?并進行比較、分析和討論。解：代碼如下%雙頻周期信號的DFT分析n=0:1000;x=cos*pi*n)+cos*pi*n);%時域的離散時間信號%N=1010點DFT%N=10;%截取的序列長度N0=10;%DFT點數(shù)N0Nxn=x(1:N);%取一個周期的長度Xk=fft

2、(xn,N0);subplot(3,1,1)stem(0:1:N0-1,abs(Xk),fill,MarkerSize,%取模畫圖title(N=1010點DFT)xlabel(k);ylabel(|X1(k)|)grid%N=10100點補零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,%取模畫圖gridtitle(N=10100點補零DFT)xlabel(k);ylabel(|X2(k)|,)%N=10

3、0100點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,%取模畫圖gridtitle(N=100100點DFT)xlabel(k);ylabel(|X3(k)|)50p*py011*1ih0102030405060708090100k圖1DFT頻譜觀察對比：1)10點DFT的|X1(k)|出現(xiàn)明顯的柵欄效應(yīng)?相鄰譜線間隔很大?而100點DFT的間隔較小?說明通過增大DFT點數(shù)N0即頻域抽樣點數(shù)?可以減小柵

4、欄效應(yīng)；2)|X2(k)|與|X3(k)|相比出現(xiàn)了明顯的頻譜泄漏?即截斷效應(yīng)。由于DFT需要對時域序列進行截斷?對x(n)乘一個N點長的矩形窗函數(shù)h(N(n)(主瓣寬度為/N),N越大主瓣越小引起的泄漏就越小。作業(yè)?設(shè) xa。)-cosQoom)+痂 i(loom)十心 0 位疝 50 瞋)o 用 FFT 分析其頻譜結(jié)構(gòu)?選擇不同的截取長度 Tp,觀察截斷效應(yīng)?試用加窗的方法降低旁瓣。解：加矩形窗的FFT代碼如下100Nxn=x(1:N);%取一個周期的長度Xk=fft(xn,N0);Xk=abs(Xk)./max(abs(Xk);subplot(3,1,1)stem(0:1:N0-1,Xk

5、,fill,MarkerSize,%取模畫圖gridtitle(Tp=);xlabel(k);ylabel(|X(k)|)axis(022000)%Tp=4*Tp=4*;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,%取模畫圖gridtitle(Tp=4*);xlabel(k);ylabel(|X(k)|)axis(022000)%Tp=16*Tp=16*;N=fs*T

6、p;%截取的序列長度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,%取模畫圖gridtitle(Tp=16*);xlabel(k);ylabel(|X(k)|)axis(022000)運行?得Tp=0.04s10.50L7/、. .02004006008001000120014001600180020002200Tp=4*0.04s1IAnA.0.5I,A.廣ffeif-0,0200400600800

7、1000120014001600180020002200kTp=16*0.04s11 1 11 10.5.Niij1j1af月.J002004006008001000120014001600180020002200k圖2DFT頻譜10.50、.A/jf模擬頻率/Hz模擬頻率/HzTp=16*0.04sTp=0.04s02004006008001000120014001600180020002200Tp=4*0.04s模擬頻率/Hz圖3DFT頻譜（模擬頻率）由模擬頻率圖可以看出?在做到Tp?附乂0XJ4S也即N取256時?出現(xiàn)了三個頻點?就是原信號的三個頻點1值k50825兀。當(dāng)Tp=0.04.

8、4-0O4S也即NB16、64時頻率分辨率較低。加Hamming窗進行分析：%單頻周期信號的DFM析cos(pi*t)fs=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=Tp=;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,1)st

9、em(0:1:N0-1,Xk,fill,MarkerSize,%取模畫圖gridtitle(Tp=);xlabel(k);ylabel(|X(k)|)axis(022000)%Tp=4*Tp=4*;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,%取模畫圖gridtitle(Tp=4*);xlabel(k);ylabel(|X(k)|)axis(0

10、22000)%Tp=16*Tp=16*;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,%取模畫圖gridtitle(Tp=16*);xlabel(k);ylabel(|X(k)|)axis(022000)002004006008001000120014001600180020002200kTp=4*0.04skTp=16*0.04s0020040

11、06008001000120014001600180020002200k圖4加Hamming窗的DFT頻譜a a1 1Tp=0.04s10.5H HII1,1,10.5對比可以明顯看出?加窗后得到的頻譜圖上非信號頻點處幅度很低。這說明旁瓣確實低了。下圖對比更明顯：右側(cè)的為加矩形窗?左側(cè)為加海明窗。也可以看出雖然海明窗降低了旁瓣幅度、但是主瓣寬度變寬?主瓣泄漏厲害?本來能夠分開的兩個頻率分量幅值很相近。這顯示了加窗的矛盾。作業(yè)三讀取一個音樂文件并播放?取其中一段數(shù)據(jù)顯示其時域波形?利用DFT進行頻譜分析并顯示其幅度譜。解：clearall%聲音文件的DF邛析%data,fs,bits=wavre

12、ad(warning);data,fs,bits=wavread(clock);sound(data,fs)subplot(2,1,1)plot(data);axis(0length(data)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,%取模畫圖title(N=,num2str(N),num2str(N0),點DFT)xlabel(k);ylabel(|X1(k)

13、|)grid%subplot(3,1,3)%Xa=abs(Xk);%Xkdb=20*log10(Xa/max(Xa)%stem(0:1:N0-1,Xkdb,fill,MarkerSize,%取模畫圖0.9810.860.840.990.9810.98t t490500510520530540550kTp=0.04sTp=4*0.04sTp=16*0.04s1480490500510520530540kTp=0.04skTp=4*0.04skTp=16*0.04s%title(N=,num2str(N),num2str(N0),點DFT)%axis(0N0+5-500)%xlabel(k);ylabel(|X1(k)|)%grid運行后