離散傅里葉變換的應(yīng)用之一——用FFT分析信號(hào)頻譜

題目如下：

欲解此題，關(guān)鍵掌握以下幾點(diǎn)：

第一，分清“截取數(shù)據(jù)長(zhǎng)度”（即窗函數(shù)長(zhǎng)度）與“DFT點(diǎn)數(shù)”二者的不同；

第二，能夠根據(jù)模擬頻率推斷出DFT譜峰處對(duì)應(yīng)的序號(hào)k的數(shù)值，方法如下：

首先，由模擬頻率轉(zhuǎn)換為數(shù)字域頻率：

圖1

然后，數(shù)字域頻率對(duì)應(yīng)到DFT的序號(hào)k

圖2

綜合以上兩式，得到：

圖3

【題目分析與解答】

我們按照DFT分析信號(hào)頻譜的三個(gè)步驟來分別求解：

第一步：采樣。所以首先，我們寫出采樣后離散時(shí)間信號(hào)的表達(dá)式：

其周期為10。

第二步：時(shí)域加窗（即截取）

截取10點(diǎn)長(zhǎng)，相當(dāng)于將該周期信號(hào)x(n)與10點(diǎn)長(zhǎng)的矩形窗相乘，得到v(n)，所以我們求v(n)的DTFT。

圖4

圖5

注意，上圖中只畫出了[-Π，Π] 區(qū)間的圖形，實(shí)際上DTFT是以2Π為周期的（所以嚇人的公式中有西格瑪求和符號(hào)）。

第三步：頻域抽樣，也就是對(duì)V(e^jw)在 [0，2Π] 區(qū)間抽取N個(gè)點(diǎn)（N為DFT點(diǎn)數(shù)，而非第二步中截取的長(zhǎng)度）

再次把狐假虎威的公式擺出來（因?yàn)橛行┩瑢W(xué)有強(qiáng)迫癥，非要看看公式長(zhǎng)什么樣，注意，要先把圖4中的V(e^jw)的西格瑪求和符號(hào)去掉，只取0~2Π區(qū)間的，所以，后面一項(xiàng)是w-2Π/5+2Π，也就是w+8Π/5，然后再把w變成2Πk/N）（在4月28日發(fā)的文章中，下面這個(gè)公式寫錯(cuò)了，沒有人給我提出來。看來你們都沒仔細(xì)看公式。）

圖6

其實(shí)V(k) 就是下圖中的紅點(diǎn)點(diǎn)啦。

圖7

好了，那么最后的問題就是，N取不同值時(shí)（也就是做不同點(diǎn)數(shù)的DFT）這些紅點(diǎn)點(diǎn)顯然也不同。對(duì)于此題來說，這三種N的取值（10、20、128）得到的結(jié)果到底是什么呢？

你看到或者聽到這里的話，暫停一下，自己算算唄。

我直接把matlab畫圖的結(jié)果給出來。

用前面的公式算一下：f=kfs/N，最后那個(gè)128點(diǎn)DFT的圖，最高的譜峰序號(hào)k是多少？

圖8

最后附上matlab程序。

clc;clear all;

f0=1;fs=5;%單位：Hz

n=0:1000;L=10;

xn=cos(2*pi*f0*n/fs);%時(shí)域離散時(shí)間信號(hào)

Xk1=fft(xn(1:L),10);

Xk2=fft(xn(1:L),20);

Xk3=fft(xn(1:L),128);

subplot(311);stem((0:length(Xk1)-1),abs(Xk1));title('10點(diǎn)DFT');

subplot(312);stem((0:length(Xk2)-1),abs(Xk2));title('20點(diǎn)DFT');

subplot(313);stem((0:length(Xk3)-1),abs(Xk3));title('128點(diǎn)DFT');

下面分析一下：

首先看N=10時(shí)，此時(shí)的結(jié)果看似最為干凈清爽，只有干干凈凈兩根線。但有的同學(xué)要問了，單頻信號(hào)，只有一個(gè)頻率成分，應(yīng)該只有一根譜線呀？為什么會(huì)有兩根？

我們先看第一根，k=2那個(gè)譜線，對(duì)應(yīng)頻率為2*fs/10=2*5/10=1Hz，與題設(shè)cos（2Πt）完全吻合。k=8那根譜線是怎么一回事呢？是負(fù)頻率周期延拓過去的，本來在-2，-2+10就等于8了。所以，（敲黑板，以下結(jié)論很重要）

