7. 濾波器設計

#Signal_Processing #DSP_Class

Posted by 高肥 on 2022-08-06

7. 濾波器設計

7.1 濾波器分類

濾波器依照濾波效果分類

  • 低通濾波器 (Lowpass Filter)

    • 使得低頻範圍的訊號通過，抑制高頻範圍的訊號，頻率的閥值稱為截止頻率(Cutoff Frequency)，定義為 $f_c$，對應的截止角頻率 $\omega_c$
    • 示意圖

  • 高通濾波器 (Highpass Filter)

    • 低通濾波器的相反，讓高頻範圍訊號通過，抑制低頻訊號
    • 示意圖

  • 帶通濾波器 (Bandpass Filter)

    • 讓頻率介於 $f_1 \sim f_2$ 或角頻率 $\omega_1 \sim \omega_2$ 之間的訊號通過，抑制其他頻率範圍的訊號
    • 示意圖

  • 帶阻濾波器 (Bandstop Filter)

    • 帶通濾波器的相反，抑制頻率範圍落在 $f_1 \sim f_2$ 或角頻率 $\omega_1 \sim \omega_2$ 之間的訊號，使得其他頻率範圍的訊號通過
    • 示意圖

  • 陷波濾波器 (Notch Filter)

    • 是帶阻濾波器的一種，僅抑制某個特定頻率，ex: 60Hz 的干擾雜訊。陷波濾波器主要是使得某特定頻率的訊號變為 0，因此也稱為歸零濾波器 Nulling Filter
  • 全通濾波器 (All-Pass Filter)
    • 讓所有頻率範圍都通過，通常會調整訊號在某些頻率範圍的相位移

7.2 頻域濾波器

到目前為止所討論到的訊號處理都是在時域做訊號處理，本小節將說明使用頻域訊號處理的手法做濾波器的設計。

處理步驟如下

  1. 輸入訊號 $x\lbrack n\rbrack$ 取傅立葉轉換，變成頻域訊號 $X(\omega)=\mathcal F{x\lbrack n\rbrack}$
  2. 定義頻域濾波器 $H(\omega)$
  3. 套用濾波器，根據卷積定理，頻率濾波為 $Y(\omega)=H(\omega)\cdot X(\omega)$，其中運算操作為點對點乘法
  4. 透過反傅立葉轉換得到濾波後的訊號：$y\lbrack n\rbrack=\mathcal F^{-1}{Y(\omega)}$

範例

若輸入訊號為：
$x(t)=\cos\left(2\pi\cdot10\cdot t\right)+\cos\left(2\pi\cdot20\cdot t\right)+\cos\left(2\pi\cdot30\cdot t\right)$
取樣頻濾為 $f_s=500Hz$
設計一個低通濾波器 (Lowpass filter)，
且截止頻率 $f_c=15Hz$
並顯示輸出其濾波結果

解答

  • 程式碼

      import numpy as np
  from numpy.fft import fft, fftshift, ifft, fftfreq
  import matplotlib.pyplot as plt

  def ideal_lowpass_filtering(x, cutoff, fs):
      # x訊號做傅立葉轉換
      X = fft(x)

      # 生成頻域濾波器
      # 在正負Cutoff Frequency的範圍內為1
      # 其餘頻帶為0
      H = np.zeros(fs)
      for i in range(-cutoff, cutoff + 1):
          H[i] = 1

      # 根據卷積定理，執行點對點乘法
      Y = H * X

      # 濾波後結果做反傅立葉轉換
      # 並取實部作為輸出
      y = ifft(Y)
      y = y.real
      return y

  def main():

      # 取樣頻率
      fs = 500
      # 濾波器截止頻率
      fc = 15
      t = np.linspace(0, 1, fs, endpoint=False)
      # 生成訊號點
      x = np.cos(2 * np.pi * 10 * t) + np.cos(2 * np.pi * 20 * t) + np.cos(2 * np.pi * 30 * t)

      # 執行低通濾波器
      y = ideal_lowpass_filtering(x, fc, fs)

      # 產生濾波結果訊號的頻譜圖
      f = fftshift(fftfreq(fs, 1 / fs))
      Xm = abs(fftshift(fft(x)))
      Ym = abs(fftshift(fft(y)))

      plt.figure(1)
      plt.plot(x)
      plt.xlabel('t (second)')
      plt.ylabel('Amplitude')

      plt.figure(2)
      plt.plot(f, Xm)
      plt.xlabel('f')
      plt.ylabel('Magnitude')

      plt.figure(3)
      plt.plot(y)
      plt.xlabel('t (second)')
      plt.ylabel('Amplitude')

      plt.figure(4)
      plt.plot(f, Ym)
      plt.xlabel('f')
      plt.ylabel('Magnitude')

      plt.show()

  main()

輸出結果

  • 原始訊號

  • 原始訊號頻譜

  • 濾波後訊號

  • 濾波後訊號頻譜

7.3 時域濾波器

