文章目錄
- 先導知識
- 時頻域互換
- 什么是濾波器?
- 如何設計濾波器?
- FT、LT、ZT之間的(映射)關系
- 零、極點分布如何影響頻率回應
- 從均值濾波看濾波器及系統函式
- 時域運算式
- 頻率運算式
- 零極點分布圖
- 頻率回應分析
- FIR和IIR濾波器的比較
- 性能上
- 結構上
- 設計工具上
- 效果上
- FIR濾波器的設計
- 基本設計思路
- 窗函式法
- 頻率采樣法
- 回應最優法
- 最小二乘法
- 切比雪夫等波紋最佳逼近法
- FIR數字濾波器設計總結
- IIR濾波器設計
- 基本設計思路
- 模擬濾波器的設計
- 巴特沃斯濾波器
- 切比雪夫I型濾波器
- 橢圓濾波器
- 高通、帶通和帶阻濾波器的設計
- 歸一化低通濾波器
- 低通到高通的頻率轉換
- 低通到帶通的頻率轉換
- 低通到帶阻的頻率轉換
- 從模擬到數字的轉換方法
- 脈沖回應不變法
- 基本思路
- 推導程序
- 雙線性變換法
- 基本思路
- 轉換性能分析
- IIR濾波器設計實體
- IIR數字濾波器設計總結
說明:本文并非全部原創,是在網上學習各位前輩的文章后寫的學習總結,文中有些圖片和代碼使用了網上的文章,若侵比刪,所參考文章的鏈接已在文中給出!
先導知識
時頻域互換
從時域我們很難看出信號包含那些頻率成分,也很難去掉信號中某些特定的成分,如果能將信號表達成不同頻率成分之和,那么不就很容易看出該信號包含那些頻率成分了嗎,此時,就需要傅里葉級數(周期信號)、傅里葉變換(非周期信號)等,對于離散信號,需要DTFT(離散時間傅里葉變換)、DFS(離散傅里葉級數、離散周期信號)、DFT(離散傅里葉變換、離散非周期信號)等,
-
連續(模擬)信號
-
FS
將周期連續信號變換到頻域,
-
FT
將非周期連續信號變換到頻域,
-
LT
變換到復平面,及s平面,
-
-
離散(數字)信號
-
DTFT
將數字信號變換到頻域,此時頻域是連續的,
-
DFS
將周期數字信號變換到頻域,
-
DFT
將數字信號變換到頻域,此時頻域也是離散的,
-
ZT
變換到z平面,
-
什么是濾波器?
濾波器就是濾掉干擾資訊(噪聲),保留有用資訊,
下面分別是按濾掉噪聲頻率的幾類濾波器的幅頻曲線圖:

設計常規濾波器的時候,我們一般采用另外一種分類,FIR(Finite Impulse Response)和IIR(Infinite Impulse Response)filter,即有限脈沖回應濾波器和無限脈沖回應濾波器,
理想脈沖信號,其傅里葉變換恒為1,也就時包含了所有頻率分量,是一個理想的測驗信號,能夠激發出所有單位頻率分量的回應,因此理想脈沖信號的回應,就代表了系統的特性,
濾波器也可以看成一個系統,如果用一個理想脈沖信號激勵,就會有輸出,我們把輸出個數有限的稱為有限脈沖回應濾波器(FIR);輸出無限多的稱為無限脈沖回應濾波器(IIR),
濾波器分為模擬濾波器和數字濾波器,模擬濾波器的理論和設計方法已經非常成熟,可以通過模擬濾波器間接設計數字濾波器,
如何設計濾波器?
就如第一幅圖一樣,我們可以通過其幅頻曲線得到其截止頻率、通帶情況,那么怎么獲得幅頻曲線圖呢?那就是通過系統函式,
H
(
j
ω
)
=
∣
H
(
j
ω
)
∣
e
θ
(
ω
)
H(j\omega)=\left|H(j\omega)\right|e^{\theta(\omega)}
H(jω)=∣H(jω)∣eθ(ω)
∣
H
(
j
ω
)
∣
|H(j\omega)|
∣H(jω)∣就是其幅頻特性曲線,
θ
(
ω
)
\theta(\omega)
θ(ω)就是其相頻特性曲線,
如何得到系統函式:
由系統的特性分析得到微分方程或差分方程,再拉普拉斯變換或z變換,
FT、LT、ZT之間的(映射)關系
LT是對連續FT的拓展,即從頻率到復頻域, j ω ? s = δ + j ω j{\omega}{\Longrightarrow}s={\delta}+j{\omega} jω?s=δ+jω,
ZT是對離散FT的拓展,也可以看是對LT的再映射,即 z = e s T s z=e^{sT_s} z=esTs?,
如何從DTFT到ZT,請看,簡單說一下:
D
F
T
:
X
(
j
ω
)
=
∑
n
=
?
∞
+
∞
x
(
n
)
e
j
ω
n
T
s
DFT:X(j\omega)=\sum_{n=-\infin}^{+\infin}x(n)e^{j{\omega}nT_s}
DFT:X(jω)=n=?∞∑+∞?x(n)ejωnTs?
同FT到LT一樣,為了滿足絕對可和條件,乘上一個衰減因子
e
δ
t
=
e
δ
n
T
s
e^{{\delta}t}=e^{{\delta}nT_s}
eδt=eδnTs?(離散化),
X
(
j
ω
)
=
∑
n
=
?
∞
+
∞
x
(
n
)
e
j
ω
n
T
s
e
δ
n
T
s
=
∑
n
=
?
∞
+
∞
x
(
n
)
e
?
(
δ
+
j
ω
)
n
T
s
X(j\omega)=\sum_{n=-\infin}^{+\infin}x(n)e^{j{\omega}nT_s}e^{{\delta}nT_s}\\ =\sum_{n=-\infin}^{+\infin}x(n)e^{-({\delta}+j{\omega})nT_s}
X(jω)=n=?∞∑+∞?x(n)ejωnTs?eδnTs?=n=?∞∑+∞?x(n)e?(δ+jω)nTs?
令
z
=
e
(
δ
+
j
ω
)
T
s
z=e^{({\delta}+j{\omega})T_s}
z=e(δ+jω)Ts?,得到DFT的式子為:
X
(
j
ω
)
=
∑
n
=
?
∞
+
∞
x
(
n
)
z
?
n
X(j{\omega})=\sum_{n=-\infin}^{+\infin}x(n)z^{-n}
X(jω)=n=?∞∑+∞?x(n)z?n
上式就是Z變換了,


