前面寫過一篇傅里葉變換

的文章:

furious:傅里葉變換學習心得?

zhuanlan.zhihu.com
圖標

但是在工程應用中,得益於數字技術的應用,絕大多數傅里葉變換的應用都是採用離散傅里葉變換(DFT),更確切的說,是它的快速演算法FFT。這篇文章再來寫寫有關離散傅里葉變換的關鍵點。

閑言少敘,直入主題。先把DFT的式子寫在這裡:

X[k]=sum_{n=0}^{N-1}{x[n]e^{-j(2pi/N)kn}} ,其中(0leq k leq N-1)

原信號 x(t) 的採樣信號 x[n] 可以用 X[k] 表示為:

x[n]=frac{1}{N}sum_{k=0}^{N-1}{X[k]e^{j(2pi/N)kn}} ,其中(0leq n leq N-1) 。這也叫DFT的逆變換式,實際的含義就是x[n]表示成X[k]為係數的不同頻率分量的和。

對離散傅里葉變換變換,我認為最重要的是搞清楚兩點:

1、實數信號變換的結果 X[k] 是一組複數,裡面一半數據和另一半是共軛的(這個好理解),意味著N點DFT,只有N/2的數據是含有有用信息的。那麼這些複數數據里的模和幅角是什麼含義?

2、用DFT的結果如何做頻譜分析,即在採樣頻率 為f_{s}的情況下, x[n] 的n只是一個離散的數值,那每一點到底代表了什麼頻率分量?

接下來就探討探討這兩個問題(推導很簡單,就是三角函數的基本運算)。

1、先來看 X[k] 的共軛性質,X[N-k]=sum_{n=0}^{N-1}{x[n]e^{-j(2pi/N)(N-k)n}} =sum_{n=0}^{N-1}{x[n]e^{j(-2npi+(2pi/N)kn)}} = sum_{n=0}^{N-1}{x[n]e^{j((2pi/N)kn)}} = X^{*}[k]

再看 X[0] 顯然沒有虛部,叫做直流分量。當N是偶數的時候, X[frac{N}{2}]=sum_{n=0}^{N-1}{x[n]e^{-j(2pi/N)nN/2}} =sum_{n=0}^{N-1}{x[n]e^{-jnpi}} ,虛部是0,也是直流分量。

再從 x[n] 出發(假設N是偶數),

N	imes x[n]=sum_{k=0}^{N-1}{X[k]e^{j(2pi/N)kn}} =X[0]+X[1]e^{j(2npi/N)}+X[2]e^{j(2(2npi/N))} +....

+ X[N/2-1]e^{j(frac{2npi}{N}(N/2-1))} + X[N/2] +X[N/2+1]e^{j(frac{2npi}{N}(N/2+1))}+...

+ X[N-2]e^{j(frac{2npi}{N}(N-2))}  +X[N-1]e^{j(frac{2npi}{N}(N-1))}

直流分量不管,只考察其中 X[k_{c}] X[N-k_{c}] 這兩項頻率分量,假設X[k_{c}] =a_{k}+jb_{k},那麼由上面的共軛特性可知X[N-k_{c}]=a_{k}-jb_{k}

那麼 X[k_{c}] X[N-k_{c}] 這兩項分量之和就 =(a_{k}+jb_{k})e^{j(2pi/N)k_{c}n} + (a_{k}-jb_{k})e^{j(2pi/N)(N-k_{c})n}

=2a_{k}cos(frac{2k_{c}npi}{N})-2b_{k}sin(frac{2k_{c}npi}{N}) ,虛數部分消去了!

=2sqrt{a_{k}^{2}+b_{k}^{2}}cos(frac{2k_{c}npi}{N}+phi_{k}) ,其中 A_{k} = 2sqrt{a_{k}^{2}+b_{k}^{2}},而 tan(phi_{k})=frac{b_{k}}{a_{k}}

所以x[n]的就可以表示成[1, 2, ..., N/2-1]這些三角函數表示的頻率分量與直流項 a_{0}= (X[0]+X[frac{N}{2}]) 的和:

N	imes x[n] = a_{0}+ A_{1}cos(n2pi/N+phi_{1}) + A_{2}cos((2n)2pi/N+phi_{2}) +

 + A_{3}cos((3n)2pi/N+phi_{3}) + ... +  A_{frac{N}{2}-1}cos(((frac{N}{2}-1)n)2pi/N+phi_{frac{N}{2}-1})

這樣,x[n]的頻率分量組成,以及每個分量的強度與初始相位與X[k]的模與幅角的關係一目了然。這也就回答了問題1。

2、假設採樣頻率是 f_{s} ,則 x[n] = x(nT_{s})=x(frac{n}{f_{s}}) ,所以事實上把 x[n] 的式子裡面的n用 f_{s}t 替換便可以得到原始的信號 x(t) ,即

N	imes x(t)= a_{0}+ A_{1}cos(f_{s}2pi t/N+phi_{1}) + A_{2}cos((2f_{s})2pi t/N+phi_{2}) +

 + A_{3}cos((3f_{s})2pi t/N+phi_{3}) + ... +  A_{frac{N}{2}-1}cos(((frac{N}{2}-1)f_{s})2pi t/N+phi_{frac{N}{2}-1})

很明顯,在 k_{c} 點的頻率分量的頻率為 2pi k_{c} f_{s}/N = 2pi f從而得出 X[k_{c}] 代表的是頻率為f= k_{c} f_{s}/N 的頻率「分量」。這就回答了問題2。

