數字信號系統?

FIR和IIR濾波器?

在數字信號處理領域中，數字濾波器占有非常重要的地位。根據其計算方式可以分為FIR(有限脈沖響應)濾波器，和IIR(無限脈沖響應)濾波器兩種。

FIR濾波器根據如下公式進行計算：

$y[m] & = b[0] x[m] + b[1] x[m-1] + \cdots + b[P] x[n-P]$

IIR濾波器根據如下公式(直接1型)進行計算：

$y[m] & = b[0] x[m] + b[1] x[m-1] + \cdots + b[P] x[m-P] \\ & - a[1] y[m-1] - a[2] y[m-2] - \cdots - a[Q] y[m-Q]$

其中x是輸入信號，數組a和b是濾波器的系數，y是濾波器的輸出。我們可以把FIR濾波器看作是IIR濾波器的一種特殊情況：當系數a都為0時就從IIR濾波器變為了FIR濾波器了。

根據FIR濾波器的計算公式我們可以知道，時刻m的輸出y[m]由時刻m的輸入x[m]以及之前的輸入x[m-1] ... x[m-P]和濾波器的系數b[0] ... b[P]求乘積和而得。而IIR濾波器只不過是再減去之前的輸出y[m-1] ... y[m-Q]和系數a[1] ... a[m-Q]的乘積和。

總之，數字濾波器的計算方法并不復雜，僅僅是數組對應元素的乘積和求和而已。然而其計算量對于Python來說是相當大的：通常FIR濾波器的系數長度都上百，而CD音質的數字聲音信號一秒鐘有44100個取樣值，假設濾波器的長度是100，那么一秒鐘需要計算4百萬次以上的乘積和加法。這對于Python這樣的動態語言來說是很困難的。

因此scipy的signal庫中提供了lfilter函數完成數字濾波器的計算工作。由于它是在C語言級別實現的，因此處理速度相當快：

signal.lfilter(b, a, x, axis=-1, zi=None)

其中的b和a是濾波器的系數，x是輸入。lfilter函數并不是直接使用上面的IIR濾波器計算公式進行計算，而是對其進行了如下的變形：

y[m] = b[0]*x[m] + z[0,m-1]                         (1)
z[0,m] = b[1]*x[m] + z[1,m-1] - a[1]*y[m]           (2)
...
z[n-3,m] = b[n-2]*x[m] + z[n-2,m-1] - a[n-2]*y[m]
z[n-2,m] = b[n-1]*x[m] - a[n-1]*y[m]

這段公式就沒有那么直白了，但是只需要仔細的觀察一下就不難發現，將式(2)的時間變為m-1，得到：

z[0,m-1] = b[1]*x[m-1] + z[1,m-2] - a[1]*y[m-1]     (3)

將其帶入到式(1)中，發現b[0]*x[m], b[1]*x[m-1], - a[1]*y[m-1]等項已經和IIR公式中的一致，依次如此代入下去最后得到的公式和IIR濾波器的公式是一致的。這個計算公式被稱為直接2型。

直接1型的公式中，為了計算m時刻的輸出y[m]，除了需要m時刻的輸入x[m]之外，還需要x[m-1]到x[m-P]和y[m-1]到y[m-Q]，這些值都需要被作為濾波器的內部狀態保存起來，因此需要保存P+Q個數。而根據直接2型的公式，只需要保存n-1個數z[0]到z[n-2]，其中n為 max(len(a), len(b)) ，即max(P, Q)。數組z就是濾波器的狀態。

濾波器的初始狀態通過關鍵字參數zi傳到lfilter函數，當zi不是None時，lfilter將返回濾波器的最終狀態zf，于是其返回值為(y, zf)，如果zi為None的話，那么只返回濾波器的輸出y。

當使用lfilter對很長的輸入進行濾波計算時，不能一次把數據都讀入到數組x中，因此需要對數據進行分段濾波，這時就需要將上一次調用lfilter時返回的數組zf，傳遞到下一次lfilter函數調用。下面的程序演示了這種分段濾波的方法：

# -*- coding: utf-8 -*-
import scipy.signal as signal
import numpy as np
import pylab as pl

# 某個均衡濾波器的參數
a = np.array([1.0, -1.947463016918843, 0.9555873701383931])
b = np.array([0.9833716591860479, -1.947463016918843, 0.9722157109523452])

# 44.1kHz， 1秒的頻率掃描波
t = np.arange(0, 0.5, 1/44100.0)
x= signal.chirp(t, f0=10, t1 = 0.5, f1=1000.0)

# 直接一次計算濾波器的輸出
y = signal.lfilter(b, a, x)

# 將輸入信號分為50個數據一組
x2 = x.reshape((-1,50))

# 濾波器的初始狀態為0， 長度是濾波器系數長度-1
z = np.zeros(max(len(a),len(b))-1, dtype=np.float)
y2 = [] # 保存輸出的列表

for tx in x2:
    # 對每段信號進行濾波，并更新濾波器的狀態z
    ty, z = signal.lfilter(b, a, tx, zi=z)
    # 將輸出添加到輸出列表中
    y2.append(ty)
    
# 將輸出y2轉換為一維數組
y2 = np.array(y2)
y2 = y2.reshape((-1,))

# 輸出y和y2之間的誤差
print np.sum((y-y2)**2)

# 繪圖
pl.plot(t, y, t, y2)
pl.show()

程序所輸出的誤差為0，經過濾波器之后的頻率掃描波形如下圖所示：

經過均衡濾波器之后的頻率掃描波形

此程序中使用的IIR濾波器的系數為二次均衡濾波器的系數，其系數的設計算法將在下節進行介紹。為了觀察濾波器的頻率特性，我們讓它對頻率掃描波進行處理。分別采用一次濾波和分段濾波兩種方式調用lfilter函數，我們看到兩個結果完全一樣。使用分段濾波結合pyaudio庫，我們很容易寫出對聲卡采集的連續的聲音信號進行濾波并輸出的實時濾波程序。

