頻域信號處理

用FFT(快速傅立葉變換)能將時域的數字信號轉換為頻域信號。轉換為頻域信號之後我們可以很方便地分析出信號的頻率成分，在頻域上進行處理，最終還可以將處理完畢的頻域信號通過IFFT(逆變換)轉換為時域信號，實現許多在時域無法完成的信號處理演算法。本章通過幾個實例，簡單地介紹有關頻域信號處理的一些基本知識。

將時域信號通過FFT轉換為頻域信號之後，將其各個頻率分量的幅值繪製成圖，可以很直觀地觀察信號的頻譜。下面的程序完成這一任務：

# -*- coding: utf-8 -*-import numpy as npimport pylab as plsampling_rate = 8000fft_size = 512t = np.arange(0, 1.0, 1.0/sampling_rate)x = np.sin(2*np.pi*156.25*t) + 2*np.sin(2*np.pi*234.375*t)xs = x[:fft_size]xf = np.fft.rfft(xs)/fft_sizefreqs = np.linspace(0, sampling_rate/2, fft_size/2+1)xfp = 20*np.log10(np.clip(np.abs(xf), 1e-20, 1e100))pl.figure(figsize=(8,4))pl.subplot(211)pl.plot(t[:fft_size], xs)pl.xlabel(u"時間(秒)")pl.title(u"156.25Hz和234.375Hz的波形和頻譜")pl.subplot(212)pl.plot(freqs, xfp)pl.xlabel(u"頻率(Hz)")pl.subplots_adjust(hspace=0.4)pl.show()

下面逐行對這個程序進行解釋：

首先定義了兩個常數：sampling_rate, fft_size，分別表示數字信號的取樣頻率和FFT的長度。

然後調用np.arange產生1秒鐘的取樣時間，t中的每個數值直接表示取樣點的時間，因此其間隔為取樣周期1/sampline_rate：

t = np.arange(0, 1.0, 1.0/sampling_rate)

用取樣時間數組t可以很方便地調用函數計算出波形數據，這裡計算的是兩個正弦波的疊加，一個頻率是156.25Hz，一個是234.375Hz：

x = np.sin(2*np.pi*156.25*t) + 2*np.sin(2*np.pi*234.375*t)

為什麼選擇這兩個奇怪的頻率呢？因為這兩個頻率的正弦波在512個取樣點中正好有整數個周期。滿足這個條件波形的FFT結果能夠精確地反映其頻譜。

N點FFT能精確計算的頻率

假設取樣頻率為fs, 取波形中的N個數據進行FFT變換。那麼這N點數據包含整數個周期的波形時，FFT所計算的結果是精確的。於是能精確計算的波形的周期是: n*fs/N。對於8kHz取樣，512點FFT來說，8000/512.0 = 15.625Hz，前面的156.25Hz和234.375Hz正好是其10倍和15倍。

下面從波形數據x中截取fft_size個點進行fft計算。np.fft庫中提供了一個rfft函數，它方便我們對實數信號進行FFT計算。根據FFT計算公式，為了正確顯示波形能量，還需要將rfft函數的結果除以fft_size：

xs = x[:fft_size]xf = np.fft.rfft(xs)/fft_size

rfft函數的返回值是N/2+1個複數，分別表示從0(Hz)到sampling_rate/2(Hz)的N/2+1點頻率的成分。於是可以通過下面的np.linspace計算出返回值中每個下標對應的真正的頻率：

freqs = np.linspace(0, sampling_rate/2, fft_size/2+1)

最後我們計算每個頻率分量的幅值，並通過 20*np.log10() 將其轉換為以db單位的值。為了防止0幅值的成分造成log10無法計算，我們調用np.clip對xf的幅值進行上下限處理：

xfp = 20*np.log10(np.clip(np.abs(xf), 1e-20, 1e100))

剩下的程序就是將時域波形和頻域波形繪製出來，這裡就不再詳細敘述了。此程序的輸出為：

圖18.1 使用FFT計算正弦波的頻譜

如果你放大其頻譜中的兩個峰值的部分的話，可以看到其值分別為：

>>> xfp[10]-6.0205999132796251>>> xfp[15]-9.6432746655328714e-16

即156.25Hz的成分為-6dB， 而234.375Hz的成分為0dB，與波形的計算公式中的各個分量的能量(振幅值/2)符合。

如果我們波形不能在fft_size個取樣中形成整數個周期的話會怎樣呢？

將波形計算公式修改為：

x = np.sin(2*np.pi*200*t) + 2*np.sin(2*np.pi*300*t)

得到的結果如下：

