代碼來源於http://bigsec.net/b52/scipydoc/frequency_process.html
文章從6個方面來寫，首先是觀察頻譜的特征，第二部分是加上窗函數之后的特征，第三部分是頻譜平均，第四部分是比較FFT與直接卷積時間效率區別，第五部分是由於FFT對輸入信號的長度有要求，因此介紹了overlap-add分段運算，最后一部分是Hilbert變換的實現。
- 觀察信號的頻譜
數據通過FFT轉換成頻域信號，對頻域信號進行分析，再通過IFFT轉換成時域信號。
import numpy as np import pylab as pl import matplotlib as mpl mpl.rcParams['font.sans-serif'] = ['KaiTi'] mpl.rcParams['font.serif'] = ['KaiTi'] mpl.rcParams['axes.unicode_minus']=False sampling_rate = 8000 #取樣頻率 fft_size = 512 #fft長度 t = np.arange(0, 1.0, 1.0/sampling_rate) #假設取樣頻率為fs, 取波形中的N個數據進行FFT變換。那么這N點數據包含整數個周期的波形時，FFT所計算的結果是精確的。於是能精確計算的波形的周期是: n*fs/N。 #對於8kHz取樣，512點FFT來說，8000/512.0 = 15.625Hz，前面的156.25Hz和234.375Hz正好是其10倍和15倍。 #選取整數倍的數據，查看當fft后的數據在頻譜中形成整數周期時的情況。 x = np.sin(2*np.pi*156.25*t) + 2*np.sin(2*np.pi*234.375*t) #選取非整數倍的數據，查看當fft后的數據在頻譜中沒有形成非整數周期時的情況。 x = np.sin(2*np.pi*200*t) + 2*np.sin(2*np.pi*300*t) xs = x[:fft_size] #取數據 xf = np.fft.rfft(xs)/fft_size #rfft：對實數信號進行FFT變換。/fft_size是為了正確顯示波形能量 freqs = np.linspace(0, sampling_rate/2, fft_size/2+1) #fft_size/2+1個點，后面的是與前面的共軛 #計算每個頻率分量的幅值，並通過 20*np.log10() 將其轉換為以db(分貝)為單位的值 xfp = 20*np.log10(np.clip(np.abs(xf), 1e-20, 1e100)) #clip將數據限制在最小值和最大值之間 pl.figure(figsize=(8,4)) #新建一個8*4英寸的圖紙 pl.subplot(211) #繪制2行1列的圖紙，這個圖形占據第一行 pl.plot(t[:fft_size], xs) pl.xlabel(u"時間(秒)") pl.title(u"156.25Hz和234.375Hz的波形和頻譜") pl.subplot(212) #繪制2行1列的圖紙，這個圖形占據第二行 pl.plot(freqs, xfp) pl.xlabel(u"頻率(Hz)") pl.subplots_adjust(hspace=0.4) pl.show()
運行結果如下圖所示：
整數周期情況：
非整數周期情況：
如上兩圖可看出非整數周期情況下，第二種情況可能會發生信號泄露狀況。由於FFT的假設前提是測量之外的信號是所測量信號的不斷重復，如代碼中那樣，取8k個樣例，從信號中取出的512個數據就是FFT的測量范圍，也就是說每個波形周期是15.625Hz，它計算的是這512個數據一直重復的波形的頻譜。顯然如果512個數據包含整數個周期的話，那么得到的結果就是原始信號的頻譜，而如果不是整數周期的話，得到的頻譜就是如下波形的頻譜，這里假設對50Hz的正弦波進行512點FFT，這種波形會發生跳變。
- 窗函數
為了減少FFT所截取的數據段前后的跳變，可以對數據先乘以一個窗函數，使得其前后數據能平滑過渡。例如常用的hann窗函數的定義如下：
hann窗的曲線如下所示：
可看到，最前和最后的值都是0，如果直接去乘的前后會出現兩個0，因此可考慮將這個波形用N+1個點表示，而取前N個點，這樣第N+1個點就是下一個波形的第一個點，也就是0,通過設置sym參數解決。這與調用linspace時指定endpoint=False類似，丟掉最后一個點。
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])
將hann窗與50hz相乘，它的曲線會更加平滑。
之前的非整數周期加了hann窗之后的結果如下圖所示：
- 頻譜平均
對於頻譜特性不隨時間變化的信號，例如引擎、壓縮機等機器噪聲，可以對其進行長時間的采樣，然后分段進行FFT計算，最后對每個頻率分量的幅值求其平均值可以准確地測量信號的頻譜，測試隨機數序列頻譜如下所示
def average_fft(x, fft_size): n = len(x) // fft_size * fft_size tmp = x[:n].reshape(-1, fft_size) #將n行轉成n列，轉為一個二維數組 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) x = np.random.rand(100000) - 0.5 xf = average_fft(x, 512) pl.plot(xf) pl.show()
結果為
這個頻譜的能量趨近於一條直線，每個窗的能量相差不大，被稱為白色噪聲。
- 快速卷積
信息處理可看作是將原始信號與一個信號進行卷積，也就需要考慮運算效率。這部分主要是比較FFT和直接卷積的運算效率
def fft_convolve(a,b): n = len(a)+len(b)-1 N = 2**(int(np.log2(n))+1) #FFT的長度：大於n的最小的2的整數次冪 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))))
結果為1.865645261656436e-12，也就是說FFT和普通卷積的結果相差很小，但速度卻快很多。
- 分段運算
對於輸入信號 x 和系統向量（eg:FIR濾波器）h而言，x的長度不固定，h的長度固定。為了加快卷積效率， 我們需要x和h的長度相當，也就是說對x進行分段處理，這種分段算法被稱為overlap-add運算。但是由於FFT在兩個數組的分段長度相當時最為有效，因此在實時性要求很強的系統中，采用直接卷積會更好一些。
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 range(int(len(x)/N)): #從輸入信號中讀取N個數據 xslice = x[i*N:(i+1)*N] #計算卷積 yslice = np.convolve(xslice, h) pl.cla() pl.plot(yslice) #將卷積的結果加入到緩沖中 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)] ) ))
- Hilbert 變換
Hilbert變換能在振幅保持不變的情況下將輸入信號的相角偏移90度，簡單地說就是能將正弦波形轉換為余弦波形。
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變換后可將直流分量變為0，正頻率成分偏移+90度，負頻率成分偏移-90度。它也可用來進行包絡檢波。
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()