如果將一個脈沖信號輸入到濾波器中，所得到的輸出被稱為濾波器的其脈沖響應。所謂脈沖信號就是在時刻0為1，其余時刻均為0的信號。根據FIR濾波器的公式，FIR濾波器的脈沖響應就是濾波器的系數。而IIR濾波器的脈沖響應就不是很直觀了，下面使用lfilter計算IIR濾波器的脈沖響應，其中的IIR濾波器的系數和前面的一樣（在IPython或者Spyder中運行上面的程序之后，再輸入下面的程序）：

>>> impulse = np.zeros(1000, dtype=np.float)
>>> impulse[0] = 1
>>> h = signal.lfilter(b, a, impulse)
>>> h[-1]
-4.2666825205952273e-12

均衡濾波器的脈沖響應

如果你觀察一下h的具體的值就會發現隨著時間的推移，h越來越小，但是始終不會為0，其脈沖響應是無限長度的，因此才被稱作無限脈沖響應濾波器。如果我們將h當作FIR濾波器的系數對信號x進行濾波的話，得到的結果

>>> y3 = signal.lfilter(h, 1, x)
>>> np.sum((y-y3)**2)
3.7835244127856444e-17

顯然由于h是逐漸衰減的，只要我們測量足夠長的脈沖響應，就可以用FIR濾波器足夠精確地模擬IIR濾波器。下圖顯示的是誤差和FIR濾波器的長度之間的關系。顯然由于IIR濾波器的脈沖響應是呈指數衰減的，因此精度隨著長度呈指數增加，請注意Y軸是對數坐標。

隨著FIR濾波器的長度的增加誤差呈指數減小

此圖的計算程序如下，注意用lfilter計算FIR濾波器時，設置參數a的值為1：

# -*- coding: utf-8 -*-
import scipy.signal as signal
import numpy as np
import pylab as pl

# 某個均衡濾波器的參數
a = np.array([1.0, -1.947463016918843, 0.9555873701383931])
b = np.array([0.9833716591860479, -1.947463016918843, 0.9722157109523452])

# 44.1kHz， 1秒的頻率掃描波
t = np.arange(0, 0.5, 1/44100.0)
x= signal.chirp(t, f0=10, t1 = 0.5, f1=1000.0)
y = signal.lfilter(b, a, x)
ns = range(10, 1100, 100)
err = []

for n in ns:
    # 計算脈沖響應
    impulse = np.zeros(n, dtype=np.float)
    impulse[0] = 1
    h = signal.lfilter(b, a, impulse)
    
    # 直接FIR濾波器的輸出
    y2 = signal.lfilter(h, 1, x)
   
    # 輸出y和y2之間的誤差
    err.append(np.sum((y-y2)**2))

# 繪圖
pl.figure(figsize=(8,4))
pl.semilogy(ns , err, "-o")
pl.xlabel(u"脈沖響應長度")
pl.ylabel(u"FIR模擬IIR的誤差")
pl.show()

FIR濾波器設計?

理想的低通濾波器頻率響應如下圖所示：

理想低通濾波器的頻率響應

其中 $f_s$ 為取樣頻率， $f_c$ 為阻帶頻率。通常為了計算方便，將取樣頻率正規化為1。于是 $f_c$ 的含義就是每個取樣點所包含的信號的周期數，例如0.1表示每個取樣點包含0.1個周期，即一個周期有10個取樣點。根據離散傅立葉變換的公式可以求出此理想低通濾波器的脈沖響應為：

$h_{ideal}(n) = \frac{\sin(2 \pi f_c)}{\pi n} = 2 f_c sinc(2 f_c n)$

其中n為負無窮到正無窮的整數。顯然此脈沖響應不但無限長，而且不滿足因果律，因為輸入信號在0時刻出現的脈沖，而輸出信號卻在0時刻之前就有值了。

這樣的脈沖響應當然無法用FIR濾波器實現，一個最直觀的近似方法就是取 $h_{ideal}$ 中 0<=n<L 的L個值當作低通FIR濾波器的系數。下面的程序計算此低通濾波器的頻率響應：

# -*- coding: utf-8 -*-
import scipy.signal as signal
import numpy as np
import pylab as pl

def h_ideal(n, fc):
    return 2*fc*np.sinc(2*fc*np.arange(-n, n, 1.0))

b = h_ideal(30, 0.25)

w, h = signal.freqz(b)

pl.figure(figsize=(8,3))
pl.plot(w/2/np.pi, 20*np.log10(np.abs(h)))
pl.xlabel(u"正規化頻率 周期/取樣")
pl.ylabel(u"幅值(dB)")
pl.show()

截取sinc函數正時間部分作為脈沖響應的低通濾波器

用freqz計算頻率響應

freqz用于計算數字濾波器的頻率響應，它的調用方式如下：

freqz(b, a=1, worN=None, whole=0, plot=None)

其中b和a是濾波器的系數，worN為所計算的頻率點數，whole為0表示計算頻率的上限為pi，whole為1表示計算頻率的上限為2*pi。

它返回一個組元 (w,h) ，其中w為所有計算了響應的頻率數組，其值為正規化的圓頻率，因此通過w/(2*pi)可以計算出對應的正規化頻率。h是一個復數數組，它表示濾波器系統在每個對應的頻率點的響應。復數的幅值表示濾波器的增益特性，相角表示濾波器的相位特性。

程序中使用freqz計算濾波器的頻率響應，并用 20*np.log10(np.abs(h)) 計算h以dB衡量的幅值。

用firwin設計濾波器?

顯然此頻率響應和理想的低通濾波器相差甚遠，并且即使增加FIR濾波器的系數也沒有作用。因為我們舍棄了n<0的那一半系數，而這些系數有著相當大的影響，因此只截取n>=0的部分是不夠的，如果我們將n<0的那一半系數也添加進濾波器的話，得到的頻率響應將會有很大的改善。如下重新定義h_ideal函數，它返回 $h_{ideal}$ 中-n到n之間的系數：

