目錄
一、一些關鍵概念的引入
1、離散傅里葉變換（DFT）
2、快速傅里葉變換（FFT）
3、採樣頻率以及採樣定理
4、如何理解採樣定理？
二、使用scipy包實現快速傅里葉變換
1、產生原始信號——原始信號是三個正弦波的疊加
2、快速傅里葉變換
3、FFT的原始頻譜
4、將振幅譜進行歸一化和取半處理
三、完整代碼
一、一些關鍵概念的引入
1、離散傅里葉變換（DFT）
離散傅里葉變換(discrete Fourier transform) 傅里葉分析方法是信號分析的最基本方法，傅里葉變換是傅里葉分析的核心，通過它把信號從時間域變換到頻率域，進而研究信號的頻譜結構和變化規律。但是它的致命缺點是：計算量太大，時間複雜度太高，當採樣點數太高的時候，計算緩慢，由此出現了DFT的快速實現，即下面的快速傅里葉變換FFT。
2、快速傅里葉變換（FFT）
計算量更小的離散傅里葉的一種實現方法。詳細細節這裡不做描述。
3、採樣頻率以及採樣定理
採樣頻率：採樣頻率，也稱為採樣速度或者採樣率，定義了每秒從連續信號中提取並組成離散信號的採樣個數，它用赫茲（Hz）來表示。採樣頻率的倒數是採樣週期或者叫作採樣時間，它是採樣之間的時間間隔。通俗的講採樣頻率是指計算機每秒鐘採集多少個信號樣本。
採樣定理：所謂採樣定理 ，又稱香農採樣定理，奈奎斯特採樣定理，是信息論，特別是通訊與信號處理學科中的一個重要基本結論。採樣定理指出，如果信號是帶限的，並且採樣頻率高於信號帶寬的兩倍，那麼，原來的連續信號可以從採樣樣本中完全重建出來。
定理的具體表述為：在進行模擬/數字信號的轉換過程中，當採樣頻率fs大於信號中最高頻率fmax的2倍時,即
fs>2*fmax
採樣之後的數字信號完整地保留了原始信號中的信息，一般實際應用中保證採樣頻率為信號最高頻率的2.56～4倍；採樣定理又稱奈奎斯特定理。
4、如何理解採樣定理？
在對連續信號進行離散化的過程中，難免會損失很多信息，就拿一個簡單地正弦波而言，如果我1秒內就選擇一個點，很顯然，損失的信號太多了，光著一個點我根本不知道這個正弦信號到底是什麼樣子的，自然也沒有辦法根據這一個採樣點進行正弦波的還原，很明顯，我採樣的點越密集，那越接近原來的正弦波原始的樣子，自然損失的信息越少，越方便還原正弦波。故而
採樣定理說明採樣頻率與信號頻率之間的關係，是連續信號離散化的基本依據。 它為採樣率建立了一個足夠的條件，該採樣率允許離散採樣序列從有限帶寬的連續時間信號中捕獲所有信息。
二、使用scipy包實現快速傅里葉變換
本節不會說明FFT的底層實現，只介紹scipy中fft的函數接口以及使用的一些細節。
1、產生原始信號——原始信號是三個正弦波的疊加
import numpy as np
from scipy.fftpack import fft,ifft
import matplotlib.pyplot as plt
from matplotlib.pylab import mpl
mpl.rcParams['font.sans-serif'] = ['SimHei'] #顯示中文
mpl.rcParams['axes.unicode_minus']=False #顯示負號
#採樣點選擇1400個，因為設置的信號頻率分量最高為600赫茲，根據採樣定理知採樣頻率要大於信號頻率2倍，所以這裡設置採樣頻率為1400赫茲（即一秒內有1400個採樣點，一樣意思的）
x=np.linspace(0,1,1400)
#設置需要採樣的信號，頻率分量有200，400和600
y=7*np.sin(2*np.pi*200*x) + 5*np.sin(2*np.pi*400*x)+3*np.sin(2*np.pi*600*x)
這裡原始信號的三個正弦波的頻率分別為，200Hz、400Hz、600Hz,最大頻率為600赫茲。根據採樣定理，fs至少是600赫茲的2倍，這裡選擇1400赫茲，即在一秒內選擇1400個點。
原始的函數圖像如下：
plt.figure()
plt.plot(x,y)
plt.title('原始波形')
plt.figure()
plt.plot(x[0:50],y[0:50])
plt.title('原始部分波形（前50組樣本）')
plt.show()
由圖可見，由於採樣點太過密集，看起來不好看，我們只顯示前面的50組數據，如下：
2、快速傅里葉變換
其實scipy和numpy一樣，實現FFT非常簡單，僅僅是一句話而已，函數接口如下：
from scipy.fftpack import fft,ifft
from numpy import fft,ifft
其中fft表示快速傅里葉變換，ifft表示其逆變換。具體實現如下：
fft_y=fft(y) #快速傅里葉變換
print(len(fft_y))
print(fft_y[0:5])
'''運行結果如下：
1400
[-4.18864943e-12+0.j 9.66210986e-05-0.04305756j 3.86508070e-04-0.08611996j
8.69732036e-04-0.12919206j 1.54641157e-03-0.17227871j]
'''
我們發現以下幾個特點：
（1）變換之後的結果數據長度和原始採樣信號是一樣的
（2）每一個變換之後的值是一個複數，為a+bj的形式，那這個複數是什麼意思呢？
我們知道，複數a+bj在座標系中表示為（a,b），故而複數具有模和角度，我們都知道快速傅里葉變換具有
“振幅譜”“相位譜”，它其實就是通過對快速傅里葉變換得到的複數結果進一步求出來的，
那這個直接變換後的結果是不是就是我需要的，當然是需要的，在FFT中，得到的結果是複數，
（3）FFT得到的複數的模（即絕對值）就是對應的“振幅譜”，複數所對應的角度，就是所對應的“相位譜”，現在可以畫圖了。
3、FFT的原始頻譜
N=1400
x = np.arange(N) # 頻率個數
abs_y=np.abs(fft_y) # 取複數的絕對值，即複數的模(雙邊頻譜)
angle_y=np.angle(fft_y) #取複數的角度
plt.figure()
plt.plot(x,abs_y)
plt.title('雙邊振幅譜（未歸一化）')
plt.figure()
plt.plot(x,angle_y)
plt.title('雙邊相位譜（未歸一化）')
plt.show()
顯示結果如下：
注意：我們在此處僅僅考慮“振幅譜”，不再考慮相位譜。
我們發現，振幅譜的縱座標很大，而且具有對稱性，這是怎麼一回事呢？
關鍵：關於振幅值很大的解釋以及解決辦法——歸一化和取一半處理
比如有一個信號如下：
Y=A1+A2*cos(2πω2+φ2）+A3*cos(2πω3+φ3）+A4*cos(2πω4+φ4）
經過FFT之後，得到的“振幅圖”中，
第一個峰值（頻率位置）的模是A1的N倍，N為採樣點，本例中為N=1400，此例中沒有，因為信號沒有常數項A1
第二個峰值（頻率位置）的模是A2的N/2倍，N為採樣點，
第三個峰值（頻率位置）的模是A3的N/2倍，N為採樣點，
第四個峰值（頻率位置）的模是A4的N/2倍，N為採樣點，
依次下去......
考慮到數量級較大，一般進行歸一化處理，既然第一個峰值是A1的N倍，那麼將每一個振幅值都除以N即可
FFT具有對稱性，一般只需要用N的一半，前半部分即可。
4、將振幅譜進行歸一化和取半處理
先進行歸一化
normalization_y=abs_y/N #歸一化處理（雙邊頻譜）
plt.figure()
plt.plot(x,normalization_y,'g')
plt.title('雙邊頻譜(歸一化)',fontsize=9,color='green')
plt.show()
結果為：
現在我們發現，振幅譜的數量級不大了，變得合理了，接下來進行取半處理：
half_x = x[range(int(N/2))] #取一半區間
normalization_half_y = normalization_y[range(int(N/2))] #由於對稱性，只取一半區間（單邊頻譜）
plt.figure()
plt.plot(half_x,normalization_half_y,'b')
plt.title('單邊頻譜(歸一化)',fontsize=9,color='blue')
plt.show()
這就是我們最終的結果，需要的
“振幅譜”。三、完整代碼
import numpy as np
from scipy.fftpack import fft,ifft
import matplotlib.pyplot as plt
from matplotlib.pylab import mpl
mpl.rcParams['font.sans-serif'] = ['SimHei'] #顯示中文
mpl.rcParams['axes.unicode_minus']=False #顯示負號
#採樣點選擇1400個，因為設置的信號頻率分量最高為600赫茲，根據採樣定理知採樣頻率要大於信號頻率2倍，所以這裡設置採樣頻率為1400赫茲（即一秒內有1400個採樣點，一樣意思的）
x=np.linspace(0,1,1400)
#設置需要採樣的信號，頻率分量有200，400和600
y=7*np.sin(2*np.pi*200*x) + 5*np.sin(2*np.pi*400*x)+3*np.sin(2*np.pi*600*x)
fft_y=fft(y) #快速傅里葉變換
N=1400
x = np.arange(N) # 頻率個數
half_x = x[range(int(N/2))] #取一半區間
abs_y=np.abs(fft_y) # 取複數的絕對值，即複數的模(雙邊頻譜)
angle_y=np.angle(fft_y) #取複數的角度
normalization_y=abs_y/N #歸一化處理（雙邊頻譜）
normalization_half_y = normalization_y[range(int(N/2))] #由於對稱性，只取一半區間（單邊頻譜）
plt.subplot(231)
plt.plot(x,y)
plt.title('原始波形')
plt.subplot(232)
plt.plot(x,fft_y,'black')
plt.title('雙邊振幅譜(未求振幅絕對值)',fontsize=9,color='black')
plt.subplot(233)
plt.plot(x,abs_y,'r')
plt.title('雙邊振幅譜(未歸一化)',fontsize=9,color='red')
plt.subplot(234)
plt.plot(x,angle_y,'violet')
plt.title('雙邊相位譜(未歸一化)',fontsize=9,color='violet')
plt.subplot(235)
plt.plot(x,normalization_y,'g')
plt.title('雙邊振幅譜(歸一化)',fontsize=9,color='green')
plt.subplot(236)
plt.plot(half_x,normalization_half_y,'blue')
plt.title('單邊振幅譜(歸一化)',fontsize=9,color='blue')
plt.show()
閱讀更多 馬士兵尚學堂 的文章
