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

17.1　FFT知識復習

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

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

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

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變換的結(jié)果可以通過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的運算結(jié)果實際上是和x相同的，由于浮點數(shù)的運算誤差，出現(xiàn)了一些非常小的虛數(shù)部分。

FFT變換和IFFT變換并沒有增加或者減少信號的數(shù)量，如果你仔細數(shù)一下的話，x中有8個實數(shù)數(shù)值，而xf中其實也只有8個有效的值。

計算FFT結(jié)果中的有用的數(shù)值

由于虛數(shù)部共軛和虛數(shù)部為0等規(guī)律，真正有用的信息保存在下標從0到N/2的N/2+1個虛數(shù)中， 又由于下標為0和N/2的值虛數(shù)部分為0，因此只有N個有效的實數(shù)值。

下面讓我們來看看FFT變換之后的那些復數(shù)都代表什么意思。

• 首先下標為0的實數(shù)表示了時域信號中的直流成分的多少
• 下標為i的復數(shù)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])

所謂直流信號，就是其值不隨時間變化，因此我們創(chuàng)建一個值全為1的數(shù)組x，我們看到它的FFT結(jié)果除了下標為0的數(shù)值不為0以外，其余的都為0。(為了計算各個成分的能量多少，需要將FFT的結(jié)果除以FFT的長度)，這表示我們的時域信號是直流的，并且其成分為1。

下面我們產(chǎn)生一個周期為8個取樣的正弦波，然后觀察其FFT的結(jié)果：

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創(chuàng)建取樣點x

要計算周期為8個取樣的正弦波，就需要把0-2*pi的區(qū)間等分為8分，如果用np.linspace(0, 2*np.pi, 8)的話，產(chǎn)生的值為：

>>> 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的話， 可以如下調(diào)用，產(chǎn)生9個點，并且設置endpoint=False，最終結(jié)果不包括最后的那個點：

>>> 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的計算結(jié)果吧。可以看到下標為1的復數(shù)的虛數(shù)部分為-0.5，而我們產(chǎn)生的正弦波的放大系數(shù)(振幅)為1，它們之間的關(guān)系是-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的復數(shù)的實數(shù)部分有有效數(shù)值0.5，和余弦波的放大系數(shù)1之間的關(guān)系是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])

上面產(chǎn)生的是周期為4個取樣(N/2)的正弦和余弦信號，其FFT的有效成分在下標為2的復數(shù)中，其中正弦波的放大系數(shù)為2,因此頻域虛數(shù)部分的值為-1；余弦波的放大系數(shù)為0.8，因此其對應的值為0.4。

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

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

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

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)

在這個例子中我們產(chǎn)生了三個不同頻率的余弦波，并且給他們不同的振幅和相位：

• 周期為128/1.0點的余弦波的相位為0， 振幅為0.3
• 周期為64/2.0點的余弦波的相位為45度， 振幅為0.5
• 周期為128/3.0點的余弦波的相位為-60度，振幅為0.8

對照yf[1], yf[2], yf[3]的復數(shù)振幅和輻角，我想你應該對FFT結(jié)果中的每個數(shù)值所表示的意思有很清楚的理解了吧。

17.2　合成時域信號

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

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計算的結(jié)果freqs中的前n項進行合成，返回合成結(jié)果，計算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 # 除去直流成分之外，其余的系數(shù)都*2data += np.real(p) * np.cos(k*index) # 余弦成分的系數(shù)為實數(shù)部data -= np.imag(p) * np.sin(k*index) # 正弦成分的系數(shù)為負的虛數(shù)部return index, data # 產(chǎn)生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項的振幅，由于不含下標為偶數(shù)的值均為0， 因此取 # log之后無窮小，無法繪圖，用np.clip函數(shù)設置數(shù)組值的上下限，保證繪圖正確 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")# 繪制原始的三角波和用正弦波逐級合成的結(jié)果，使用取樣點為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函數(shù)產(chǎn)生一個周期的三角波形，注意我們使用np.where函數(shù)計算區(qū)間函數(shù)的值。triangle函數(shù)返回兩個數(shù)組，分別表示x軸和y軸的值。注意后面的計算和繪圖不使用x軸坐標，而是直接用取樣次數(shù)作為x軸坐標。

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

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

合成方波

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

計算方波的波形可以采用如下的函數(shù)：

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?正弦波合成方波在跳變處出現(xiàn)都抖動

17.3　三角波FFT演示程序

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

?

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