def h_ideal(n, fc):
    return 2*fc*np.sinc(2*fc*np.arange(-n, n, 1.0))

下面是添加n<0系數之后的頻率響應：

對稱截取sinc函數的低通濾波器

這樣做雖然改善了頻率響應，但是給系統帶來了許多延時，為了頻率響應更好必須增加濾波器的點數，然而為了減少延時，必須減少點數，為了解決這個矛盾，我們給系數乘上一個窗函數，讓它快速收斂。

SciPy提供了firwin用窗函數設計低通濾波器，firwin的調用形式如下：

firwin(N, cutoff, width=None, window='hamming')

其中N為濾波器的長度；cutoff為以 $f_s/2$ 正規化的頻率；window為所使用的窗函數。

下面的程序用firwin設計低通濾波器，并且和上面的結果進行比較，注意由于firwin的cutoff頻率是以取樣頻率/2正規化的，因此它是前面所介紹的 $f_c$ 的兩倍。

# -*- coding: utf-8 -*-
import scipy.signal as signal
import numpy as np
import pylab as pl

def h_ideal(n, fc):
    return 2*fc*np.sinc(2*fc*np.arange(-n, n, 1.0))

b = h_ideal(30, 0.25) # 以fs正規化的頻率
b2 = signal.firwin(len(b), 0.5) # 以fs/2正規化的頻率

w, h = signal.freqz(b)
w2, h2 = signal.freqz(b2)

pl.figure(figsize=(8,6))
pl.subplot(211)
pl.plot(w/2/np.pi, 20*np.log10(np.abs(h)), label=u"h_ideal")
pl.plot(w2/2/np.pi, 20*np.log10(np.abs(h2)), label=u"firwin")
pl.xlabel(u"正規化頻率 周期/取樣")
pl.ylabel(u"幅值(dB)")
pl.legend()
pl.subplot(212)
pl.plot(b, label=u"h_ideal")
pl.plot(b2, label=u"firwin")
pl.legend()
pl.show()

firwin使用窗函數設計的低通濾波器的頻率響應和脈沖響應

使用firwin函數設計的濾波器并不是最優化的，為了實現同樣效果頻率響應，還存在長度更短的FIR濾波器。

用remez設計濾波器?

remez函數能夠幫助我們找到更優的濾波器系數。remez的調用形式如下：

remez(numtaps, bands, desired,
    weight=None, Hz=1, type='bandpass', maxiter=25, grid_density=16)

其中：

numtaps : 所設計的FIR濾波器的長度
bands ：一個遞增序列，它包括頻率響應中的所有頻帶的邊界，其值在0到Hz/2之間，如果參數Hz為缺省值1的話，那么可以把它當作是以取樣頻率正規化的頻率
desired : 長度為bands的一半的增益序列，它給出頻率響應在bands中的每個頻帶的增益值
weight : 長度和desired一樣的權重序列，它給出desired中的每個增益所占的權重，即給出desired中的每個增益的重要性，值越大表示其越重要
type : 'bandpass'或者'differentiator'，本書只介紹type為'bandpass'的情況

remez算法

remez是一種迭代算法，它能夠找到一個n階多項式，使得在指定的區間中此多項式和指定函數之間的最大誤差最小化。由于FIR濾波器的頻率響應實際上是一個多項式函數（請參考下節內容），因此可以用remez算法進行FIR濾波器系數設計。

remez返回經過remez算法最優化之后的FIR濾波器的系數。此系數和用firwin所設計的結果一樣是對稱的。當numtaps為偶數時，所設計的濾波器對于取樣頻率/2的響應為0，因此無法設計出長度為偶數的高通濾波器。

下面的程序演示通過remez設計高通濾波器：

# -*- coding: utf-8 -*-
import scipy.signal as signal
import numpy as np
import pylab as pl

for length in [11, 31, 51, 101, 201]:
    b = signal.remez(length, (0, 0.18,  0.2,  0.50), (0.01, 1))
    w, h = signal.freqz(b, 1)
    pl.plot(w/2/np.pi, 20*np.log10(np.abs(h)), label=str(length))
pl.legend()
pl.xlabel(u"正規化頻率 周期/取樣")
pl.ylabel(u"幅值(dB)")
pl.title(u"remez設計高通濾波器 - 濾波器長度和頻響的關系")
pl.show()

程序中，remez函數的bands參數給出兩個頻帶(以取樣頻率正規化)：0到0.18和0.2到0.5，而desired給出兩個頻帶的增益分別為0.01和1，因此它所設計的是一個通帶頻率為0.2、阻帶增益為-40dB的高通濾波器。

此程序顯示出濾波器長度和頻率響應之間存在如下關系，可以看出濾波器越長，頻率響應越接近設計值：

remez設計的高通濾波器，長度越長頻率響應越接近設計值

下圖顯示權值和頻率響應之間的關系，圖中的濾波器長度為101。我們注意到，當權值為1, 0.01(紅色曲線)時，兩個頻帶的增益抖動量相同，這個權值正好和增益desired的設置相反。這時因為缺省情況下，增益越大的頻帶的頻率響應要求越精確，而當權值和增益的乘積相等時，頻率響應的誤差也就相同了。

remez設計濾波器時權值影響頻率響應

濾波器級聯?

假設有兩個濾波器h1和h2，我們將h1的輸出輸入到h2，這樣得到的濾波器稱為h1和h2的級聯。級聯后的濾波器的脈沖響應為h1和h2的脈沖響應的卷積，而其頻率響應為兩個濾波器的頻率響應的乘積。

下面的程序先用remez分別設計一個高通濾波器h1和一個低通濾波器h2，然后通過卷積計算出它們的級聯濾波器h3的系數：

# -*- coding: utf-8 -*-
import scipy.signal as signal
import numpy as np
import pylab as pl

h1 = signal.remez(201, (0, 0.18,  0.2,  0.50), (0.01, 1))
h2 = signal.remez(201, (0, 0.38,  0.4,  0.50), (1, 0.01))
h3 = np.convolve(h1, h2)

