頻域信號處理

賢人好客 2010-09-08

展開全文

18　頻域信號處理¶

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

18.1　觀察信號的頻譜¶

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

# -*- 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()

下面逐行對這個程序進(jìn)行解釋：

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

然后調(diào)用np.arange產(chǎn)生1秒鐘的取樣時間，t中的每個數(shù)值直接表示取樣點的時間，因此其間隔為取樣周期1/sampline_rate：

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

用取樣時間數(shù)組t可以很方便地調(diào)用函數(shù)計算出波形數(shù)據(jù)，這里計算的是兩個正弦波的疊加，一個頻率是156.25Hz，一個是234.375Hz：

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

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

N點FFT能精確計算的頻率

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

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

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

rfft函數(shù)的返回值是N/2+1個復(fù)數(shù)，分別表示從0(Hz)到sampling_rate/2(Hz)的N/2+1點頻率的成分。于是可以通過下面的np.linspace計算出返回值中每個下標(biāo)對應(yīng)的真正的頻率：

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

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

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

剩下的程序就是將時域波形和頻域波形繪制出來，這里就不再詳細(xì)敘述了。此程序的輸出為：

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

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

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

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

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

將波形計算公式修改為：

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

得到的結(jié)果如下：

圖18.2 非完整周期的正弦波經(jīng)過FFT變換之后出現(xiàn)頻譜泄漏

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

頻譜泄漏的解釋

我們只能在有限的時間段中對信號進(jìn)行測量，無法知道在測量范圍之外的信號是怎樣的。因此只能對測量范圍之外的信號進(jìn)行假設(shè)。而傅立葉變換的假設(shè)很簡單：測量范圍之外的信號是所測量到的信號的重復(fù)。

現(xiàn)在考慮512點FFT，從信號中取出的512個數(shù)據(jù)就是FFT的測量范圍，它計算的是這512個數(shù)據(jù)一直重復(fù)的波形的頻譜。顯然如果512個數(shù)據(jù)包含整數(shù)個周期的話，那么得到的結(jié)果就是原始信號的頻譜，而如果不是整數(shù)周期的話，得到的頻譜就是如下波形的頻譜，這里假設(shè)對50Hz的正弦波進(jìn)行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所計算的頻譜的實際波形

由于這個波形的前后不是連續(xù)的，出現(xiàn)波形跳變，而跳變處的有著非常廣泛的頻譜，因此FFT的結(jié)果中出現(xiàn)頻譜泄漏。

18.1.1　窗函數(shù)¶

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

$w(n)= 0.5\; \left(1 - \cos \left ( \frac{2 \pi n}{N-1} \right) \right)$

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

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

圖18.4 hann窗函數(shù)

窗函數(shù)都在scipy.signal庫中定義，它們的第一個參數(shù)為點數(shù)N?？梢钥闯鰄ann窗函數(shù)是完全對稱的，也就是說第0點和第511點的值完全相同，都為0。在這樣的函數(shù)和信號數(shù)據(jù)相乘的話，結(jié)果中會出現(xiàn)前后兩個連續(xù)的0，這樣FFT的結(jié)果所表示的周期信號中有兩個連續(xù)的0值，會對信號的周期性有一定的影響。

計算周期信號的一個周期的數(shù)據(jù)

考慮對一個正弦波取樣10個點，那么第一個點的值為0，而最后一個點的值不應(yīng)該是0，這樣這10個數(shù)據(jù)的重復(fù)才能是精確的正弦波，下面的兩種計算中，前者是正確的：

>>> 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])

為了解決連續(xù)0值的問題，hann函數(shù)提供了一個sym關(guān)鍵字參數(shù)，如果設(shè)置其為0的話，那么將產(chǎn)生一個N+1點的hann窗函數(shù)，然后取其前N個數(shù)，這樣得到的窗函數(shù)適合于周期信號：

>>> 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正弦波與窗函數(shù)乘積之后的重復(fù)波形如下：

>>> 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))/512
0.5

因此如果需要嚴(yán)格保持信號的能量的話，還需要在乘以hann窗之后再乘以2。

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

18.1.2　頻譜平均¶

對于頻譜特性不隨時間變化的信號，例如引擎、壓縮機等機器噪聲，可以對其進(jìn)行長時間的采樣，然后分段進(jìn)行FFT計算，最后對每個頻率分量的幅值求其平均值可以準(zhǔn)確地測量信號的頻譜。

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

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)對數(shù)組x進(jìn)行fft_size點FFT運算，以dB為單位返回其平均后的幅值。由于x的長度可能不是fft_size的整數(shù)倍，因此首先將其縮短為fft_size的整數(shù)倍，然后用reshape函數(shù)將其轉(zhuǎn)換為一個二維數(shù)組tmp。tmp的第1軸的長度為fft_size：

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

