主頁 >  其他 > 數字濾波器的設計

數字濾波器的設計

2020-11-03 04:48:44 其他

文章目錄

  • 先導知識
    • 時頻域互換
    • 什么是濾波器?
    • 如何設計濾波器?
    • 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)=Ar=1N?(1?dr?z?1)r=1M?(1?cr?z?1)?H(ejω)=Ar=1N?Pr?r=1M?zr???(ω)=ω(N?M)+r=1N?αr??r=1M?β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=0N?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=0N?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濾波器,其基本步驟如下:

  1. 確定頻域的回應函式 H d ( e j ω ) H_d(e^{j\omega}) Hd?(ejω),低通、高通、帶通或者其它;
  2. 確定頻域的回應函式 H d ( e j ω ) H_d(e^{j\omega}) Hd?(ejω)的傅里葉逆變換,找到連續脈沖回應函式 h d ( t ) h_d(t) hd?(t)
  3. 對連續脈沖回應函式 h d ( t ) h_d(t) hd?(t)按照一定的采樣頻率進行采樣,獲得離散信號 h d ( m ) h_d(m) hd?(m)
  4. 選擇合適的窗函式,對離散信號 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=0N?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=0N?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=0N?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=1N?Wi?[H(fi?)?Hd?(fi?)]2min
其中 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 maxWi?(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數字低通數字濾波器是常用的方法,設計程序是:

  1. 按照數字濾波器技術指標要求設計一個過渡模擬低通濾波器 H a ( s ) H_a(s) Ha?(s)
  2. 再按照一定的轉換關系將 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),這種轉換關系需要滿足一定的條件:

  1. 因果穩定的模擬濾波器轉換成數字濾波器后,仍是因果穩定的;
  2. s平面的左半平面映射到z平面的單位圓內部;
  3. 數字濾波器的頻率回應應模仿模擬濾波器的頻響特性,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} ssω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} sBss2+ω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} ss2+ω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=1N?s?si?Ai??ha?(t)=L[Ha?(s)]=i=1N?Ai?esi?tu(t)h(n)=ha?(nT)=i=1N?Ai?esi?nTu(nT)H(z)=i=1N?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=1N?s?si?Ai??(1)H(z)=i=1N?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?nz=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濾波器了,基本步驟如下:

  1. 選擇一個歸一化的模擬濾波器 H ( s ) H(s) H(s)
  2. 確定數字濾波器的截止頻率,也就是回應為-3dB(幅值衰減為 1 / 2 1/\sqrt{2} 1/2 ?)時對應的頻率;
  3. 利用公式 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??)計算對應的模擬截止頻率
  4. 選擇合適的變換,得到 H ′ ( s ) H'(s) H(s);比如對于低通濾波器: s → s / ω c s{\rightarrow}s/\omega_c ss/ωc?,其中 ω c = 2 π f a \omega_c=2{\pi}f_a ωc?=2πfa?
  5. H ′ ( s ) H'(s) H(s)中,用 s → 2 f s ( z ? 1 z + 1 ) s{\rightarrow}2f_s(\frac{z-1}{z+1}) s2fs?(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

  1. 選擇歸一化的濾波器 H ( s ) = 1 s + 1 H(s)=\frac{1}{s+1} H(s)=s+11?
  2. 數字截止頻率為10Hz;
  3. 對應的模擬截止頻率為: 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
  4. 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 ss/ω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
  5. 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平面的變換,該變換需要滿足兩個基本要求 :

    1. 頻率回應保持一致
    2. 因果穩定性保持一致

    因此產生了兩種映射方法,即脈沖回應不變法和雙線性變換法,這里的映射包含四個:

    1. 模擬低通—>數字低通
    2. 模擬高通—>數字高通
    3. 模擬帶通—>數字帶通
    4. 模擬帶阻—>數字帶阻
  • 脈沖回應不變法的優點是設計簡單,缺點是可能會產生頻譜混疊失真,因此只能用于帶限的模擬濾波器,如低通和帶通濾波器,雙線性變換法是一一映射,優點是不會產生頻譜混疊,缺點是變換是非線性的,會產生頻率畸形,但是可以通過頻率的預畸變來加以校正,

  • 除了上面的兩種方法外,還有其它的方法,具體參見(點我)

在這里插入圖片描述
在這里插入圖片描述

轉載請註明出處,本文鏈接:https://www.uj5u.com/qita/200453.html

標籤:其他

上一篇:數字信號處理實驗一 T3

下一篇:ENVI5.4完美實作MODIS NDVI資料格式轉換和投影變換