上述的頻域濾波器明顯是理想濾波器，但是在時域做濾波器是無法做到 "理想” 濾波器的。

7.3.1 時域理想濾波器無法實現原因說明

現在考慮一個理想低通濾波器 (Ideal Lowpass Filter)
在頻域的定義為：


試推導該濾波器的時域脈衝響應

透過反傅立葉轉換，推導理想低通濾波器的時域脈衝響應


當 $n=0$ 時


理想低通濾波器的時域脈衝響應為：

  • 理想低通濾波器脈衝響應模擬程式

      import numpy as np
  import matplotlib.pyplot as plt

  filter_size = 501
  filter_half = int( filter_size / 2 )
  # 定義 Cutoff Frequency
  wc = np.pi / 16

  na = np.arange( -filter_half, filter_half + 1 ) # 定義 n 陣列 
  h = np.zeros( filter_size )                     # 計算脈衝響應
  for n in na:
      if n == 0: 
          h[n+filter_half] = 0.0625
      else:
          h[n+filter_half] = np.sin( wc * n ) / ( np.pi * n )

  # plt.stem( na, h )
  plt.plot( na, h )
  plt.xlabel( 'n' )
  plt.ylabel( 'h[n]' )

  plt.show( )

  • 理想低通濾波器脈衝響應圖

討論

由上述的理想低通濾波器的脈衝響應可以歸納出下列兩點，使得理想低通濾波器無法在時域實現

  • 在時域要應用 $h_{LP}\lbrack n\rbrack$ 做卷積時需要的資料量是 $-\infty<n<\infty$ ，
    需要有無限多的輸入點，這在真實世界是不可能發生的
  • 因為 $-\infty<n<\infty$ ，代表所需要的資料的時間點有未來的資料，
    這就不存在因果關係 (Causality)，因此在真實世界是不可實現的

7.3.2 時域濾波器的設計方法列表

在上面的章節說明了頻域的理想濾波器，轉換到時域時就無法實現。

但是聰明的數學家們想了一些辦法讓時域濾波器還是有辦法被設計出來，只是是近似法，無法做到一個真正的理想濾波器。

  • 濾波器設計方法

在時域設計濾波器的方法主要有下列幾種

  1. Window design method (FIR)
  2. Frequency sampling method (FIR)
  3. Least MSE (mean square error) method (FIR)
  4. Parks–McClellan method (FIR)
  5. Impulse Invariance (IIR)
  6. Step Invariance (IIR)
  7. Bilinear Transform (IIR)

上面列出的各種濾波器設計方法其實我也沒有真的實現過。

畢竟在實務應用上一定是透過Python或是Matlab等等的軟體程式工具做設計，而這些軟體背後的方法使用的就是上述的演算法。

在我們的實際應用面主要還是了解濾波器的頻率響應圖如何看，然後從各種的濾波器中選取符合使用的濾波器。

下面的章節主要講解 FIR 與 IIR 兩種濾波器有甚麼定義與差別，並且說明他們在Python上面要怎麼實作，並給出一些範例。

7.3.3 FIR濾波器

假設輸入訊號為 $x\lbrack n\rbrack$，輸出訊號為 $y\lbrack n\rbrack$
則 有限脈衝響應 (Finite Impulse Response, FIR)可以定義為
$y\lbrack n\rbrack=\sum_{k=0}^Mb_k\cdot x\lbrack n-k\rbrack$
其中，$b_k,\;k=0,\;1,\;\dots,\;M$ 稱為FIR濾波器的係數 (Coefficients)
$M$ 稱為濾波器的階數 (Order)

由上面的濾波器數學式定義，可以看出濾波器的輸出只與輸入有關。
因此，當輸入式有限的，輸出結果也一定是有限的，故而得名。

濾波器設計參數 — 窗函數 (主要用在FIR濾波器設計上)

由於是理想濾波器，在時域的脈衝響應是是無限長序列。但是我們所要設計的 FIR 濾波器，其 $h\lbrack n\rbrack$是有限長的。為了能用 FIR 濾波器近似理想濾波器，需透過窗函數將理想濾波器的無限長單位脈衝響應分別從左右進行截斷。

