FFT頻譜分析原理

取樣定理：取樣頻率要大於訊號頻率的兩倍。

N個取樣點經過FFT變換後得到N個點的以複數形式記錄的FFT結果。

假設取樣頻率為Fs，取樣點數為N。那麼FFT運算的結果就是N個複數（或N個點），每一個複數就對應著一個頻率值以及該頻率訊號的幅值和相位。第一個點對應的頻率為0Hz（即直流分量），最後一個點N的下一個點對應取樣頻率Fs。其中任意一個取樣點n所代表的訊號頻率：

Fn=(n-1)*Fs/N。

這表明，頻譜分析得到的訊號頻率最大為(N-1)*Fs/N,對頻率的分辨能力是Fs/N。取樣頻率和取樣時間制約著通過FFT運算能分析得到的訊號頻率上限，同時也限定了分析得到的訊號頻率的解析度。

每一個複數的模值對應該點所對應的頻率值的幅度特性，具體的定量關係如下：

假設訊號由以下週期的原始訊號疊加而成：

Y = A1 + A2 Cos (2*PI*ω2*t + φ2 * PI/180) + A3 Cos (2*PI*ω3*t + φ3 * PI/180)

那麼，在經過FFT分析後得到的第一個點的模值是A1的N倍，而且只有在FFT結果點對應的頻率在ω2,ω3時，其模值才明顯放大，在其他頻率點，模值接近於0。在這些模值明顯放大的點中，除第一個點之外的其它點模值是相應訊號幅值的N/2倍。

每個複數的相位就是在該頻率值下訊號的相位：φ2，φ3。

FFT結果有對稱性，通常我們只是用前半部分的結果，也就是小於取樣頻率一半的結果。同時也只有取樣頻率一半以內、具有一定幅值的訊號頻率才是真正的訊號頻率。

Python實踐FFT頻譜分析

假如訊號S是有1個直流訊號和4個週期訊號疊加而成，如下公式所列（t為自變數，pi為圓周率值）現要對其進行FFT分析並繪製頻譜圖。

S =  2.0  + 3.0 * cos(2.0 * pi * 50 * t - pi * 30/180)
      + 1.5 * cos(2.0 * pi * 75 * t + pi * 90/180)
      +  1.0 * cos(2.0 * pi * 150 * t + pi * 120/180)
      +  2.0 * cos(2.0 * pi * 220 * t + pi * 30/180)

我們先使用Python繪製其1秒內的波形圖：

import numpy as np
　　import pylab as pl
　　import math
　　# 取樣步長
　　t = [x/1048.0 for x in range(1048)]
　　# 設計的取樣值
　　y = [2.0 + 3.0 * math.cos(2.0 * math.pi * 50 * t0 - math.pi * 30/180)
      　　　　+ 1.5 * math.cos(2.0 * math.pi * 75 * t0 + math.pi * 90/180)
      　　　　+  1.0 * math.cos(2.0 * math.pi * 150 * t0 + math.pi * 120/180)
      　　　　+  2.0 * math.cos(2.0 * math.pi * 220 * t0 + math.pi * 30/180)
      　　　　for t0 in t ]
　　pl.plot(t,y)
　　pl.show()

波形圖如圖1所示。

圖 1 訊號S在1秒內的波形

現在我們要對該訊號在0-1秒時間內進行頻譜分析，我們在這1秒時間內取樣1048點，則我們的取樣頻率為1048Hz。這兩個設定決定了我們頻譜分析所能得到的最高訊號頻率是1047Hz (1048Hz*(1048-1)/1048)，但只有小於524Hz的頻率才可能是真實的訊號頻率。我們能分辨的最小頻率差為1Hz（1048Hz/1048)。輸入如下的python程式碼：

N=len(t)    # 取樣點數
　　fs=1048.0     # 取樣頻率
　　df = fs/(N-1)   # 解析度
　　f = [df*n for n in range(0,N)]   # 構建頻率陣列

使用下面的程式碼進行快速傅立葉變換並對結果模值進行縮放：

Y = np.fft.fft(y)*2/N  #*2/N 反映了FFT變換的結果與實際訊號幅值之間的關係
　　absY = [np.abs(x) for x in Y]      #求傅立葉變換結果的模

此時我們得到了訊號的FFT分析結果，它儲存在列表Y內，同時我們也獲取了Y內每一個元素（是複數）的模組成的列表absY。根據理論分析，absY內的元素大部分都極其接近0，只有極少數峰值提示的是訊號頻率的幅度。我們寫一小段程式碼顯示幅度較大的訊號的資訊。程式碼如下：

i=0
　　while i < len(absY):
       　　 if absY[i]>0.01:
              　　   print("freq:M, Y: %5.2f + %5.2f j, A:%3.2f, phi: %6.1f"

　　　　　%(i, Y[i].real, Y[i].imag, absY[i],math.atan2(Y[i].imag,Y[i].real)*180/math.pi))
    　　　i+=1

該段程式碼的執行結果如下：

freq:   0, Y:  4.00 +  0.00 j, A:4.00, phi:    0.0

freq:  50, Y:  2.60 + -1.50 j, A:3.00, phi:  -30.0

freq:  75, Y:  0.00 +  1.50 j, A:1.50, phi:   90.0

freq: 150, Y: -0.50 +  0.87 j, A:1.00, phi:  120.0

freq: 220, Y:  1.73 +  1.00 j, A:2.00, phi:   30.0

freq: 828, Y:  1.73 + -1.00 j, A:2.00, phi:  -30.0

freq: 898, Y: -0.50 + -0.87 j, A:1.00, phi: -120.0

freq: 973, Y:  0.00 + -1.50 j, A:1.50, phi:  -90.0

freq: 998, Y:  2.60 +  1.50 j, A:3.00, phi:   30.0

該段結果很好的反映了原始訊號的頻率、幅度和相位，我們在程式碼中已經對FFT結果與訊號幅值做了一個換算2/N，但在0Hz處，實際的訊號應該是FFT結果的1/N，也就是A：4.00再除以2。

使用如下程式碼繪製頻譜圖，得到圖2：

pl.plot(f,absY)
　　pl.show()

 
圖 2 訊號S的頻譜分析圖

　由圖2可見，頻譜分析的結果具有對稱性，同220Hz對稱的頻率為828Hz，可以計算出對稱軸為524Hz。但只有小於524Hz的頻率才是訊號頻率。我們通過修改取樣頻率和取樣點數，我們發現作為對稱軸的頻率會變化，同樣的對稱軸右側的訊號頻率也會發生變化。提示右側的頻率沒有實際意義，這也從一個方面印證了取樣定理。

本文大量參考了博文：利用matlab怎樣進行頻譜分析 在此對博主表示感謝，原文連結：http://blog.sina.com.cn/s/blog_a07f4fe301013gj3.html.