如何理解上面的規律?
- 需要明白幅頻曲線自變數是頻率 ω \omega ω,而不是 s = δ + j ω s={\delta}+j{\omega} s=δ+jω,
- H ( j ω ) = ∣ H ( j ω ) ∣ e θ ( j ω ) = N u m ( j ω ) / D e n ( j ω ) H(j{\omega})=|H(j{\omega})|e^{{\theta}(j{\omega})}=Num(j{\omega})/Den(j{\omega}) H(jω)=∣H(jω)∣eθ(jω)=Num(jω)/Den(jω),其中的 ∣ H ( j ω ) ∣ |H(j\omega)| ∣H(jω)∣就是幅頻曲線、 θ ( j ω ) \theta{(j\omega)} θ(jω)就是相頻回應,
- s平面的虛軸(頻率軸)被ZT映射到z平面的單位圓上,
因此,分析幅頻回應時,頻率軸在單位圓上,這也是為什么上面提到的規律都在說單位圓上的點怎么怎么樣,
- 對于規律一,因為零點就在單位圓上,因此當頻率為零點對應的頻率時, ∣ H ( j ω ) ∣ = 0 |H(j\omega)|=0 ∣H(jω)∣=0,因此幅值回應為零,
- 對于規律二,不在單位圓上的零點,單位圓上越靠近零點的頻率,其 ∣ H ( j ω ) ∣ |H(j\omega)| ∣H(jω)∣也就越小,
- 對于規律三,因為是極點,極點為0時即系統函式的分母為0, ∣ H ( j ω ) ∣ |H(j\omega)| ∣H(jω)∣很大,(至于為什么要指定是單位圓內部的極點呢?我猜是單位圓外部的極點使系統不穩定,因此不必討論)

零、極點分布如何影響頻率回應
H ( z ) = A ∏ r = 1 M ( 1 ? c r z ? 1 ) ∏ r = 1 N ( 1 ? d r z ? 1 ) ∣ H ( e j ω ) ∣ = ∣ A ∣ ∏ r = 1 M ∣ z r ∣ ∏ r = 1 N ∣ P r ∣ ? ( ω ) = ω ( N ? M ) + ∑ r = 1 N α r ? ∑ r = 1 M β r H(z)=A\frac{\prod_{r=1}^M(1-c_rz^{-1})}{\prod_{r=1}^{N}(1-d_rz^{-1})} \\ |H(e^{j\omega})|=|A|\frac{\prod_{r=1}^M|z_r|}{\prod_{r=1}^N{|P_r|}} \\ \phi(\omega)=\omega(N-M)+\sum_{r=1}^{N}\alpha_r-\sum_{r=1}^{M}\beta_r H(z)=A∏r=1N?(1?dr?z?1)∏r=1M?(1?cr?z?1)?∣H(ejω)∣=∣A∣∏r=1N?∣Pr?∣∏r=1M?∣zr?∣??(ω)=ω(N?M)+r=1∑N?αr??r=1∑M?βr?
即幅頻回應等于所有的零點到 ω \omega ω的模長的乘積除以所有的極點到 ω \omega ω的模長的乘積;相頻回應等于所有的零點與 ω \omega ω形成的夾角之和減去所有的極點與 ω \omega ω形成的夾角之和,
從均值濾波看濾波器及系統函式
時域運算式
y ( n ) = 1 N ∑ k = 0 N x ( n ? k ) = 1 N [ x ( n ) + . . . + x ( n ? N + 1 ) ] y(n)=\frac{1}{N}\sum_{k=0}^Nx(n-k)\\=\frac{1}{N}\left[x(n)+...+x(n-N+1)\right] y(n)=N1?k=0∑N?x(n?k)=N1?[x(n)+...+x(n?N+1)]
x ( n ) x(n) x(n)是當前測量值, x ( n ? N + 1 ) x(n-N+1) x(n?N+1)是前第N-1次測量值,
頻率運算式
對均值濾波的時域運算式進行Z變換得:
H
(
z
)
=
1
N
∑
k
=
0
N
z
?
N
=
1
N
1
?
z
?
N
1
?
z
?
1
H(z)=\frac{1}{N}\sum_{k=0}^{N}z^{-N}\\=\frac{1}{N}\frac{1-z^{-N}}{1-z^{-1}}
H(z)=N1?k=0∑N?z?N=N1?1?z?11?z?N?
令
z
=
e
j
ω
z=e^{j\omega}
z=ejω得:
H
(
e
j
ω
)
=
1
N
1
?
e
?
j
ω
N
1
?
e
?
j
ω
=
1
N
e
?
j
ω
(
N
?
1
)
2
sin
?
(
ω
N
)
/
2
sin
?
(
ω
)
/
2
H(e^{j\omega})=\frac{1}{N}\frac{1-e^{-j{\omega}N}}{1-e^{-j{\omega}}}\\=\frac{1}{N}e^{-j\frac{\omega(N-1)}{2}}\frac{\sin({\omega}N)/2}{\sin(\omega)/2}
H(ejω)=N1?1?e?jω1?e?jωN?=N1?e?j2ω(N?1)?sin(ω)/2sin(ωN)/2?
零極點分布圖
H ( z ) = 1 N 1 ? z ? N 1 ? z ? 1 = 1 N z N ? 1 z N ? z N ? 1 H(z)=\frac{1}{N}\frac{1-z^{-N}}{1-z^{-1}}\\=\frac{1}{N}\frac{z^{N}-1}{z^{N}-z^{N-1}} H(z)=N1?1?z?11?z?N?=N1?zN?zN?1zN?1?
matlab代碼:
NN = [3 6 12];
for i = 1:3
N = NN(i);
B = ones(1, N);
B = [B -1];
A = zeros(1, N-1);
A = [1 -1 A];
subplot(3, 1, i);
zplane(B, A);
xlabel(['N=', num2str(N)]);
end

頻率回應分析
python代碼分析均值濾波頻率回應,采樣率為200Hz,階數為7,
from scipy.optimize import newton
from scipy.signal import freqz, dimpulse, dstep
from math import sin, cos, sqrt, pi
import numpy as np
import matplotlib.pyplot as plt
import sys
#函式用于計算移動平均濾波器的截止頻率
def get_filter_cutoff(N, **kwargs):
func = lambda w: sin(N*w/2) - N/sqrt(2) * sin(w/2)
deriv = lambda w: cos(N*w/2) * N/2 - N/sqrt(2) * cos(w/2) / 2
omega_0 = pi/N # Starting condition: halfway the first period of sin
return newton(func, omega_0, deriv, **kwargs)
#設定采樣率
sample_rate = 200 #Hz
N = 7
# 計算截止頻率
w_c = get_filter_cutoff(N)
cutoff_freq = w_c * sample_rate / (2 * pi)
# 濾波器引數
b = np.ones(N)
a = np.array([N] + [0]*(N-1))
#頻率回應
w, h = freqz(b, a, worN=4096)
#轉為頻率
w *= sample_rate / (2 * pi)
#繪制波特圖
plt.subplot(2, 1, 1)
plt.suptitle("Bode")
#轉換為分貝
plt.plot(w, 20 * np.log10(abs(h)))
plt.ylabel('Magnitude [dB]')
plt.xlim(0, sample_rate / 2)
plt.ylim(-60, 10)
plt.axvline(cutoff_freq, color='red')
plt.axhline(-3.01, linewidth=0.8, color='black', linestyle=':')
# 相頻回應
plt.subplot(2, 1, 2)
plt.plot(w, 180 * np.angle(h) / pi)
plt.xlabel('Frequency [Hz]')
plt.ylabel('Phase [°]')
plt.xlim(0, sample_rate / 2)
plt.ylim(-180, 90)
plt.yticks([-180, -135, -90, -45, 0, 45, 90])
plt.axvline(cutoff_freq, color='red')
plt.show()

改變濾波器階數,觀察幅頻回應(N=3,7,18):