實務上所謂的窗函數截斷，其實就是將窗函數與脈衝響應做點對點乘法，就可以將無限長度的脈衝響應截斷。而窗函數的設計方式有很多種，下面列出一些找到有說明其特點的窗函數

  • 矩形窗 (Rectangle or Boxcar)

    • 數學式
      $
      \omega\lbrack k\rbrack =
      \begin{cases}
      1, & 0\leq k\leq M \\
      0, & otherwise
      \end{cases}
      $

    • 特點

      矩形窗使用最多，習慣上不加窗就是使信號通過了矩形窗。

      這種窗的優點是主瓣比較集中，缺點是旁瓣較高，並有負旁瓣，導致變換中帶進了高頻干擾和洩漏，甚至出現負譜現象。頻率識別精度最高，幅值識別精度最低，所以矩形窗不是一個理想的窗。

    • 應用

      如果僅要求精確讀出主瓣頻率，而不考慮幅值精度，則可選用矩形窗，例如測量物體的自振頻率等，也可以用在階次分析中。

    • 範例程式

        from scipy import signal
  from scipy.fft import fft, fftshift
  import matplotlib.pyplot as plt
  import numpy as np

  window = signal.windows.boxcar(51)
  plt.plot(window)
  plt.title("Boxcar window")
  plt.ylabel("Amplitude")
  plt.xlabel("Sample")
  plt.show()
  plt.savefig('01.png', bbox_inches='tight')

  plt.figure()
  A = fft(window, 2048) / (len(window)/2.0)
  freq = np.linspace(-0.5, 0.5, len(A))
  response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
  plt.plot(freq, response)
  plt.axis([-0.5, 0.5, -120, 0])
  plt.title("Frequency response of the boxcar window")
  plt.ylabel("Normalized magnitude [dB]")
  plt.xlabel("Normalized frequency [cycles per sample]")
  plt.show()
  plt.savefig('02.png', bbox_inches='tight')

    • 結果圖

  • 漢寧窗 (Hanning)

    • 數學式

      $\omega\lbrack k\rbrack=0.5-0.5\cdot\cos\left(\frac{2\pi n}{M-1}\right),\;\;0\leq n\leq M-1$

    • 特點

      主瓣加寬並降低，旁瓣則顯著減小，從減小洩漏觀點出發，漢寧窗優於矩形窗。但漢甯窗主瓣加寬，相當於分析頻寬加寬，頻率分辨力下降。它與矩形窗相比，洩漏、波動都減小了,並且選擇性也提高。

    • 應用

      漢甯窗是很有用的窗函數。如果測試信號有多個頻率分量，頻譜表現的十分複雜，且測試的目的更多關注頻率點而非能量的大小，需要選擇漢寧窗。如果被測信號是隨機或者未知的，選擇漢寧窗。

    • 範例程式

        from scipy import signal
  from scipy.fft import fft, fftshift
  import matplotlib.pyplot as plt
  import numpy as np

  window = signal.windows.hann(51)
  plt.plot(window)
  plt.title("Hann window")
  plt.ylabel("Amplitude")
  plt.xlabel("Sample")
  plt.savefig('01.png', bbox_inches='tight')

  plt.figure()
  A = fft(window, 2048) / (len(window)/2.0)
  freq = np.linspace(-0.5, 0.5, len(A))
  response = np.abs(fftshift(A / abs(A).max()))
  response = 20 * np.log10(np.maximum(response, 1e-10))
  plt.plot(freq, response)
  plt.axis([-0.5, 0.5, -120, 0])
  plt.title("Frequency response of the Hann window")
  plt.ylabel("Normalized magnitude [dB]")
  plt.xlabel("Normalized frequency [cycles per sample]")
  plt.savefig('02.png', bbox_inches='tight')

    • 結果圖

  • 漢明窗 (Hamming)

    • 數學式

      $\omega\lbrack k\rbrack=0.54-0.46\cdot\cos\left(\frac{2\pi n}{M-1}\right),\;\;0\leq n\leq M-1$

    • 特點

      與漢寧窗都是余弦窗，又稱改進的升余弦窗，只是加權係數不同，使旁瓣達到更小。但其旁瓣衰減速度比漢寧窗衰減速度慢。

    • 應用

      與漢明窗類似，也是很有用的窗函數。

    • 範例程式

        from scipy import signal
  from scipy.fft import fft, fftshift
  import matplotlib.pyplot as plt
  import numpy as np

  window = signal.windows.hamming(51)
  plt.plot(window)
  plt.title("Hamming window")
  plt.ylabel("Amplitude")
  plt.xlabel("Sample")
  plt.savefig('01.png', bbox_inches='tight')

  plt.figure()
  A = fft(window, 2048) / (len(window)/2.0)
  freq = np.linspace(-0.5, 0.5, len(A))
  response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
  plt.plot(freq, response)
  plt.axis([-0.5, 0.5, -120, 0])
  plt.title("Frequency response of the Hamming window")
  plt.ylabel("Normalized magnitude [dB]")
  plt.xlabel("Normalized frequency [cycles per sample]")
  plt.savefig('02.png', bbox_inches='tight')

    • 結果圖

  • 高斯窗 (Gaussian)

    • 數學式

      $\omega\lbrack k\rbrack=e^{-\frac12{(\frac n\sigma)}^2}$

    • 特點

      是一種指數窗。主瓣較寬，故而頻率分辨力低；無負的旁瓣，第一旁瓣衰減達一55dB。常被用來截短一些非週期信號，如指數衰減信號等。

    • 應用

      對於隨時間按指數衰減的函數，可採用指數窗來提高信噪比。

    • 範例程式

        from scipy import signal
  from scipy.fft import fft, fftshift
  import matplotlib.pyplot as plt
  import numpy as np

  window = signal.windows.gaussian(51, std=7)
  plt.plot(window)
  plt.title(r"Gaussian window ($\sigma$=7)")
  plt.ylabel("Amplitude")
  plt.xlabel("Sample")
  plt.savefig('01.png', bbox_inches='tight')

  plt.figure()
  A = fft(window, 2048) / (len(window)/2.0)
  freq = np.linspace(-0.5, 0.5, len(A))
  response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
  plt.plot(freq, response)
  plt.axis([-0.5, 0.5, -120, 0])
  plt.title(r"Frequency response of the Gaussian window ($\sigma$=7)")
  plt.ylabel("Normalized magnitude [dB]")
  plt.xlabel("Normalized frequency [cycles per sample]")
  plt.savefig('02.png', bbox_inches='tight')

    • 結果圖

  • 布萊克曼窗 (Blackman)

