數(shù)字信號處理研討_DFT近似計算信號頻譜

1、 Beijing Jiaotong University 數(shù)字信號處理研討DFT近似計算信號頻譜學(xué) 院：電子信息工程學(xué)院小組成員：指導(dǎo)教師：時 間：2013.06.24DFT近似計算信號頻譜專題研討【目的】(1) 掌握利用DFT近似計算不同類型信號頻譜的原理和方法。(2) 理解誤差產(chǎn)生的原因及減小誤差的方法。(3) 培養(yǎng)學(xué)生自主學(xué)習(xí)能力，以及發(fā)現(xiàn)問題、分析問題和解決問題的能力?！狙杏戭}目】 基本題 1. 已知一離散序列為 (1)用L=32點DFT計算該序列的頻譜，求出頻譜中譜峰的頻率；(2)對序列進行補零，然后分別用L=64、128、256、512點DFT計算該序列的頻譜，求出頻譜中譜峰的頻率

2、；(3)討論所獲得的結(jié)果，給出你的結(jié)論。該結(jié)論對序列的頻譜計算有何指導(dǎo)意義？【題目分析】本題討論補零對離散序列頻譜計算的影響?！緶仨嗵崾尽吭谟嬎汶x散非周期序列頻譜時常用W/p作為橫坐標，稱W/p為歸一化頻率(normalized frequency)。在畫頻譜時需給出橫坐標。每幅圖下都需給出簡要的文字說明。由于離散非周期序列頻譜是周期的，所以在計算時不必用fftshift 函數(shù)對fft計算的結(jié)果進行重新排列?！拘蛄蓄l譜計算的基本方法】連續(xù)信號，通過其抽樣的離散信號，和離散信號的DFT變換存在如下關(guān)系：通過如上關(guān)系，我們就可以通過DFT來求信號的頻譜?！痉抡娼Y(jié)果】fm = 0.9375fm =

3、0.9219fm = 0.9141fm = 0.9102fm = 0.9082【結(jié)果分析】增加DFT的點數(shù)可以使頻譜更容易觀察，即減輕了柵欄效應(yīng)帶來的影響。頻譜的橫坐標為歸一化頻率，所以原信號的峰值第一次應(yīng)該出現(xiàn)在0.2處，隨著DFT點數(shù)的增大，頻譜表示也越來越精確?！咀灾鲗W(xué)習(xí)內(nèi)容】1. 歸一化頻率相關(guān)知識。2. 通過matlab計算DFT和matlab的繪圖操作。【閱讀文獻】1數(shù)字信號處理2補零對有限長序列頻譜及DFT的影響。【發(fā)現(xiàn)問題】 (專題研討或相關(guān)知識點學(xué)習(xí)中發(fā)現(xiàn)的問題)：DFT點數(shù)的增多是否能提高頻譜分辨率？【問題探究】雖然DFT點數(shù)增加使圖像更加細致，但是因為不論DFT點數(shù)是多少

4、，抽樣點數(shù)都是相同的，所以每個頻譜所包含的信息相同，頻譜分辨率只與抽樣點數(shù)相關(guān)，與DFT點數(shù)無關(guān)，頻譜分辨率相同。所以不能通過增大DFT點數(shù)而減少信息損失。DFT點數(shù)的增多不能提過頻率分辨率?！痉抡娉绦颉縦 = 0:31;N = 16; x = sin(0.2*pi*k);for i = 1:5 N = 2*N; X = fft(x,N); k = 0:N-1; subplot (5,1,i); plot (k/N,abs(X); xlabel ('f/HZ'); ylabel ('Magnitude'); fm = find( X = max(X)/N; fm

5、end2 已知一離散序列為 x k=AcosW0k+Bcos ( (W0+DW)k)。用長度N=64的哈明窗對信號截短后近似計算其頻譜。試用不同的A和B的值(如 A和B近似相等，A和B差距較大)，確定用哈明窗能分辯的最小的譜峰間隔中c的值?！绢}目分析】本題討論用哈明窗計算序列頻譜時的頻率分辨率?！痉抡娼Y(jié)果】【結(jié)果分析】第一排是c=4，第二排c=3，遠大于教材定義的標準c=2。無論A/B的比值是多少，區(qū)分譜峰毫無壓力。第3排中，c=2是教材定義的標準，A與B接近的時候，區(qū)分較為容易，當A/B為1.6的時候，幾乎是無法區(qū)分，而A/B=4的時候，完全無法區(qū)分。所以說在設(shè)計使用哈明窗的時候，如果不同頻

6、譜的信號所占比例相差較大，哈明窗的c要比2大。從第一列可以看出，如果A/B象近的時候，即使c稍小于2，也是可以區(qū)分的。每一列中，A/B的比值是相同的，隨著c的減小，譜峰的區(qū)別能力很容易看出是降低的?！咀灾鲗W(xué)習(xí)內(nèi)容】不同窗函數(shù)的屏幕分辨率?！鹃喿x文獻】數(shù)字信號處理【發(fā)現(xiàn)問題】 (專題研討或相關(guān)知識點學(xué)習(xí)中發(fā)現(xiàn)的問題)：哈明窗所起的作用？哈明窗與矩形窗誰的分辨率高？在仿真的時候c該選用怎樣的值？【問題探究】哈明窗主瓣較寬，分辨率比矩形窗差，但是哈明窗可以使能量集中在主瓣處。所以哈明窗可以起到減小泄漏的作用。在仿真的時候c該選用怎樣的值？在開始的時候選用的c為 1 2 4 8 16，對比效果不明顯，

7、發(fā)現(xiàn)c的取值過于大，反復(fù)嘗試確定在2附近多取值，其他值也不可取得過于大，因為c稍大實驗現(xiàn)象基本相同，所以沒有意義。【仿真程序】N = 64;k = 0:N-1;c_set = 4 3 2.2 2 1.9 1.8 1.5;A_set = 1 1.2 1.6 4; B = 1;w0 = 0.5*pi;for c_index = 1:length(c_set) c = c_set(c_index); dw = c*2*pi/N; for A_index = 1:length(A_set) A = A_set(A_index); x = A*cos(w0*k)+B*cos(w0+dw)*k); Xh

