IIR數字濾波器的設計

Mrakdown常見數學符號

Markdown基礎語法

沖激響應不變法

• 沖激響應不變法：就是用其單位沖激響應序列模仿模擬濾波器的單位沖激響應的抽樣值

• 設計的具體步驟及方法

  ? 首先要設計一個響應的模擬濾波器。
  ? 模擬濾波器的單位沖激響應非周期 \(\Rightarrow\) 頻率離散（- \(\infty\)-+\(\infty\) ),頻譜范圍可能大于折疊頻率，即奈奎斯特頻率 \(\Rightarrow\) 取樣，頻率周期延拓 \(\Rightarrow\) 頻譜可能產生混疊。
  ? 數字域中極點在單位圓內\(\Leftrightarrow\)模擬域中極點在左半平面。
  可以通過提高抽樣頻率來減少混疊，但設計指標若以數字域頻率給定時，不能通過提高抽樣頻率改善混疊現象。因為\(\Omega=\frac{\omega}{T}\) ，\(\omega\)是數字域頻率，\(\Omega\)是模擬域頻率，增大抽樣頻率\(\Rightarrow\)減小采樣周期\(\Rightarrow\)增大截止頻率。

  ? 模擬域 ：H_a(s) = \(\sum_{k=1}^N \frac{A_k}{s-s_k}\) 數字域：H(z) = \(\sum_{k=1}^N \frac{A_k}{1-e^{s_kT}z^{-1}}\)
  ? s平面：s=s_k z平面：z=\(e^{s_kT}\) 兩者系數相同

  • 當采樣周期很小時，數字濾波器增益很大，容易溢出，可以修正，令h(n)=T\(h_a(nT)\)，因為H(e^jw)=\(\frac{1}{T}H_a(j\frac{\omega}{T})\)
  • 該方法對應的MATLAB程序：[BZ, AZ] = impinvar(B,A,Fs),其中BZ，AZ為數字濾波器單位沖擊響應的分子和分母系數；B，A為模擬濾波器的分子和分母的系數，Fs為抽樣頻率。

• 優點

  • 時域逼近良好
  • 保持線性關系：\(\omega\)=\(\Omega\)T

• 缺點

  • 頻率響應混疊，只適用低通、帶通濾波器。

階躍響應不變法

• 變換原理：數字濾波器的階躍響應模仿模擬濾波器的階躍響應，即g(n)=g_a(t)|_t=nt=g_a(nT).
• g(n)=u(n)*h(n)\(\rightarrow\) G(z)=\(\frac{z}{z-1}H(z)\) ,g_a(t)=u(t)*h_a(t)\(\rightarrow\)G_a(s)=\(\frac{1}{s}H_a(s)\)
  變換步驟：
  \(H_a(s)\rightarrow G_a(s)\rightarrow g_a(t)\rightarrow g(n)\rightarrow G(z)\rightarrow H(z)\)
• 優缺點同沖激響應不變法

雙線性不變法

• 為什么需要雙線性變換法？--避免頻譜混疊

? 在推導過程中有一個中間域，最終的關系式為z=\(\frac{1+s}{1-s}\) ,為了使模擬濾波器與數字濾波器的某一頻率有對應關系，引入一個常數c，則有s=c\(\frac{1-\frac{1}{z}}{1+\frac{1}{z}}\) ,即z=\(\frac{c+s}{c-s}\)。若低頻處有確切的對應關系應為c=\(\frac{2}{T}\) 。
? s平面虛軸與z平面單位圓的對應關系\(\Omega\) =ctan\(\frac{w}{2}\) \(\Omega\) 為s域,w為數字域

• 優點：無頻率混疊，避免了多值映射。

• 缺點：除零頻率附近，嚴重非線性關系，在臨界頻率點（通帶截止頻率和阻帶截止頻率--這個是技術指標)產生畸變。