更多的Python 窗函數可參考下列連結

FIR濾波器設計 (Python)

FIR濾波器設計最直接的方法就是窗函數法 (Window Method)，主要步驟如下

  1. 選取濾波器種類

    首先選取濾波器的種類: Lowpass、Highpass、Bandpass or Bandstop

  2. 定義濾波器規格

    1. 若是Lowpass或Highpass，需定義Cutoff Frequency ($f_c$)，以 $Hz$ 為單位
    2. 若是Bandpass或Bandstop，需定義兩個Cutoff Frequency ($f_1$ 與 $f_2$)，以 $Hz$ 為單位

  3. 選取取樣頻率

    以 $Hz$ 為單位

  4. 套用窗函數

    Python中可套用的窗函數 ⇒ Window functions (scipy.signal.windows)

  5. 計算 FIR濾波器係數

    此部分透過 scipy.signal.firwin 實作

  6. 頻率響應觀察

    根據頻率響應的結果觀察所設計的濾波器性能是否符合需求

範例

使用Hamming Window設計下列濾波器

  1. 低通濾波器：$f_c=250 Hz$
  2. 高通濾波器：$f_c=250 Hz$
  3. 帶通濾波器：$f_1=200 Hz、f_2=300 Hz$
  4. 帶阻濾波器：$f_1=200 Hz、f_2=300 Hz$
    其中，取樣頻率 $f_s=1000 Hz$
    顯示其頻率響應的強度頻譜 (Magnitude Spectrum)

  • Lowpass程式碼

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

  # Cutoff Frequency
  cutoff = 250
  # Sampling Rate
  freq = 1000
  # Number of Impulse Response
  n = 31

  # 計算產生Lowpass FIR Filter係數
  h = signal.firwin( n, cutoff, window = 'hamming', pass_zero = 'lowpass', fs = freq )

  w, H = signal.freqz( h )
  w = w*freq/(2*np.pi)
  magnitude = abs( H )
  phase = np.angle( H )

  plt.figure( 1 )
  plt.plot( w, magnitude )
  plt.xlabel( r'f (Hz)' )
  plt.ylabel( 'Magnitude' )

  plt.figure( 2 )
  plt.plot( w, phase )
  plt.xlabel( r'$\omega$' )
  plt.ylabel( 'Phase' )

  plt.show( )

  • Highpass程式碼

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

  # Cutoff Frequency
  cutoff = 250
  # Sampling Rate
  freq = 1000
  # Number of Impulse Response
  n = 31

  # 計算產生Lowpass FIR Filter係數
  h = signal.firwin( n, cutoff, window = 'hamming', pass_zero = 'highpass', fs = freq )

  w, H = signal.freqz( h )
  w = w*freq/(2*np.pi)
  magnitude = abs( H )
  phase = np.angle( H )

  plt.figure( 1 )
  plt.plot( w, magnitude )
  plt.xlabel( r'f (Hz)' )
  plt.ylabel( 'Magnitude' )

  plt.figure( 2 )
  plt.plot( w, phase )
  plt.xlabel( r'$\omega$' )
  plt.ylabel( 'Phase' )

  plt.show( )

  • Bandpass程式碼

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

  # Stop Frequency
  f1 = 200
  f2 = 300
  # Sampling Rate
  freq = 1000
  # Number of Impulse Response
  n = 31

  # 計算產生Lowpass FIR Filter係數
  h = signal.firwin( n, [f1, f2], window = 'hamming', pass_zero = 'bandpass', fs = freq )

  w, H = signal.freqz( h )
  w = w*freq/(2*np.pi)
  magnitude = abs( H )
  phase = np.angle( H )

  plt.figure( 1 )
  plt.plot( w, magnitude )
  plt.xlabel( r'f (Hz)' )
  plt.ylabel( 'Magnitude' )

  plt.figure( 2 )
  plt.plot( w, phase )
  plt.xlabel( r'$\omega$' )
  plt.ylabel( 'Phase' )

  plt.show( )

  • Bandstop程式碼

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

  # Stop Frequency
  f1 = 200
  f2 = 300
  # Sampling Rate
  freq = 1000
  # Number of Impulse Response
  n = 31

  # 計算產生Lowpass FIR Filter係數
  h = signal.firwin( n, [f1, f2], window = 'hamming', pass_zero = 'bandstop', fs = freq )

  w, H = signal.freqz( h )
  w = w*freq/(2*np.pi)
  magnitude = abs( H )
  phase = np.angle( H )

  plt.figure( 1 )
  plt.plot( w, magnitude )
  plt.xlabel( r'f (Hz)' )
  plt.ylabel( 'Magnitude' )

  plt.figure( 2 )
  plt.plot( w, phase )
  plt.xlabel( r'$\omega$' )
  plt.ylabel( 'Phase' )

  plt.show( )