然后將tmp的第1軸上的數(shù)據(jù)和窗函數(shù)相乘，這里選用的是hann窗：

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

調(diào)用rfft對tmp每的行數(shù)據(jù)進(jìn)行FFT計算，并求其幅值：

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

接下來調(diào)用average函數(shù)對xf沿著第0軸進(jìn)行平均，這樣就得到每個頻率分量的平均幅值：

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

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

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

圖18.7 白色噪聲的頻譜接近水平直線(注意Y軸的范圍)

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

如果我們利用scipy.signal庫中的濾波器設(shè)計函數(shù)，設(shè)計一個IIR低通濾波器，將白色噪聲輸入到此低通濾波器，繪制其輸出數(shù)據(jù)的平均頻譜的話，就能夠觀察到IIR濾波器的頻率響應(yīng)特性，下面的程序利用iirdesign設(shè)計一個8kHz取樣的1kHz的 Chebyshev I 型低通濾波器，iirdesign函數(shù)需要用正規(guī)化的頻率(取值范圍為0-1)，然后調(diào)用filtfilt對白色噪聲信號x進(jìn)行低通濾波：

>>> 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 經(jīng)過低通濾波器的白色噪聲的頻譜

18.2　快速卷積¶

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

信號系統(tǒng)理論中有這樣一個規(guī)律：時域的卷積等于頻域的乘積，因此要計算時域的卷積，可以將時域信號轉(zhuǎn)換為頻域信號，進(jìn)行乘積運算之后再將結(jié)果轉(zhuǎn)換為時域信號，實現(xiàn)快速卷積。

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

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

例如，如果我們要計算數(shù)組a和b的卷積，a和b的長度都為128，那么它們的卷積結(jié)果的長度為 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快速卷積的結(jié)果之間的誤差，大約為5e-12左右。

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

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

在調(diào)用fft函數(shù)對其進(jìn)行變換時，傳遞第二個參數(shù)為N(FFT的長度)，這樣fft函數(shù)將自動對a,b進(jìn)行補零。最后通過ifft得到的卷積結(jié)果c2的長度為N，比實際的卷積結(jié)果c要多出一個數(shù)，這個多出來的元素應(yīng)該接近于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

顯然計算兩個很長的數(shù)組的卷積，F(xiàn)FT快速卷積要比直接卷積快很多。但是對于較短的數(shù)組，直接卷積運算將更快一些。下圖顯示了直接卷積和快速卷積的每點的平均計算時間和長度之間的關(guān)系：

圖18.9 用FFT計算卷積和直接卷積的時間復(fù)雜度比較

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

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

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

18.2.1　分段運算¶

現(xiàn)在考慮對于輸入信號x和系統(tǒng)響應(yīng)h的卷積運算，通常x是非常長的，例如要對某段錄音進(jìn)行濾波處理，假設(shè)取樣頻率為8kHz，錄音長度為1分鐘的話，那么x的長度為480000。而且x的長度也可能不是固定的，例如我們可能需要對麥克風(fēng)的連續(xù)輸入信號進(jìn)行濾波處理。而h的長度通常都是固定的，例如它是某個房間的沖擊響應(yīng)，或者是某種FIR濾波器。

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

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

圖18.10 分段卷積的過程演示

原始信號x長度為300，將它分為三段，分別與濾波器系數(shù)h進(jìn)行卷積計算，h的長度為101，因此每段輸出200個數(shù)據(jù)，圖中用綠色標(biāo)出每段輸出的200個數(shù)據(jù)。這3段數(shù)據(jù)按照時間順序進(jìn)行求和之后得到結(jié)果和原始信號的卷積是相同的。

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

建立一個緩存，其大小為N+M-1，初始值為0
每次從x中讀取N個數(shù)據(jù)，和h進(jìn)行卷積，得到N+M-1個數(shù)據(jù)，和緩存中的數(shù)據(jù)進(jìn)行求和，并放進(jìn)緩存中，然后輸出緩存前N個數(shù)據(jù)
將緩存中的數(shù)據(jù)向左移動N個元素，也就是讓緩存中的第N個元素成為第0個元素，后面的N個元素全部設(shè)置為0
跳轉(zhuǎn)到2重復(fù)運行