圖18.2 非完整周期的正弦波經過FFT變換之後出現頻譜泄漏

這次得到的頻譜不再是兩個完美的峰值，而是兩個峰值頻率周圍的頻率都有能量。這顯然和兩個正弦波的疊加波形的頻譜有區別。本來應該屬於200Hz和300Hz的能量分散到了周圍的頻率中，這個現象被稱為頻譜泄漏。出現頻譜泄漏的原因在於fft_size個取樣點無法放下整數個200Hz和300Hz的波形。

頻譜泄漏的解釋

我們只能在有限的時間段中對信號進行測量，無法知道在測量範圍之外的信號是怎樣的。因此只能對測量範圍之外的信號進行假設。而傅立葉變換的假設很簡單：測量範圍之外的信號是所測量到的信號的重複。

現在考慮512點FFT，從信號中取出的512個數據就是FFT的測量範圍，它計算的是這512個數據一直重複的波形的頻譜。顯然如果512個數據包含整數個周期的話，那麼得到的結果就是原始信號的頻譜，而如果不是整數周期的話，得到的頻譜就是如下波形的頻譜，這裡假設對50Hz的正弦波進行512點FFT：

>>> t = np.arange(0, 1.0, 1.0/8000)>>> x = np.sin(2*np.pi*50*t)[:512]>>> pl.plot(np.hstack([x,x,x]))>>> pl.show()

圖18.3 50Hz正弦波的512點FFT所計算的頻譜的實際波形

由於這個波形的前後不是連續的，出現波形跳變，而跳變處的有著非常廣泛的頻譜，因此FFT的結果中出現頻譜泄漏。

為了減少FFT所截取的數據段前後的跳變，可以對數據先乘以一個窗函數，使得其前後數據能平滑過渡。例如常用的hann窗函數的定義如下：

其中N為窗函數的點數，下面是一個512點hann窗的曲線：

>>> import pylab as pl>>> import scipy.signal as signal>>> pl.figure(figsize=(8,3))>>> pl.plot(signal.hann(512))

圖18.4 hann窗函數

窗函數都在scipy.signal庫中定義，它們的第一個參數為點數N。可以看出hann窗函數是完全對稱的，也就是說第0點和第511點的值完全相同，都為0。在這樣的函數和信號數據相乘的話，結果中會出現前後兩個連續的0，這樣FFT的結果所表示的周期信號中有兩個連續的0值，會對信號的周期性有一定的影響。

計算周期信號的一個周期的數據

考慮對一個正弦波取樣10個點，那麼第一個點的值為0，而最後一個點的值不應該是0，這樣這10個數據的重複才能是精確的正弦波，下面的兩種計算中，前者是正確的：

>>> np.sin(np.arange(0, 2*np.pi, 2*np.pi/10))array([ 0.00000000e+00, 5.87785252e-01, 9.51056516e-01, 9.51056516e-01, 5.87785252e-01, 1.22464680e-16, -5.87785252e-01, -9.51056516e-01, -9.51056516e-01, -5.87785252e-01])>>> np.sin(np.linspace(0, 2*np.pi, 10))array([ 0.00000000e+00, 6.42787610e-01, 9.84807753e-01, 8.66025404e-01, 3.42020143e-01, -3.42020143e-01, -8.66025404e-01, -9.84807753e-01, -6.42787610e-01, -2.44929360e-16])

為了解決連續0值的問題，hann函數提供了一個sym關鍵字參數，如果設置其為0的話，那麼將產生一個N+1點的hann窗函數，然後取其前N個數，這樣得到的窗函數適合於周期信號：

>>> signal.hann(8)array([ 0. , 0.1882551 , 0.61126047, 0.95048443, 0.95048443, 0.61126047, 0.1882551 , 0. ])>>> signal.hann(8, sym=0)array([ 0. , 0.14644661, 0.5 , 0.85355339, 1. , 0.85355339, 0.5 , 0.14644661])

50Hz正弦波與窗函數乘積之後的重複波形如下：

>>> t = np.arange(0, 1.0, 1.0/8000)>>> x = np.sin(2*np.pi*50*t)[:512] * signal.hann(512, sym=0)>>> pl.plot(np.hstack([x,x,x]))>>> pl.show()

圖18.5 加hann窗的50Hz正弦波的512點FFT所計算的實際波形

回到前面的例子，將200Hz, 300Hz的疊加波形與hann窗乘積之後再計算其頻譜，得到如下頻譜圖：

圖18.6 加hann窗前後的頻譜，hann窗能降低頻譜泄漏

可以看到與hann窗乘積之後的信號的頻譜能量更加集中於200Hz和300Hz，但是其能量有所降低。這是因為hann窗本身有一定的能量衰減：

