FFT演示程序 · 用 Python 做科學計算

# FFT演示程序本章詳細介紹如何綜合利用之前所學習的numpy，traits，traitsUI和Chaco等多個庫，編寫一個FFT演示程序。此程序可以幫助你理解FFT是如何將時域信號轉換為頻域信號的，在開始正式的程序編寫之前，讓我們來復習一下有關FFT變換的知識，如果你沒有學習過FFT，現在就是一個不錯的學習機會。 ## FFT知識復習 FFT變換是針對一組數值進行運算的，這組數的長度N必須是2的整數次冪，例如64, 128, 256等等；數值可以是實數也可以是復數，通常我們的時域信號都是實數，因此下面都以實數為例。我們可以把這一組實數想像成對某個連續信號按照一定取樣周期進行取樣而得來，如果對這組N個實數值進行FFT變換，將得到一個有N個復數的數組，我們稱此復數數組為頻域信號，此復數數組符合如下規律： * 下標為0和N/2的兩個復數的虛數部分為0， * 下標為i和N-i的兩個復數共軛，也就是其虛數部分數值相同、符號相反。下面的例子演示了這一個規律，先以rand隨機產生有8個元素的實數數組x，然后用fft對其運算之后，觀察其結果為8個復數，并且滿足上面兩個條件： ``` >>> x = np.random.rand(8) >>> x array([ 0.15562099, 0.56862756, 0.54371949, 0.06354358, 0.60678158, 0.78360968, 0.90116887, 0.1588846 ]) >>> xf = np.fft.fft(x) >>> xf array([ 3.78195634+0.j , -0.53575962+0.57688097j, -0.68248579-1.12980906j, -0.36656155-0.13801778j, 0.63262552+0.j , -0.36656155+0.13801778j, -0.68248579+1.12980906j, -0.53575962-0.57688097j]) ``` FFT變換的結果可以通過IFFT變換（逆FFT變換）還原為原來的值： ``` >>> np.fft.ifft(xf) array([ 0.15562099 +0.00000000e+00j, 0.56862756 +1.91940002e-16j, 0.54371949 +1.24900090e-16j, 0.06354358 -2.33573365e-16j, 0.60678158 +0.00000000e+00j, 0.78360968 +2.75206729e-16j, 0.90116887 -1.24900090e-16j, 0.15888460 -2.33573365e-16j]) ``` 注意ifft的運算結果實際上是和x相同的，由于浮點數的運算誤差，出現了一些非常小的虛數部分。 FFT變換和IFFT變換并沒有增加或者減少信號的數量，如果你仔細數一下的話，x中有8個實數數值，而xf中其實也只有8個有效的值。計算FFT結果中的有用的數值由于虛數部共軛和虛數部為0等規律，真正有用的信息保存在下標從0到N/2的N/2+1個虛數中，又由于下標為0和N/2的值虛數部分為0，因此只有N個有效的實數值。下面讓我們來看看FFT變換之后的那些復數都代表什么意思。 * 首先下標為0的實數表示了時域信號中的直流成分的多少 * 下標為i的復數a+b*j表示時域信號中周期為N/i個取樣值的正弦波和余弦波的成分的多少，其中a表示cos波形的成分，b表示sin波形的成分讓我們通過幾個例子來說明一下，下面是對一個直流信號進行FFT變換： ``` >>> x = np.ones(8) >>> x array([ 1., 1., 1., 1., 1., 1., 1., 1.]) >>> np.fft.fft(x)/len(x) array([ 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j]) ``` 所謂直流信號，就是其值不隨時間變化，因此我們創建一個值全為1的數組x，我們看到它的FFT結果除了下標為0的數值不為0以外，其余的都為0。(為了計算各個成分的能量多少，需要將FFT的結果除以FFT的長度)，這表示我們的時域信號是直流的，并且其成分為1。下面我們產生一個周期為8個取樣的正弦波，然后觀察其FFT的結果： ``` >>> x = np.arange(0, 2*np.pi, 2*np.pi/8) >>> y = np.sin(x) >>> np.fft.fft(y)/len(y) array([ 1.42979161e-18 +0.00000000e+00j, -4.44089210e-16 -5.00000000e-01j, 1.53075794e-17 -1.38777878e-17j, 3.87737802e-17 -1.11022302e-16j, 2.91853672e-17 +0.00000000e+00j, 0.00000000e+00 -9.71445147e-17j, 1.53075794e-17 +1.38777878e-17j, 3.44085112e-16 +5.00000000e-01j]) ``` 如何用linspace創建取樣點x 要計算周期為8個取樣的正弦波，就需要把0-2*pi的區間等分為8分，如果用np.linspace(0, 2*np.pi, 8)的話，產生的值為： ``` >>> np.linspace(0, 2*np.pi, 8) array([ 0\. , 0.8975979 , 1.7951958 , 2.6927937 , 3.5903916 , 4.48798951, 5.38558741, 6.28318531]) >>> 2*np.pi / 0.8975979 7.0000000079986666 ``` 可以看出上面用linspace只等分為7份，如果要正確使用np.linspace的話，可以如下調用，產生9個點，并且設置endpoint=False，最終結果不包括最后的那個點： ``` >>> np.linspace(0, 2*np.pi, 9, endpoint=False) array([ 0\. , 0.6981317 , 1.3962634 , 2.0943951 , 2.7925268 , 3.4906585 , 4.1887902 , 4.88692191, 5.58505361]) ``` 讓我們再來看看對正弦波的FFT的計算結果吧。可以看到下標為1的復數的虛數部分為-0.5，而我們產生的正弦波的放大系數(振幅)為1，它們之間的關系是-0.5*（-2）=1。再來看一下余弦波形: ``` >>> np.fft.fft(np.cos(x))/len(x) array([ -4.30631550e-17 +0.00000000e+00j, 5.00000000e-01 -2.52659764e-16j, 1.53075794e-17 +0.00000000e+00j, 1.11022302e-16 +1.97148613e-16j, 1.24479962e-17 +0.00000000e+00j, -1.11022302e-16 +1.91429446e-16j, 1.53075794e-17 +0.00000000e+00j, 5.00000000e-01 -1.35918295e-16j]) ``` 只有下標為1的復數的實數部分有有效數值0.5，和余弦波的放大系數1之間的關系是0.5*2=1。再來看2個例子： ``` >>> np.fft.fft(2*np.sin(2*x))/len(x) array([ 6.12303177e-17 +0.00000000e+00j, 6.12303177e-17 +6.12303177e-17j, -1.83690953e-16 -1.00000000e+00j, 6.12303177e-17 -6.12303177e-17j, 6.12303177e-17 +0.00000000e+00j, 6.12303177e-17 +6.12303177e-17j, -1.83690953e-16 +1.00000000e+00j, 6.12303177e-17 -6.12303177e-17j]) >>> np.fft.fft(0.8*np.cos(2*x))/len(x) array([ -2.44921271e-17 +0.00000000e+00j, -3.46370983e-17 +2.46519033e-32j, 4.00000000e-01 -9.79685083e-17j, 3.46370983e-17 -3.08148791e-32j, 2.44921271e-17 +0.00000000e+00j, 3.46370983e-17 -2.46519033e-32j, 4.00000000e-01 +9.79685083e-17j, -3.46370983e-17 +3.08148791e-32j]) ``` 上面產生的是周期為4個取樣(N/2)的正弦和余弦信號，其FFT的有效成分在下標為2的復數中，其中正弦波的放大系數為2,因此頻域虛數部分的值為-1；余弦波的放大系數為0.8，因此其對應的值為0.4。同頻率的正弦波和余弦波通過不同的系數疊加，可以產生同頻率的各種相位的余弦波，因此我們可以這樣來理解頻域中的復數: * 復數的模（絕對值）代表了此頻率的余弦波的振幅 * 復數的輻角代表了此頻率的余弦波的相位讓我們來看最后一個例子： ``` >>> x = np.arange(0, 2*np.pi, 2*np.pi/128) >>> y = 0.3*np.cos(x) + 0.5*np.cos(2*x+np.pi/4) + 0.8*np.cos(3*x-np.pi/3) >>> yf = np.fft.fft(y)/len(y) >>> yf[:4] array([ 1.00830802e-17 +0.00000000e+00j, 1.50000000e-01 +6.27820821e-18j, 1.76776695e-01 +1.76776695e-01j, 2.00000000e-01 -3.46410162e-01j]) >>> np.angle(yf[1]) 4.1854721366992471e-017 >>> np.abs(yf[1]), np.rad2deg(np.angle(yf[1])) (0.15000000000000008, 2.3980988870246962e-015) >>> np.abs(yf[2]), np.rad2deg(np.angle(yf[2])) (0.25000000000000011, 44.999999999999993) >>> np.abs(yf[3]), np.rad2deg(np.angle(yf[3])) (0.39999999999999991, -60.000000000000085) ``` 在這個例子中我們產生了三個不同頻率的余弦波，并且給他們不同的振幅和相位： * 周期為128/1.0點的余弦波的相位為0，振幅為0.3 * 周期為64/2.0點的余弦波的相位為45度，振幅為0.5 * 周期為128/3.0點的余弦波的相位為-60度，振幅為0.8 對照yf[1], yf[2], yf[3]的復數振幅和輻角，我想你應該對FFT結果中的每個數值所表示的意思有很清楚的理解了吧。 ## 合成時域信號前面說過通過ifft函數可以將頻域信號轉換回時域信號，這種轉換是精確的。下面我們要寫一個小程序，完成類似的事情，不過可以由用戶選擇只轉換一部分頻率回到時域信號，這樣轉換的結果和原始的時域信號會有誤差，我們通過觀察可以發現使用的頻率信息越多，則此誤差越小，直觀地看到如何通過多個余弦波逐步逼近任意的曲線信號的: ``` # -*- coding: utf-8 -*- # 本程序演示如何用多個正弦波合成三角波 import numpy as np import pylab as pl # 取FFT計算的結果freqs中的前n項進行合成，返回合成結果，計算loops個周期的波形 def fft_combine(freqs, n, loops=1): length = len(freqs) * loops data = np.zeros(length) index = loops * np.arange(0, length, 1.0) / length * (2 * np.pi) for k, p in enumerate(freqs[:n]): if k != 0: p *= 2 # 除去直流成分之外，其余的系數都*2 data += np.real(p) * np.cos(k*index) # 余弦成分的系數為實數部 data -= np.imag(p) * np.sin(k*index) # 正弦成分的系數為負的虛數部 return index, data # 產生size點取樣的三角波，其周期為1 def triangle_wave(size): x = np.arange(0, 1, 1.0/size) y = np.where(x<0.5, x, 0) y = np.where(x>=0.5, 1-x, y) return x, y fft_size = 256 # 計算三角波和其FFT x, y = triangle_wave(fft_size) fy = np.fft.fft(y) / fft_size # 繪制三角波的FFT的前20項的振幅，由于不含下標為偶數的值均為0，因此取 # log之后無窮小，無法繪圖，用np.clip函數設置數組值的上下限，保證繪圖正確 pl.figure() pl.plot(np.clip(20*np.log10(np.abs(fy[:20])), -120, 120), "o") pl.xlabel("frequency bin") pl.ylabel("power(dB)") pl.title("FFT result of triangle wave") # 繪制原始的三角波和用正弦波逐級合成的結果，使用取樣點為x軸坐標 pl.figure() pl.plot(y, label="original triangle", linewidth=2) for i in [0,1,3,5,7,9]: index, data = fft_combine(fy, i+1, 2) # 計算兩個周期的合成波形 pl.plot(data, label = "N=%s" % i) pl.legend() pl.title("partial Fourier series of triangle wave") pl.show() ``` ![](https://box.kancloud.cn/2016-03-19_56ed1bb4616aa.png) 三角波的頻譜 ![](https://box.kancloud.cn/2016-03-19_56ed1bb471362.png) 部分頻譜重建的三角波第18行的triangle_wave函數產生一個周期的三角波形，注意我們使用np.where函數計算區間函數的值。triangle函數返回兩個數組，分別表示x軸和y軸的值。注意后面的計算和繪圖不使用x軸坐標，而是直接用取樣次數作為x軸坐標。第7行的fft_combine的函數使用fft的結果freqs中的前n個數據重新合成時域信號，由于合成所使用的信號都是正弦波這樣的周期信號，所以我們可以通過第三個參數loops指定計算幾個周期。通過這個例子，我們可以看出使用的頻率越多，最終合成的波形越接近原始的三角波。合成方波由于方波的波形中存在跳變，因此用有限個正弦波合成的方波在跳變出現抖動現象，如[下圖](#fig-fftstudy03)所示，用正弦波合成的方波的收斂速度比三角波慢得多：計算方波的波形可以采用如下的函數： ``` def square_wave(size): x = np.arange(0, 1, 1.0/size) y = np.where(x<0.5, 1.0, 0) return x, y ``` ![](https://box.kancloud.cn/2016-03-19_56ed1bb48254e.png) 正弦波合成方波在跳變處出現都抖動 ## 三角波FFT演示程序我們希望制作一個用戶友好的界面，交互式地觀察各種三角波的頻譜以及其正弦合成的近似波形。制作界面是一件很費工的事情，幸好我們有TraitsUI庫的幫忙，不到200行程序就可以制作出如下的效果了: <object classid="clsid:D27CDB6E-AE6D-11cf-96B8-444553540000" width="589" height="447" codebase="http://active.macromedia.com/flash5/cabs/swflash.cab#version=7,0,0,0"><param name="movie" value="img/fft_study_04.swf"> <param name="play" value="true"> <param name="loop" value="false"> <param name="wmode" value="transparent"> <param name="quality" value="high"> <embed src="img/fft_study_04.swf" width="589" height="447" quality="high" loop="false" wmode="transparent" type="application/x-shockwave-flash" pluginspage="http://www.macromedia.com/shockwave/download/index.cgi?P1_Prod_Version=ShockwaveFlash"> </object> ![](https://box.kancloud.cn/2016-03-19_56ed1bb49f6ec.swf) 程序中已經給出了詳細的注釋，相信大家能夠讀懂并掌握這類程序的寫法，其中需要注意的幾點: * 16行，用ScrubberEditor創建一個自定義樣式的拖動調整值的控件，77-80行設置Item的editor = scrubber，這樣就用我們自定義的控件修改trait屬性了，如果不指定editor的話，Range類型的trait屬性將以一個滾動條做為編輯器。 * 用Range traits可以指定一個帶范圍的屬性，它可以設置上限下限，上下限可以用整數或者浮點數直接指定，也可以用另外一個trait屬性指定(用字符串指定trait屬性名)，但是幾種類型不能混用，因此程序中專門設計了兩個常數的trait屬性(36, 37行)，它們的作用只是用來指定其它trait屬性的上下限。 ``` low = Float(0.02) hi = Float(1.0) ``` * 190行的triangle_func函數返回一個用frompyfunc創建的ufunc函數，150行用此函數計算三角波，但是用frompyfunc創建的ufunc函數返回的數組的元素的類型(dtype)為object，因此需要用cast["float64"]函數強制將其轉換為類型為float64的數組。 * 處理trait屬性的改變事件最簡單的方式就是用固定的函數名: _trait屬性名_changed，但是當多個trait屬性需要共用某一個處理函數時，用@on_trait_change更加簡潔。下面是完整的程序: ``` # -*- coding: utf-8 -*- from enthought.traits.api import \ Str, Float, HasTraits, Property, cached_property, Range, Instance, on_trait_change, Enum from enthought.chaco.api import Plot, AbstractPlotData, ArrayPlotData, VPlotContainer from enthought.traits.ui.api import \ Item, View, VGroup, HSplit, ScrubberEditor, VSplit from enthought.enable.api import Component, ComponentEditor from enthought.chaco.tools.api import PanTool, ZoomTool import numpy as np # 鼠標拖動修改值的控件的樣式 scrubber = ScrubberEditor( hover_color = 0xFFFFFF, active_color = 0xA0CD9E, border_color = 0x808080 ) # 取FFT計算的結果freqs中的前n項進行合成，返回合成結果，計算loops個周期的波形 def fft_combine(freqs, n, loops=1): length = len(freqs) * loops data = np.zeros(length) index = loops * np.arange(0, length, 1.0) / length * (2 * np.pi) for k, p in enumerate(freqs[:n]): if k != 0: p *= 2 # 除去直流成分之外，其余的系數都*2 data += np.real(p) * np.cos(k*index) # 余弦成分的系數為實數部 data -= np.imag(p) * np.sin(k*index) # 正弦成分的系數為負的虛數部 return index, data class TriangleWave(HasTraits): # 指定三角波的最窄和最寬范圍，由于Range似乎不能將常數和traits名混用 # 所以定義這兩個不變的trait屬性 low = Float(0.02) hi = Float(1.0) # 三角波形的寬度 wave_width = Range("low", "hi", 0.5) # 三角波的頂點C的x軸坐標 length_c = Range("low", "wave_width", 0.5) # 三角波的定點的y軸坐標 height_c = Float(1.0) # FFT計算所使用的取樣點數，這里用一個Enum類型的屬性以供用戶從列表中選擇 fftsize = Enum( [(2**x) for x in range(6, 12)]) # FFT頻譜圖的x軸上限值 fft_graph_up_limit = Range(0, 400, 20) # 用于顯示FFT的結果 peak_list = Str # 采用多少個頻率合成三角波 N = Range(1, 40, 4) # 保存繪圖數據的對象 plot_data = Instance(AbstractPlotData) # 繪制波形圖的容器 plot_wave = Instance(Component) # 繪制FFT頻譜圖的容器 plot_fft = Instance(Component) # 包括兩個繪圖的容器 container = Instance(Component) # 設置用戶界面的視圖，注意一定要指定窗口的大小，這樣繪圖容器才能正常初始化 view = View( HSplit( VSplit( VGroup( Item("wave_width", editor = scrubber, label=u"波形寬度"), Item("length_c", editor = scrubber, label=u"最高點x坐標"), Item("height_c", editor = scrubber, label=u"最高點y坐標"), Item("fft_graph_up_limit", editor = scrubber, label=u"頻譜圖范圍"), Item("fftsize", label=u"FFT點數"), Item("N", label=u"合成波頻率數") ), Item("peak_list", style="custom", show_label=False, width=100, height=250) ), VGroup( Item("container", editor=ComponentEditor(size=(600,300)), show_label = False), orientation = "vertical" ) ), resizable = True, width = 800, height = 600, title = u"三角波FFT演示" ) # 創建繪圖的輔助函數，創建波形圖和頻譜圖有很多類似的地方，因此單獨用一個函數以 # 減少重復代碼 def _create_plot(self, data, name, type="line"): p = Plot(self.plot_data) p.plot(data, name=name, title=name, type=type) p.tools.append(PanTool(p)) zoom = ZoomTool(component=p, tool_mode="box", always_on=False) p.overlays.append(zoom) p.title = name return p def __init__(self): # 首先需要調用父類的初始化函數 super(TriangleWave, self).__init__() # 創建繪圖數據集，暫時沒有數據因此都賦值為空，只是創建幾個名字，以供Plot引用 self.plot_data = ArrayPlotData(x=[], y=[], f=[], p=[], x2=[], y2=[]) # 創建一個垂直排列的繪圖容器，它將頻譜圖和波形圖上下排列 self.container = VPlotContainer() # 創建波形圖，波形圖繪制兩條曲線：原始波形(x,y)和合成波形(x2,y2) self.plot_wave = self._create_plot(("x","y"), "Triangle Wave") self.plot_wave.plot(("x2","y2"), color="red") # 創建頻譜圖，使用數據集中的f和p self.plot_fft = self._create_plot(("f","p"), "FFT", type="scatter") # 將兩個繪圖容器添加到垂直容器中 self.container.add( self.plot_wave ) self.container.add( self.plot_fft ) # 設置 self.plot_wave.x_axis.title = "Samples" self.plot_fft.x_axis.title = "Frequency pins" self.plot_fft.y_axis.title = "(dB)" # 改變fftsize為1024，因為Enum的默認缺省值為枚舉列表中的第一個值 self.fftsize = 1024 # FFT頻譜圖的x軸上限值的改變事件處理函數，將最新的值賦值給頻譜圖的響應屬性 def _fft_graph_up_limit_changed(self): self.plot_fft.x_axis.mapper.range.high = self.fft_graph_up_limit def _N_changed(self): self.plot_sin_combine() # 多個trait屬性的改變事件處理函數相同時，可以用@on_trait_change指定 @on_trait_change("wave_width, length_c, height_c, fftsize") def update_plot(self): # 計算三角波 global y_data x_data = np.arange(0, 1.0, 1.0/self.fftsize) func = self.triangle_func() # 將func函數的返回值強制轉換成float64 y_data = np.cast["float64"](func(x_data)) # 計算頻譜 fft_parameters = np.fft.fft(y_data) / len(y_data) # 計算各個頻率的振幅 fft_data = np.clip(20*np.log10(np.abs(fft_parameters))[:self.fftsize/2+1], -120, 120) # 將計算的結果寫進數據集 self.plot_data.set_data("x", np.arange(0, self.fftsize)) # x坐標為取樣點 self.plot_data.set_data("y", y_data) self.plot_data.set_data("f", np.arange(0, len(fft_data))) # x坐標為頻率編號 self.plot_data.set_data("p", fft_data) # 合成波的x坐標為取樣點，顯示2個周期 self.plot_data.set_data("x2", np.arange(0, 2*self.fftsize)) # 更新頻譜圖x軸上限 self._fft_graph_up_limit_changed() # 將振幅大于-80dB的頻率輸出 peak_index = (fft_data > -80) peak_value = fft_data[peak_index][:20] result = [] for f, v in zip(np.flatnonzero(peak_index), peak_value): result.append("%s : %s" %(f, v) ) self.peak_list = "\n".join(result) # 保存現在的fft計算結果，并計算正弦合成波 self.fft_parameters = fft_parameters self.plot_sin_combine() # 計算正弦合成波，計算2個周期 def plot_sin_combine(self): index, data = fft_combine(self.fft_parameters, self.N, 2) self.plot_data.set_data("y2", data) # 返回一個ufunc計算指定參數的三角波 def triangle_func(self): c = self.wave_width c0 = self.length_c hc = self.height_c def trifunc(x): x = x - int(x) # 三角波的周期為1，因此只取x坐標的小數部分進行計算 if x >= c: r = 0.0 elif x < c0: r = x / c0 * hc else: r = (c-x) / (c-c0) * hc return r # 用trifunc函數創建一個ufunc函數，可以直接對數組進行計算, 不過通過此函數 # 計算得到的是一個Object數組，需要進行類型轉換 return np.frompyfunc(trifunc, 1, 1) if __name__ == "__main__": triangle = TriangleWave() triangle.configure_traits() ```