下面是實現(xiàn)這一算法的演示程序：

# -*- 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個數(shù)據(jù)
xslice = x[i*N:(i+1)*N]
#計算卷積
yslice = np.convolve(xslice, h)
#將卷積的結(jié)果加入到緩沖中
buffer += yslice
#輸出緩存中的前N個數(shù)據(jù)，注意使用copy，否則輸出的是buffer的一個視圖
output.append( buffer[:N].copy() )
#緩存中的數(shù)據(jù)左移動N個元素
buffer[0:M-1] = buffer[N:]
#后面的補0
buffer[M-1:] = 0
#將輸出的數(shù)據(jù)組合為數(shù)組
y2 = np.hstack(output)
#計算和直接卷積的結(jié)果之間的誤差
print np.sum(np.abs( y2 - y[:len(x)] ) )

注意第23行需要輸出緩存前N個數(shù)據(jù)的拷貝，否則輸出的是數(shù)組的一個視圖，當(dāng)此后buffer更新時，視圖中的數(shù)據(jù)會一起更新。

將FFT快速卷積和overlap-add相結(jié)合，可以制作出一些快速的實時數(shù)據(jù)濾波算法。但是由于FFT卷積對于兩個長度相當(dāng)?shù)臄?shù)組時最為有效，因此在分段時也會有所限制：例如如果濾波器的長度為2048，那么理想的分段長度也為2048，如果將分段長度設(shè)置得過低，反而會增加運算量。因此在實時性要求很強的系統(tǒng)中，只能采用直接卷積。

18.3　Hilbert變換¶

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

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

關(guān)于程序的幾點說明：

hilbert轉(zhuǎn)換函數(shù)在scipy.fftpack函數(shù)庫中
為了生成完美的正弦波，需要計算整數(shù)個周期，因此調(diào)用linspace時指定endpoint=False，這樣就不會包括區(qū)間的結(jié)束點：8*np.pi

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

圖18.11 Hilbert變換將正弦波變?yōu)橛嘞也?/p>

Hilbert正變換的相角偏移符號

本書中將相角偏移+90度成為Hilbert正變換。有的文獻(xiàn)書籍正好將定義倒轉(zhuǎn)過來：將偏移+90度成為Hilbert負(fù)變換，而偏移-90度成為Hilbert正變換。

相角偏移90度相當(dāng)于復(fù)數(shù)平面上的點與虛數(shù)單位1j相乘，因此Hilbert變換的頻率響應(yīng)可以用如下公式給出：

$H(\omega) = j \cdot sgn(\omega)$

其中 $\omega$ 為圓頻率，sgn函數(shù)為符號函數(shù)，即：

$sgn(\omega) =\begin{cases} \ \ 1, & \mbox{for } \omega > 0,\\ \ \ 0, & \mbox{for } \omega = 0,\\ -1, & \mbox{for } \omega <0, \end{cases}$

我們可以將其頻率響應(yīng)理解為：

直流分量為0
正頻率成分偏移+90度
負(fù)頻率成分偏移-90度

由于對于實數(shù)信號來說，正負(fù)頻率成分共軛，因此對實數(shù)信號進(jìn)行Hilbert變換之后仍然是實數(shù)信號。下面的程序驗證Hilbert變換的頻率響應(yīng)：

>>> 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.])

對信號進(jìn)行N點FFT變換之后：

下標(biāo)為0的頻率分量表示直流分量
下標(biāo)為N/2的的頻率分量為取樣頻率/2的頻率分量
1到N/2-1為正頻率分量
N/2+1到N為負(fù)頻率分量

對照Y/X的虛數(shù)部分，不難看出它是符合Hilbert的頻率響應(yīng)的。如果你用np.real(Y/X)觀察實數(shù)部分的話，它們?nèi)拷咏?。

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

$envelope = \sqrt{H(x)^2 + x^2}$

其中x為原始載波波形，H(x)為x的Hilbert變換之后的波形，envelope為信號x的包絡(luò)。其原理很容易理解：假設(shè)x為正弦波，那么H(x)為余弦波，根據(jù)公式：

$\sin^2(t) + \cos^2(t) = 1$

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

# -*- 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"檢出的包絡(luò)信號")
pl.title(u"使用Hilbert變換進(jìn)行包絡(luò)檢波")
pl.legend()
pl.show()

圖18.12 使用Hilbert變換對載波信號進(jìn)行包絡(luò)檢波

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

>>> 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"頻率掃描波的包絡(luò)")
>>> pl.show()