>>> np.sum(signal.hann(512, sym=0))/5120.5

因此如果需要嚴格保持信號的能量的話，還需要在乘以hann窗之後再乘以2。

上面完整繪圖程序請參照： 頻譜泄漏和hann窗

對於頻譜特性不隨時間變化的信號，例如引擎、壓縮機等機器雜訊，可以對其進行長時間的採樣，然後分段進行FFT計算，最後對每個頻率分量的幅值求其平均值可以準確地測量信號的頻譜。

下面的程序完成這一計算：

import numpy as npimport scipy.signal as signalimport pylab as pldef average_fft(x, fft_size):n = len(x) // fft_size * fft_sizetmp = x[:n].reshape(-1, fft_size)tmp *= signal.hann(fft_size, sym=0)xf = np.abs(np.fft.rfft(tmp)/fft_size)avgf = np.average(xf, axis=0)return 20*np.log10(avgf)

average_fft(x, fft_size)對數組x進行fft_size點FFT運算，以dB為單位返回其平均後的幅值。由於x的長度可能不是fft_size的整數倍，因此首先將其縮短為fft_size的整數倍，然後用reshape函數將其轉換為一個二維數組tmp。tmp的第1軸的長度為fft_size：

n = len(x) // fft_size * fft_sizetmp = x[:n].reshape(-1, fft_size)

然後將tmp的第1軸上的數據和窗函數相乘，這裡選用的是hann窗：

tmp *= signal.hann(fft_size, sym=0)

調用rfft對tmp每的行數據進行FFT計算，並求其幅值：

xf = np.abs(np.fft.rfft(tmp)/fft_size)

接下來調用average函數對xf沿著第0軸進行平均，這樣就得到每個頻率分量的平均幅值：

avgf = np.average(xf, axis=0)

下面是利用averagge_fft函數計算隨機數序列頻譜的例子：

>>> x = np.random.rand(100000) - 0.5>>> xf = average_fft(x, 512)>>> pl.plot(xf)>>> pl.show()

圖18.7 白色雜訊的頻譜接近水平直線(注意Y軸的範圍)

我們可以看到隨機雜訊的頻譜接近一條水平的直線，也就是說每個頻率窗口的能量都相同，這種雜訊我們稱之為白色雜訊。

如果我們利用scipy.signal庫中的濾波器設計函數，設計一個IIR低通濾波器，將白色雜訊輸入到此低通濾波器，繪製其輸出數據的平均頻譜的話，就能夠觀察到IIR濾波器的頻率響應特性，下面的程序利用iirdesign設計一個8kHz取樣的1kHz的 Chebyshev I 型低通濾波器，iirdesign函數需要用正規化的頻率(取值範圍為0-1)，然後調用filtfilt對白色雜訊信號x進行低通濾波：

>>> b,a=signal.iirdesign(1000/4000.0, 1100/4000.0, 1, -40, 0, "cheby1")>>> x = np.random.rand(100000) - 0.5>>> y = signal.filtfilt(b, a, x)

如果用average_fft計算輸出信號y的平均頻譜，得到如下頻譜圖：

圖18.8 經過低通濾波器的白色雜訊的頻譜

我們知道，信號x經過系統h之後的輸出y是x和h的卷積，雖然卷積的計算方法很簡單，但是當x和h都很長的時候，卷積計算是非常耗費時間的。因此對於比較長的系統h，需要找到比直接計算卷積更快的方法。

信號系統理論中有這樣一個規律：時域的卷積等於頻域的乘積，因此要計算時域的卷積，可以將時域信號轉換為頻域信號，進行乘積運算之後再將結果轉換為時域信號，實現快速卷積。

由於FFT運算可以高效地將時域信號轉換為頻域信號，其運算的複雜度為 O(N*log(N))，因此三次FFT運算加一次乘積運算的總複雜度仍然為 O(N*log(N)) 級別，而卷積運算的複雜度為 O(N*N)，顯然通過FFT計算卷積要比直接計算快速得多。這裡假設需要卷積的兩個信號的長度都為N。

但是有一個問題：FFT運算假設其所計算的信號為周期信號，因此通過上述方法計算出的結果實際上是兩個信號的循環卷積，而不是線性卷積。為了用FFT計算線性卷積，需要對信號進行補零擴展，使得其長度長於線性卷積結果的長度。

例如，如果我們要計算數組a和b的卷積，a和b的長度都為128，那麼它們的卷積結果的長度為 len(a) + len(b) - 1 = 257。為了用FFT能夠計算其線性卷積，需要將a和b都擴展到256。下面的程序演示這個計算過程：