如前所述,X[k]只有[0, 1, 2, ..., N/2-1]這些項的頻率就夠了,後續的點在頻率信息上是冗餘的。那麼從另一個角度去考慮,採樣定理決定了f_{s}最多只能不失真地採樣最大f_{s}/2的頻率, N/2-1 這個點剛好到了上限了,神奇吧?

最後結合一個實例來加深上述對離散傅里葉變換的結果的理解。

假設原始信號:

x(t)=sin(2pi	imes40t)+0.5cos(2pi	imes90t) ,即40Hz和90Hz兩個正餘弦信號的疊加。來看看不同採樣過程得到的DFT(FFT)結果。

1、在0~0.5s內採樣500點得到 x[n] ,即採樣頻率為 500/0.5s=1000Hz

x(t)
FFT結果

Value at index 0: (4.618527782440651e-14+0j)
Value at index 1: (0.00038048349283981153-0.06055503176189958j)
Value at index 499: (0.00038048349283897887+0.06055503176189887j)
Value at index 2: (0.0015317714831372398-0.12188808528069522j)
Value at index 498: (0.001531771483137656+0.12188808528069539j)
Value at index 3: (0.003484103499379887-0.1848155396034572j)
Value at index 497: (0.003484103499377167+0.18481553960345828j)
Value at index 19: (1.1433174456370772-9.531545524428282j)
Value at index 481: (1.1433174456370745+9.531545524428282j)
Value at index 20: (31.252728294296478-247.3908181827295j)
Value at index 480: (31.252728294296475+247.3908181827295j)
Value at index 21: (-1.363915207723784+10.276792004132425j)
Value at index 479: (-1.3639152077237846-10.276792004132425j)
Value at index 44: (2.5907927302365366-9.131330368139245j)
Value at index 456: (2.5907927302365366+9.131330368139245j)
Value at index 45: (34.21663413671753-117.77442719621841j)
Value at index 455: (34.21663413671753+117.7744271962184j)
Value at index 46: (-3.6756991040418767+12.36140199986314j)
Value at index 454: (-3.675699104041877-12.36140199986314j)

根據 f= k_{c} f_{s}/N ,40Hz對應的點應該是40x500/1000 = 20,90Hz對應的點是90x500/1000 = 45。仔細看上面的結果,峰值就出現在20和45這兩個X[k]上,與上面的結果是相符的。

2、在0~0.5s內採樣200點。即採樣頻率為 200/0.5s=400Hz

Value at index 0: (6.372680161348399e-14+0j)
Value at index 1: (0.0009006884640155866-0.05733489242258602j)
Value at index 199: (0.000900688464015087+0.057334892422586464j)
Value at index 2: (0.0036263425132948557-0.11539208979600157j)
Value at index 198: (0.0036263425132945226+0.11539208979600174j)
Value at index 3: (0.008249466986120169-0.1749295249833991j)
Value at index 197: (0.008249466986122056+0.17492952498339767j)
Value at index 19: (2.624895912836278-8.532357512798926j)
Value at index 181: (2.6248959128362785+8.532357512798928j)
Value at index 20: (30.380956546694442-93.50296980739842j)
Value at index 180: (30.38095654669445+93.50296980739842j)
Value at index 21: (-3.5621007673958807+10.40404496136136j)
Value at index 179: (-3.5621007673958793-10.404044961361361j)
Value at index 44: (4.903695343956267-5.927549420513547j)
Value at index 156: (4.903695343956271+5.927549420513546j)
Value at index 45: (29.319402903322846-34.32861016803424j)
Value at index 155: (29.319402903322846+34.328610168034245j)
Value at index 46: (-9.357297882864374+10.613771038768512j)
Value at index 154: (-9.357297882864373-10.613771038768512j)

同樣根據 f= k_{c} f_{s}/N ,40Hz對應的點應該是40x200/400 = 20,90Hz對應的點是90x200/400 = 45。與上面的結果是相符的。

3、在0~0.2s內採樣200點。即採樣頻率還是 200/0.2s=1000Hz

Value at index 0: (-6.483702463810914e-14+0j)
Value at index 1: (0.0023975099150431722-0.15261766809813948j)
Value at index 199: (0.0023975099150432833+0.1526176680981392j)
Value at index 2: (0.009996346441271342-0.3180889014147228j)
Value at index 198: (0.009996346441271842+0.3180889014147237j)
Value at index 3: (0.0242274834823393-0.5137425465454714j)
Value at index 197: (0.024227483482339385+0.5137425465454717j)
Value at index 7: (0.41780330667369625-3.7844150089860715j)
Value at index 193: (0.41780330667369603+3.784415008986071j)
Value at index 8: (12.501978196875637-98.96334764449756j)
Value at index 192: (12.50197819687564+98.96334764449756j)
Value at index 9: (-0.5766958686580562+4.052076375966517j)
Value at index 191: (-0.5766958686580559-4.052076375966516j)
Value at index 17: (0.8826413246935703-3.226393289054478j)
Value at index 183: (0.8826413246935703+3.2263932890544775j)
Value at index 18: (13.569220126648087-46.70556202371618j)
Value at index 182: (13.56922012664809+46.70556202371618j)
Value at index 19: (-1.633027099955663+5.308237548305262j)
Value at index 181: (-1.633027099955663-5.308237548305263j)

同樣根據 f= k_{c} f_{s}/N ,40Hz對應的點應該是40x200/1000 = 8,90Hz對應的點是90x200/1000 = 18。與上面的結果是相符的。

文章就到這裡。限於水平,如果有紕漏之處,請不吝指出。


推薦閱讀:
相关文章