說明：本文適合信號處理方面有一定的基礎的人閱讀琅捏，能夠理解什么時候傅里葉級數(shù)和傅里葉變換，能夠理解他們的核心思想以及基本原理递雀，能夠理解到底什么是“頻率域”柄延，能夠從頻率的角度分析信號。

一缀程、一些關鍵概念的引入

1搜吧、離散傅里葉變換（DFT)

離散傅里葉變換(discrete Fourier transform) 傅里葉分析方法是信號分析的最基本方法，傅里葉變換是傅里葉分析的核心杨凑，通過它把信號從時間域變換到頻率域滤奈，進而研究信號的頻譜結(jié)構(gòu)和變化規(guī)律。但是它的致命缺點是：計算量太大撩满，時間復雜度太高蜒程，當采樣點數(shù)太高的時候绅你，計算緩慢，由此出現(xiàn)了DFT的快速實現(xiàn)昭躺，即下面的快速傅里葉變換FFT忌锯。

2、快速傅里葉變換（FFT）

計算量更小的離散傅里葉的一種實現(xiàn)方法领炫。詳細細節(jié)這里不做描述偶垮。

3、采樣頻率以及采樣定理

采樣頻率：采樣頻率帝洪，也稱為采樣速度或者采樣率似舵，定義了每秒從連續(xù)信號中提取并組成離散信號的采樣個數(shù)，它用赫茲（Hz）來表示碟狞。采樣頻率的倒數(shù)是采樣周期或者叫作采樣時間啄枕，它是采樣之間的時間間隔婚陪。通俗的講采樣頻率是指計算機每秒鐘采集多少個信號樣本族沃。

采樣定理：所謂采樣定理 ，又稱香農(nóng)采樣定理泌参，奈奎斯特采樣定理脆淹，是信息論，特別是通訊與信號處理學科中的一個重要基本結(jié)論沽一。采樣定理指出盖溺，如果信號是帶限的，并且采樣頻率高于信號帶寬的兩倍铣缠，那么烘嘱，原來的連續(xù)信號可以從采樣樣本中完全重建出來。

定理的具體表述為：在進行模擬/數(shù)字信號的轉(zhuǎn)換過程中蝗蛙，當采樣頻率fs大于信號中最高頻率fmax的2倍時,即

fs>2*fmax

采樣之后的數(shù)字信號完整地保留了原始信號中的信息蝇庭，一般實際應用中保證采樣頻率為信號最高頻率的2.56～4倍；采樣定理又稱奈奎斯特定理捡硅。

4哮内、如何理解采樣定理？

在對連續(xù)信號進行離散化的過程中壮韭，難免會損失很多信息北发，就拿一個簡單地正弦波而言，如果我1秒內(nèi)就選擇一個點喷屋，很顯然琳拨，損失的信號太多了，光著一個點我根本不知道這個正弦信號到底是什么樣子的屯曹，自然也沒有辦法根據(jù)這一個采樣點進行正弦波的還原狱庇，很明顯寄疏，我采樣的點越密集，那越接近原來的正弦波原始的樣子僵井，自然損失的信息越少陕截，越方便還原正弦波。故而

采樣定理說明采樣頻率與信號頻率之間的關系批什，是連續(xù)信號離散化的基本依據(jù)农曲。 它為采樣率建立了一個足夠的條件，該采樣率允許離散采樣序列從有限帶寬的連續(xù)時間信號中捕獲所有信息驻债。

二乳规、使用scipy包實現(xiàn)快速傅里葉變換

本節(jié)不會說明FFT的底層實現(xiàn)，只介紹scipy中fft的函數(shù)接口以及使用的一些細節(jié)合呐。

1暮的、產(chǎn)生原始信號——原始信號是三個正弦波的疊加

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赫茲淌实，根據(jù)采樣定理知采樣頻率要大于信號頻率2倍冻辩，所以這里設置采樣頻率為1400赫茲（即一秒內(nèi)有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赫茲。根據(jù)采樣定理淤年，fs至少是600赫茲的2倍钧敞，這里選擇1400赫茲，即在一秒內(nèi)選擇1400個點麸粮。

原始的函數(shù)圖像如下：

plt.figure()
plt.plot(x,y)   
plt.title('原始波形') 
plt.figure()
plt.plot(x[0:50],y[0:50])   
plt.title('原始部分波形（前50組樣本）')
plt.show()

由圖可見溉苛，由于采樣點太過密集，看起來不好看豹休，我們只顯示前面的50組數(shù)據(jù)炊昆，如下：

2、快速傅里葉變換

其實scipy和numpy一樣威根，實現(xiàn)FFT非常簡單凤巨，僅僅是一句話而已，函數(shù)接口如下：

from scipy.fftpack import fft,ifft

from numpy import fft,ifft

其中fft表示快速傅里葉變換洛搀，ifft表示其逆變換敢茁。具體實現(xiàn)如下：

fft_y=fft(y)                          #快速傅里葉變換
print(len(fft_y))
print(fft_y[0:5])
'''
運行結(jié)果如下：
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]
'''

