傅里葉變換通俗解釋及快速傅里葉變換的python實現

　　通俗理解傅里葉變換，先看這篇文章傅里葉變換的通俗理解！

　　接下來即是使用python進行傅里葉FFT-頻譜分析：

1、一些關鍵概念的引入

一、離散傅里葉變換（DFT）

       離散傅里葉變換(discrete Fourier transform) 傅里葉分析方法是信號分析的最基本方法，傅里葉變換是傅里葉分析的核心，經過它把信號從時間域變換到頻率域，進而研究信號的頻譜結構和變化規律。可是它的致命缺點是：計算量太大，時間複雜度過高，當採樣點數過高的時候，計算緩慢，由此出現了DFT的快速實現，即下面的快速傅里葉變換FFT。

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

       計算量更小的離散傅里葉的一種實現方法。詳細細節這裏不作描述。

三、採樣頻率以及採樣定理

採樣頻率：採樣頻率，也稱爲採樣速度或者採樣率，定義了每秒從連續信號中提取並組成離散信號的採樣個數，它用赫茲（Hz）來表示。採樣頻率的倒數是採樣週期或者叫做採樣時間，它是採樣之間的時間間隔。通俗的講採樣頻率是指計算機每秒鐘採集多少個信號樣本。

採樣定理：所謂採樣定理 ，又稱香農採樣定理，奈奎斯特採樣定理，是信息論，特別是通信與信號處理學科中的一個重要基本結論。採樣定理指出，若是信號是帶限的，而且採樣頻率高於信號帶寬的兩倍，那麼，原來的連續信號能夠從採樣樣本中徹底重建出來。

定理的具體表述爲：在進行模擬/數字信號的轉換過程當中，當採樣頻率fs大於信號中最高頻率fmax的2倍時,即

　　fs>2*fmax

採樣以後的數字信號完整地保留了原始信號中的信息，通常實際應用中保證採樣頻率爲信號最高頻率的2.56～4倍；採樣定理又稱奈奎斯特定理。

四、如何理解採樣定理？

      在對連續信號進行離散化的過程當中，不免會損失不少信息，就拿一個簡單地正弦波而言，若是我1秒內就選擇一個點，很顯然，損失的信號太多了，光着一個點我根本不知道這個正弦信號究竟是什麼樣子的，天然也沒有辦法根據這一個採樣點進行正弦波的還原，很明顯，我採樣的點越密集，那越接近原來的正弦波原始的樣子，天然損失的信息越少，越方便還原正弦波。故而

       採樣定理說明採樣頻率與信號頻率之間的關係，是連續信號離散化的基本依據。 它爲採樣率創建了一個足夠的條件，該採樣率容許離散採樣序列從有限帶寬的連續時間信號中捕獲全部信息。

 

2、使用scipy包實現快速傅里葉變換

      本節不會說明FFT的底層實現，只介紹scipy中fft的函數接口以及使用的一些細節。

一、產生原始信號——原始信號是三個正弦波的疊加

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組數據，以下：

 

 

二、快速傅里葉變換

其實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獲得的複數的模（即絕對值）就是對應的「振幅譜」，複數所對應的角度，就是所對應的「相位譜」，如今能夠畫圖了。

 

三、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的一半，前半部分便可。

 

四、將振幅譜進行歸一化和取半處理

先進行歸一化

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()

這就是咱們最終的結果，須要的「振幅譜」。

 

3、完整代碼

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()

 

疑問解答，關於歸一化和取一半處理需看快速傅里葉變換在信號處理中的應用，具體爲：