w, h = signal.freqz(h3, 1)
pl.plot(w/2/np.pi, 20*np.log10(np.abs(h)))

pl.legend()
pl.xlabel(u"正規化頻率 周期/取樣")
pl.ylabel(u"幅值(dB)")
pl.title(u"低通和高通級聯為帶通濾波器")
pl.show()

最后使用freqz函數計算h3的頻率響應：

低通和高通濾波器級聯之后是帶通濾波器

可以看出，所得到的是一個帶通濾波器。

我們也可以直接用remez設計帶通濾波器：

>>> h4 = signal.remez(201, (0, 0.18, 0.2, 0.38, 0.4, 0.50), (0.01, 1, 0.01))

如果你觀察此濾波器的頻率響應的話，發現它和h3的基本一致，如果比較h3和h4的話，我們得到如下結果：

級聯的濾波器和remz設計的帶通濾波器的脈沖響應近似

可以看出雖然h3的長度幾乎是h4的兩倍，但是由于它的許多系數都接近于0，因此h3和h4的頻率響應近似相同。

IIR濾波器設計?

通常在設計數字IIR濾波器時，都會先設計一個對應的模擬濾波器，然后通過雙線性變換將模擬濾波器轉換為數字濾波器。這意味著我們需要在s復平面上設計濾波器的傳遞函數H(s)。當H(s)的所有的極點都在s的左半平面時，濾波器的響應是穩定的。下面以巴特沃斯濾波器為例，說明這一設計過程。

巴特沃斯低通濾波器?

巴特沃斯低通濾波器的振幅的平方和頻率之間的關系可以用如下公式表示：

$|H(j\omega)|^2 = \frac {1}{1+\left(\frac{\omega}{\omega_c}\right)^{2n}}$

其中n為濾波器的階數， $\omega_c$ 為振幅下降3dB時的截止頻率。這個公式很容易理解：

當 $\omega<\omega_c$ 時， $\omega$ 越小，振幅越接近于1
當 $\omega>\omega_c$ 時， $\omega$ 越大，振幅越接近于0
隨著n的增大，振幅接近于1或者0的速度將變快，即n越大，低通濾波器在阻頻帶的衰減速度將越快
當 $\omega=\omega_c$ 時，振幅的平方為1/2，即-3dB

下面我們推導出巴特沃斯低通濾波器的傳遞函數H(s)，其中 $s = \delta + j \omega$ ，為復數平面上的點。

由于當 $s =j \omega$ 時， $H(s)H(-s) = |H(j\omega)|^2$ ，因此將 $\omega =s/j$ 帶入到巴特沃斯低通濾波器的振幅平方公式中可以得到：

$H(s)H(-s) = \frac {{G_0}^2}{1+\left (\frac{-s^2}{\omega_c^2}\right)^n}$

此公式有2n個極點，其中n個在左半平面，n個在右半平面，由于H(s)必須是穩定的，因此左半平面的n的極點屬于H(s)。

最后得到的傳遞函數為：

$H(s)=\frac{1}{\prod_{k=1}^n (s-s_k)/\omega_c}$

其中 $s_k$ 為左半平面上的極點：

$s_k = \omega_c e^{\frac{j(2k+n - 1)\pi}{2n}}\qquad\mathrm{k = 1,2,3, \ldots, n}$

下面的程序繪制6、7階巴特沃斯低通濾波器的S復平面上的極點：

# -*- coding: utf-8 -*-
from scipy import signal
import numpy as np
import matplotlib.pyplot as pl

pl.figure(figsize=(5,5))

b, a = signal.butter(6, 1.0, analog=1)
z,p,k = signal.tf2zpk(b, a)
pl.plot(np.real(p), np.imag(p), '^', label=u"6階巴特沃斯極點")

b, a = signal.butter(7, 1.0, analog=1)
z,p,k = signal.tf2zpk(b, a)
pl.plot(np.real(p), np.imag(p), 's', label=u"7階巴特沃斯極點")

pl.axis("equal")
pl.legend(loc="center right")
pl.show()

程序中，使用butter函數設計巴特沃斯濾波器，缺省情況下它設計的是數字濾波器，為了設計模擬濾波器，需要傳遞關鍵字參數analog=1。獲得傳遞函數的b和a的系數之后，通過tf2zpk函數將它們轉換為零點和極點：

巴特沃斯低通濾波器在S復平面上的極點分布

雙線性變換?

有了連續時間的傳遞函數H(s)，下一步就是如何將它轉換為離散時間的傳遞函數H(z)。轉換的方法有幾種，其中最常用的是雙線性變換，其變換公式為：

$s \leftarrow \frac{2}{T} \frac{z - 1}{z + 1}$

其中T為離散時間的取樣周期。雙線性變換公式的推導過程請參考下面的鏈接：

雙線性變換公式推導： http://en.wikipedia.org/wiki/Bilinear_transform

雙線性變換實際上是s復平面和z復平面上的點的映射變換，他將s復平面上的豎線變換成z復平面上的圓，而s復平面上的Y軸對應于z復平面上的單位圓。下面的程序演示了這一對應關系：

# -*- coding: utf-8 -*-
import numpy as np
import pylab as pl

def stoz(s):
    """
    將s復平面映射到z復平面
    為了方便起見，假設取樣周期T=1
    """
    return (2+s)/(2-s)
    
def make_vline(x):
    return x + 1j*np.linspace(-100.0,100.0,20000)
    
fig = pl.figure(figsize=(7,3))    
axs = pl.subplot(121)
axz = pl.subplot(122)
for x in np.arange(-3, 4, 1):
    s = make_vline(x)
    z = stoz(s)
    axs.plot(np.real(s), np.imag(s))
    axz.plot(np.real(z), np.imag(z))

axs.set_xlim(-4,4)
axz.axis("equal")
axz.set_ylim(-3,3)

pl.show()