# -*- coding: utf-8 -*-import numpy as npdef fft_convolve(a,b):n = len(a)+len(b)-1N = 2**(int(np.log2(n))+1)A = np.fft.fft(a, N)B = np.fft.fft(b, N)return np.fft.ifft(A*B)[:n]if __name__ == "__main__":a = np.random.rand(128)b = np.random.rand(128)c = np.convolve(a,b)print np.sum(np.abs(c - fft_convolve(a,b)))

此程序的輸出為直接卷積和FFT快速卷積的結果之間的誤差，大約為5e-12左右。

在這段程序中，a,b的長度為128，其卷積結果c的長度為n=255，我們通過下面的算式找到大於n的最小的2的整數次冪：

N = 2**(int(np.log2(n))+1)

在調用fft函數對其進行變換時，傳遞第二個參數為N(FFT的長度)，這樣fft函數將自動對a,b進行補零。最後通過ifft得到的卷積結果c2的長度為N，比實際的卷積結果c要多出一個數，這個多出來的元素應該接近於0，請讀者自行驗證。

下面測試一下速度：

>>> import timeit>>> setup="""import numpy as npa=np.random.rand(10000)b=np.random.rand(10000)from spectrum_fft_convolve import fft_convolve""">>> timeit.timeit("np.convolve(a,b)",setup, number=10)1.852900578146091>>> timeit.timeit("fft_convolve(a,b)",setup, number=10)0.19475575806416145

顯然計算兩個很長的數組的卷積，FFT快速卷積要比直接卷積快很多。但是對於較短的數組，直接卷積運算將更快一些。下圖顯示了直接卷積和快速卷積的每點的平均計算時間和長度之間的關係：

圖18.9 用FFT計算卷積和直接卷積的時間複雜度比較

由於圖中的Y軸表示每點的計算時間，因此對於直接卷積它是線性的：O(N*N)/N。我們看到對於1024點以上的計算，快速卷積顯示出明顯的優勢。

具體的程序請參照 FFT卷積的速度比較

由於FFT卷積很常用，因此scipy.signal庫中提供了fftconvolve函數，此函數採用FFT運算可以計算多維數組的卷積。讀者也可以參照此函數的源代碼幫助理解。

現在考慮對於輸入信號x和系統響應h的卷積運算，通常x是非常長的，例如要對某段錄音進行濾波處理，假設取樣頻率為8kHz，錄音長度為1分鐘的話，那麼x的長度為480000。而且x的長度也可能不是固定的，例如我們可能需要對麥克風的連續輸入信號進行濾波處理。而h的長度通常都是固定的，例如它是某個房間的衝擊響應，或者是某種FIR濾波器。

根據前面的介紹，為了有效地利用FFT計算卷積，我們希望它的兩個輸入長度相當，於是就需要對信號x進行分段處理。對卷積的分段運算被稱作：overlap-add運算。

overlap-add的計算方法如下圖所示：

圖18.10 分段卷積的過程演示

原始信號x長度為300，將它分為三段，分別與濾波器係數h進行卷積計算，h的長度為101，因此每段輸出200個數據，圖中用綠色標出每段輸出的200個數據。這3段數據按照時間順序進行求和之後得到結果和原始信號的卷積是相同的。

因此將持續的輸入信號x和濾波器h進行卷積的運算可以按照如下步驟進行，假設h的長度為M：

1. 建立一個緩存，其大小為N+M-1，初始值為0
2. 每次從x中讀取N個數據，和h進行卷積，得到N+M-1個數據，和緩存中的數據進行求和，並放進緩存中，然後輸出緩存前N個數據
3. 將緩存中的數據向左移動N個元素，也就是讓緩存中的第N個元素成為第0個元素，後面的N個元素全部設置為0
4. 跳轉到2重複運行

下面是實現這一演算法的演示程序：

# -*- coding: utf-8 -*-import numpy as npx = np.random.rand(1000)h = np.random.rand(101)y = np.convolve(x, h)N = 50 # 分段大小M = len(h) # 濾波器長度output = []#緩存初始化為0buffer = np.zeros(M+N-1,dtype=np.float64)for i in xrange(len(x)/N):#從輸入信號中讀取N個數據xslice = x[i*N:(i+1)*N]#計算卷積yslice = np.convolve(xslice, h)#將卷積的結果加入到緩衝中buffer += yslice#輸出緩存中的前N個數據，注意使用copy，否則輸出的是buffer的一個視圖output.append( buffer[:N].copy() )#緩存中的數據左移動N個元素buffer[0:M-1] = buffer[N:]#後面的補0buffer[M-1:] = 0#將輸出的數據組合為數組y2 = np.hstack(output)#計算和直接卷積的結果之間的誤差print np.sum(np.abs( y2 - y[:len(x)] ) )