8、= fft(x.*hamming(N)'); subplot(length(c_set),length(A_set),length(A_set)*(c_index-1)+A_index); plot(k,abs(Xh); hold on; title('c=' num2str(c) ' ;A/B=' num2str(A); axis(12 25 0 25); endend3 已知一離散序列為 xk=cos(W0k)+0.75cos(W1k), 0£ k £ 63 其中W0=0.4p, W1=W0+p/64(1) 對xk做64點FFT,

9、 畫出此時信號的頻譜。(2) 如果(1)中顯示的譜不能分辨兩個譜峰，是否可對(1)中的64點信號補零而分辨出兩個譜峰。通過編程進行證實，并解釋其原因 。(3) 給出一種能分辨出信號中兩個譜峰的計算方案，并進行仿真實驗。【題目分析】分析影響譜峰分辨率的主要因數(shù)，進一步認識補零在在頻譜計算中的作用?！痉抡娼Y(jié)果】【結(jié)果分析】頻譜的橫坐標是歸一化頻率，圖示只觀察要區(qū)別的部分。每一個列的抽樣點相同，每一行DFT點相同。但看第一列，雖然增加DFT點數(shù)，甚至是增加到32768，也不能夠區(qū)分兩個譜峰，所以要想增加頻率的分辨率，可以增加抽樣點數(shù)，第二列和第三列抽樣點數(shù)分別是128和256，都可以區(qū)分兩個譜峰。【

10、自主學(xué)習(xí)內(nèi)容】如何增加頻譜分辨率。【閱讀文獻】數(shù)字信號處理。【發(fā)現(xiàn)問題】 (專題研討或相關(guān)知識點學(xué)習(xí)中發(fā)現(xiàn)的問題)：如何計算所需抽樣點數(shù)？可以通過信號頻率的差值計算。【問題探究】1、2、3題討論的是離散信號頻譜的計算問題。與連續(xù)信號頻譜計算問題相比較，其計算誤差有何不同？連續(xù)信號的頻率是非周期的，離散信號的頻譜是連續(xù)信號頻譜的周期話，可能會有混疊誤差?！痉抡娉绦颉縲0 = 0.4*pi;dw = pi/64;w1 = w0+dw;N_set = 64 128 256;L_set = 64 128 256 512 32768;for N_index = 1:length(N_set); N =

11、N_set(N_index); k = 0:N-1; x = cos(w0*k)+0.75*cos(w1*k); for L_index =1:length(L_set) L = L_set(L_index); if L<N continue end X = fft(x,L); m = (0:L-1)*2/L; subplot(length(L_set),length(N_set),(L_index-1)*length(N_set)+N_index); plot(m,abs(X); axis(0.38 0.44 0 40); title('N=' num2str(N) &

12、#39; L=' num2str(L); endend4 試用DFT近似計算高斯信號的頻譜抽樣值。高斯信號頻譜的理論值為通過與理論值比較，討論信號的時域截取長度和抽樣頻率對計算誤差的影響?！绢}目分析】連續(xù)非周期信號頻譜計算的基本方法。計算中出現(xiàn)誤差的主要原因及減小誤差的方法?！痉抡娼Y(jié)果】【結(jié)果分析】如圖所示，橫軸坐標為角頻率，左斜線上的圖截取時間相同，垂線上的抽樣頻率相同，很容易能否看出，增加截取時間和抽樣頻率都可以減小誤差?！咀灾鲗W(xué)習(xí)內(nèi)容】連續(xù)非周期信號的頻譜計算方法?！痉抡娉绦颉縟 = pi;fsam_set = 1 2 4 8 16;N_set = 1 2 4 8;L = 64;

13、for fsam_index = 1:length(fsam_set); fsam = fsam_set(fsam_index); Tsam = 1/fsam; for N_index = 1:length(N_set); N = N_set(N_index); cuttime = N*Tsam; k = -N/2:N/2-1; t = k*Tsam; x = exp(-d*t.2); X = Tsam*fftshift(fft(x,L); w = linspace(-pi*fsam,pi*fsam,L); subplot(length(N_set),length(fsam_set),fsam

14、_index+length(fsam_set)*(N_index-1); plot(w,abs(X); hold on; G = exp(-w.*w/(4*d); plot(w,G,'r'); title('cuttime=' num2str(cuttime) ' fsam=' num2str(fsam); endend擴展題5 本題研究連續(xù)周期信號頻譜的近似計算問題。 周期為T0的連續(xù)時間周期信號x(t)可用Fourier級數(shù)表示為其中X(nw0)稱為連續(xù)時間周期信號x(t)的頻譜函數(shù)。稱為信號的基頻（基波），稱為信號的諧波。如果信號x(t)函

15、數(shù)表達式已知，則可由積分得出信號的頻譜。如果信號x(t)函數(shù)表達式未知，或者x(t)函數(shù)表達式非常復(fù)雜，則很難由積分得信號的頻譜。本題的目的就是研究如何利用DFT近似計算連續(xù)時間周期信號的頻譜。(1)若在信號x(t)的一個周期T0內(nèi)抽樣N個點，即， T為抽樣周期(間隔)，可獲得序列xk試分析序列xk的DFT與連續(xù)時間周期信號x(t)的頻譜X(nw0)的關(guān)系；(2)由(1)的結(jié)論，給出由DFT近似計算周期信號頻譜X(nw0)的方案；(3)周期信號x(t)的周期T0=1，x(t)在區(qū)間0，1的表達式為x(t)=20t2(1-t)4cos(12pt)(a)試畫出信號x(t)在區(qū)間0，1的波形；(b)

16、若要用6次以內(nèi)的諧波近似表示x(t)，試給出計算方案，并計算出近似表示的誤差。討論出現(xiàn)誤差的原因及減小誤差的方法。 【理論推導(dǎo)】DFT計算所得結(jié)果Xm與連續(xù)周期信號頻譜X(nw0)的關(guān)系。 由于,于是得令n=m+rN;m=0,1,2,，為整數(shù)，上式可化為化簡可得對比IDFT 可得 【計算方案】根據(jù)理論推導(dǎo)結(jié)果設(shè)計近似計算方案。分析產(chǎn)生誤差的主要原因。DFT是連續(xù)信號頻譜的周期化，如果頻譜信號高頻信號較多，就會產(chǎn)生混疊，所以有誤差。如果主要是低頻信號，周期話后互補產(chǎn)生影響，就沒有誤差?！緮U展分析】如果周期信號x(t)是帶限信號，即信號的最高頻率分量為Mw0(是正整數(shù))，試確定在一個周期內(nèi)的最少抽

17、樣點N，使得在頻譜的計算過程當中不存在混疊誤差。與抽樣定理給出的結(jié)論比較，發(fā)表你的看法。很容易看出，要想不產(chǎn)生混疊，必須N ³ 2M+1。該結(jié)論與抽樣定理相似。因為周期信號是一種特殊的連續(xù)信號，所以抽樣定理對周期信號也是適用的?！痉抡娼Y(jié)果】【結(jié)果分析】在N=16,k=6時，誤差最小，最大誤差為0.1229，平均誤差為0.07187。無論是減小N和減小k都很大的增大了誤差。為了分別測試N和k對誤差的影響，分別設(shè)k=6，求不同N時候的誤差，N=16，不同k對誤差的影響，分別得到了下圖（在原程序上稍加改動，所以在仿真程序中不再另寫了）此圖可以看出，在N等于13的時候誤差有明顯的減小。在k=

18、6的時候，誤差有明顯的減小?？梢娫瘮?shù)最大的高頻分量為6*w0，根據(jù)前面推出的理論，N取6*2+1=13點的時侯誤差減小。實驗結(jié)果符合結(jié)論【自主學(xué)習(xí)內(nèi)容】周期信號的頻譜與DFT的關(guān)系?！景l(fā)現(xiàn)問題】 ：原函數(shù)最高頻分量大約是多少？【問題探究】根據(jù)前面推理，最高頻分量為6*w0?！痉抡娉绦颉縩 = 1024;t = linspace(0,1,n);N_set = 12 16;H_set = 4 6;T0 = 1;count = 0;for N_index = 1:length(N_set) N = N_set(N_index); dT = T0/N; k = (0:N-1)*dT; x = 20*k.2.*(1-k).4.*cos(12*pi.*k); X = fft(x); for H_index = 1:length(H_set) H = H_set(H_index); sum = 0; X1 = X(N-H+1:N) X(1:H+1)/N; for k = -H:H;