傅裏葉變換通俗解釋及快速傅裏葉變換的python實現
通俗理解傅裏葉變換,先看這篇文章傅裏葉變換的通俗理解!
接下來便是使用python進行傅裏葉FFT-頻譜分析:
一、一些關鍵概念的引入
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()
疑問解答,關於歸一化和取一半處理需看快速傅裏葉變換在信號處理中的應用,具體為:
傅裏葉變換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,就能在一定程度上提高頻率分辨率。
傅裏葉變換通俗解釋及快速傅裏葉變換的python實現