7.3.4 IIR濾波器

差分方程式

常係數差分方程式 (Constant Coefficient Difference Equation)可定義為：


其中，${a_k},\;k=0,\;1,\;\dots,\;N$，${b_k},\;k=0,\;1,\;\dots,\;M$ 皆為常數係數

IIR濾波器

假設輸入訊號為 $x\lbrack n\rbrack$，輸出訊號為 $y\lbrack n\rbrack$
則 無限脈衝響應 (Infinite Impulse Response, IIR)可以定義為


其中，${a_k},\;k=1,\;\dots,\;N$，${b_k},\;k=0,\;1,\;\dots,\;M$ 稱為IIR濾波器的係數
$N$ 為濾波器的階數

由IIR濾波器的數學式可觀察到，濾波器的輸出結果不僅僅與輸入有關，還與之前的輸出結果有關。
因此，當輸入中止時，之前的輸出結果還是有可能造成現在的輸出，進而造成一個類似於無止盡的輸出的效果，因此命名為 “無限" 脈衝響應。

濾波器設計參數 — 誤差容許範圍 (主要用在IIR濾波器設計上)

實際低通濾波器設計時，通常在通過頻帶或是抑制頻帶內，容許頻率響應有少許的波動(或誤差)。

$\delta_P$： 通過頻帶波紋 (Passband Ripple)

$\delta_S$： 抑制頻帶波紋 (Stopband Ripple)

$\omega_P$：通過頻帶邊緣 (Passband Edge)

$\omega_S$：抑制頻帶邊緣 (Stopband Edge)

在 $\omega_P\sim\omega_S$之間，則是容許頻率響應從1逐漸降為0，稱為過度頻帶 (Transition Band)。

IIR濾波器設計 (Python)

IIR濾波器設計主要步驟與FIR非常接近，如下

  1. 選取濾波器種類

    首先選取濾波器的種類: Lowpass、Highpass、Bandpass or Bandstop

  2. 定義濾波器規格

    1. 若是Lowpass或Highpass，需定義通過頻帶邊緣 (Passband Edge，$\omega_P$)與抑制頻帶邊緣 (Stopband Edge，$\omega_S$)，以 $Hz$ 為單位
    2. 若是Bandpass或Bandstop，需定義兩組通過頻帶邊緣 (Passband Edge，$\omega_{P1}$與$\omega_{P2}$)與抑制頻帶邊緣 (Stopband Edge，$\omega_{S1}$與$\omega_{S2}$)，，以 $Hz$ 為單位
    3. 此外，需要定義通過頻帶波紋與抑制頻帶波紋，單位為 $dB$

    4. 階數設定

      這部分目前來說是實際做實驗測試，階數多寡影響到Ripple的嚴重程度，就看實際使用能接受到多嚴重的誤差

  3. 選取取樣頻率

    以 $Hz$ 為單位

  4. 計算 IIR濾波器係數

    此部分透過 scipy.signal.butter 等function實作

  5. 頻率響應觀察

    根據頻率響應的結果觀察所設計的濾波器性能是否符合需求

IIR濾波器四大類型與實作

在DSP的書籍中講到IIR濾波器基本上都會講下列四種典型IIR濾波器

  1. Butterworth Filter
  2. Chebyshev Type I Filter
  3. Chebyshev Type II Filter
  4. Elliptic Filter

但是DSP並不是只有這四種濾波器的設計，只是發展較久、特性已經明確才被納入大部分DSP教材中

以下將透過一個範例說明如何在Python中呼叫使用上述四種IIR濾波器與其相對應的頻率響應分析

範例

設計IIR濾波器，相關規格如下

  1. Lowpass & Highpass Filter
    Passband Edge: 50 Hz
    Stopband Edge: 100 Hz
    Passband Ripple: 1 dB
    Stopband Ripple: 20 dB
  2. Bandpass & Bandstop FIlter
    Passband Edge: 100, 400 Hz
    Stopband Edge: 100, 400 Hz
    Passband Ripple: 1 dB
    Stopband Ripple: 20 dB
    通過Butterworth、Chebyshev Type I & II、Elliptic Filter等四種濾波器
    實現上述的四個濾波器的規格