標籤雲
其他(157675) Python(38076) JavaScript(25376) Java(17977) C(15215) 區塊鏈(8255) C#(7972) AI(7469) 爪哇(7425) MySQL(7132) html(6777) 基礎類(6313) sql(6102) 熊猫(6058) PHP(5869) 数组(5741) R(5409) Linux(5327) 反应(5209) 腳本語言(PerlPython)(5129) 非技術區(4971) Android(4554) 数据框(4311) css(4259) 节点.js(4032) C語言(3288) json(3245) 列表(3129) 扑(3119) C++語言(3117) 安卓(2998) 打字稿(2995) VBA(2789) Java相關(2746) 疑難問題(2699) 细绳(2522) 單片機工控(2479) iOS(2429) ASP.NET(2402) MongoDB(2323) 麻木的(2285) 正则表达式(2254) 字典(2211) 循环(2198) 迅速(2185) 擅长(2169) 镖(2155) 功能(1967) .NET技术(1958) Web開發(1951) python-3.x(1918) HtmlCss(1915) 弹簧靴(1913) C++(1909) xml(1889) PostgreSQL(1872) .NETCore(1853) 谷歌表格(1846) Unity3D(1843) for循环(1842)

熱門瀏覽
  • 網閘典型架構簡述

    網閘架構一般分為兩種:三主機的三系統架構網閘和雙主機的2+1架構網閘。 三主機架構分別為內端機、外端機和仲裁機。三機無論從軟體和硬體上均各自獨立。首先從硬體上來看,三機都用各自獨立的主板、記憶體及存盤設備。從軟體上來看,三機有各自獨立的作業系統。這樣能達到完全的三機獨立。對于“2+1”系統,“2”分為 ......

    uj5u.com 2020-09-10 02:00:44 more
  • 如何從xshell上傳檔案到centos linux虛擬機里

    如何從xshell上傳檔案到centos linux虛擬機里及:虛擬機CentOs下執行 yum -y install lrzsz命令,出現錯誤:鏡像無法找到軟體包 前言 一、安裝lrzsz步驟 二、上傳檔案 三、遇到的問題及解決方案 總結 前言 提示:其實很簡單,往虛擬機上安裝一個上傳檔案的工具 ......

    uj5u.com 2020-09-10 02:00:47 more
  • 一、SQLMAP入門

    一、SQLMAP入門 1、判斷是否存在注入 sqlmap.py -u 網址/id=1 id=1不可缺少。當注入點后面的引數大于兩個時。需要加雙引號, sqlmap.py -u "網址/id=1&uid=1" 2、判斷文本中的請求是否存在注入 從文本中加載http請求,SQLMAP可以從一個文本檔案中 ......

    uj5u.com 2020-09-10 02:00:50 more
  • Metasploit 簡單使用教程

    metasploit 簡單使用教程 浩先生, 2020-08-28 16:18:25 分類專欄: kail 網路安全 linux 文章標簽: linux資訊安全 編輯 著作權 metasploit 使用教程 前言 一、Metasploit是什么? 二、準備作業 三、具體步驟 前言 Msfconsole ......

    uj5u.com 2020-09-10 02:00:53 more
  • 游戲逆向之驅動層與用戶層通訊

    驅動層代碼: #pragma once #include <ntifs.h> #define add_code CTL_CODE(FILE_DEVICE_UNKNOWN,0x800,METHOD_BUFFERED,FILE_ANY_ACCESS) /* 更多游戲逆向視頻www.yxfzedu.com ......

    uj5u.com 2020-09-10 02:00:56 more
  • 北斗電力時鐘(北斗授時服務器)讓網路資料更精準

    北斗電力時鐘(北斗授時服務器)讓網路資料更精準 北斗電力時鐘(北斗授時服務器)讓網路資料更精準 京準電子科技官微——ahjzsz 近幾年,資訊技術的得了快速發展,互聯網在逐漸普及,其在人們生活和生產中都得到了廣泛應用,并且取得了不錯的應用效果。計算機網路資訊在電力系統中的應用,一方面使電力系統的運行 ......

    uj5u.com 2020-09-10 02:01:03 more
  • 【CTF】CTFHub 技能樹 彩蛋 writeup

    ?碎碎念 CTFHub:https://www.ctfhub.com/ 筆者入門CTF時時剛開始刷的是bugku的舊平臺,后來才有了CTFHub。 感覺不論是網頁UI設計,還是題目質量,賽事跟蹤,工具軟體都做得很不錯。 而且因為獨到的金幣制度的確讓人有一種想去刷題賺金幣的感覺。 個人還是非常喜歡這個 ......

    uj5u.com 2020-09-10 02:04:05 more
  • 02windows基礎操作

    我學到了一下幾點 Windows系統目錄結構與滲透的作用 常見Windows的服務詳解 Windows埠詳解 常用的Windows注冊表詳解 hacker DOS命令詳解(net user / type /md /rd/ dir /cd /net use copy、批處理 等) 利用dos命令制作 ......

    uj5u.com 2020-09-10 02:04:18 more
  • 03.Linux基礎操作

    我學到了以下幾點 01Linux系統介紹02系統安裝,密碼啊破解03Linux常用命令04LAMP 01LINUX windows: win03 8 12 16 19 配置不繁瑣 Linux:redhat,centos(紅帽社區版),Ubuntu server,suse unix:金融機構,證券,銀 ......

    uj5u.com 2020-09-10 02:04:30 more
  • 05HTML

    01HTML介紹 02頭部標簽講解03基礎標簽講解04表單標簽講解 HTML前段語言 js1.了解代碼2.根據代碼 懂得挖掘漏洞 (POST注入/XSS漏洞上傳)3.黑帽seo 白帽seo 客戶網站被黑帽植入劫持代碼如何處理4.熟悉html表單 <html><head><title>TDK標題,描述 ......

    uj5u.com 2020-09-10 02:04:36 more