注意第23行需要輸出緩存前N個數據的拷貝，否則輸出的是數組的一個視圖，當此後buffer更新時，視圖中的數據會一起更新。

將FFT快速卷積和overlap-add相結合，可以製作出一些快速的實時數據濾波演算法。但是由於FFT卷積對於兩個長度相當的數組時最為有效，因此在分段時也會有所限制：例如如果濾波器的長度為2048，那麼理想的分段長度也為2048，如果將分段長度設置得過低，反而會增加運算量。因此在實時性要求很強的系統中，只能採用直接卷積。

18.3　Hilbert變換?

Hilbert變換能在振幅保持不變的情況下將輸入信號的相角偏移90度，簡單地說就是能將正弦波形轉換為餘弦波形，下面的程序驗證這一特性：

# -*- coding: utf-8 -*-from scipy import fftpackimport numpy as npimport matplotlib.pyplot as pl# 產生1024點4個周期的正弦波t = np.linspace(0, 8*np.pi, 1024, endpoint=False)x = np.sin(t)# 進行Hilbert變換y = fftpack.hilbert(x)pl.plot(x, label=u"原始波形")pl.plot(y, label=u"Hilbert轉換後的波形")pl.legend()pl.show()

關於程序的幾點說明：

此程序的輸出圖表如下，我們可以很清楚地看出hilbert將正弦波變換為了餘弦波形。

圖18.11 Hilbert變換將正弦波變為餘弦波

Hilbert正變換的相角偏移符號

本書中將相角偏移+90度成為Hilbert正變換。有的文獻書籍正好將定義倒轉過來：將偏移+90度成為Hilbert負變換，而偏移-90度成為Hilbert正變換。

相角偏移90度相當於複數平面上的點與虛數單位1j相乘，因此Hilbert變換的頻率響應可以用如下公式給出：

其中

為圓頻率，sgn函數為符號函數，即：

我們可以將其頻率響應理解為：

由於對於實數信號來說，正負頻率成分共軛，因此對實數信號進行Hilbert變換之後仍然是實數信號。下面的程序驗證Hilbert變換的頻率響應：

>>> x = np.random.rand(16)>>> y = fftpack.hilbert(x)>>> X = np.fft.fft(x)>>> Y = np.fft.fft(y)>>> np.imag(Y/X)array([ 0., 1., 1., 1., 1., 1., 1., 1., 0., -1., -1., -1., -1., -1., -1., -1.])

對信號進行N點FFT變換之後：

對照Y/X的虛數部分，不難看出它是符合Hilbert的頻率響應的。如果你用np.real(Y/X)觀察實數部分的話，它們全部接近於0。

Hilbert變換可以用作包絡檢波。具體演算法如下所示：

其中x為原始載波波形，H(x)為x的Hilbert變換之後的波形，envelope為信號x的包絡。其原理很容易理解：假設x為正弦波，那麼H(x)為餘弦波，根據公式：

可知envelope恆等於1，為sin(t)信號的包絡。下面的程序驗證這一演算法：

# -*- coding: utf-8 -*-import numpy as npimport pylab as plfrom scipy import fftpackt = np.arange(0, 0.3, 1/20000.0)x = np.sin(2*np.pi*1000*t) * (np.sin(2*np.pi*10*t) + np.sin(2*np.pi*7*t) + 3.0)hx = fftpack.hilbert(x)pl.plot(x, label=u"載波信號")pl.plot(np.sqrt(x**2 + hx**2), "r", linewidth_=2, label=u"檢出的包絡信號")pl.title(u"使用Hilbert變換進行包絡檢波")pl.legend()pl.show()

圖18.12 使用Hilbert變換對載波信號進行包絡檢波

前面介紹過可以使用頻率掃描波形測量濾波器的頻率響應，我們可以使用這個演算法計算出掃描波的包絡：

>>> run filter_lfilter_example01.py # 運行濾波器的例子>>> hy = fftpack.hilbert(y)>>> pl.plot( np.sqrt(y**2 + hy**2),"r", linewidth_=2)>>> pl.plot(y)>>> pl.title(u"頻率掃描波的包絡")>>> pl.show()

得到的包絡波形如下圖所示：

圖18.13 使用Hilbert變換對頻率掃描波進行包絡檢波

可以看出在高頻和低頻處包絡計算出現較大的誤差。而在中頻部分能很好地計算出包絡的形狀。