本章詳細介紹如何綜合利用之前所學習的numpy，traits，traitsUI和Chaco等多個庫，編寫一個FFT演示程序。此程序可以幫助你理解FFT是如何將時域信號轉換為頻域信號的，在開始正式的程序編寫之前，讓我們來復習一下有關FFT變換的知識，如果你沒有學習過FFT，現在就是一個不錯的學習機會。

17.1　FFT知識復習

FFT變換是針對一組數值進行運算的，這組數的長度N必須是2的整數次冪，例如64, 128, 256等等； 數值可以是實數也可以是復數，通常我們的時域信號都是實數，因此下面都以實數為例。我們可以把這一組實數想像成對某個連續信號按照一定取樣周期進行取樣而得來，如果對這組N個實數值進行FFT變換，將得到一個有N個復數的數組，我們稱此復數數組為頻域信號，此復數數組符合如下規律：

• 下標為0和N/2的兩個復數的虛數部分為0，
• 下標為i和N-i的兩個復數共軛，也就是其虛數部分數值相同、符號相反。

下面的例子演示了這一個規律，先以rand隨機產生有8個元素的實數數組x，然后用fft對其運算之后，觀察其結果為8個復數，并且滿足上面兩個條件：

123456789 10

>>> 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的結果：

123456789 10

>>> 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個例子：

123456789 10 11 12 13 14 15 16

>>> 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。

同頻率的正弦波和余弦波通過不同的系數疊加，可以產生同頻率的各種相位的余弦波，因此我們可以這樣來理解頻域中的復數:

• 復數的模（絕對值）代表了此頻率的余弦波的振幅
• 復數的輻角代表了此頻率的余弦波的相位

讓我們來看最后一個例子：

123456789 10 11 12 13 14 15

>>> 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結果中的每個數值所表示的意思有很清楚的理解了吧。

17.2　合成時域信號

前面說過通過ifft函數可以將頻域信號轉換回時域信號，這種轉換是精確的。下面我們要寫一個小程序，完成類似的事情，不過可以由用戶選擇只轉換一部分頻率回到時域信號，這樣轉換的結果和原始的時域信號會有誤差，我們通過觀察可以發現使用的頻率信息越多，則此誤差越小，直觀地看到如何通過多個余弦波逐步逼近任意的曲線信號的:

123456789 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46

# -*- 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) * loopsdata = 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 # 除去直流成分之外，其余的系數都*2data += 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, yfft_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()

圖17.1?三角波的頻譜

圖17.2?部分頻譜重建的三角波

第18行的triangle_wave函數產生一個周期的三角波形，注意我們使用np.where函數計算區間函數的值。triangle函數返回兩個數組，分別表示x軸和y軸的值。注意后面的計算和繪圖不使用x軸坐標，而是直接用取樣次數作為x軸坐標。

第7行的fft_combine的函數使用fft的結果freqs中的前n個數據重新合成時域信號，由于合成所使用的信號都是正弦波這樣的周期信號，所以我們可以通過第三個參數loops指定計算幾個周期。

通過這個例子，我們可以看出使用的頻率越多，最終合成的波形越接近原始的三角波。

合成方波

由于方波的波形中存在跳變，因此用有限個正弦波合成的方波在跳變出現抖動現象，如圖17.3所示，用正弦波合成的方波的收斂速度比三角波慢得多：

計算方波的波形可以采用如下的函數：

def square_wave(size):x = np.arange(0, 1, 1.0/size)y = np.where(x<0.5, 1.0, 0)return x, y

圖17.3?正弦波合成方波在跳變處出現都抖動

17.3　三角波FFT演示程序

我們希望制作一個用戶友好的界面，交互式地觀察各種三角波的頻譜以及其正弦合成的近似波形。 制作界面是一件很費工的事情，幸好我們有TraitsUI庫的幫忙，不到200行程序就可以制作出如下的效果了:

?

程序中已經給出了詳細的注釋，相信大家能夠讀懂并掌握這類程序的寫法，其中需要注意的幾點:

• 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更加簡潔。

下面是完整的程序:

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208

# -*- coding: utf-8 -*- from enthought.traits.api import \Str, Float, HasTraits, Property, cached_property, Range, Instance, on_trait_change, Enumfrom enthought.chaco.api import Plot, AbstractPlotData, ArrayPlotData, VPlotContainerfrom enthought.traits.ui.api import \Item, View, VGroup, HSplit, ScrubberEditor, VSplitfrom 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) * loopsdata = 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 # 除去直流成分之外，其余的系數都*2data += 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 = namereturn pdef __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和pself.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_limitdef _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_datax_data = np.arange(0, 1.0, 1.0/self.fftsize)func = self.triangle_func()# 將func函數的返回值強制轉換成float64y_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_parametersself.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_widthc0 = self.length_chc = self.height_cdef trifunc(x):x = x - int(x) # 三角波的周期為1，因此只取x坐標的小數部分進行計算if x >= c: r = 0.0elif x < c0: r = x / c0 * hcelse: r = (c-x) / (c-c0) * hcreturn r# 用trifunc函數創建一個ufunc函數，可以直接對數組進行計算, 不過通過此函數# 計算得到的是一個Object數組，需要進行類型轉換return np.frompyfunc(trifunc, 1, 1) if __name__ == "__main__":triangle = TriangleWave()triangle.configure_traits()

總結

以上是生活随笔為你收集整理的基于python的FFT演示程序的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。