結果

  • 低通濾波器 (Lowpass Fileter)

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

  # 取樣頻率
  fs = 1000
  # 濾波器階數
  N1 = 10
  # Passband Frequency (Hz，需轉換成Normalized Nyquist Frequency)
  wp = 50 / ( fs / 2)
  # Stopband Frequency (Hz，需轉換成Normalized Nyquist Frequency)
  ws = 100 / ( fs / 2)
  # Passband Ripple (dB)
  rp = 1
  # Stopband Ripple (dB)
  rs = 20

  # Butterworth Lowpass Filter
  b, a = signal.butter(N1, ws, 'lowpass')
  w, H = signal.freqz( b, a )
  magnitude_butter = abs( H )

  # Chebyshev Type I Lowpass Filter
  b, a = signal.cheby1( N1, rp, ws, 'lowpass' )
  w, H = signal.freqz( b, a )
  magnitude_ch1 = abs( H )

  # Chebyshev Type II Lowpass Filter
  b, a = signal.cheby2( N1, rs, ws, 'lowpass' )
  w, H = signal.freqz( b, a )
  magnitude_ch2 = abs( H )

  # Elliptic Lowpass Filter
  b, a = signal.ellip( N1, rp, rs, ws, 'lowpass' )
  w, H = signal.freqz( b, a )
  magnitude_ell = abs( H )

  # 計算頻率響應圖的x座標，需轉換成Hz
  w = w*fs/(2*np.pi)

  fig, ax = plt.subplots(2, 2)
  ax[0, 0].plot( w, magnitude_butter, color='#CC0000')
  ax[0, 0].set_title("Butterworth")
  ax[0, 0].set_xlabel(r'f (Hz)')
  ax[0, 0].set_ylabel('Magnitude')
  ax[0, 1].plot( w, magnitude_ch1, color='#00CC00')
  ax[0, 1].set_title("Chebyshev Type I")
  ax[0, 1].set_xlabel(r'f (Hz)')
  ax[0, 1].set_ylabel('Magnitude')
  ax[1, 0].plot( w, magnitude_ch2, color='#0000CC')
  ax[1, 0].set_title("Chebyshev Type II")
  ax[1, 0].set_xlabel(r'f (Hz)')
  ax[1, 0].set_ylabel('Magnitude')
  ax[1, 1].plot( w, magnitude_ell, color='#FFA042')
  ax[1, 1].set_title("Elliptic")
  ax[1, 1].set_xlabel(r'f (Hz)')
  ax[1, 1].set_ylabel('Magnitude')
  fig.tight_layout()
  plt.show()

  • 高通濾波器 (Highpass Fileter)

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

  # 取樣頻率
  fs = 1000
  # 濾波器階數
  N1 = 10
  # Passband Frequency (Hz，需轉換成Normalized Nyquist Frequency)
  wp = 450 / ( fs / 2)
  # Stopband Frequency (Hz，需轉換成Normalized Nyquist Frequency)
  ws = 400 / ( fs / 2)
  # Passband Ripple (dB)
  rp = 1
  # Stopband Ripple (dB)
  rs = 20

  # Butterworth Highpass Filter
  b, a = signal.butter(N1, ws, 'highpass')
  w, H = signal.freqz( b, a )
  magnitude_butter = abs( H )

  # Chebyshev Type I Highpass Filter
  b, a = signal.cheby1( N1, rp, ws, 'highpass' )
  w, H = signal.freqz( b, a )
  magnitude_ch1 = abs( H )

  # Chebyshev Type II Highpass Filter
  b, a = signal.cheby2( N1, rs, ws, 'highpass' )
  w, H = signal.freqz( b, a )
  magnitude_ch2 = abs( H )

  # Elliptic Highpass Filter
  b, a = signal.ellip( N1, rp, rs, ws, 'highpass' )
  w, H = signal.freqz( b, a )
  magnitude_ell = abs( H )

  # 計算頻率響應圖的x座標，需轉換成Hz
  w = w*fs/(2*np.pi)

  fig, ax = plt.subplots(2, 2)
  ax[0, 0].plot( w, magnitude_butter, color='#CC0000')
  ax[0, 0].set_title("Butterworth")
  ax[0, 0].set_xlabel(r'f (Hz)')
  ax[0, 0].set_ylabel('Magnitude')
  ax[0, 1].plot( w, magnitude_ch1, color='#00CC00')
  ax[0, 1].set_title("Chebyshev Type I")
  ax[0, 1].set_xlabel(r'f (Hz)')
  ax[0, 1].set_ylabel('Magnitude')
  ax[1, 0].plot( w, magnitude_ch2, color='#0000CC')
  ax[1, 0].set_title("Chebyshev Type II")
  ax[1, 0].set_xlabel(r'f (Hz)')
  ax[1, 0].set_ylabel('Magnitude')
  ax[1, 1].plot( w, magnitude_ell, color='#FFA042')
  ax[1, 1].set_title("Elliptic")
  ax[1, 1].set_xlabel(r'f (Hz)')
  ax[1, 1].set_ylabel('Magnitude')
  fig.tight_layout()
  plt.show()

  • 帶通濾波器 (Bandpass Fileter)

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

  # 取樣頻率
  fs = 1000
  # 濾波器階數
  N1 = 10
  # Passband Frequency (Hz，需轉換成Normalized Nyquist Frequency)
  wp = np.array([ 100 , 400 ]) / ( fs / 2)
  # Stopband Frequency (Hz，需轉換成Normalized Nyquist Frequency)
  ws = np.array([ 100 , 400 ]) / ( fs / 2)
  # Passband Ripple (dB)
  rp = 1
  # Stopband Ripple (dB)
  rs = 20

  # Butterworth Bandpass Filter
  b, a = signal.butter( N1, ws, 'bandpass' )
  w, H = signal.freqz( b, a )
  magnitude_butter = abs( H )

  # Chebyshev Type I Bandpass Filter
  b, a = signal.cheby1( N1, rp, ws, 'bandpass' )
  w, H = signal.freqz( b, a )
  magnitude_ch1 = abs( H )

  # Chebyshev Type II Bandpass Filter
  b, a = signal.cheby2( N1, rs, ws, 'bandpass' )
  w, H = signal.freqz( b, a )
  magnitude_ch2 = abs( H )

  # Elliptic Bandpass Filter
  b, a = signal.ellip( N1, rp, rs, ws, 'bandpass' )
  w, H = signal.freqz( b, a )
  magnitude_ell = abs( H )

  # 計算頻率響應圖的x座標，需轉換成Hz
  w = w*fs/(2*np.pi)

  fig, ax = plt.subplots(2, 2)
  ax[0, 0].plot( w, magnitude_butter, color='#CC0000')
  ax[0, 0].set_title("Butterworth")
  ax[0, 0].set_xlabel(r'f (Hz)')
  ax[0, 0].set_ylabel('Magnitude')
  ax[0, 1].plot( w, magnitude_ch1, color='#00CC00')
  ax[0, 1].set_title("Chebyshev Type I")
  ax[0, 1].set_xlabel(r'f (Hz)')
  ax[0, 1].set_ylabel('Magnitude')
  ax[1, 0].plot( w, magnitude_ch2, color='#0000CC')
  ax[1, 0].set_title("Chebyshev Type II")
  ax[1, 0].set_xlabel(r'f (Hz)')
  ax[1, 0].set_ylabel('Magnitude')
  ax[1, 1].plot( w, magnitude_ell, color='#FFA042')
  ax[1, 1].set_title("Elliptic")
  ax[1, 1].set_xlabel(r'f (Hz)')
  ax[1, 1].set_ylabel('Magnitude')
  fig.tight_layout()
  plt.show()

  • 帶阻濾波器 (Bandstop Fileter)

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

  # 取樣頻率
  fs = 1000
  # 濾波器階數
  N1 = 10
  # Passband Frequency (Hz，需轉換成Normalized Nyquist Frequency)
  wp = np.array([ 100 , 400 ]) / ( fs / 2)
  # Stopband Frequency (Hz，需轉換成Normalized Nyquist Frequency)
  ws = np.array([ 100 , 400 ]) / ( fs / 2)
  # Passband Ripple (dB)
  rp = 1
  # Stopband Ripple (dB)
  rs = 20

  # Butterworth Bandstop Filter
  b, a = signal.butter( N1, ws, 'bandstop' )
  w, H = signal.freqz( b, a )
  magnitude_butter = abs( H )

  # Chebyshev Type I Bandstop Filter
  b, a = signal.cheby1( N1, rp, ws, 'bandstop' )
  w, H = signal.freqz( b, a )
  magnitude_ch1 = abs( H )

  # Chebyshev Type II Bandstop Filter
  b, a = signal.cheby2( N1, rs, ws, 'bandstop' )
  w, H = signal.freqz( b, a )
  magnitude_ch2 = abs( H )

  # Elliptic Bandstop Filter
  b, a = signal.ellip( N1, rp, rs, ws, 'bandstop' )
  w, H = signal.freqz( b, a )
  magnitude_ell = abs( H )

  # 計算頻率響應圖的x座標，需轉換成Hz
  w = w*fs/(2*np.pi)

  fig, ax = plt.subplots(2, 2)
  ax[0, 0].plot( w, magnitude_butter, color='#CC0000')
  ax[0, 0].set_title("Butterworth")
  ax[0, 0].set_xlabel(r'f (Hz)')
  ax[0, 0].set_ylabel('Magnitude')
  ax[0, 1].plot( w, magnitude_ch1, color='#00CC00')
  ax[0, 1].set_title("Chebyshev Type I")
  ax[0, 1].set_xlabel(r'f (Hz)')
  ax[0, 1].set_ylabel('Magnitude')
  ax[1, 0].plot( w, magnitude_ch2, color='#0000CC')
  ax[1, 0].set_title("Chebyshev Type II")
  ax[1, 0].set_xlabel(r'f (Hz)')
  ax[1, 0].set_ylabel('Magnitude')
  ax[1, 1].plot( w, magnitude_ell, color='#FFA042')
  ax[1, 1].set_title("Elliptic")
  ax[1, 1].set_xlabel(r'f (Hz)')
  ax[1, 1].set_ylabel('Magnitude')
  fig.tight_layout()
  plt.show()