程序中的stoz函數是將s變換為z的變換公式，只需要對上述的雙線性變換公式稍作變形即可得到，這里為了方便起見，假設取樣周期T=1。

雙線性變換將s平面(左圖)上的豎線變換為z平面上的圓(右圖)

通過雙線性變換之后，濾波器的頻率響應會發生變化。在下一節中我們會介紹，離散時間的濾波器的頻率響應是將其傳遞函數H(z)用 $z=e^{j\omega T}$ 進行替換。將其帶入到雙線性變換公式得到：

$s = \frac{2}{T} \frac{z - 1}{z + 1} = \frac{2}{T} \frac{e^{j\omega T} - 1}{e^{j\omega T} + 1} = \frac{2}{T} \frac{e^{j \omega T/2} - e^{-j \omega T/2}}{e^{j \omega T/2} + e^{-j \omega T/2}} = j \frac{2}{T} \tan{(\omega T / 2)}$

對于s平面上的點來說，過原點的豎線 $s=j \omega$ 就是連續時間傳遞函數H(s)的頻率響應，因此雙線性變換將離散時間的頻率 $\omega$ 通過如下的公式轉換為連續時間的頻率 $\omega_a$ :

$\omega_a = \frac{2}{T} \tan \left( \omega \frac{T}{2} \right)$

而其反函數為：

$\omega = \frac{2}{T} \arctan \left( \omega_a \frac{T}{2} \right)$

下面讓我們用程序來驗證這個頻率轉換公式。首先我們載入所需要的庫，并且定義離散時間的取樣頻率fs為8kHz，設計的巴特沃斯低通濾波器通帶截至頻率為1kHz：

>>> from scipy import signal
>>> from numpy import *
>>> fs = 8000.0
>>> f = 1000.0

下面使用butter函數設計一個3階的巴特沃斯濾波器，注意關鍵字參數analog=1，表示設計連續時間傳遞函數H(s)的系數，由于通帶頻率參數為圓頻率，因此需要乘以2*pi：

>>> b, a = signal.butter(3, 2*pi*f, analog=1)

模擬濾波器的系數b和a和H(s)的關系

假設b和a的長度分別為M和N，模擬濾波器的系數中b[0]為分子中 $s^{M-1}$ 項的系數，a[0]為分母中 $s^{N-1}$ ，b[-1]和a[-1]為分子分母的常數項的系數，即：

$H(s) = \frac{b_0 s^{M-1} + b_1 s^{M-2} + \ldots + b_{M-1}}{a_0 s^{N-1} + a_1 s^{N-2} + \ldots + a_{N-1}}$

然后調用雙線性變換函數bilinear，將系數轉換為離散時間的傳遞函數系數，通過關鍵字參數fs指定取樣頻率：

>>> b2, a2 = signal.bilinear(b,a,fs=fs)

接下來調用freqz函數得到此數字濾波器的頻率響應，為了得到盡可能精確的值，我們通過worN關鍵字參數讓它計算10000點的頻率響應：

>>> w2, h2 = signal.freqz(b2,a2,worN=1000)

接下來將h2轉換為增益，并且找到增益為-3dB(精確值為10*log10(0.5))時所對應的正規化圓頻率w的下標idx，w/(2*pi)*fs就是其對應的實際頻率值：

>>> p2 = 20*log10(abs(h2))
>>> idx = argmin(abs(p2-10*log10(0.5)))
>>> w2[idx]/2/pi*8000
952.8

通過頻率轉換公式得到的頻率為：

>>> 2*fs*arctan(2*pi*f/2/fs) /2/pi
952.8840223

實際使用scipy.signal庫設計IIR濾波器沒有這么麻煩，因為它所提供的濾波器設計函數缺省都是直接設計數字濾波器。這些函數設計數字濾波器時采用的取樣頻率為2，即以香農頻率fs/2為1進行正規化。因此要設計取樣頻率為fs、通帶頻率為f的濾波器需要將通帶頻率正規化為f/(fs/2)，下面調用butter函數設計數字低通濾波器，這里使用上述計算所得的通帶頻率：

>>>  b3,a3 = signal.butter(3, 952.8840223/(fs/2))
>>> sum(abs(b3-b2))
1.3226225670237568e-13
>>> sum(abs(a3-a2))
7.0876637892069994e-13

數字濾波器的系數b和a和H(z)的關系

假設b和a的長度分別為M和N，數字濾波器的系數中b[0]和a[0]分別為分子分母中常數項的系數，a[-1]為分母中 $z^{-(N-1)}$ ，b[-1]為分子中 $z^{-(M-1)}$ 的系數，即：

$H(s) = \frac{b_0 + b_1 z^{-1} + \ldots + b_{M-1} z^{-(M-1)}}{a_0 + a_1 z^{-1} + \ldots + a_{N-1} z^{-(N-1)}}$

所得的濾波器的系數b3和a3與手工通過bilinear函數計算的系數b2和a2是一致的。在signal庫設計數字濾波器時，其內部會先通過頻率轉換公式對頻率進行轉換，然后設計連續時間的傳遞函數系數，最后通過bilinear函數進行系數轉換。有興趣的讀者可以查看signal.iirfilter函數的源代碼。

濾波器的頻帶轉換?

只要知道了低通濾波器的傳遞函數H(s)，就很容易利用變量替換設計出同樣階數的高通、帶通或者其它通帶頻率的低通濾波器。讓我們來看看低通濾波器的變換。

假設我們使用巴特沃斯低通濾波器設計公式，設計出通帶頻率為1弧度/秒的標準低通濾波器：

>>> b, a = signal.butter(2, 1.0, analog=1)
>>> np.real(b)
array([ 1.])
>>> np.real(a)
array([ 1.        ,  1.41421356,  1.        ])

$H(s) = \frac{1}{s^2 + 1.4142 s + 1}$

為了讓它變為通帶頻率為 $\omega_c$ 的低通濾波器，只需要進行如下替換：

$s \to \frac{s}{\omega_c}$