最新发布
  • 2023年最新微信小程式抓包教程

    01 開門見山 隔一個月發一篇文章,不過分。 首先回顧一下《微信系結手機號資料庫被脫庫事件》,我也是第一時間得知了這個訊息,然后跟蹤了整件事情的經過。下面是這起事件的相關截圖以及近日流出的一萬條資料樣本: 個人認為這件事也沒什么,還不如關注一下之前45億快遞資料查詢渠道疑似在近日復活的訊息。 訊息是 ......

    uj5u.com 2023-04-20 08:48:24 more
  • web3 產品介紹:metamask 錢包 使用最多的瀏覽器插件錢包

    Metamask錢包是一種基于區塊鏈技術的數字貨幣錢包,它允許用戶在安全、便捷的環境下管理自己的加密資產。Metamask錢包是以太坊生態系統中最流行的錢包之一,它具有易于使用、安全性高和功能強大等優點。 本文將詳細介紹Metamask錢包的功能和使用方法。 一、 Metamask錢包的功能 數字資 ......

    uj5u.com 2023-04-20 08:47:46 more
  • vulnhub_Earth

    前言 靶機地址->>>vulnhub_Earth 攻擊機ip:192.168.20.121 靶機ip:192.168.20.122 參考文章 https://www.cnblogs.com/Jing-X/archive/2022/04/03/16097695.html https://www.cnb ......

    uj5u.com 2023-04-20 07:46:20 more
  • 從4k到42k,軟體測驗工程師的漲薪史,給我看哭了

    清明節一過,盲猜大家已經無心上班,在數著日子準備過五一,但一想到銀行卡里的余額……瞬間心情就不美麗了。最近,2023年高校畢業生就業調查顯示,本科畢業月平均起薪為5825元。調查一出,便有很多同學表示自己又被平均了。看著這一資料,不免讓人想到前不久中國青年報的一項調查:近六成大學生認為畢業10年內會 ......

    uj5u.com 2023-04-20 07:44:00 more
  • 最新版本 Stable Diffusion 開源 AI 繪畫工具之中文自動提詞篇

    🎈 標簽生成器 由于輸入正向提示詞 prompt 和反向提示詞 negative prompt 都是使用英文,所以對學習母語的我們非常不友好 使用網址:https://tinygeeker.github.io/p/ai-prompt-generator 這個網址是為了讓大家在使用 AI 繪畫的時候 ......

    uj5u.com 2023-04-20 07:43:36 more
  • 漫談前端自動化測驗演進之路及測驗工具分析

    隨著前端技術的不斷發展和應用程式的日益復雜,前端自動化測驗也在不斷演進。隨著 Web 應用程式變得越來越復雜,自動化測驗的需求也越來越高。如今,自動化測驗已經成為 Web 應用程式開發程序中不可或缺的一部分,它們可以幫助開發人員更快地發現和修復錯誤,提高應用程式的性能和可靠性。 ......

    uj5u.com 2023-04-20 07:43:16 more
  • CANN開發實踐:4個DVPP記憶體問題的典型案例解讀

    摘要:由于DVPP媒體資料處理功能對存放輸入、輸出資料的記憶體有更高的要求(例如,記憶體首地址128位元組對齊),因此需呼叫專用的記憶體申請介面,那么本期就分享幾個關于DVPP記憶體問題的典型案例,并給出原因分析及解決方法。 本文分享自華為云社區《FAQ_DVPP記憶體問題案例》,作者:昇騰CANN。 DVPP ......

    uj5u.com 2023-04-20 07:43:03 more
  • msf學習

    msf學習 以kali自帶的msf為例 一、msf核心模塊與功能 msf模塊都放在/usr/share/metasploit-framework/modules目錄下 1、auxiliary 輔助模塊,輔助滲透(埠掃描、登錄密碼爆破、漏洞驗證等) 2、encoders 編碼器模塊,主要包含各種編碼 ......

    uj5u.com 2023-04-20 07:42:59 more
  • Halcon軟體安裝與界面簡介

    1. 下載Halcon17版本到到本地 2. 雙擊安裝包后 3. 步驟如下 1.2 Halcon軟體安裝 界面分為四大塊 1. Halcon的五個助手 1) 影像采集助手:與相機連接,設定相機引數,采集影像 2) 標定助手:九點標定或是其它的標定,生成標定檔案及內參外參,可以將像素單位轉換為長度單位 ......

    uj5u.com 2023-04-20 07:42:17 more
  • 在MacOS下使用Unity3D開發游戲

    第一次發博客,先發一下我的游戲開發環境吧。 去年2月份買了一臺MacBookPro2021 M1pro(以下簡稱mbp),這一年來一直在用mbp開發游戲。我大致分享一下我的開發工具以及使用體驗。 1、Unity 官網鏈接: https://unity.cn/releases 我一般使用的Apple ......

    uj5u.com 2023-04-20 07:40:19 more