• 解決臨界頻率點的畸變--預畸變
  \(\Omega\) =\(c\tan\frac{w}{2}\) ,c一般取2/T，在臨界頻率點根據w去設計\(\Omega\) ,得到模擬濾波器對應的技術指標，再去設計模擬濾波器，再用雙線性變換法設計數字濾波器。

• 該方法對應的MATLAB方法：[BZ, AZ] = bilinear(B,A,Fs)

常用模擬濾波器設計

• 由幅度平方函數確定模擬濾波器的系統函數
  \(|H_a(j\Omega)|^2=H_a(j\Omega)H_a^*(j\Omega)=H_a(j\Omega)H_a(-j\Omega)=H_a(s)H_a(-s)|_{(s=j\Omega)}\) ,\(H_a(j\Omega)\)可以看做模擬濾波器單位沖激響應的傅里葉變化。
  \(H_a(s)\)和\(H_a(-s)\)的零極點對稱，成對出現，且其零極點分別分布在左半平面（只有分布在左半平面系統才是穩定的）和右半平面。

巴特沃斯低通逼近

幅度平方函數為\(|H_a(j\Omega)|^2 = \frac{1}{1+(\frac{\Omega}{\Omega_c})^{2N}}\) ,其中N為濾波器階數，\(\Omega_c\)為3dB截止頻率，因為\(\Omega=\Omega_c\)時，幅度為原來的\(\frac{1}{\sqrt{2}}\)，

巴特沃斯幅度特性及其與N的關系

• 極點分布

\(|H_a(j\Omega)|_{\Omega=\frac{s}{j}}^2 = H_a(s)H_a(-s) = \frac{1}{1+(\frac{s}{j\Omega})^{2N}}\) ，令\(1+(\frac{s}{j\Omega})^{2N}=0\)即可解得極點分布。

極點在s平面呈象限對稱，個數為2N，極點間角度間隔為\(\frac{\pi}{N}\)，極點不落在虛軸上，若N為奇數，實軸上有極點；若為偶數，實軸上無極點。

• 濾波器的系統函數

\(H_a(s)=\frac{\Omega_c^N}{\prod_{k=1}^N(s-s_k)}\)，s_k為左半平面的極點。若令\(\Omega_c=1rad/s\),得到的各個階次的系統函數為歸一化的系統函數\(H_{an}(s)\)，其與原系統函數之間的關系為\(H_a(s)=H_{an}(s)|_{s=\frac{s}{\Omega_c}}\)（去歸一化）。

• 如何確定N和\(\Omega_c\)

已知技術指標：\(\Omega_c(通帶截止頻率)、\delta_1(通帶內允許的最大衰減)、\Omega_s(阻帶截止頻率)、\delta_2(阻帶允許的最小衰減)\).
\(\delta_1=-20lg|H_a(j\Omega)|,|H_a(j\Omega)|^2=\frac{1}{1+(\frac{\Omega_p}{\Omega_c})^{2N}}\rightarrow1+(\frac{\Omega_p}{\Omega_c})^{2N}=10^{0.1\delta_1}\rightarrow N\) ①
同理,有\(1+(\frac{\Omega_s}{\Omega_c})^{2N}=10^{0.1\delta_2}\) ② 由這兩個方程可以解得N和\(\Omega_c\) 。

\(\Omega_c\)既可以若由第一個方程解得，那么等同于將通帶的技術指標卡死，表示阻帶指標有富余，若由第二個方程解得，則與前一種情況相反。

• MATLAB仿真

[N, Wc]=buttord(wp,ws,Rp,Rs,'s') %根據性能指標得到巴特沃斯濾波器的階數及截止頻率
[z,p,k]=buttap(N) %返回N階歸一化原型模擬濾波器的零極點以及增益
[B,A]=zp2tf(z,p,k) %根據零極點增益得到原型多項式的系數（歸一化）
[Bt,At]=lp2lp(B,A,Wo) %得到模擬低通濾波器的多項式系數（去歸一化）
[BZ,AZ]=bilinear(Bt,At,Fs) %使用雙線性變換法求得數字低通濾波器額的多項式系數

