# 頻域信號處理 用FFT(快速傅立葉變換)能將時域的數字信號轉換為頻域信號。轉換為頻域信號之后我們可以很方便地分析出信號的頻率成分，在頻域上進行處理，最終還可以將處理完畢的頻域信號通過IFFT(逆變換)轉換為時域信號，實現許多在時域無法完成的信號處理算法。本章通過幾個實例，簡單地介紹有關頻域信號處理的一些基本知識。 ## 觀察信號的頻譜 將時域信號通過FFT轉換為頻域信號之后，將其各個頻率分量的幅值繪制成圖，可以很直觀地觀察信號的頻譜。下面的程序完成這一任務： ``` # -*- coding: utf-8 -*- import numpy as np import pylab as pl sampling_rate = 8000 fft_size = 512 t = 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_size freqs = 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)) ``` 剩下的程序就是將時域波形和頻域波形繪制出來，這里就不再詳細敘述了。此程序的輸出為：  使用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) ``` 得到的結果如下：  非完整周期的正弦波經過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() ```  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)) ```  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() ```  加hann窗的50Hz正弦波的512點FFT所計算的實際波形 回到前面的例子，將200Hz, 300Hz的疊加波形與hann窗乘積之后再計算其頻譜，得到如下頻譜圖：  加hann窗前后的頻譜，hann窗能降低頻譜泄漏 可以看到與hann窗乘積之后的信號的頻譜能量更加集中于200Hz和300Hz，但是其能量有所降低。這是因為hann窗本身有一定的能量衰減： ``` >>> np.sum(signal.hann(512, sym=0))/512 0.5 ``` 因此如果需要嚴格保持信號的能量的話，還需要在乘以hann窗之后再乘以2。 上面完整繪圖程序請參照： [_頻譜泄漏和hann窗_](example_spectrum_example_hann.html) ### 頻譜平均 對于頻譜特性不隨時間變化的信號，例如引擎、壓縮機等機器噪聲，可以對其進行長時間的采樣，然后分段進行FFT計算，最后對每個頻率分量的幅值求其平均值可以準確地測量信號的頻譜。 下面的程序完成這一計算： ``` import numpy as np import scipy.signal as signal import pylab as pl def average_fft(x, fft_size): n = len(x) // fft_size * fft_size tmp = 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_size tmp = 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() ```  白色噪聲的頻譜接近水平直線(注意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的平均頻譜，得到如下頻譜圖：  經過低通濾波器的白色噪聲的頻譜 ## 快速卷積 我們知道，信號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 np def fft_convolve(a,b): n = len(a)+len(b)-1 N = 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 np a=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快速卷積要比直接卷積快很多。但是對于較短的數組，直接卷積運算將更快一些。下圖顯示了直接卷積和快速卷積的每點的平均計算時間和長度之間的關系：  用FFT計算卷積和直接卷積的時間復雜度比較 由于圖中的Y軸表示每點的計算時間，因此對于直接卷積它是線性的：O(N*N)/N。我們看到對于1024點以上的計算，快速卷積顯示出明顯的優勢。 具體的程序請參照 [_FFT卷積的速度比較_](example_spectrum_fft_convolve_timeit.html) 由于FFT卷積很常用，因此scipy.signal庫中提供了fftconvolve函數，此函數采用FFT運算可以計算多維數組的卷積。讀者也可以參照此函數的源代碼幫助理解。 ### 分段運算 現在考慮對于輸入信號x和系統響應h的卷積運算，通常x是非常長的，例如要對某段錄音進行濾波處理，假設取樣頻率為8kHz，錄音長度為1分鐘的話，那么x的長度為480000。而且x的長度也可能不是固定的，例如我們可能需要對麥克風的連續輸入信號進行濾波處理。而h的長度通常都是固定的，例如它是某個房間的沖擊響應，或者是某種FIR濾波器。 根據前面的介紹，為了有效地利用FFT計算卷積，我們希望它的兩個輸入長度相當，于是就需要對信號x進行分段處理。對卷積的分段運算被稱作：overlap-add運算。 overlap-add的計算方法如下圖所示：  分段卷積的過程演示 原始信號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 np x = np.random.rand(1000) h = np.random.rand(101) y = np.convolve(x, h) N = 50 # 分段大小 M = len(h) # 濾波器長度 output = [] #緩存初始化為0 buffer = 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:] #后面的補0 buffer[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，如果將分段長度設置得過低，反而會增加運算量。因此在實時性要求很強的系統中，只能采用直接卷積。 ## Hilbert變換 Hilbert變換能在振幅保持不變的情況下將輸入信號的相角偏移90度，簡單地說就是能將正弦波形轉換為余弦波形，下面的程序驗證這一特性： ``` # -*- coding: utf-8 -*- from scipy import fftpack import numpy as np import 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轉換函數在scipy.fftpack函數庫中 * 為了生成完美的正弦波，需要計算整數個周期，因此調用linspace時指定endpoint=False，這樣就不會包括區間的結束點：8*np.pi 此程序的輸出圖表如下，我們可以很清楚地看出hilbert將正弦波變換為了余弦波形。  Hilbert變換將正弦波變為余弦波 Hilbert正變換的相角偏移符號 本書中將相角偏移+90度成為Hilbert正變換。有的文獻書籍正好將定義倒轉過來：將偏移+90度成為Hilbert負變換，而偏移-90度成為Hilbert正變換。 相角偏移90度相當于復數平面上的點與虛數單位1j相乘，因此Hilbert變換的頻率響應可以用如下公式給出：  其中  為圓頻率，sgn函數為符號函數，即：  我們可以將其頻率響應理解為： * 直流分量為0 * 正頻率成分偏移+90度 * 負頻率成分偏移-90度 由于對于實數信號來說，正負頻率成分共軛，因此對實數信號進行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變換之后： * 下標為0的頻率分量表示直流分量 * 下標為N/2的的頻率分量為取樣頻率/2的頻率分量 * 1到N/2-1為正頻率分量 * N/2+1到N為負頻率分量 對照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 np import pylab as pl from scipy import fftpack t = 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() ```  使用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() ``` 得到的包絡波形如下圖所示：  使用Hilbert變換對頻率掃描波進行包絡檢波 可以看出在高頻和低頻處包絡計算出現較大的誤差。而在中頻部分能很好地計算出包絡的形狀。