　　傅里葉變換FT（Fourier Transform）是一種將信號從時域變換到頻域的變換形式。它在聲學、信號處理等領域有普遍的應用。計算機處理信號的要求是：在時域和頻域都應該是離散的，並且都應該是有限長的。而傅里葉變換僅能處理連續信號，離散傅里葉變換DFT（Discrete Fourier Transform）就是應這種須要而誕生的。它是傅里葉變換在離散域的表示形式。可是通常來講，DFT的運算量是很是大的。在1965年首次提出快速傅里葉變換算法FFT（Fast Fourier Transform）以前，其應用領域一直難以拓展，是FFT的提出使DFT的實現變得接近實時。DFT的應用領域也得以迅速拓展。除了一些速度要求很是高的場合以外，FFT算法基本上能夠知足工業應用的要求。因爲數字信號處理的其它運算均可以由DFT來實現，所以FFT算法是數字信號處理的重要基石。

　　傅立葉原理代表：任何連續測量的時序或信號，均可以表示爲不一樣頻率的正弦波信號的無限疊加。而根據該原理創立的傅立葉變換算法利用直接測量到的原始信號，以累加方式來計算該信號中不一樣正弦波信號的頻率、振幅和相位。如圖1所示，即爲時域信號與不一樣頻率的正弦波信號的關係。圖中最右側展現的是時域中的一個信號，這是一個近似於矩形的波，而圖的正中間則是組成該信號的各個頻率的正弦波。從圖中咱們能夠看出，即便角度幾乎爲直角的正弦波，其實也是由衆多的弧度圓滑的正弦波來組成的。在時域圖像中，咱們看到的只有一個矩形波，咱們無從得知他是由這些正弦波組成。但當咱們經過傅里葉變換將該矩形波轉換到頻域以後，咱們可以很清楚的看到許多脈衝，其中頻域圖中的橫軸爲頻率，縱軸爲振幅。所以能夠經過這個頻域圖像得知，時域中的矩形波是由這麼多頻率的正弦波疊加而成的。

　　這就是傅里葉變換的最基本最簡單的應用，固然這是從數學的角度去看傅立葉變換。在信號分析過程當中，傅里葉變換的做用就是將組成這個回波信號的全部輸入源在頻域中按照頻率的大小來表示出來。傅里葉變換以後，信號的幅度譜可表示對應頻率的能量，而相位譜可表示對應頻率的相位特徵。通過傅立葉變換能夠在頻率中很容易的找出雜亂信號中各頻率份量的幅度譜和相位譜，而後根據需求，進行高通或者低通濾波處理，最終獲得所須要頻率域的回波。

　　傅里葉變換在圖像處理過程當中也有很是重要的做用，設信號f是一個能量有限的模擬信號，則其傅里葉變換就表示信號f的頻譜。從純粹的數學意義上看，傅里葉變換是將一個函數轉換爲一系列周期函數來處理的。從物理效果看，傅里葉變換是將圖像從空間域轉換到頻率域，其逆變換是將圖像從頻率域轉換到空間域。換句話說，傅里葉變換的物理意義是將圖像的灰度分佈函數變換爲圖像的頻率分佈函數。傅里葉逆變換是將圖像的頻率分佈函數變換爲灰度分佈函數。傅里葉頻譜圖上咱們看到的明暗不一的亮點，其意義是指圖像上某一點與鄰域點差別的強弱，即梯度的大小，也即該點的頻率的大小。通常來說，梯度大則該點的亮度強，不然該點亮度弱。這樣經過觀察傅里葉變換後的頻譜圖，也叫功率圖，咱們就能夠直觀地看出圖像的能量分佈：若是頻譜圖中暗的點數更多，那麼實際圖像是比較柔和的，這是由於各點與鄰域差別都不大，梯度相對較小；反之，若是頻譜圖中亮的點數多，那麼實際圖像必定是尖銳的、邊界分明且邊界兩邊像素差別較大的。