切比雪夫低通逼近

• 類型
  I型：通帶范圍內波動；II型：阻帶范圍內波動。
  相比于巴特沃斯：巴特沃斯在整個頻帶內都是下降的，卡死通帶或阻帶中其中一個，都將引起另一個的快速變化。

• 幅度平方函數（以I型為例）
  \(|H_a(j\Omega)|^2=\frac{1}{1+\varepsilon^2C_N^2(\frac{\Omega}{\Omega_c})}\)
  \(\varepsilon：\)通帶波紋大小，\(0<\varepsilon<1，\)值越大，波紋越大
  \(\Omega_c：\)通帶截止頻率，不一定為3dB帶寬
  \(N：\)濾波器階數
  \(C_N(x)：\)N階切比雪夫多項式

  \[C_N(x)=\begin{cases} \cos(N\cos^{-1}x)& \text{|x|<=1}\\ ch(Nch^{-1}x)&|x|>1 \end{cases}$$,該多項式的圖像如下： <div align=center> <img src="https://ss0.bdstatic.com/70cFuHSh_Q1YnxGkpoWK1HF6hhy/it/u=1536110891,3134924895&fm=26&gp=0.jpg"/> <p> 切比雪夫多項式圖像 </p> </div> 從圖像可以看出，n為奇數時，$|C_N(0)|=0$；n為偶數時，$|C_N(0)|=1$。 $\Omega=\Omega_c時，|C_N(x)|=1，|H_a(j\Omega)|=\frac{1}{\sqrt{1+\varepsilon^2}}$ $在|x|<=1時，|C_N(x)|在0-1$間反復變化，導致$|H_a(j\Omega)|在\frac{1}{\sqrt{1+\varepsilon^2}}$之間反復變化，從而形成波紋$ 。 $而$\delta_1$表示通帶內允許的最大衰減，故應有$\delta_1=20lg\sqrt{1+\varepsilon^2}，$從而可以解得參數$\varepsilon。$ 再將阻帶允許的最小衰減和阻帶截止頻率帶入即可得到階數N。 \]

  極點個數為2N個，分布在橢圓上。

• 設計步驟
  確定技術指標：\(\Omega_p(=\Omega_c),\delta_1,\Omega_s,\delta_2\)
  根據技術指標求出濾波器階數N以及\(\varepsilon\)：\(N>=\frac{ch^{-1}(k_1^{-1})}{ch^{-1}\lambda_s}\)，\(\varepsilon^2=10^{0.1\delta}-1\)，其中\(k_1^{-1}=\sqrt{\frac{10^{0.1\delta_2-1}}{10^{0.1\delta_1-1}}}\) 。
  由N和\(\delta_1\)，直接查表得到歸一化的系統函數\(H_{an}(s)\)。
  去歸一化。

• MATLAB仿真

  [N,Wc]=cheblord(wp,ws,Rp,Rs,'s') %根據性能指標求出切比雪夫濾波器的階數和截止頻率
  [z,p,k]=cheb1ap(N,rp) %返回N階歸一化原型模擬濾波器的零極點及增益
  [B,A]=zp2tf(z,p,k) %求出歸一化低通原型多項式系數
  [Bt,At]=lp2lp(B,A,Wo) %求出模擬低通濾波器的系數
  [BZ,AZ]=impinvar(Bt,At,Fs) %沖激響應不變法求數字低通濾波器的系數
  

頻率變換法

• 為什么需要頻率變換法？
  之間的設計方法都是基于低通濾波器進行設計的，為了得到其他類型的濾波器，需要進行頻率變換

模擬域頻帶變換法