我們發(fā)現(xiàn)以下幾個特點：

（1）變換之后的結(jié)果數(shù)據(jù)長度和原始采樣信號是一樣的

（2）每一個變換之后的值是一個復數(shù)，為a+bj的形式留美，那這個復數(shù)是什么意思呢彰檬？

我們知道伸刃，復數(shù)a+bj在坐標系中表示為（a,b），故而復數(shù)具有模和角度逢倍，我們都知道快速傅里葉變換具有

“振幅譜”“相位譜”捧颅，它其實就是通過對快速傅里葉變換得到的復數(shù)結(jié)果進一步求出來的，

那這個直接變換后的結(jié)果是不是就是我需要的较雕，當然是需要的碉哑，在FFT中，得到的結(jié)果是復數(shù)亮蒋，

（3）FFT得到的復數(shù)的模（即絕對值）就是對應的“振幅譜”扣典，復數(shù)所對應的角度，就是所對應的“相位譜”慎玖，現(xiàn)在可以畫圖了贮尖。

3、FFT的原始頻譜

N=1400x = np.arange(N)           # 頻率個數(shù) 
abs_y=np.abs(fft_y)                # 取復數(shù)的絕對值趁怔，即復數(shù)的模(雙邊頻譜)
angle_y=np.angle(fft_y)              #取復數(shù)的角度 
plt.figure()
plt.plot(x,abs_y)   
plt.title('雙邊振幅譜（未歸一化）') 
plt.figure()
plt.plot(x,angle_y)   
plt.title('雙邊相位譜（未歸一化）')
plt.show()

顯示結(jié)果如下：

雙邊振幅譜（未歸一化）

雙邊相位譜（未歸一化）

注意：我們在此處僅僅考慮“振幅譜”湿硝，不再考慮相位譜。

我們發(fā)現(xiàn)痕钢，振幅譜的縱坐標很大图柏，而且具有對稱性，這是怎么一回事呢任连？

關鍵：關于振幅值很大的解釋以及解決辦法——歸一化和取一半處理

比如有一個信號如下：

Y=A1+A2cos(2πω2+φ2）+A3cos(2πω3+φ3）+A4*cos(2πω4+φ4）

經(jīng)過FFT之后，得到的“振幅圖”中例诀，

第一個峰值（頻率位置）的模是A1的N倍随抠，N為采樣點，本例中為N=1400繁涂，此例中沒有拱她，因為信號沒有常數(shù)項A1

第二個峰值（頻率位置）的模是A2的N/2倍，N為采樣點扔罪，

第三個峰值（頻率位置）的模是A3的N/2倍秉沼，N為采樣點，

第四個峰值（頻率位置）的模是A4的N/2倍矿酵，N為采樣點唬复，

依次下去......

考慮到數(shù)量級較大，一般進行歸一化處理全肮，既然第一個峰值是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()

結(jié)果為：

雙邊頻譜(歸一化)

現(xiàn)在我們發(fā)現(xiàn)，振幅譜的數(shù)量級不大了测砂，變得合理了茵烈，

接下來進行取半處理：

half_x = x[range(int(N/2))]                                  #取一半?yún)^(qū)間
normalization_half_y = normalization_y[range(int(N/2))]      #由于對稱性，只取一半?yún)^(qū)間（單邊頻譜）
plt.figure()
plt.plot(half_x,normalization_half_y,'b')
plt.title('單邊頻譜(歸一化)',fontsize=9,color='blue')
plt.show()

單邊頻譜(歸一化)'

這就是我們最終的結(jié)果砌些，需要的“振幅譜”瞧毙。

PS:還需要注意fft的第0個數(shù)據(jù)，是表示直流分量寄症，本例中沒有宙彪，大部分情況不考慮第0個數(shù)據(jù)

三、完整代碼

import numpy as npfrom scipy.fftpack import fft,ifftimport matplotlib.pyplot as pltfrom matplotlib.pylab import mpl mpl.rcParams['font.sans-serif'] = ['SimHei']   #顯示中文mpl.rcParams['axes.unicode_minus']=False       #顯示負號 #采樣點選擇1400個有巧，因為設置的信號頻率分量最高為600赫茲释漆，根據(jù)采樣定理知采樣頻率要大于信號頻率2倍，所以這里設置采樣頻率為1400赫茲（即一秒內(nèi)有1400個采樣點篮迎，一樣意思的）x=np.linspace(0,1,1400)       #設置需要采樣的信號男图，頻率分量有200，400和600y=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=1400x = np.arange(N)             # 頻率個數(shù)half_x = x[range(int(N/2))]  #取一半?yún)^(qū)間 abs_y=np.abs(fft_y)                # 取復數(shù)的絕對值甜橱，即復數(shù)的模(雙邊頻譜)angle_y=np.angle(fft_y)            #取復數(shù)的角度normalization_y=abs_y/N            #歸一化處理（雙邊頻譜）                              normalization_half_y = normalization_y[range(int(N/2))]      #由于對稱性逊笆，只取一半?yún)^(qū)間（單邊頻譜） 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() 

匯總