四大類型IIR濾波器特性

  • Butterworth Filter特點
    • Butterworth 特點是Bandpass內的頻率響應曲線最大限度平坦，沒有起伏，而在Bandstop則逐漸下降爲零。
    • Butterworth 在Bandpass內外都有平穩的幅頻特性，但有較長的Transition Band，在Transition Band上很容易造成失真
  • Chebyshev Type I Filter特點
    • Chebyshev在Transition Band比Butterworth 的衰減快，但頻率響應的幅頻特性不如Butterworth 平坦。
    • Chebyshev和理想濾波器的頻率響應曲線之間的誤差最小，但在Bandpass存在Ripple。
    • Ⅰ型在Bandpass範圍內是等波紋的，在Bandstop範圍內是單調的。
  • Chebyshev Type II Filter特點
    • Ⅱ型在Bandpass範圍內是單調的，在Bandstop範圍內是等波紋的。
  • Elliptic Filter特點
    • Elliptic的Bandpass和Bandstop都具有等波紋特性，因此通帶，阻帶逼近特性良好。
    • 對於同樣的性能要求，比前幾種濾波器所需用的階數都低
    • Elliptic的Transition Band比較窄、比較筆直

7.4 時域濾波器 v.s. 頻率濾波器

優點 缺點
時域濾波器 資料需求量少、計算量較少 無法做到理想的濾波器
頻域濾波器 頻率濾波可做到理想濾波 資料需求量大、傅立葉轉換計算量高