• 變換步驟
  歸一化模擬低通濾波器\(\Rightarrow\)模擬域頻帶變換\(\Rightarrow\)模擬低通、高通、帶通、帶阻濾波器\(\Rightarrow\)雙線性變換\(\Rightarrow\)數字低通、高通、帶通、帶阻濾波器。
  總結來講，就是在模擬域就進行頻帶變換。
• 具體方法
  歸一化低通濾波器的系統函數：\(H_a(s),s=\sigma+j\Omega\)
  頻帶變換后的系統函數：\(\overline{H}_s(\overline{s}),\overline{s}=\overline{\sigma}+j\overline{\Omega}\)
  • 模擬低通\(\rightarrow\)模擬低通
    \(s=\frac{\Omega_c\overline{s}}{\overline{\Omega_c}}\)，當\(\Omega_c=1\)時，相當于去歸一化。
    對應的MATLAB函數為[Bt,At]=lp2lp(B,A,Wo)
  • 模擬低通\(\rightarrow\)模擬帶通
    \(s=\frac{\overline{s}^2+\overline{\Omega}_0^2}{\overline{s}}\)，\(\overline{\Omega}_0\)為帶通濾波器的中心頻率，\(\overline{\Omega}_0=\overline{\Omega}_1\overline{\Omega}_2\)。\(\overline{\Omega}_1、\overline{\Omega}_2\)分別為帶通濾波器的上通帶截止頻率和下通帶截止頻率。為了得到歸一化的帶通濾波器，可令\(s=\frac{\overline{s}^2+\overline{\Omega}_0^2}{\overline{s}}\)兩邊同時除以\(\Omega_c(B)\).
    對應的MATLAB函數為[Bt,At]=lp2bp(B,A,Wo,Bw)，Bw為帶寬。
    數字帶通濾波器性能指標\(\rightarrow\)模擬帶通濾波器性能能指標\(\rightarrow\)歸一化\(\rightarrow\)模擬低通濾波器歸一化性能指標\(\rightarrow\)模擬低通濾波器。
    帶通濾波器歸一化頻率計算方法：\(\eta=\frac{\overline{\Omega}}{B}\)
    低通濾波器歸一化頻率計算方法：\(\lambda=\frac{\Omega}{\Omega_c},\Omega_c\)為通帶截止頻率。
    兩者間關系：\(\lambda=\frac{\eta^2-\eta_0^2}{\eta}\)。

數字域頻帶變換法

• 變換步驟
  歸一化模擬低通濾波器\(\Rightarrow\)數字低通濾波器\(\Rightarrow\)數字域頻帶變換\(\Rightarrow\)數字低通、高通、帶通、帶阻濾波器。
  總結來講，就是在數字域進行頻帶變換。

• 具體方法
  低通數字濾波器的系統函數：\(H_L(z)\)
  頻帶變化后的數字濾波器的系統函數：\(Hd(Z)\)
  要求：\(\tag{1}\)兩個系統函數的單位圓要對應；\(\tag{2}\) 都是因果穩定；
  關系：\(H_d(Z)=H_L(z)|_{z^{-1}=G(Z^{-1})}\) (\(G(Z^{-1})=\pm\prod_{k=1}^N\frac{Z^{-1}-a^*}{1-a_kZ^{-1}}\)（全通函數）應為有理函數)

  • 數字低通\(\rightarrow\)數字低通
    \(N=1,z^{-1}=G(Z^{-1})=\frac{Z^{-1}-\alpha}{1-\alpha Z^{-1}}\),根據兩個低通濾波器的對應關系（技術指標）,最終得到\(\alpha=\frac{\sin\frac{\theta_c-\omega_c}{2}}{\sin\frac{\theta_c+\omega_c}{2}}\)

  • 數字低通\(\rightarrow\)數字高通
    低通頻率響應在單位圓上旋轉\(180^\circ\),即得到高通頻率響應：Z=-Z。

    \(G(Z^{-1})=\frac{-Z^{-1}-\alpha}{1+\alpha Z^{-1}}\)