對(duì)實(shí)信號(hào)做N點(diǎn)DFT，我們只需要看前N/2根譜線就行了，不用關(guān)注N/2~N-1之間的。

再看N=20和N=128的DFT結(jié)果，怎么出來那么多根譜線呢？

回過頭去看一看，做DFT之前的截取L點(diǎn)長(zhǎng)的序列cos(2Πn/5)的頻譜到底是什么樣子呢？是圖7中的虛線所示。而N點(diǎn)DFT，是對(duì)V(e^jw)在 [0，2Π] 區(qū)間抽取N個(gè)點(diǎn)。相當(dāng)于把連續(xù)的頻譜圖（V(e^jw)，如圖7中的虛線所示），用一張不透明的紙蓋住，紙上以2Π/N為間隔開了一些縫，露出來的點(diǎn)才是我們得到的DFT的結(jié)果。這就是頻域抽樣產(chǎn)生的”柵欄效應(yīng)“。

所以，不管是10點(diǎn)DFT干干凈凈的兩根線，還是128點(diǎn)DFT密密麻麻的那么多根線，背后隱藏的，都是連續(xù)的頻譜函數(shù)。之所以10點(diǎn)DFT的結(jié)果看起來更順眼，無非是因?yàn)橐驗(yàn)?Π/N也恰好是旁瓣的寬度（因?yàn)樾盘?hào)的周期和截取長(zhǎng)度也是10）

在DFT譜分析中，當(dāng)DFT點(diǎn)數(shù)N大于數(shù)據(jù)本身的實(shí)際點(diǎn)數(shù)L時(shí)，相當(dāng)于在數(shù)據(jù)后面補(bǔ)上了L-N個(gè)0再做DFT，稱為“補(bǔ)零DFT“。補(bǔ)零經(jīng)常是必要的，補(bǔ)零相當(dāng)于對(duì)信號(hào)頻譜以更小的間隔采樣，得到更多的頻譜的信息。而且有時(shí)數(shù)據(jù)長(zhǎng)度不是2的整數(shù)次冪，如果我們想采用基2FFT算法，就必須進(jìn)行補(bǔ)零。

此題以單頻周期信號(hào)為例，展示了不同點(diǎn)數(shù)DFT時(shí)結(jié)果的不同。給人一種錯(cuò)覺：補(bǔ)零，似乎沒帶來任何好處，反倒是點(diǎn)數(shù)最少的10點(diǎn)DFT的結(jié)果最好看。

是不是這樣呢？當(dāng)然了，如果你知道這個(gè)周期信號(hào)的周期是多少，毫無疑問，就截取一個(gè)周期的數(shù)據(jù)，做同樣點(diǎn)數(shù)的DFT，結(jié)果最好看。但在分析實(shí)際問題時(shí)，要么信號(hào)根本沒有周期性；要么雖然有周期性，但你不知道；而且實(shí)際信號(hào)也不會(huì)是簡(jiǎn)單的單頻信號(hào)，會(huì)包含多個(gè)頻率分量。這個(gè)時(shí)候，在系統(tǒng)性能、實(shí)時(shí)性、存儲(chǔ)量等等容許的范圍內(nèi)，截取盡量長(zhǎng)的數(shù)據(jù)（即獲取更多的信息，得到更高的頻率分辨率，下一篇會(huì)專門講頻率分辨率），做盡量多點(diǎn)數(shù)的DFT（即對(duì)頻譜進(jìn)行更為精細(xì)的采樣）。

而且，實(shí)際應(yīng)用中，由于DFT的點(diǎn)數(shù)一般都比較大，我們一般不會(huì)以離散的形式畫頻譜圖，而是直接將頻譜圖化成連續(xù)的曲線。例如，上例中，我們截取128點(diǎn)長(zhǎng)的數(shù)據(jù)，做128點(diǎn)DFT，用連續(xù)曲線形式畫圖（matlab中為plot函數(shù)），并且只畫出前一半（即0~N/2-1）的點(diǎn)，并且把橫軸直接轉(zhuǎn)換為Hz，如下圖所示：

Matlab代碼如下：

上圖中，橫軸單位為：Hz。采用如下公式，將序號(hào)k轉(zhuǎn)換為模擬頻率Hz：

DFT分析信號(hào)頻譜，是實(shí)際中應(yīng)用最廣泛的數(shù)字信號(hào)處理算法，還有很多種題目可以出。還是那句話，題目無窮無盡，而原理就那么多，大家只有掌握了其真正含義，才能以不變應(yīng)萬變。