2021/03/15

濾波器設計

首先介紹濾波器的基本概念：窗函數。濾波器的設計分成 FIR, IIR 兩種，其中 FIR Filter 使用 Window Method，IIR Filter 則包含 Butterworth Filter, Chebyshev Type I Filter, Chebyshev Type II Filter, 及 Elliptic Filter

濾波器在 DSP 系統中，可用來選取或抑制某些頻率範圍的輸入訊號，藉以產生符合需求的輸出訊號。

由於有限長度的脈衝響應無法實現 ideal filter，在設計實際的濾波器時，通常是考慮濾波器的頻率響應是否符合事先定義的規格，同時容許頻率響應有些微誤差。

實際的 lowpass filter 其頻率響應應如下圖，通過頻帶 passband 的範圍介於 \(0 ～ ω_p\) 之間，理想的頻率響應為 1；抑制頻帶 stopband 的範圍介於 \(ω_s ～ \pi\) 之間，理想的頻率響應為 0。其中 \(ω_p\) 稱為通過頻帶邊緣頻率 Passband Edge Frequency，\(ω_s\) 稱為抑制頻帶邊緣頻率 Stopband Edge Frequency。

實際的 lowpass filter 的設計，通常在通過頻帶或抑制頻帶內，容許頻率響應有些許波動（誤差），其中\(δ_p\) 稱為通過頻帶波紋 (Passband Ripple)，\(δ_s\) 稱為抑制頻帶波紋 (Stopband Ripple)。在 \(ω_p ～ ω_s\) 之間，則是容許頻率響應從 1 逐漸降為 0，稱為過渡頻帶 (Transition Band)

filter 的設計，主要包含以下步驟：

  1. 定義濾波器規格：根據 DSP 應用的需求，選取濾波器的種類，例如：lowpass, highpass, bandpass, bandstop 等等。然後定義濾波器的規格，包含：passbadn 與 stopband 的截止頻率、取樣頻率等等參數

  2. 選取濾波器的種類：選取 FIR 或 IIR 濾波器作為設計目標。通常 FIR Filter 為穩定系統，若選取 IIR Filter ，則需要進一步檢驗其穩定性。

  3. 求轉換函式：根據濾波器規格，求濾波器的轉換函式 H(z)，例如 FIR Filter 的轉換函式為

    \(H(z) = b_0+b_1z^{-1}+...+b_Mz^{-M}\)

    IIR 濾波器的轉換函式為：

    \(H(z) = \frac{b_0+b_1z^{-1}+...+b_Mz^{-M}}{a_0+a_1z^{-1}+...+a_Nz^{-N}}\)

    這個步驟的目的是求濾波器的係數：\(\{a_k\} k=0,1,...,N\) 與 \(\{b_k\} k=0,1,...,M\)

  4. 實現濾波器：在完成設計後，可採用硬體或軟體方式實現。

Window Function

idea lowpass filter 在時間域的脈衝響應為無限序列，因此，為了實現有限的脈衝響應，同時近似理想的 lowpass filter，最直接的方法是 Window Method。也就是定義 Window Function，用來擷取有限的脈衝響應。

Window Function 定義為：介於某特定區間的數學函示，區間外的數值為 0

任意函數與 window function 相乘後，結果就像是透過窗戶觀察該函數一樣。在 DSP，window function 經常用來進行頻譜分析、濾波器設計、音頻壓縮等等。