由于當 $s =j \omega$ 時，H(s)就是濾波器的頻率響應。因此所設計的標準低通濾波器H(s)在 $s =j$ 時振幅下降3dB的，而 $H(s/\omega_c)$ 則在 $s =j\omega_c$ 處下降3dB。下面是通帶頻率為2弧度/秒的2階低通濾波器的系數：

>>> b2, a2 = signal.butter(2, 2.0, analog=1)
>>> np.real(b2)
array([ 4.])
>>> np.real(a2)
array([ 1.        ,  2.82842712,  4.        ])

可以看出將 $s \to \frac{s}{2}$ 代入到前面的H(s)中即可得到這些系數。

低通濾波器轉高通濾波器的替代公式為：

$s \to \frac{\omega_c}{s}$

此替代公式很容易理解：

若 $\omega$ 為0，則替代之后的頻率為無窮大，而低通濾波器無窮大處的頻率響應為0，即轉換之后的濾波器在0處的頻率響應為0；
若 $\omega$ 為無窮大，則替代之后的頻率為0，因此轉換之后的濾波器在無窮大處頻率響應為1。

下面設計通帶頻率為1弧度/秒的高通濾波器：

>>> b3,a3 = signal.butter(2,1.0,btype="high",analog=1)
>>> np.real(b3)
array([ 1.,  0.,  0.])
>>> np.real(a3)
array([ 1.        ,  1.41421356,  1.        ])

可以看出這些系數是將 $s \to \frac{1}{s}$ 帶入到H(s)中之后，上下分母乘以 $s^2$ 之后得到的。

低通濾波器還可以轉換為帶通濾波器，這可能有點難以理解，讓我們先來看看替代公式，假設帶通濾波器的高低通帶頻率為 $\omega_2$ 和 $\omega_1$ ：

$s \to \frac{\omega_0}{\Delta\omega} \left( \frac{s}{\omega_0} + \frac{\omega_0}{s} \right)$

其中 $\Delta\omega=\omega_2-\omega_1$ ， $\omega_0=\sqrt{\omega_1\omega_2}$ 。 $\Delta\omega$ 被稱為通帶帶寬，而 $\omega_0$ 則是通帶的中心頻率。

讓我們通過下面的程序來研究一下為何這種替代能夠將低通映射為帶通濾波器：

# -*- coding: utf-8 -*-
import numpy as np
from scipy import signal
import pylab as pl

b, a = signal.butter(2, 1.0, analog=1)

# 低通->帶通的頻率變換函數
w1 = 1.0 # 低通帶頻率
w2 = 2.0 # 高通帶頻率
dw = w2 - w1 # 通帶寬度
w0 = np.sqrt(w1*w2) # 通帶中心頻率

# 產生10**-2到10**2的頻率點
w = np.logspace(-2, 2, 1000)

# 使用頻率變換公式計算出轉換之后的頻率
nw = np.imag(w0/dw*(1j*w/w0 + w0/(1j*w)))

_, h = signal.freqs(b, a, worN=nw)
h = 20*np.log10(np.abs(h))

pl.figure(figsize=(8,5))

pl.subplot(221)
pl.semilogx(w, nw) # X軸使用log坐標繪圖
pl.xlabel(u"變換前圓頻率(弧度/秒)")
pl.ylabel(u"變換后圓頻率(弧度/秒)")

pl.subplot(222)
pl.plot(h, nw)
pl.xlabel(u"低通濾波器的頻率響應(dB)")

pl.subplot(212)
pl.semilogx(w, h)
pl.xlabel(u"變換前圓頻率(弧度/秒)")
pl.ylabel(u"帶通濾波器的頻率響應(dB)")

pl.subplots_adjust(wspace=0.3, hspace=0.3, top=0.95, bottom=0.14)

print "center:", w[np.argmin(np.abs(nw))]
pl.show()

程序中先使用butter函數設計一個模擬的二階標準低通濾波器：

b, a = signal.butter(2, 1.0, analog=1)

我們要將其轉換為通帶頻率為w1=1到w2=2的帶通濾波器：

# 低通->帶通的頻率變換函數
w1 = 1.0 # 低通帶頻率
w2 = 2.0 # 高通帶頻率
dw = w2 - w1 # 通帶帶寬
w0 = np.sqrt(w1*w2) # 通帶中心頻率

假設我們關心的頻率響應的頻率段為0.01到100，使用logspace函數產生一個這個區間的等比數列w：

w = np.logspace(-2, 2, 1000)

通過帶通頻率轉換公式將其轉換為新的頻率序列nw：

nw = np.imag(w0/dw*(1j*w/w0 + w0/(1j*w)))

使用此新的頻率序列nw計算出每個頻率點對應的低通濾波器的頻率響應，注意我們通過worN關鍵字傳輸讓freqs函數計算指定頻率的頻率響應：

_, h = signal.freqs(b, a, worN=nw)
h = 20*np.log10(np.abs(h))

下面的語句就繪制出w和h的關系，也就是帶通濾波器的頻率響應：

pl.semilogx(w, h)

下圖是程序的輸出。其中左上圖繪制的是頻率轉換函數，右上圖繪制的是低通濾波器的頻率響應(X軸為響應，Y軸為頻率)，最下面繪制的是最終的帶通濾波器的頻率響應。

使用頻率轉換公式將低通變為帶通濾波器

由于帶通的頻率轉換公式將0到無窮大映射到正無窮到負無窮，而低通濾波器在正負無窮處的頻率響應都為0，因此可以想象轉換后的濾波器是一個帶通濾波器。而轉換之后的頻率nw為0的點所對應的原始頻率就是帶通濾波器的中心頻率，此處的頻率響應為1，下面的程序找到nw絕對值最小的下標，并輸出其對應的轉換前的頻率，我們看到它和w0是一致的：

>>> print w[np.argmin(np.abs(nw))]
1.4130259906

事實上，scipy.signal庫已經為我們提供了頻帶轉換的函數：