可以看出,均值濾波器:
- 是一個低通濾波器
- 隨著濾波器得階數(長度)增加(求和的項越來越多),截止頻率變小,通帶變窄,濾波器的回應變慢,延遲變大,
- 是FIR濾波器
FIR和IIR濾波器的比較
穩定和線性相位特性是FIR濾波器最突出的優點,但階數較高;IIR濾波器相位非線性,但階數低,
性能上
- IIR濾波器系統函式的極點可位于單位圓內的任何地方,因此零點和極點相結合,可用較低的階數獲得較高的選擇性,所用的存盤單元少,計算量小,經濟效應高,但是這個高效率是以相位的非線性為代價的,
- 相反,FIR濾波器可以得到嚴格的線性相位,然而由于FIR濾波器系統函式的極點固定在原點,因而只能用較高的階數達到較高的選擇性,
- 對于同樣的濾波器幅頻特性指標,FIR濾波器所要求的階數一般比IIR濾波器高5~10倍,信號延時也比較大,
結構上
- IIR濾波器必須采用遞回結構,極點位置必須在單位圓內,否則系統不穩定,另外,這種結構中由于運算程序中對序列的舍入處理,會引起寄生振蕩,
- FIR濾波器主要采用非遞回結構,有限精度運算中不存在穩定性問題,運算誤差引起的輸出信號噪聲功率也較小,此外,FIR濾波器可以采用FFT演算法實作,速度大大提高,
設計工具上
- IIR濾波器可以借助成熟的模擬濾波器設計理論及成果,
- FIR濾波器計算帶通和阻帶衰減等仍無顯示運算式,邊界頻率也不易精確控制,
效果上
- IIR濾波器雖然設計簡單,但主要是用于設計具有片段常數特性的選頻型濾波器,比如低通、高通、帶通及帶阻等,往往脫離不了幾種典型模擬濾波器的頻響特性的約束,
- 而FIR濾波器要靈活得多,易于適應某些特殊的應用,比如只想濾掉特定的不連續的頻率,
從上面可以看出,IIR和FIR濾波器各有所長,在實際運用當中應全面考慮再加以選擇,比如,對于相位要求不高的場合,如語音通訊等可以選用IIR濾波器;對于圖象信號處理等相位要求高的場合應采用FIR濾波器,
FIR濾波器的設計
如何快速設計一個FIR濾波器(一)
如何快速設計一個FIR濾波器(二)
基本設計思路
FIR濾波器的設計方法和IIR濾波器的設計方法有很大差別,FIR濾波器的設計任務是選擇有限長度的 h ( n ) h(n) h(n),使頻率回應函式 H ( e j ω ) H(e^{j\omega}) H(ejω)滿足技術指標要求,
下面介紹幾種設計FIR濾波器的方法:窗函式法、頻率采樣法、回應最優法(切比雪夫等波紋逼近法、最小二乘法),
窗函式法
窗函式的作用
加窗原理和頻譜泄露 深入理解
FFT變換只能對有限長度的時域資料進行變換,因此,需要對時域信號進行信號截斷,即使是周期信號,如果截斷的時間長度不是周期的整數倍(周期截斷),那么,截取后的信號將會存在泄漏,為了將這個泄漏誤差減少到最小程度(注意我說是的減少,而不是消除),我們需要使用加權函式,也叫窗函式,加窗主要是為了使時域信號似乎更好地滿足FFT處理的周期性要求,減少泄漏,
如下圖所示:
- 若周期截斷,則FFT頻譜為單一譜線;
- 若為非周期截斷,則頻譜出現拖尾,如圖中部所示,可以看出泄漏很嚴重,

為了減少泄漏,給信號施加一個窗函式(如圖中上部紅色曲線所示),原始截斷后的信號與這個窗函式相乘之后得到的信號為上面右側的信號,可以看出,此時,信號的起始時刻和結束時刻幅值都為0,也就是說在這個時間長度內,信號為周期信號,但是只有一個周期,對這個信號做FFT分析,得到的頻譜如下部右側所示,相比較之前未加窗的頻譜,可以看出,泄漏已明顯改善,但并沒有完全消除,因此,窗函式只能減少泄漏,不能消除泄漏,
基于窗函式法設計FIR濾波器,其基本步驟如下:
- 確定頻域的回應函式 H d ( e j ω ) H_d(e^{j\omega}) Hd?(ejω),低通、高通、帶通或者其它;
- 確定頻域的回應函式 H d ( e j ω ) H_d(e^{j\omega}) Hd?(ejω)的傅里葉逆變換,找到連續脈沖回應函式 h d ( t ) h_d(t) hd?(t);
- 對連續脈沖回應函式 h d ( t ) h_d(t) hd?(t)按照一定的采樣頻率進行采樣,獲得離散信號 h d ( m ) h_d(m) hd?(m);
- 選擇合適的窗函式,對離散信號 h d ( m ) h_d(m) hd?(m)加窗,獲得有限長度的脈沖回應 h ( m ) h(m) h(m),
MATLAB中基于窗函式法設計濾波器:
% h = fir1(n, Wn, ftype, window);
% n表示濾波器階數;Wn表示歸一化后的濾波器截止頻率,可表示成[fl fh]的形式;ftype表示濾波器型別;window表示窗函式型別,
% 這是一個24階的帶通濾波器,歸一化截止頻率為[0.2 0.4]
h = fir1(24, [0.2 0.4]);
freqz(h, 1, 512);

頻率采樣法
設希望逼近的濾波器的頻響函式用
H
d
(
e
j
ω
)
H_d(e^{j\omega})
Hd?(ejω)表示,對
H
d
(
e
j
ω
)
H_d(e^{j\omega})
Hd?(ejω)在
ω
=
0
\omega=0
ω=0到
2
π
2\pi
2π之間等間隔采樣N點,得到
H
d
(
k
)
H_d(k)
Hd?(k):
H
d
(
k
)
=
H
d
(
e
j
ω
)
∣
ω
=
2
π
N
k
k
=
0
,
1
,
.
.
.
,
N
?
1
H_d(k)=H_d(e^{j\omega})|_{\omega=\frac{2\pi}{N}k}{\quad}k=0,1,...,N-1
Hd?(k)=Hd?(ejω)∣ω=N2π?k?k=0,1,...,N?1
再對
H
d
(
k
)
H_d(k)
Hd?(k)進行N點IDFT,得到
h
(
n
)
h(n)
h(n):
h
(
n
)
=
I
D
F
T
[
H
d
(
k
)
]
=
1
N
∑
k
=
0
N
?
1
H
d
(
k
)
W
N
?
k
n
n
=
0
,
1
,
2
,
.
.
.
,
N
?
1
h(n)=IDFT[H_d(k)]=\frac{1}{N}\sum_{k=0}^{N-1}H_d(k)W_N^{-kn}{\quad}n=0,1,2,...,N-1
h(n)=IDFT[Hd?(k)]=N1?k=0∑N?1?Hd?(k)WN?kn?n=0,1,2,...,N?1
將
h
(
n
)
h(n)
h(n)作為所設計的FIR濾波器的單位沖激回應,其系統函式
H
(
z
)
H(z)
H(z)為:
H
(
z
)
=
∑
n
=
0
N
?
1
h
(
n
)
z
?
n
H(z)=\sum_{n=0}^{N-1}h(n)z^{-n}
H(z)=n=0∑N?1?h(n)z?n
或者根據頻率域采樣理論,利用頻率域采樣值
H
d
(
k
)
H_d(k)
Hd?(k)恢復原信號的Z變換,得到
H
(
z
)
H(z)
H(z)的內插表示形式為:
H
(
z
)
=
1
?
z
?
1
N
∑
k
=
0
N
?
1
H
d
(
k
)
1
?
W
N
?
k
z
?
1
H(z)=\frac{1-z^{-1}}{N}\sum_{k=0}^{N-1}\frac{H_d(k)}{1-W_N^{-k}z^{-1}}
H(z)=N1?z?1?k=0∑N?1?1?WN?k?z?1Hd?(k)?
上述兩種方法如下:

MATLAB中使用頻率采樣法設計FIR濾波器:
% h = fir2(n, f, m);
% n表示階數;f表示分立點頻率矢量;m表示分立點對應的幅值回應矢量
f = [0 0.7 0.7 1];
m = [10 1 0 0];
h = fir2(24, f, m);
freqz(h, 1);

回應最優法
主要思路就是找到一組脈沖回應,讓它的頻域回應 H ( e j ω ) H(e^{j\omega}) H(ejω)與期望的濾波器的頻域回應 H d ( e j ω ) H_d(e^{j\omega}) Hd?(ejω)盡可能一致,主要通過兩種方法來實作,一個是最小二乘法、另一個是(切比雪夫)等波紋逼近法,
最小二乘法
與頻率采樣法近似,期望的頻率回應用一組分立的點
(
f
i
,
A
(
f
i
)
)
i
=
1
,
2
,
.
.
.
,
N
(f_i,\ A(f_i)){\quad}i=1,2,...,N
(fi?, A(fi?))i=1,2,...,N來表示,優化目標如下:
∑
i
=
1
N
W
i
[
H
(
f
i
)
?
H
d
(
f
i
)
]
2
→
m
i
n
\sum_{i=1}^{N}W_i\left[H(f_i)-H_d(f_i)\right]^2{\rightarrow}min
i=1∑N?Wi?[H(fi?)?Hd?(fi?)]2→min
其中
W
i
W_i
Wi?為不同頻率下的權重,
MATLAB指令為:
% h = firls(n, f, m);
% n為濾波器階數,f表示分立點的矢量,m表示分立點對應的幅值回應矢量,
% 下面以低通濾波器為例:
f = [0 0.3 0.3 1];
m = [10 1 0 0];
h = firls(24, f, m);
freqz(h, 1)

切比雪夫等波紋最佳逼近法
書上的:
用
H
d
(
ω
)
H_d(\omega)
Hd?(ω)表示希望逼近的幅度特性函式,要求設計線性相位FIR數字濾波器時,
H
d
(
ω
)
H_d(\omega)
Hd?(ω)必須滿足線性相位約束條件,用
H
g
(
ω
)
H_g(\omega)
Hg?(ω)表示實際設計的濾波器幅度特性函式,定義加權誤差函式
E
(
ω
)
E(\omega)
E(ω)為:
E
(
ω
)
=
W
(
ω
)
[
H
d
(
ω
)
?
H
g
(
ω
)
]
m
a
x
E
(
ω
)
→
m
i
n
E(\omega)=W(\omega)[H_d(\omega)-H_g(\omega)]\\ max\ E(\omega){\rightarrow}min
E(ω)=W(ω)[Hd?(ω)?Hg?(ω)]max E(ω)→min
另外:
與最小二乘法不同(方差最小),切比雪夫法采用的方案是最大誤差最小,即:
max
?
∣
W
i
(
H
(
f
i
)
?
H
d
(
f
i
)
)
∣
→
min
?
\max{\left|W_i(H(f_i)-H_d(f_i))\right|}{\rightarrow}\min
max∣Wi?(H(fi?)?Hd?(fi?))∣→min
MATLAB指令:
h = firpm(n, f, m);
% 對于切比雪夫法,頻率回應不能直接從1降到0
% f = [0 0.1 0.1 1]這樣是不行的
f = [0 0.1 0.101 1];
m = [100 1 0 0];
h = firpm(24, f, m);
freqz(h, 1)

FIR數字濾波器設計總結
對于回應最優法更像是擬合,下圖是窗函式法和頻率采樣法的流程圖:

IIR濾波器設計
如何快速設計一個IIR濾波器
基本設計思路
利用模擬濾波器成熟的理論及其設計方法來設計IIR數字低通數字濾波器是常用的方法,設計程序是:
- 按照數字濾波器技術指標要求設計一個過渡模擬低通濾波器 H a ( s ) H_a(s) Ha?(s);
- 再按照一定的轉換關系將 H a ( s ) H_a(s) Ha?(s)轉換成數字低通濾波器的系統函式 H ( z ) H(z) H(z);
所以可見,設計的關鍵就是找到某種轉換關系,將s平面上的 H a ( s ) H_a(s) Ha?(s)轉換到z平面上的 H ( z ) H(z) H(z),這種轉換關系需要滿足一定的條件:
- 因果穩定的模擬濾波器轉換成數字濾波器后,仍是因果穩定的;
- s平面的左半平面映射到z平面的單位圓內部;
- 數字濾波器的頻率回應應模仿模擬濾波器的頻響特性,s平面的虛軸映射為z平面的單位圓,相應的頻率之間呈線性關系,
基于以上要求,目前常用的轉換方法有脈沖回應不變法、雙線性變換法,脈沖回應不變法是線性變換,而雙線性變換法是非線性變換,
注:我們一般是設計相應的低通濾波器,然后根據頻率變換公式,得到相應的高通、帶通、帶阻濾波器,
模擬濾波器的設計
模擬濾波器的技術指標給定后,需要設計一個系統函式
H
a
(
s
)
H_a(s)
Ha?(s),希望其幅度平方函式滿足給定的指標,一般濾波器的單位沖激回應為實函式,因此:
∣
H
a
(
j
Ω
)
∣
2
=
H
a
(
s
)
H
a
(
?
s
)
∣
s
=
j
Ω
=
H
a
(
j
Ω
)
H
a
?
(
j
Ω
)
|H_a(j\Omega)|^2=H_a(s)H_a(-s)|_{s=j\Omega}=H_a(j\Omega)H_a^*(j\Omega)
∣Ha?(jΩ)∣2=Ha?(s)Ha?(?s)∣s=jΩ?=Ha?(jΩ)Ha??(jΩ)
如果能由
α
p
\alpha_p
αp?,
Ω
p
\Omega_p
Ωp?,
α
s
\alpha_s
αs?,
Ω
s
\Omega_s
Ωs?求出
∣
H
a
(
j
Ω
)
∣
2
|H_a(j\Omega)|^2
∣Ha?(jΩ)∣2,那么就可以求出
H
a
(
s
)
H
a
(
?
s
)
H_a(s)H_a(-s)
Ha?(s)Ha?(?s),由此可以求出需要的
H
a
(
s
)
H_a(s)
Ha?(s),其極點必須落在s平面的左半平面;相應地,
H
a
(
?
s
)
H_a(-s)
Ha?(?s)的極點就應落在s平面的右半平面,這就是模擬低通濾波器的逼近方法,
巴特沃斯濾波器
幅度平方函式
∣
H
a
(
j
Ω
)
∣
2
|H_a(j\Omega)|^2
∣Ha?(jΩ)∣2為:
∣
H
a
(
j
Ω
)
∣
2
=
1
1
+
(
Ω
Ω
c
)
2
N
=
1
1
+
(
s
j
Ω
c
)
2
N
∣
s
=
j
Ω
|H_a(j\Omega)|^2=\frac{1}{1+(\frac{\Omega}{\Omega_c})^{2N}}\\ =\frac{1}{1+(\frac{s}{j\Omega_c})^{2N}}|_{s=j\Omega}
∣Ha?(jΩ)∣2=1+(Ωc?Ω?)2N1?=1+(jΩc?s?)2N1?∣s=jΩ?
N為濾波器階數,該濾波器是一個低通濾波器,原因如下:
- Ω = 0 \Omega=0 Ω=0時, ∣ H a ( j Ω ) ∣ 2 = 1 |H_a(j\Omega)|^2=1 ∣Ha?(jΩ)∣2=1;
- Ω = Ω c \Omega=\Omega_c Ω=Ωc?時, ∣ H a ( j Ω ) ∣ 2 = 1 2 |H_a(j\Omega)|^2=\frac{1}{\sqrt{2}} ∣Ha?(jΩ)∣2=2 ?1?,即3dB截至頻率;
- Ω > > Ω c \Omega>>\Omega_c Ω>>Ωc?時, ∣ H a ( j Ω ) ∣ 2 = 0 |H_a(j\Omega)|^2=0 ∣Ha?(jΩ)∣2=0;
該系統函式有2N個極點,極點為:
s
k
=
(
?
1
)
1
2
N
(
j
Ω
c
)
=
Ω
c
e
j
π
(
1
2
+
2
k
+
1
2
N
)
,
k
=
0
,
1
,
.
.
.
,
2
N
?
1
s_k=(-1)^{\frac{1}{2N}}(j\Omega_c)={\Omega_c}e^{j{\pi}(\frac{1}{2}+\frac{2k+1}{2N})},k=0,1,...,2N-1
sk?=(?1)2N1?(jΩc?)=Ωc?ejπ(21?+2N2k+1?),k=0,1,...,2N?1
這2N個極點等間隔分布在半徑為
Ω
c
\Omega_c
Ωc?的圓上(此圓稱為巴特沃斯圓),間隔是
π
/
N
\pi/N
π/Nrad,
N = 3;
k = 0:(2*N-1);
k = k';
Z = [];
P = exp(1i*pi*(0.5 + (2*k+1) ./ (2*N)));
% 零極點
zplane(Z, P);
title('零極點分布圖');