　　以信號處理過程當中的一個例子來詳細說明FFT的效果：假設採樣頻率爲Fs，信號頻率爲F，採樣點數爲N。那麼FFT處理以後的結果就是一個點數爲N點的複數。每個點就對應着一個頻率點，而每一個點的模值，就是該頻率值下的幅度特性。假設原始信號的峯值爲A，那麼在處理後除第一個點以外的其餘點的模值就是A的N/2倍。而第一個點就是直流份量，它的模值就是直流份量的N倍。而每一個點的相位呢，就是在該頻率下的信號的相位。第一個點表示直流份量(即頻率爲0Hz)，而最後一個點N的再下一個點(實際上這個點是不存在的，這裏是假設的第N+1個點，也能夠看作是將第一個點分作兩半分，另外一半移到最後)則表示採樣頻率Fs，這中間被N-1個點平均分紅N等份，每一個點的頻率依次增長。例如某點n所表示的頻率爲：Fn=(n-1)*Fs/N。由上面的公式能夠看出，Fn 所能分辨到的最小頻率爲Fs/N，若是採樣頻率Fs爲1024Hz，採樣點數N爲1024點，則最小分辨率能夠精確到1Hz。1024Hz的採樣率採樣1024點，恰好是1秒，也就是說，採樣1秒時間的信號並作FFT處理，則結果能夠分析到1Hz；若是採樣2秒時間的信號並作FFT處理，則結果能夠精確到0.5Hz。

　　假設如今咱們有一個輸入信號，該信號總共包含3種成分信號，其一是5V的直流份量；其二是頻率爲50Hz、相位爲-60度、幅度爲10V的交流信號；第三個成分信號是頻率爲100Hz、相位爲90度、幅度爲5V的交流信號。該輸入信號用數學表達式表示以下：

 　　　　S=5+10*cos(2*pi*50*t-pi*60/180)+5*cos(2*pi*100*t+pi*90/180)

　　圖2即爲S信號的圖像表示。如今，咱們以256Hz的採樣率Fs對這個信號進行採樣，採樣點數N一樣爲256點。根據公式咱們能夠算出其頻譜圖中的頻率精度爲1Hz。所以對於輸入信號頻率包含0Hz、50Hz和100hz的複合信號，在其通過FFT處理以後，應該會在頻譜圖中出現3個峯值，並且頻率分別爲0Hz、50Hz和100Hz，處理結果如圖3所示：

　　結果正如咱們所預料的，對輸入信號’S’作FFT處理以後，圖3中出現了5個峯值，這是由於對輸入信號作256點的FFT處理以後並無第257個頻點信息，這也是前文中所提到的第一個點的模值是N倍的緣由。所以，信號的 FFT結果具備必定的對稱性。通常狀況下，咱們只使用前半部分的結果，即小於採樣頻率一半的結果。對於圖像進行簡單處理後，咱們的前半部分的FFT結果如圖4所示：

　　從圖4中能夠看出，三個輸入信號頻點的幅值依次爲1280、1280、640；其餘頻率所對應的幅值均爲0。按照公式，能夠計算直流份量（頻率爲0Hz）的幅值爲：1280/N= 1280/256=5;頻率爲50Hz的交流信號的幅值爲：1280/(N/2)= 1280/(256/2)=10;而75Hz的交流信號的幅值爲640/(N/2)=640/(256/2)=5。這也正是咱們輸入信號中的三個份量的直流份量值，因而可知，從頻譜分析出來的幅值是正確的。

　　經過上面的例子能夠看出，對於一個輸入信號，假如咱們不能肯定該輸入信號的頻率組成，咱們對其進行FFT處理以後，便能夠很輕鬆的看出其頻率份量，而且能夠經過簡單的計算來獲知該信號的幅值信息等。另外，若是想要提升頻率分辨率，咱們根據計算公式首先想到的就是須要增長採樣點數，但增長採樣點數也就意味着計算量增長，這在工程應用中增長了工程難度。解決這個問題的方法有頻率細分法，比較簡單的方法是採樣較短期的信號，而後在後面補充必定數量的0，使其長度達到須要的點數（通常爲2的冪次方的點數），而後再作FFT，就能在必定程度上提升頻率分辨率。