lp2lp : 低通轉低通
lp2hp : 低通轉高通
lp2bp : 低通轉帶通
lp2bs : 低通轉帶阻，轉帶阻的公式留給讀者思考

下面以lp2bp為例簡要說明一下函數的用法，假設b,a為二階標準低通濾波器，下面的語句將轉換為通帶為1到2弧度的帶通濾波器，前兩個參數為濾波器的系數，后兩個參數分別為中心頻率和通帶帶寬：

>>> b, a = signal.butter(2, 1.0, analog=1)
>>> b3, a3 = signal.lp2bp(b,a,np.sqrt(2), 1)

我們也可以直接調用butter設計一個低通濾波器：

>>> b4, a4 = signal.butter(2, [1,2], btype='bandpass', analog=1)

兩個結果是完全一致的：

>>> np.all(b3==b4)
True
>>> np.all(a3==a4)
True

濾波器的頻率響應?

前面的許多例子中都使用函數freqz計算濾波器的頻率響應，在這一節中，讓我們來深入研究一下freqz是如何計算頻率響應的。在IPython中輸入：

>>> import scipy.signal
>>> signal.freqz??

即可看到freqz的完整的實現函數，下面是其完整的源程序(除去文檔說明):

def freqz(b, a=1, worN=None, whole=0, plot=None):
    b, a = map(atleast_1d, (b,a))
    if whole:
        lastpoint = 2*pi
    else:
        lastpoint = pi
    if worN is None:
        N = 512
        w = numpy.arange(0,lastpoint,lastpoint/N)
    elif isinstance(worN, types.IntType):
        N = worN
        w = numpy.arange(0,lastpoint,lastpoint/N)
    else:
        w = worN
    w = atleast_1d(w)
    zm1 = exp(-1j*w)
    h = polyval(b[::-1], zm1) / polyval(a[::-1], zm1)
    if not plot is None:
        plot(w, h)
    return w, h

研究一下這段代碼，不難發現真正的計算頻率響應的代碼可以用如下3行程序概括：

w = numpy.arange(0,pi,pi/N)
zm1 = exp(-1j*w)
h = polyval(b[::-1], zm1) / polyval(a[::-1], zm1)

為了弄清楚為什么這3行代碼能夠計算濾波器的頻率響應，讓我們先來學習一下相關的理論知識。

濾波器的頻率響應由濾波器的傳遞函數給出，IIR濾波器的計算公式如下：

$y[m] & = b[0] x[m] + b[1] x[m-1] + \cdots + b[P] x[m-P] \\ & - a[1] y[m-1] - a[2] y[m-2] - \cdots - a[Q] y[m-Q]$

根據Z變換的相關公式，容易求得其傳遞函數為：

$H(z) = \frac{Y(z)}{X(z)} = \frac{b[0] + z^{-1}b[1] + z^{-2}b[2] + \cdots + z^{-M}b[P]}{1 + z^{-1}a[1] + z^{-2}a[2] + \cdots + z^{-N}a[Q]}$

其中z為復數平面上的任意一點。當z為單位圓上的點，即 $z=e^{j\omega}$ 時， $H(\omega)$ 就是濾波器的頻率響應。 $\omega$ 被稱之為圓頻率，當其取值從0到2*pi變化時， $e^{j\omega}$ 正好繞復數平面單位圓轉一圈。由于復數平面上下兩個半平面的復數存在共軛關系，因此通常只需要求上半圓的頻率響應，因此下面的語句將上半圓等分為N份：

w = numpy.arange(0,pi,pi/N)

然后計算w中每點對應的復數值 $z^{-1}$ ，注意這里將負號帶入，于是傳遞函數的分子分母部分就都變成了zm1的多項式函數：

zm1 = exp(-1j*w)

最后帶入傳遞函數的公式中計算出頻率響應h：

h = polyval(b[::-1], zm1) / polyval(a[::-1], zm1)

polyval(p, x)函數對于數組x中的每個元素計算多項式p的值，其計算公式如下：

p[0]*(x**N-1) + p[1]*(x**N-2) + ... + p[N-2]*x + p[N-1]

由于濾波器系數b和a的順序正好和polyval的多項式系數p的順序相反，因此通過數組切片運算b[::-1]將濾波器的系數反轉。由于數組zm1中的值都為復數，因此所得到的頻率響應h的值也都是復數。復數的幅值對應于頻率響應中的增益特性，而其相角對應于頻率響應中相位特性。

freqz中在Z平面單位圓上所取的點是等距線性的，然而我們經常需要在繪制頻率響應圖表時要求頻率坐標為對數坐標，對于對數坐標，等距的頻率點會造成低頻過疏，高頻過密的問題，因此我們可以如下改造freqz函數，使其更適合計算對數頻率坐標的頻率響應。

# -*- coding: utf-8 -*-
import numpy as np
import pylab as pl
import scipy.signal as signal

def logfreqz(b, a, f0, f1, fs, N):
    """
    以對數頻率坐標計算濾波器b,a的頻率響應
    f0, f1: 計算頻率響應的開始頻率和結束頻率
    fs: 取樣頻率
    """
    w0, w1 = np.log10(f0/fs*2*np.pi), np.log10(f1/fs*2*np.pi)
    # 不包括結束頻率
    w = np.logspace(w0, w1, N, endpoint=False)
    zm1 = np.exp(-1j*w)
    h = np.polyval(b[::-1], zm1) / np.polyval(a[::-1], zm1)
    return w/2/np.pi*fs, h
    
for n in range(1, 6):    
    # 設計n階的通頻為0.1*4000 = 400Hz的高通濾波器
    b, a = signal.iirfilter(n, [0.1, 1])    
    f, h = logfreqz(b, a, 10.0, 4000.0, 8000.0, 400)
    gain = 20*np.log10(np.abs(h))
    pl.semilogx(f, gain, label="N=%s" % n)
    slope = (gain[100]-gain[10]) / (np.log2(f[100]) - np.log2(f[10]))
    print "N=%s, slope=%s dB" % (n, slope)