以下是幾個典型的 window function，這幾種 window function 都具有對稱性，SciPy 的Signal 有支援這些 functions。在濾波器設計過程中，窗函數通常是先經過平移，平移後以原點為中心，且對原點對稱，再與無限脈衝響應進行運算，藉以擷取有限的離散序列。

  • 矩形窗 Rectangular Window

    \( W(n) = \left\{ \begin{array}{ll} 1 & \mbox{if 0 ≤ n ≤ M-1} \\ 0 & \mbox{otherwise} \end{array} \right.\)

    其中 window 的大小為 M。由於矩形窗的形狀如同箱型車 boxcar，因此也被稱為 boxcar function

  • Hamming 窗

    \( W(n) = \left\{ \begin{array}{ll} 0.54-0.46cos(\frac{2 \pi n}{M-1}) & \mbox{if 0 ≤ n ≤ M-1} \\ 0 & \mbox{otherwise} \end{array} \right.\)

    這是根據 Richard W. Hamming 而命名的

  • Hanning 窗

    \( W(n) = \left\{ \begin{array}{ll} 0.5-0.5cos(\frac{2 \pi n}{M-1}) & \mbox{if 0 ≤ n ≤ M-1} \\ 0 & \mbox{otherwise} \end{array} \right.\)

    這是根據 Julus von Hann 而命名的，形狀跟 Hamming 窗相似

  • Bartlet 窗

    \( W(n) = \left\{ \begin{array}{ll} \frac{2}{M-1}(\frac{M-1}{2} - |n-\frac{M-1}{2}|) & \mbox{if 0 ≤ n ≤ M-1} \\ 0 & \mbox{otherwise} \end{array} \right.\)

    Bartlet window 的形狀是三角形

  • Blackman 窗

    \( W(n) = \left\{ \begin{array}{ll} 0.42-0.5cos(\frac{2 \pi n}{M-1})+0.08cos(\frac{4 \pi n}{M-1}) & \mbox{if 0 ≤ n ≤ M-1} \\ 0 & \mbox{otherwise} \end{array} \right.\)

    這是根據 Ralph B. Blackmand 而命名的

  • Kaiser 窗

    \(W(n) = I_0(β\sqrt{1-\frac{4n^2}{(M-1)^2}})/I_0(β), -\frac{M-1}{2} ≤ n ≤ \frac{M-1}{2}\)

    這是根據 Jim Kaiser 命名的。

    \(I_0\) 是修正的零階 Bessel 函式，參數 β 可用來近似其他窗函式。例如：β=0 為矩形窗、β=5 近似 Hamming 窗、β=6 近似 Hanning 窗、β=8.6 近似 Blackman 窗、β=14 則是 Kaiser 窗的建議初始值

統一選取 M=2*32+1 = 65，為有限脈衝響應的長度。window function 為對稱圖形，可平移到以原點為中心。

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt

M = 65
w1 = signal.boxcar( M )
w2 = signal.hamming( M )
w3 = signal.hann( M )
w4 = signal.bartlett( M )
w5 = signal.barthann( M )
w6 = signal.kaiser( M, 14 )

plt.subplot(231)
plt.plot( w1 )
plt.xlabel( 'n (rectangular, boxcar)' )
plt.ylabel( 'Amplitude' )

plt.subplot(232)
plt.plot( w2 )
plt.xlabel( 'n (hamming)' )
plt.ylabel( 'Amplitude' )

plt.subplot(233)
plt.plot( w3 )
plt.xlabel( 'n (hanning)' )
plt.ylabel( 'Amplitude' )

plt.subplot(234)
plt.plot( w4 )
plt.xlabel( 'n (bartlett)' )
plt.ylabel( 'Amplitude' )

plt.subplot(235)
plt.plot( w5 )
plt.xlabel( 'n (blackman)' )
plt.ylabel( 'Amplitude' )

plt.subplot(236)
plt.plot( w6 )
plt.xlabel( 'n (kaiser)' )
plt.ylabel( 'Amplitude' )

plt.show( )

FIR 濾波器設計

最簡單的方法是 window method，步驟如下

  1. 選取濾波器的種類：包含 lowpass, highpass, bandpass, bandstop 濾波器
  2. 定義濾波器的規格：如果是 lowpass, highpass，需定義截止頻率 (cutoff frequency)，以 Hz 為單位，如果是 bandpass, bandstop，需定義兩個截止頻率
  3. 選取取樣頻率：sampling frequency，以 Hz 為單位。選取原則以輸入數位訊號為主，且對應的 Nyquist 頻率需大於上述的截止頻率
  4. 套用窗法：選取 window function。在此需選取濾波器的長度，也就是離散時間域的脈衝響應樣本數 (number of taps)
  5. FIR 濾波器係數：根據濾波器的種類，計算 FIR filter 的係數，並套用 window function，藉以擷取有限長度的脈衝響應。理論上，樣本數越多，越接近 ideal filter
  6. 頻率響應：根據脈衝響應，求其在頻率域的頻率響應，包含強度頻譜與相位頻譜，用來觀察濾波器的設計，是否符合原先地定義的規格需求

ex: 試用 window method 設計濾波器規格為

(1) lowpass filter：截止頻率 \(f_c=250Hz\)

(2) highpass filter：截止頻率 \(f_c=250Hz\)

(3) bandpass filter：截止頻率 \(f_1=200Hz, f_2=300Hz\)

(4) bandstop filter：截止頻率 \(f_1=200Hz, f_2=300Hz\)

並顯示強度頻率與相位頻譜

以下是用不同 window method 實作 lowpass filter 的結果

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt

def fir_filter(n, cutoff, win, pass_zero, freq, plotpos1, plotpos2):
    h = signal.firwin( n, cutoff, window = win, pass_zero = pass_zero, fs = freq )

    w, H = signal.freqz( h )
    magnitude = abs( H )
    phase = np.angle( H )

    plt.subplot(plotpos1)
    plt.plot( w, magnitude )
    plt.xlabel( r'$\omega$ '+'強度頻譜 ('+str(win)+')' )
    plt.ylabel( 'Magnitude' )

    plt.subplot(plotpos2)
    plt.plot( w, phase )
    plt.xlabel( r'$\omega$ '+'相位頻譜 ('+str(win)+')' )
    plt.ylabel( 'Phase' )

# cutoff frequency(Hz)
cutoff = 250
# numeber of taps
n = 31
# pass_zero = {True, False, 'bandpass', 'lowpass', 'highpass', 'bandstop'}
pass_zero = 'lowpass'
# sampling frequency (Hz)
freq = 1000


# 替換中文字型，處理 matplotlib 中文 label 問題
# https://blog.csdn.net/gmr2453929471/article/details/78655834
plt.rcParams['font.sans-serif'] = ['Heiti TC'] # 步驟一（替換sans-serif字型）
plt.rcParams['axes.unicode_minus'] = False  # 步驟二（解決座標軸負數的負號顯示問題）

plt.figure( 1 )

# rectangular (boxcar) window method
win='boxcar'
fir_filter(n, cutoff, win, pass_zero, freq, 231, 234)

# hamming window method
win='hamming'
fir_filter(n, cutoff, win, pass_zero, freq, 232, 235)

# hanning window method
win='hanning'
fir_filter(n, cutoff, win, pass_zero, freq, 233, 236)

plt.figure( 2 )

# bartlett window method
win='bartlett'
fir_filter(n, cutoff, win, pass_zero, freq, 231, 234)

# blackman window method
win='blackman'
fir_filter(n, cutoff, win, pass_zero, freq, 232, 235)

# kaiser window method
win=('kaiser',14)
fir_filter(n, cutoff, win, pass_zero, freq, 233, 236)

plt.show( )

使用 rectangular window method 設計的 lowpass filter，其頻率響應有 Gibbs 現象，結果較不理想。使用 blackman 或kaiser 設計的 lowpass filter，結果較理想。

截止角頻率經過正規化為 \(2 \pi f_c/f_s = 2 \pi (250) / 1000 = \pi/2\)，角頻率介於 \(0～\pi\) 之間

以下是用 hamming window method 實作四種 filter 的結果

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt

def fir_filter(n, cutoff, win, pass_zero, freq, plotpos1):
    h = signal.firwin( n, cutoff, window = win, pass_zero = pass_zero, fs = freq )

    w, H = signal.freqz( h )
    magnitude = abs( H )
    phase = np.angle( H )

    plt.subplot(plotpos1)
    plt.plot( w, magnitude )
    plt.xlabel( r'$\omega$ '+'('+str(pass_zero)+')' )
    plt.ylabel( 'Magnitude' )

    # plt.subplot(plotpos2)
    # plt.plot( w, phase )
    # plt.xlabel( r'$\omega$ '+'相位頻譜 ('+str(win)+')' )
    # plt.ylabel( 'Phase' )

# cutoff frequency(Hz)
cutoff = 250
f1 = 200
f2 = 300

# numeber of taps
n = 31
# sampling frequency (Hz)
freq = 1000

# hamming window method
win='hamming'

plt.figure( 1 )

# lowpass
# pass_zero = {True, False, 'bandpass', 'lowpass', 'highpass', 'bandstop'}
pass_zero='lowpass'
fir_filter(n, cutoff, win, pass_zero, freq, 221)

pass_zero='highpass'
fir_filter(n, cutoff, win, pass_zero, freq, 222)

pass_zero='bandpass'
fir_filter(n, [f1, f2], win, pass_zero, freq, 223)

pass_zero='bandstop'
fir_filter(n, [f1, f2], win, pass_zero, freq, 224)

plt.show( )

IIR 濾波器設計

IIR filters 分為以下幾種

  1. Butterworth Filter
  2. Chebyshev Type I Filter
  3. Chebyshev Type II Filter
  4. Elliptic Filter 橢圓濾波器

IIR filter 的設計步驟，跟 FIR filter 的步驟類似

  1. 選取濾波器的種類：包含 lowpass, highpass, bandpass, bandstop 濾波器

  2. 定義濾波器的規格：如果是 lowpass, highpass，需定義截止頻率 (cutoff frequency)，以 Hz 為單位，如果是 bandpass, bandstop，需定義兩個截止頻率。另外需要定義通過頻帶波紋與抑制頻帶波紋，單位為 dB

  3. 選取取樣頻率：sampling frequency，以 Hz 為單位。選取原則以輸入數位訊號為主，且對應的 Nyquist 頻率需大於上述的截止頻率

  4. 濾波器參數：根據濾波器的規格，計算濾波器的階數 (Order) 與 -3dB 截止頻率

  5. IIR 濾波器係數：根據濾波器的種類及相關參數，計算轉換函式

    \(H(z)=\frac{Y(z)}{X(z)} = \frac{b_0+b_1z^{-1}+...+b_Mz^{-M}}{a_0+z_1z^{-1}+...+a_Nz^{-N}}\)

    也就是計算濾波器的係數 \(\{a_k\}, k=1,2,...,N\) \(\{b_k\}, k=0,1,...M\)

  6. 頻率響應：根據脈衝響應，求其在頻率域的頻率響應，包含強度頻譜與相位頻譜，用來觀察濾波器的設計，是否符合原先地定義的規格需求

Butterworth 濾波器

是一種通過頻帶的頻率響應非常平坦的 DSP，由英國工程師 Stephen Butterworth 提出

根據下列規格，設計 Butterworth filters

(1) lowpass filter & highpass filter

通過頻帶邊緣頻率 200 Hz

抑制頻帶邊緣頻率 300 Hz

通過頻帶波紋 0.5 dB

抑制頻帶波紋 50 dB

取樣頻率 1000Hz

(2) bandpass & bandstop filter

通過頻帶邊緣頻率 200, 300 Hz

抑制頻帶邊緣頻率 100, 400 Hz

通過頻帶波紋 0.5 dB

抑制頻帶波紋 50 dB

取樣頻率 1000Hz

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt

def plot(a, b, plotname, plotpos1, plotpos2):
    w, H = signal.freqz( b, a )
    magnitude = abs( H )
    phase = np.angle( H )

    plt.subplot(plotpos1)
    plt.plot( w, magnitude )
    plt.xlabel( r'$\omega$ 強度頻譜('+plotname+')' )
    plt.ylabel( 'Magnitude' )

    plt.subplot(plotpos2)
    plt.plot( w, phase )
    plt.xlabel( r'$\omega$ 相位頻譜('+plotname+')' )
    plt.ylabel( 'Phase' )

# 替換中文字型，處理 matplotlib 中文 label 問題
# https://blog.csdn.net/gmr2453929471/article/details/78655834
plt.rcParams['font.sans-serif'] = ['Heiti TC'] # 步驟一（替換sans-serif字型）
plt.rcParams['axes.unicode_minus'] = False  # 步驟二（解決座標軸負數的負號顯示問題）


# lowpass filter
# passband edge frequency(Hz)
fp = 200
# stopband edge frequency(Hz)
fs = 300
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp = 2 * fp / Fs
ws = 2 * fs / Fs

n, wn = signal.buttord( wp, ws, rp, rs )
b, a = signal.butter( n, wn, 'lowpass' )

plt.figure( 1 )
plot(a, b, 'lowpass', '121', '122')


# highpass filter
# passband edge frequency(Hz)
fp = 200
# stopband edge frequency(Hz)
fs = 300
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp = 2 * fp / Fs
ws = 2 * fs / Fs

n, wn = signal.buttord( wp, ws, rp, rs )
b, a = signal.butter( n, wn, 'highpass' )

plt.figure( 2 )
plot(a, b, 'highpass', '121', '122')


# bandpass filter
# 1st passband edge frequency(Hz)
fp1 = 200
# 2nd passband edge frequency(Hz)
fp2 = 300
# 1st stopband edge frequency(Hz)
fs1 = 100
# 2nd stopband edge frequency(Hz)
fs2 = 400
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp1 = 2 * fp1 / Fs
wp2 = 2 * fp2 / Fs
ws1 = 2 * fs1 / Fs
ws2 = 2 * fs2 / Fs

n, wn = signal.buttord( [wp1, wp2], [ws1, ws2], rp, rs )
b, a = signal.butter( n, wn, 'bandpass' )

plt.figure( 3 )
plot(a, b, 'bandpass', '121', '122')


# bandstop filter
# 1st passband edge frequency(Hz)
fp1 = 200
# 2nd passband edge frequency(Hz)
fp2 = 300
# 1st stopband edge frequency(Hz)
fs1 = 100
# 2nd stopband edge frequency(Hz)
fs2 = 400
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp1 = 2 * fp1 / Fs
wp2 = 2 * fp2 / Fs
ws1 = 2 * fs1 / Fs
ws2 = 2 * fs2 / Fs

n, wn = signal.buttord( [wp1, wp2], [ws1, ws2], rp, rs )
b, a = signal.butter( n, wn, 'bandstop' )

plt.figure( 4 )
plot(a, b, 'bandstop', '121', '122')

plt.show( )

Chebyshev Type I Filter

Chebyshev Filter 是一種在通過頻帶或抑制頻帶上頻率響應等波動的濾波器。在通過頻帶波動的濾波器為 Chebyshev Type I Filter，在抑制頻帶上波動的濾波器為 Chebyshev Type II Filter

根據下列規格，設計 Chebyshev Type I Filters

(1) lowpass filter & highpass filter

通過頻帶邊緣頻率 200 Hz

抑制頻帶邊緣頻率 300 Hz

通過頻帶波紋 0.5 dB

抑制頻帶波紋 50 dB

取樣頻率 1000Hz

(2) bandpass & bandstop filter

通過頻帶邊緣頻率 200, 300 Hz

抑制頻帶邊緣頻率 100, 400 Hz

通過頻帶波紋 0.5 dB

抑制頻帶波紋 50 dB

取樣頻率 1000Hz

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt

def plot(a, b, plotname, plotpos1, plotpos2):
    w, H = signal.freqz( b, a )
    magnitude = abs( H )
    phase = np.angle( H )

    plt.subplot(plotpos1)
    plt.plot( w, magnitude )
    plt.xlabel( r'$\omega$ 強度頻譜('+plotname+')' )
    plt.ylabel( 'Magnitude' )

    plt.subplot(plotpos2)
    plt.plot( w, phase )
    plt.xlabel( r'$\omega$ 相位頻譜('+plotname+')' )
    plt.ylabel( 'Phase' )

# 替換中文字型，處理 matplotlib 中文 label 問題
# https://blog.csdn.net/gmr2453929471/article/details/78655834
plt.rcParams['font.sans-serif'] = ['Heiti TC'] # 步驟一（替換sans-serif字型）
plt.rcParams['axes.unicode_minus'] = False  # 步驟二（解決座標軸負數的負號顯示問題）

# lowpass filter
# passband edge frequency(Hz)
fp = 200
# stopband edge frequency(Hz)
fs = 300
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp = 2 * fp / Fs
ws = 2 * fs / Fs

n, wn = signal.cheb1ord( wp, ws, rp, rs )
b, a = signal.cheby1( n, rp, wn, 'lowpass' )

plt.figure( 1 )
plot(a, b, 'lowpass', '221', '222')


# highpass filter
# passband edge frequency(Hz)
fp = 200
# stopband edge frequency(Hz)
fs = 300
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp = 2 * fp / Fs
ws = 2 * fs / Fs

n, wn = signal.cheb1ord( wp, ws, rp, rs )
b, a = signal.cheby1( n, rp, wn, 'highpass' )

plot(a, b, 'highpass', '223', '224')


# bandpass filter
# 1st passband edge frequency(Hz)
fp1 = 200
# 2nd passband edge frequency(Hz)
fp2 = 300
# 1st stopband edge frequency(Hz)
fs1 = 100
# 2nd stopband edge frequency(Hz)
fs2 = 400
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp1 = 2 * fp1 / Fs
wp2 = 2 * fp2 / Fs
ws1 = 2 * fs1 / Fs
ws2 = 2 * fs2 / Fs

n, wn = signal.cheb1ord( [ wp1, wp2 ], [ ws1, ws2 ], rp, rs )
b, a = signal.cheby1( n, rp, wn, 'bandpass' )

plt.figure( 2 )
plot(a, b, 'bandpass', '221', '222')


# bandstop filter
# 1st passband edge frequency(Hz)
fp1 = 200
# 2nd passband edge frequency(Hz)
fp2 = 300
# 1st stopband edge frequency(Hz)
fs1 = 100
# 2nd stopband edge frequency(Hz)
fs2 = 400
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp1 = 2 * fp1 / Fs
wp2 = 2 * fp2 / Fs
ws1 = 2 * fs1 / Fs
ws2 = 2 * fs2 / Fs

n, wn = signal.cheb1ord( [ wp1, wp2 ], [ ws1, ws2 ], rp, rs )
b, a = signal.cheby1( n, rp, wn, 'bandstop' )

plot(a, b, 'bandstop', '223', '224')

plt.show( )

Chebyshev Type II Filter

在抑制頻帶上波動的濾波器為 Chebyshev Type II Filter

根據下列規格，設計 Chebyshev Type II Filters

(1) lowpass filter & highpass filter

通過頻帶邊緣頻率 200 Hz

抑制頻帶邊緣頻率 300 Hz

通過頻帶波紋 0.5 dB

抑制頻帶波紋 50 dB

取樣頻率 1000Hz

(2) bandpass & bandstop filter

通過頻帶邊緣頻率 200, 300 Hz

抑制頻帶邊緣頻率 100, 400 Hz

通過頻帶波紋 0.5 dB

抑制頻帶波紋 50 dB

取樣頻率 1000Hz

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt

def plot(a, b, plotname, plotpos1, plotpos2):
    w, H = signal.freqz( b, a )
    magnitude = abs( H )
    phase = np.angle( H )

    plt.subplot(plotpos1)
    plt.plot( w, magnitude )
    plt.xlabel( r'$\omega$ 強度頻譜('+plotname+')' )
    plt.ylabel( 'Magnitude' )

    plt.subplot(plotpos2)
    plt.plot( w, phase )
    plt.xlabel( r'$\omega$ 相位頻譜('+plotname+')' )
    plt.ylabel( 'Phase' )

# 替換中文字型，處理 matplotlib 中文 label 問題
# https://blog.csdn.net/gmr2453929471/article/details/78655834
plt.rcParams['font.sans-serif'] = ['Heiti TC'] # 步驟一（替換sans-serif字型）
plt.rcParams['axes.unicode_minus'] = False  # 步驟二（解決座標軸負數的負號顯示問題）

# lowpass filter
# passband edge frequency(Hz)
fp = 200
# stopband edge frequency(Hz)
fs = 300
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp = 2 * fp / Fs
ws = 2 * fs / Fs

n, wn = signal.cheb2ord( wp, ws, rp, rs )
b, a = signal.cheby2( n, rp, wn, 'lowpass' )

plt.figure( 1 )
plot(a, b, 'lowpass', '221', '222')


# highpass filter
# passband edge frequency(Hz)
fp = 200
# stopband edge frequency(Hz)
fs = 300
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp = 2 * fp / Fs
ws = 2 * fs / Fs

n, wn = signal.cheb2ord( wp, ws, rp, rs )
b, a = signal.cheby2( n, rp, wn, 'highpass' )

plot(a, b, 'highpass', '223', '224')


# bandpass filter
# 1st passband edge frequency(Hz)
fp1 = 200
# 2nd passband edge frequency(Hz)
fp2 = 300
# 1st stopband edge frequency(Hz)
fs1 = 100
# 2nd stopband edge frequency(Hz)
fs2 = 400
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp1 = 2 * fp1 / Fs
wp2 = 2 * fp2 / Fs
ws1 = 2 * fs1 / Fs
ws2 = 2 * fs2 / Fs

n, wn = signal.cheb2ord( [ wp1, wp2 ], [ ws1, ws2 ], rp, rs )
b, a = signal.cheby2( n, rp, wn, 'bandpass' )

plt.figure( 2 )
plot(a, b, 'bandpass', '221', '222')


# bandstop filter
# 1st passband edge frequency(Hz)
fp1 = 200
# 2nd passband edge frequency(Hz)
fp2 = 300
# 1st stopband edge frequency(Hz)
fs1 = 100
# 2nd stopband edge frequency(Hz)
fs2 = 400
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp1 = 2 * fp1 / Fs
wp2 = 2 * fp2 / Fs
ws1 = 2 * fs1 / Fs
ws2 = 2 * fs2 / Fs

n, wn = signal.cheb2ord( [ wp1, wp2 ], [ ws1, ws2 ], rp, rs )
b, a = signal.cheby2( n, rp, wn, 'bandstop' )

plot(a, b, 'bandstop', '223', '224')

plt.show( )

Elliptic Filter

是一種在通過頻帶與抑制頻帶波紋的濾波器。橢圓濾波器比其他類型的濾波器，在階數相同的條件下，有比較小的通過頻帶與抑制頻帶波動。

根據下列規格，設計 Elliptic Filters

(1) lowpass filter & highpass filter

通過頻帶邊緣頻率 200 Hz

抑制頻帶邊緣頻率 300 Hz

通過頻帶波紋 0.5 dB

抑制頻帶波紋 50 dB

取樣頻率 1000Hz

(2) bandpass & bandstop filter

通過頻帶邊緣頻率 200, 300 Hz

抑制頻帶邊緣頻率 100, 400 Hz

通過頻帶波紋 0.5 dB

抑制頻帶波紋 50 dB

取樣頻率 1000Hz

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt

def plot(a, b, plotname, plotpos1, plotpos2):
    w, H = signal.freqz( b, a )
    magnitude = abs( H )
    phase = np.angle( H )

    plt.subplot(plotpos1)
    plt.plot( w, magnitude )
    plt.xlabel( r'$\omega$ 強度頻譜('+plotname+')' )
    plt.ylabel( 'Magnitude' )

    plt.subplot(plotpos2)
    plt.plot( w, phase )
    plt.xlabel( r'$\omega$ 相位頻譜('+plotname+')' )
    plt.ylabel( 'Phase' )

# 替換中文字型，處理 matplotlib 中文 label 問題
# https://blog.csdn.net/gmr2453929471/article/details/78655834
plt.rcParams['font.sans-serif'] = ['Heiti TC'] # 步驟一（替換sans-serif字型）
plt.rcParams['axes.unicode_minus'] = False  # 步驟二（解決座標軸負數的負號顯示問題）

# lowpass filter
# passband edge frequency(Hz)
fp = 200
# stopband edge frequency(Hz)
fs = 300
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp = 2 * fp / Fs
ws = 2 * fs / Fs

n, wn = signal.ellipord( wp, ws, rp, rs )
b, a = signal.ellip( n, rp, rs, wn, 'lowpass' )

plt.figure( 1 )
plot(a, b, 'lowpass', '221', '222')


# highpass filter
# passband edge frequency(Hz)
fp = 200
# stopband edge frequency(Hz)
fs = 300
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp = 2 * fp / Fs
ws = 2 * fs / Fs

n, wn = signal.ellipord( wp, ws, rp, rs )
b, a = signal.ellip( n, rp, rs, wn, 'highpass' )

plot(a, b, 'highpass', '223', '224')


# bandpass filter
# 1st passband edge frequency(Hz)
fp1 = 200
# 2nd passband edge frequency(Hz)
fp2 = 300
# 1st stopband edge frequency(Hz)
fs1 = 100
# 2nd stopband edge frequency(Hz)
fs2 = 400
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp1 = 2 * fp1 / Fs
wp2 = 2 * fp2 / Fs
ws1 = 2 * fs1 / Fs
ws2 = 2 * fs2 / Fs

n, wn = signal.ellipord( [ wp1, wp2 ], [ ws1, ws2 ], rp, rs )
b, a = signal.ellip( n, rp, rs, wn, 'bandpass' )

plt.figure( 2 )
plot(a, b, 'bandpass', '221', '222')


# bandstop filter
# 1st passband edge frequency(Hz)
fp1 = 200
# 2nd passband edge frequency(Hz)
fp2 = 300
# 1st stopband edge frequency(Hz)
fs1 = 100
# 2nd stopband edge frequency(Hz)
fs2 = 400
# passband ripple(dB)
rp = 0.5
# stopband ripple(dB)
rs = 50
# sampling frequency
Fs = 1000

wp1 = 2 * fp1 / Fs
wp2 = 2 * fp2 / Fs
ws1 = 2 * fs1 / Fs
ws2 = 2 * fs2 / Fs

n, wn = signal.ellipord( [ wp1, wp2 ], [ ws1, ws2 ], rp, rs )
b, a = signal.ellip( n, rp, rs, wn, 'bandstop' )

plot(a, b, 'bandstop', '223', '224')

plt.show( )

References

數位訊號處理：Python程式實作（附範例光碟）（第二版）

張貼者： 於
標籤：

沒有留言:

張貼留言

訂閱： 張貼留言 (Atom)