取左半平面的N個極點構成系統函式
H
a
(
s
)
H_a(s)
Ha?(s),并用3dB截止頻率
Ω
c
\Omega_c
Ωc?歸一化得:
G
a
(
s
Ω
c
)
=
1
∏
k
=
0
N
?
1
(
s
Ω
c
?
s
k
Ω
c
)
=
1
∏
k
=
0
N
?
1
(
p
?
p
k
)
(
1
)
G_a(\frac{s}{\Omega_c})=\frac{1}{\prod_{k=0}^{N-1}(\frac{s}{\Omega_c}-\frac{s_k}{\Omega_c})}\\ ={\color{blue}\frac{1}{\prod_{k=0}^{N-1}(p-p_k)}}{\quad}(1)
Ga?(Ωc?s?)=∏k=0N?1?(Ωc?s??Ωc?sk??)1?=∏k=0N?1?(p?pk?)1?(1)
式中
p
k
=
s
k
/
Ω
c
=
e
j
π
(
1
2
+
2
K
+
1
2
N
)
,
k
=
0
,
1
,
.
.
.
,
N
?
1
(
2
)
p_k=s_k/\Omega_c=e^{j{\pi}(\frac{1}{2}+\frac{2K+1}{2N})},k=0,1,...,N-1{\qquad}(2)
pk?=sk?/Ωc?=ejπ(21?+2N2K+1?),k=0,1,...,N?1(2)
為歸一化極點,
這樣,只要根據技術指標求出階數N,按照式(2)求出N個極點,再按照式(1)得到歸一化的系統函式 G a ( s ) G_a(s) Ga?(s),如果知道截止頻率 Ω c \Omega_c Ωc?,即可知道系統函式 H a ( s ) H_a(s) Ha?(s),
注:設計出滿足要求的低通濾波器后,通過頻率變換公式即可得到高通濾波器、帶通濾波器、帶阻濾波器等,
用MATLAB設計巴特沃斯濾波器(成電DSP書上P160有詳細說明):
[Z, P, K] = buttap(N);
[N, wc] = buttord(wp, ws, Rp, As);
[N, wc] = buttord(wp, ws, Rp, As, 's');
[B, A] = butter(N, wc, 'ftype');
[B, A] = butter(N, wc, 'ftype', 's')
切比雪夫I型濾波器
其幅度平方函式為:
∣
H
a
(
j
Ω
)
∣
2
=
1
1
+
ε
2
C
N
2
(
Ω
Ω
p
)
|H_a(j\Omega)|^2=\frac{1}{1+\varepsilon^2C_N^2(\frac{\Omega}{\Omega_p})}
∣Ha?(jΩ)∣2=1+ε2CN2?(Ωp?Ω?)1?
式中,
ε
\varepsilon
ε為小于1的正數,
橢圓濾波器
會用MATLAB提供的函式即可,
高通、帶通和帶阻濾波器的設計
通過某些方法得到了低通濾波器后,再通過頻率變換就可以得到高通、帶通和帶阻濾波器,
歸一化低通濾波器
令截止頻率
ω
c
=
1
\omega_c=1
ωc?=1:
H
(
s
)
=
ω
c
s
+
ω
c
=
1
s
+
1
H(s)=\frac{\omega_c}{s+\omega_c}\\=\frac{1}{s+1}
H(s)=s+ωc?ωc??=s+11?
低通到高通的頻率轉換
我們的目標是:
- 當 ω → 0 \omega{\rightarrow}0 ω→0時, ∣ H ′ ( ω ) ∣ → 0 |H'(\omega)|{\rightarrow}0 ∣H′(ω)∣→0;
- 當 ω → ω c \omega{\rightarrow}\omega_c ω→ωc?時, ∣ H ′ ( ω ) ∣ → 1 2 |H'(\omega)|{\rightarrow}\frac{1}{\sqrt{2}} ∣H′(ω)∣→2 ?1?;
- 當 ω → ∞ \omega{\rightarrow}\infin ω→∞時, ∣ H ′ ( ω ) ∣ → 1 |H'(\omega)|{\rightarrow}1 ∣H′(ω)∣→1;
所以,令
s
→
ω
c
s
s{\rightarrow}\frac{\omega_c}{s}
s→sωc??即可,則傳遞函式變為:
H
′
(
s
)
=
s
s
+
ω
c
H'(s)=\frac{s}{s+\omega_c}
H′(s)=s+ωc?s?
低通到帶通的頻率轉換
我們的目標是:
- 當 ω → 0 \omega{\rightarrow}0 ω→0時, ∣ H ′ ( ω ) ∣ → 0 |H'(\omega)|{\rightarrow}0 ∣H′(ω)∣→0;
- 當 ω → ω l \omega{\rightarrow}\omega_l ω→ωl?時, ∣ H ′ ( ω ) ∣ → 1 2 |H'(\omega)|{\rightarrow}\frac{1}{\sqrt{2}} ∣H′(ω)∣→2 ?1?;
- 當 ω ? ( ω l , ω h ) \omega{\subset}(\omega_l,\omega_h) ω?(ωl?,ωh?)時, ∣ H ′ ( ω ) ∣ ≥ 1 2 |H'(\omega)|{\ge}\frac{1}{\sqrt{2}} ∣H′(ω)∣≥2 ?1?
- 當 ω → ω h \omega{\rightarrow}\omega_h ω→ωh?時, ∣ H ′ ( ω ) ∣ → 1 2 |H'(\omega)|{\rightarrow}\frac{1}{\sqrt{2}} ∣H′(ω)∣→2 ?1?;
- 當 ω → ∞ \omega{\rightarrow}\infin ω→∞時, ∣ H ′ ( ω ) ∣ → 0 |H'(\omega)|{\rightarrow}0 ∣H′(ω)∣→0;
所以,令
s
→
s
2
+
ω
0
2
B
s
s{\rightarrow}\frac{s^2+\omega_0^2}{Bs}
s→Bss2+ω02??即可,其中
ω
0
2
=
ω
l
ω
h
,
B
=
ω
h
?
ω
l
\omega_0^2=\omega_l\omega_h,\ B=\omega_h-\omega_l
ω02?=ωl?ωh?, B=ωh??ωl?,則傳遞函式變為:
H
′
(
s
)
=
B
s
s
2
+
B
s
+
ω
0
2
H'(s)=\frac{Bs}{s^2+Bs+\omega_0^2}
H′(s)=s2+Bs+ω02?Bs?
低通到帶阻的頻率轉換
我們的目標是:
- 當 ω → 0 \omega{\rightarrow}0 ω→0時, ∣ H ′ ( ω ) ∣ → 1 |H'(\omega)|{\rightarrow}1 ∣H′(ω)∣→1;
- 當 ω → ω l \omega{\rightarrow}\omega_l ω→ωl?時, ∣ H ′ ( ω ) ∣ → 1 2 |H'(\omega)|{\rightarrow}\frac{1}{\sqrt{2}} ∣H′(ω)∣→2 ?1?;
- 當 ω ? ( ω l , ω h ) \omega{\subset}(\omega_l,\omega_h) ω?(ωl?,ωh?)時, ∣ H ′ ( ω ) ∣ ≤ 1 2 |H'(\omega)|{\le}\frac{1}{\sqrt{2}} ∣H′(ω)∣≤2 ?1?
- 當 ω → ω h \omega{\rightarrow}\omega_h ω→ωh?時, ∣ H ′ ( ω ) ∣ → 1 2 |H'(\omega)|{\rightarrow}\frac{1}{\sqrt{2}} ∣H′(ω)∣→2 ?1?;
- 當 ω → ∞ \omega{\rightarrow}\infin ω→∞時, ∣ H ′ ( ω ) ∣ → 1 |H'(\omega)|{\rightarrow}1 ∣H′(ω)∣→1;
所以,令
s
→
B
s
s
2
+
ω
0
2
s{\rightarrow}\frac{Bs}{s^2+\omega_0^2}
s→s2+ω02?Bs?即可,其中
ω
0
2
=
ω
l
ω
h
,
B
=
ω
h
?
ω
l
\omega_0^2=\omega_l\omega_h,\ B=\omega_h-\omega_l
ω02?=ωl?ωh?, B=ωh??ωl?,則傳遞函式變為:
H
′
(
s
)
=
s
2
+
ω
0
2
s
2
+
B
s
+
ω
0
2
H'(s)=\frac{s^2+\omega_0^2}{s^2+Bs+\omega_0^2}
H′(s)=s2+Bs+ω02?s2+ω02??
從模擬到數字的轉換方法
脈沖回應不變法
基本思路
設模擬濾波器的系統函式為 H a ( s ) H_a(s) Ha?(s),單位沖激回應是 h a ( t ) h_a(t) ha?(t), H a ( s ) = L T [ h a ( t ) ] H_a(s)=LT[h_a(t)] Ha?(s)=LT[ha?(t)],對 h a ( t ) h_a(t) ha?(t)進行等間隔采樣,采樣間隔為T,得到 h a ( n T ) h_a(nT) ha?(nT),將 h ( n ) = h a ( n T ) h(n)=h_a(nT) h(n)=ha?(nT)作為數字濾波器的單位沖激回應,那么數字濾波器的系統函式 H ( z ) H(z) H(z)便是 h ( n ) h(n) h(n)的Z變換,
因此脈沖回應不變法是一種時域逼近方法,它使 h ( n ) h(n) h(n)在采樣點上等于 h a ( t ) h_a(t) ha?(t),由于模擬濾波器的設計結果是 H a ( s ) H_a(s) Ha?(s),所以下面基于脈沖回應不變法的思想,匯出直接從 H a ( s ) H_a(s) Ha?(s)到 H ( z ) H(z) H(z)的轉換公式,
推導程序
設模擬濾波器
H
a
(
s
)
H_a(s)
Ha?(s)只有單階極點,且分母多項式的階次高于分子多項式的階次,將
H
a
(
s
)
H_a(s)
Ha?(s)用部分分式表示并進行推導:
H
a
(
s
)
=
∑
i
=
1
N
A
i
s
?
s
i
h
a
(
t
)
=
L
[
H
a
(
s
)
]
=
∑
i
=
1
N
A
i
e
s
i
t
u
(
t
)
h
(
n
)
=
h
a
(
n
T
)
=
∑
i
=
1
N
A
i
e
s
i
n
T
u
(
n
T
)
H
(
z
)
=
∑
i
=
1
N
A
i
T
1
?
e
s
i
T
z
?
1
H_a(s)=\sum_{i=1}^{N}\frac{A_i}{s-s_i}\\ h_a(t)=\mathscr{L}[H_a(s)]=\sum_{i=1}^{N}A_ie^{s_it}u(t)\\ h(n)=h_a(nT)=\sum_{i=1}^{N}A_ie^{s_inT}u(nT)\\ H(z)=\sum_{i=1}^{N}\frac{A_iT}{1-e^{s_iT}z^{-1}}\\
Ha?(s)=i=1∑N?s?si?Ai??ha?(t)=L[Ha?(s)]=i=1∑N?Ai?esi?tu(t)h(n)=ha?(nT)=i=1∑N?Ai?esi?nTu(nT)H(z)=i=1∑N?1?esi?Tz?1Ai?T?
即:
H
a
(
s
)
=
∑
i
=
1
N
A
i
s
?
s
i
(
1
)
H
(
z
)
=
∑
i
=
1
N
A
i
T
1
?
e
s
i
T
z
?
1
(
2
)
H_a(s)=\sum_{i=1}^{N}\frac{A_i}{s-s_i}{\quad}(1)\\ H(z)=\sum_{i=1}^{N}\frac{A_iT}{1-e^{s_iT}z^{-1}}{\quad}(2)
Ha?(s)=i=1∑N?s?si?Ai??(1)H(z)=i=1∑N?1?esi?Tz?1Ai?T?(2)
我們已經得到直接從
H
a
(
s
)
H_a(s)
Ha?(s)到
H
(
z
)
H(z)
H(z)的轉換公式了,下面分析從模擬濾波器
H
a
(
s
)
H_a(s)
Ha?(s)轉換到數字濾波器
H
(
z
)
H(z)
H(z),s平面和z平面之間的映射關系,從而找到這種轉換方式(即式(1)到式(2))的優缺點,這里以理想采樣信號
h
^
a
(
t
)
\hat{h}_a(t)
h^a?(t)作為橋梁,推導其映射關系,即推導
h
^
a
(
t
)
\hat{h}_a(t)
h^a?(t)與
H
(
z
)
H(z)
H(z)之間的映射關系,
∵
h
^
a
(
t
)
=
∑
n
=
?
∞
+
∞
h
a
(
t
)
δ
(
t
?
n
T
)
∴
H
^
a
(
s
)
=
L
[
h
^
a
(
t
)
]
=
∫
?
∞
+
∞
h
^
a
(
t
)
e
?
s
t
d
t
=
∫
?
∞
+
∞
[
∑
n
=
?
∞
+
∞
h
a
(
t
)
δ
(
t
?
n
T
)
]
e
?
s
t
d
t
=
∫
?
∞
+
∞
[
∑
n
=
?
∞
+
∞
h
a
(
n
T
)
δ
(
t
?
n
T
)
]
e
?
s
n
T
d
t
=
∑
n
=
?
∞
+
∞
h
a
(
n
T
)
e
?
s
n
T
∫
?
∞
+
∞
δ
(
t
?
n
T
)
d
t
=
∑
n
=
?
∞
+
∞
h
a
(
n
T
)
e
?
s
n
T
{\because}{\quad}\hat{h}_a(t)=\sum_{n=-\infin}^{+\infin}h_a(t)\delta(t-nT)\\ {\therefore}{\quad}\hat{H}_a(s)=\mathscr{L}[\hat{h}_a(t)]=\int_{-\infin}^{+\infin}\hat{h}_a(t)e^{-st}d_t\\ =\int_{-\infin}^{+\infin}[\sum_{n=-\infin}^{+\infin}h_a(t)\delta(t-nT)]e^{-st}d_t\\ =\int_{-\infin}^{+\infin}[\sum_{n=-\infin}^{+\infin}h_a({\color{blue}nT})\delta(t-nT)]e^{-s{\color{blue}nT}}d_t\\ =\sum_{n=-\infin}^{+\infin}h_a(nT)e^{-snT}{\color{red}\int_{-\infin}^{+\infin}\delta(t-nT)d_t}\\ =\sum_{n=-\infin}^{+\infin}h_a(nT)e^{-snT}
∵h^a?(t)=n=?∞∑+∞?ha?(t)δ(t?nT)∴H^a?(s)=L[h^a?(t)]=∫?∞+∞?h^a?(t)e?stdt?=∫?∞+∞?[n=?∞∑+∞?ha?(t)δ(t?nT)]e?stdt?=∫?∞+∞?[n=?∞∑+∞?ha?(nT)δ(t?nT)]e?snTdt?=n=?∞∑+∞?ha?(nT)e?snT∫?∞+∞?δ(t?nT)dt?=n=?∞∑+∞?ha?(nT)e?snT
式中,
h
a
(
n
T
)
h_a(nT)
ha?(nT)是
h
a
(
t
)
h_a(t)
ha?(t)在采樣點
t
=
n
T
t=nT
t=nT時的幅度值,它與序列
h
(
n
)
h(n)
h(n)的幅度值相等,即
h
(
n
)
=
h
a
(
n
T
)
h(n)=h_a(nT)
h(n)=ha?(nT),因此得到:
H
^
a
(
s
)
=
∑
n
=
?
∞
+
∞
h
(
n
)
e
?
s
n
T
=
∑
n
=
?
∞
+
∞
h
(
n
)
z
?
n
∣
z
=
e
s
T
=
H
(
z
)
∣
z
=
e
s
T
\hat{H}_a(s)=\sum_{n=-\infin}^{+\infin}h(n)e^{-snT}\\ =\sum_{n=-\infin}^{+\infin}h(n)z^{-n}|_{z=e^{sT}}\\ =H(z)|_{z=e^{sT}}
H^a?(s)=n=?∞∑+∞?h(n)e?snT=n=?∞∑+∞?h(n)z?n∣z=esT?=H(z)∣z=esT?
上式表明理想采樣信號
h
^
a
(
t
)
\hat{h}_a(t)
h^a?(t)的拉氏變換與相應采樣序列
h
(
n
)
h(n)
h(n)的Z變換之間的映射關系為:
z
=
e
s
T
l
e
t
s
=
δ
+
j
Ω
,
z
=
r
e
j
ω
{
r
=
e
δ
T
ω
=
Ω
T
z=e^{sT}\\ let{\quad}s={\delta}+j\Omega,\ z=re^{j\omega}\\ \begin{cases} r=e^{{\delta}T} \\ \omega={\Omega}T \end{cases}
z=esTlets=δ+jΩ, z=rejω{r=eδTω=ΩT?
因此:
{
δ
=
0
,
r
=
1
δ
<
0
,
r
<
1
δ
>
0
,
r
>
1
\begin{cases} \delta=0,\ r=1\\ \delta<0,\ r<1\\ \delta>0,\ r>1 \end{cases}
??????δ=0, r=1δ<0, r<1δ>0, r>1?
上式說明,s平面的虛軸映射為z平面的單位圓;s平面的左半平面映射為z平面單位圓內;s平面的右半平面映射為z平面單位圓外,這說明,如果
H
a
(
s
)
H_a(s)
Ha?(s)因果穩定,轉換后得到的
H
(
z
)
H(z)
H(z)仍是穩定的,
待補:討論數字濾波器的頻響特性與模擬濾波器的頻響特性之間的關系(P179),
雙線性變換法
脈沖回應不變法的主要缺點是會產生頻譜混疊現象,使數字濾波器的頻響偏離模擬濾波器的頻響特性,產生的原因是模擬低通濾波器不是帶限于折疊頻率 π / T {\pi}/T π/T,在離散化(采樣)后產生了頻譜混疊,再通過映射關系 z = e s T z=e^{sT} z=esT,使數字濾波器在 ω = π \omega=\pi ω=π附近形成了頻譜混疊,
基本思路
為了克服這一缺點,可以采用非線性頻率壓縮方法,將整個模擬頻率軸壓縮到
±
π
/
T
\pm\pi/T
±π/T之間,再用
z
=
e
s
T
z=e^{sT}
z=esT轉換到z平面上,這里用正切變換實作頻率壓縮:
Ω
=
2
T
tan
?
(
1
2
Ω
1
T
)
\Omega=\frac{2}{T}\tan{\left(\frac{1}{2}{\Omega_1}T\right)}
Ω=T2?tan(21?Ω1?T)
式中,T仍是采樣間隔,當
Ω
1
\Omega_1
Ω1?從
?
π
/
T
-\pi/T
?π/T經過0變化到
?
π
/
T
-\pi/T
?π/T時,
Ω
\Omega
Ω則由
?
∞
-\infin
?∞經過0變換到
+
∞
+\infin
+∞,實作了s平面上整個虛軸完全壓縮到s1平面上虛軸的
±
π
/
T
\pm\pi/T
±π/T之間的轉換,
s
=
2
T
1
?
z
?
1
1
+
z
?
1
(
1
)
z
=
2
T
+
s
2
T
?
s
s=\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}{\quad}(1)\\ z=\frac{\frac{2}{T}+s}{\frac{2}{T}-s}
s=T2?1+z?11?z?1?(1)z=T2??sT2?+s?

轉換性能分析
先分析模擬頻率
Ω
\Omega
Ω和數字頻率
ω
\omega
ω之間的關系,令
s
=
j
Ω
,
z
=
e
j
ω
s=j\Omega,\ z=e^{j\omega}
s=jΩ, z=ejω,并帶入式(1)得到:
j
Ω
=
2
T
1
?
e
?
j
ω
1
+
e
?
j
ω
Ω
=
2
T
tan
?
(
1
2
ω
)
j\Omega=\frac{2}{T}\frac{1-e^{-j\omega}}{1+e^{-j\omega}}\\ \Omega=\frac{2}{T}\tan{(\frac{1}{2}\omega)}
jΩ=T2?1+e?jω1?e?jω?Ω=T2?tan(21?ω)
上式說明,s平面上的
Ω
\Omega
Ω與z平面的
ω
\omega
ω成非線性正切關系,

- 在 ω = 0 \omega=0 ω=0附近接近線性關系;
- 當 ω \omega ω增加時, Ω \Omega Ω加得愈來愈快;
- 當 ω \omega ω趨近于 π \pi π時, Ω \Omega Ω趨近于 ∞ \infin ∞,
正是這種非線性關系,消除了頻譜混疊現象,
可見,對于雙線性變換而言,在s域的頻率和z域對應的頻率不同,發生了一定的彎曲,也就意味著截止頻率在s域和在z域是不一樣的,現在我們需要找到這種關系(省略掉推導),在z域,假設截止頻率為
f
d
f_d
fd?,則:
f
a
=
f
s
π
tan
?
(
π
f
d
f
s
)
f_a=\frac{f_s}{\pi}\tan{\left(\frac{{\pi}f_d}{f_s}\right)}
fa?=πfs??tan(fs?πfd??)
也就是說對于雙線性變而言,模擬的截止頻率和數字的截止頻率是不同的,不同的原因是因為雙線性變換是近似變換,不是準確換算,
但是,當
f
s
>
>
f
d
f_s>>f_d
fs?>>fd?時,即采樣頻率遠大于截止頻率時,可以得到:
f
a
=
f
s
π
tan
?
(
π
f
d
f
s
)
≈
f
s
π
(
π
f
d
f
s
)
=
f
d
f_a=\frac{f_s}{\pi}\tan{\left(\frac{{\pi}f_d}{f_s}\right)}{\approx}\frac{f_s}{\pi}\left(\frac{{\pi}f_d}{f_s}\right)=f_d
fa?=πfs??tan(fs?πfd??)≈πfs??(fs?πfd??)=fd?
有了前面的說明,就可以設計IIR濾波器了,基本步驟如下:
- 選擇一個歸一化的模擬濾波器 H ( s ) H(s) H(s);
- 確定數字濾波器的截止頻率,也就是回應為-3dB(幅值衰減為 1 / 2 1/\sqrt{2} 1/2 ?)時對應的頻率;
- 利用公式 f a = f s π tan ? ( π f d f s ) f_a=\frac{f_s}{\pi}\tan{\left(\frac{{\pi}f_d}{f_s}\right)} fa?=πfs??tan(fs?πfd??)計算對應的模擬截止頻率;
- 選擇合適的變換,得到 H ′ ( s ) H'(s) H′(s);比如對于低通濾波器: s → s / ω c s{\rightarrow}s/\omega_c s→s/ωc?,其中 ω c = 2 π f a \omega_c=2{\pi}f_a ωc?=2πfa?,
- 在 H ′ ( s ) H'(s) H′(s)中,用 s → 2 f s ( z ? 1 z + 1 ) s{\rightarrow}2f_s(\frac{z-1}{z+1}) s→2fs?(z+1z?1?)進行替換,得到數字化的濾波器 H ( z ) H(z) H(z),
IIR濾波器設計實體
比如:現在假設要設計一個低通濾波器,截止頻率為 f d = 10 H z f_d=10Hz fd?=10Hz,采樣頻率為 f s = 50 H z f_s=50Hz fs?=50Hz,
- 選擇歸一化的濾波器 H ( s ) = 1 s + 1 H(s)=\frac{1}{s+1} H(s)=s+11?;
- 數字截止頻率為10Hz;
- 對應的模擬截止頻率為: f a = f s π tan ? ( π f d f s ) = 11.56 H z f_a=\frac{f_s}{\pi}\tan{\left(\frac{{\pi}f_d}{f_s}\right)}=11.56Hz fa?=πfs??tan(fs?πfd??)=11.56Hz;
- 將 H ( s ) H(s) H(s)中的s進行替換: s → s / ω c , ω c = 2 π f a s{\rightarrow}s/\omega_c,\ \omega_c=2{\pi}f_a s→s/ωc?, ωc?=2πfa?,得到: H ′ ( s ) = ω c s + ω c H'(s)=\frac{\omega_c}{s+\omega_c} H′(s)=s+ωc?ωc??,其中 ω c = 2 π ( 11.56 ) = 72.6 r a d / s \omega_c=2{\pi}(11.56)=72.6rad/s ωc?=2π(11.56)=72.6rad/s;
- 用 s = 2 f s ( z ? 1 z + 1 ) s=2f_s\left(\frac{z-1}{z+1}\right) s=2fs?(z+1z?1?)進行替換,得到 H ( z ) = 0.4206 ( z + 1 ) z ? 0.1587 H(z)=\frac{0.4206(z+1)}{z-0.1587} H(z)=z?0.15870.4206(z+1)?,
IIR數字濾波器設計總結
IIR數字濾波器設計總結
-
無限沖激回應濾波器的設計,首先從三個模擬原型低通濾波器的設計開始,即巴特沃斯濾波器、切比雪夫濾波器和橢圓濾波器,
由這三種模擬低通濾波器,通過變數變換的方法,可以得到其他三種模擬濾波器,即模擬高通、模擬帶通和模擬帶阻濾波器,
-
利用模擬濾波器設計數字濾波器,即從s平面向z平面的變換,該變換需要滿足兩個基本要求 :
- 頻率回應保持一致
- 因果穩定性保持一致
因此產生了兩種映射方法,即脈沖回應不變法和雙線性變換法,這里的映射包含四個:
- 模擬低通—>數字低通
- 模擬高通—>數字高通
- 模擬帶通—>數字帶通
- 模擬帶阻—>數字帶阻
-
脈沖回應不變法的優點是設計簡單,缺點是可能會產生頻譜混疊失真,因此只能用于帶限的模擬濾波器,如低通和帶通濾波器,雙線性變換法是一一映射,優點是不會產生頻譜混疊,缺點是變換是非線性的,會產生頻率畸形,但是可以通過頻率的預畸變來加以校正,
-
除了上面的兩種方法外,還有其它的方法,具體參見(點我)


轉載請註明出處,本文鏈接:https://www.uj5u.com/qita/200453.html
標籤:其他
上一篇:數字信號處理實驗一 T3