pl.ylim(-100, 20)
pl.xlabel(u"頻率(Hz)")
pl.ylabel(u"增益(dB)")
pl.legend()
pl.show()

程序中的logfreqz函數計算系數為b和a的濾波器在f0到f1之間的頻率響應，fs為取樣頻率，N為計算的點數。首先通過f/fs*2*pi將實際頻率轉換為與之對應的圓頻率。然后通過logspace函數計算頻率點的等比數列。最后和freqz一樣通過調用polyval計算頻率響應，返回值為實際頻率點和對應的頻率響應。

接下來通過調用iirfilter設計5個不同階的IIR高通濾波器，通頻為0.1，如果取樣頻率為8kHz的話，那么實際的通頻為0.1*4kHz=400Hz。5個IIR濾波器的增益特性如下圖所示：

iirfilter設計5個不同階的IIR高通濾波器

由此圖可知，隨著IIR濾波器的階數的增加，增益的下降速度增加，程序中第25行計算出下降處兩個倍頻之間的增益差值，其結果如下：

N=1, slope=5.9955865774 dB
N=2, slope=12.0417201051 dB
N=3, slope=18.0630802032 dB
N=4, slope=24.0841135443 dB
N=5, slope=30.1051375912 dB

即IIR濾波器的階數每增加1，其增益的下降速度增加 6dB/oct (6dB每倍頻)。并且所有曲線的相交于一點：此處的頻率正好是400Hz，增益為-3dB。

二次均衡器設計工具?

也許很少有人知道，無論是古老的盒式錄音機還是現代的流行數碼音響設備，以及眾多的音樂播放軟件，其中絕大多數的均衡器都只是一系列簡單的二次IIR濾波器組合而成。

二次IIR濾波器的傳遞函數如下：

$H(z) = \frac{b0 + b1 \cdot z^{-1} + b2 \cdot z^{-2}}{a0 + a1 \cdot z^{-1} + a2 \cdot z^{-2}}$

當a0的值不為1時，可以將所有系數b0, b1, b2, a0, a1, a2都除以a0，這樣就能得到a0=1的五個系數： b0, b1, b2, a1, a2。即二次均衡濾波器的頻率響應曲線由這5個獨立的參數決定。其頻率響應如下圖所示：

二次均衡濾波器的頻率響應

從圖中可以看出，根據頻率響應曲線的兩個參數：中心頻率f0，振幅峰值peak：

頻率為0時頻率響應為1
頻率為取樣頻率一半時(圓頻率為pi)，頻率響應為1
頻率為中心頻率f0時，振幅達到峰值peak
頻率響應在中心頻率f0處的導數為0

這樣就有了4個方程，再加上一個Q值決定振幅駝峰的寬度，因此一共5個方程決定5個系數。

Audio EQ Cookbook ( http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt ) 中提供了二次均衡器系數的設計公式，根據這個設計手冊，很容易寫出如下的設計均衡器參數的函數：

# -*- coding: utf-8 -*-

import scipy.signal as signal
import pylab as pl
import math
import numpy as np

def design_equalizer(freq, Q, gain, Fs):
    '''設計二次均衡濾波器的系數'''
    A = 10**(gain/40.0)
    w0 = 2*math.pi*freq/Fs
    alpha = math.sin(w0) / 2 / Q
    
    b0 = 1 + alpha * A
    b1 = -2*math.cos(w0)
    b2 = 1 - alpha * A
    a0 = 1 + alpha / A
    a1 = -2*math.cos(w0)
    a2 = 1 - alpha / A
    return [b0/a0,b1/a0,b2/a0], [1.0, a1/a0, a2/a0]    
    
pl.figure(figsize=(8,4))
for freq in [1000, 2000, 4000]:  
    for q in [0.5, 1.0]:
        for p in [5, -5, -10]:
            b,a = design_equalizer(freq, q, p, 44100)
            w, h = signal.freqz(b, a)
            pl.semilogx(w/np.pi*44100, 20*np.log10(np.abs(h)))
pl.xlim(100, 44100)      
pl.xlabel(u"頻率(Hz)")
pl.ylabel(u"振幅(dB)")
pl.subplots_adjust(bottom=0.15)
pl.show()

使用上節介紹的對數頻率響應的求法以及TraitsUI和Chaco等界面庫，我們可以設計如下界面的二次均衡器設計程序：

二次均衡器設計工具的界面

用戶可以使用此程序添加、刪除和編輯二次均衡器，并且即時查看均衡器級聯之后的頻率響應。完整的程序請參照： 二次均衡器設計

ScrubberEditor的BUG

如果界面中的ScrubberEditor無法用鼠標拖動修改的話，那么你需要修改site-packages目錄下的 scrubber_editor.py 文件：

TraitsBackendWX-3.2.0-py2.6.egg\enthought\traits\ui\wx\scrubber_editor.py

按照如下打星號的行進行修改：

# Establish the slider increment:
increment = self.factory.increment
if increment <= 0.0:
    if (low is None) or (high is None) or isinstance( low, int ):
        increment = 1.0
    else:
        increment = pow( 10, round( log10( (high - low) / 100.0 ) ) )

self.increment = increment # ** 將此行移出if作用域

數字信號系統?

FIR和IIR濾波器?

FIR濾波器設計?

用firwin設計濾波器?

用remez設計濾波器?

濾波器級聯?

IIR濾波器設計?

巴特沃斯低通濾波器?

雙線性變換?

濾波器的頻帶轉換?

濾波器的頻率響應?

二次均衡器設計工具?

首頁目錄

上一篇文章

下一篇文章

Navigation

數字信號系統?

FIR和IIR濾波器?

FIR濾波器設計?

用firwin設計濾波器?

用remez設計濾波器?

濾波器級聯?

IIR濾波器設計?

巴特沃斯低通濾波器?

雙線性變換?

濾波器的頻帶轉換?

濾波器的頻率響應?

二次均衡器設計工具?

首頁目錄

上一篇文章

下一篇文章

快速搜索

Navigation