備註: 針對資料量與計算量討論

  • 時域濾波器其實就是一個差分方程式，可能因為階數會影響到需要保存的歷史輸入、輸出，但是資料量畢竟較低
  • 頻域濾波器的實作方式是把原始資料整段作傅立葉轉換後，刪除不需要的頻段。此種作法如果資料量太少，會造成頻譜的解析度變差，而導致誤刪除資料。相反的如果需要做到頻譜解析度夠高，需要累積很大的資料量，才可做濾波計算，就不即時。

7.5 濾波器套用

當前述的濾波器設計已經完成後，所設計出的濾波器參數透過下列兩個 function執行濾波

  • scipy.signal.filtfilt
  • scipy.signal.lfilter

選用哪個function都可以正確濾波，差別在於 "lfilter" 有相位偏移，而 "filtfilt” 沒有相位偏移。

filtfilt的作法是法原始訊號經過lfilter運作後，將濾波後的訊號左右顛倒，然後再用lfilter濾波一次，經過反向與二次濾波可以將相位偏移校正。

據一個例子說明實際應用效果

有一個隨機訊號
假設 $f_s=1000Hz$
有1000筆資料
透過一個Lowpass Filter
將雜訊點濾除，僅保留下訊號的趨勢

  • 程式碼

      from __future__ import division, print_function
  import numpy as np
  from numpy.random import randn
  from numpy.fft import rfft
  from scipy import signal
  import matplotlib.pyplot as plt

  # 濾波器設計
  b, a = signal.butter(4, 0.07, analog=False)
  # 分析濾波器頻率響應
  w, H = signal.freqz( b, a )
  magnitude_butter = abs( H )

  # 產生有雜訊的隨機訊號
  # 取樣頻率有1000 Hz
  # 共1000筆訊號
  sig = np.cumsum(randn(1000))  # Brownian noise
  # 套用濾波器參數
  # signal.filtfilt 不會有相位偏移
  sig_ff = signal.filtfilt(b, a, sig)
  # 套用濾波器參數
  # signal.lfilter 有相位偏移
  sig_lf = signal.lfilter(b, a, sig)

  plt.subplot(2, 1, 1)
  plt.plot( w*1000/(2*np.pi), magnitude_butter, label = "Butterworth Frequency Reponse" )
  plt.ylabel( 'Magnitude' )
  plt.title("Lowpass Filter")
  plt.legend()

  plt.subplot(2, 1, 2)
  plt.plot(sig, color='silver', label='Original')
  plt.plot(sig_ff, color='#3465a4', label='filtfilt')
  plt.plot(sig_lf, color='#cc0000', label='lfilter')
  plt.grid(True, which='both')
  plt.legend(loc="best")
  plt.show( )

  • 結果

Ref


#Signal_Processing #DSP_Class

高肥
軟體工程師 @NCKU
紙上得來終覺淺，絕知此事要躬行




Related Posts

Day05 state 及 props 的應用

Day05 state 及 props 的應用

 ronntsai
[day-9] DOM-選取網頁元素

[day-9] DOM-選取網頁元素

 Xiang
秘密：十週年紀念版

秘密：十週年紀念版

 JingTeng

Comments