• 16行，用ScrubberEditor創(chuàng)建一個自定義樣式的拖動調(diào)整值的控件，77-80行設置Item的editor = scrubber，這樣就用我們自定義的控件修改trait屬性了，如果不指定editor的話，Range類型的trait屬性將以一個滾動條做為編輯器。

• 用Range traits可以指定一個帶范圍的屬性，它可以設置上限下限，上下限可以用整數(shù)或者浮點數(shù)直接指定，也可以用另外一個trait屬性指定(用字符串指定trait屬性名)，但是幾種類型不能混用，因此程序中專門設計了兩個常數(shù)的trait屬性(36, 37行)，它們的作用只是用來指定其它trait屬性的上下限。

  low = Float(0.02) hi = Float(1.0)

• 190行的triangle_func函數(shù)返回一個用frompyfunc創(chuàng)建的ufunc函數(shù)，150行用此函數(shù)計算三角波，但是用frompyfunc創(chuàng)建的ufunc函數(shù)返回的數(shù)組的元素的類型(dtype)為object，因此需要用cast["float64"]函數(shù)強制將其轉(zhuǎn)換為類型為float64的數(shù)組。

• 處理trait屬性的改變事件最簡單的方式就是用固定的函數(shù)名: _trait屬性名_changed，但是當多個trait屬性需要共用某一個處理函數(shù)時，用@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計算的結(jié)果freqs中的前n項進行合成，返回合成結(jié)果，計算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 # 除去直流成分之外，其余的系數(shù)都*2data += np.real(p) * np.cos(k*index) # 余弦成分的系數(shù)為實數(shù)部data -= np.imag(p) * np.sin(k*index) # 正弦成分的系數(shù)為負的虛數(shù)部return index, data class TriangleWave(HasTraits):# 指定三角波的最窄和最寬范圍，由于Range似乎不能將常數(shù)和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計算所使用的取樣點數(shù)，這里用一個Enum類型的屬性以供用戶從列表中選擇fftsize = Enum( [(2**x) for x in range(6, 12)])# FFT頻譜圖的x軸上限值fft_graph_up_limit = Range(0, 400, 20)# 用于顯示FFT的結(jié)果peak_list = Str# 采用多少個頻率合成三角波N = Range(1, 40, 4)# 保存繪圖數(shù)據(jù)的對象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點數(shù)"),Item("N", label=u"合成波頻率數(shù)")),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演示")# 創(chuàng)建繪圖的輔助函數(shù)，創(chuàng)建波形圖和頻譜圖有很多類似的地方，因此單獨用一個函數(shù)以# 減少重復代碼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):# 首先需要調(diào)用父類的初始化函數(shù)super(TriangleWave, self).__init__()# 創(chuàng)建繪圖數(shù)據(jù)集，暫時沒有數(shù)據(jù)因此都賦值為空，只是創(chuàng)建幾個名字，以供Plot引用self.plot_data = ArrayPlotData(x=[], y=[], f=[], p=[], x2=[], y2=[]) # 創(chuàng)建一個垂直排列的繪圖容器，它將頻譜圖和波形圖上下排列self.container = VPlotContainer()# 創(chuàng)建波形圖，波形圖繪制兩條曲線： 原始波形(x,y)和合成波形(x2,y2)self.plot_wave = self._create_plot(("x","y"), "Triangle Wave")self.plot_wave.plot(("x2","y2"), color="red")# 創(chuàng)建頻譜圖，使用數(shù)據(jù)集中的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軸上限值的改變事件處理函數(shù)，將最新的值賦值給頻譜圖的響應屬性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屬性的改變事件處理函數(shù)相同時，可以用@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函數(shù)的返回值強制轉(zhuǎn)換成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)# 將計算的結(jié)果寫進數(shù)據(jù)集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)# 保存現(xiàn)在的fft計算結(jié)果，并計算正弦合成波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計算指定參數(shù)的三角波def triangle_func(self):c = self.wave_widthc0 = self.length_chc = self.height_cdef trifunc(x):x = x - int(x) # 三角波的周期為1，因此只取x坐標的小數(shù)部分進行計算if x >= c: r = 0.0elif x < c0: r = x / c0 * hcelse: r = (c-x) / (c-c0) * hcreturn r# 用trifunc函數(shù)創(chuàng)建一個ufunc函數(shù)，可以直接對數(shù)組進行計算, 不過通過此函數(shù)# 計算得到的是一個Object數(shù)組，需要進行類型轉(zhuǎn)換return np.frompyfunc(trifunc, 1, 1) if __name__ == "__main__":triangle = TriangleWave()triangle.configure_traits()

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。