2021/1/25

IIR 濾波器

無限脈衝響應 (Infinite Impulse Response) IIR filter 是線性時間不變性 LTI 系統,而且是因果系統,適合用在即時性 DSP 應用

IIR Filter

通常因果的 LTI 系統可用常係數差異方程式(Constant Coefficient Difference Equation)定義

常係數差異方程式(Constant Coefficient Difference Equation) 定義:

\(\sum_{k=0}^{N}a_ky[n-k] = \sum_{k=0}^{M}b_kx[n-k]\) ,其中 \(\{a_k\}, k=1,2,...,N\) \(\{b_k\}, k=0,1,...M\) 都是常係數

IIR Filter 定義:

如果輸入訊號 x[n], 輸出訊號 y[n],IIR 定義為:

\(y[n]=-\sum_{k=1}^{N}a_ky[n-k]+\sum_{k=0}^{M}b_kx[n-k]\) ,其中 \(\{a_k\}, k=0,1,...,N\) \(\{b_k\}, k=0,1,...M\) 都是IIR 的常係數

根據定義,IIR 也可以表示為

\(y[n]=-a_1y[n-1]-a_2y[n-2]-...-a_Ny[n-N] +b_0x[n]+b_1x[n-1]+...+b_Mx[n-M]\)

因此 IIR 在產生出輸出訊號時,同時牽涉 x[n] 及 yn,所以 IIR 是一種回饋系統 Feedback System。IIR 也常被稱為 Recursive Filter

IIR 的係數包含 回饋Feedback 係數 \(\{a_k\}\) 及前饋 Feed-Forward 係數 \(\{b_k\}\),共需要 M+N+1 個係數。如果將 回饋Feedback 係數 \(\{a_k\}\) 設定為 0,則 IIR 就變成一個 FIR Filter

在 FIR filter 定義中,M 稱為階數 Order,但在 IIR filter 中,M與N都是計算量,為避免混淆,改用 N 代表 IIR filter 的階數Order

根據差異方程式的定義:\(\sum_{k=0}^{N}a_ky[n-k] = \sum_{k=0}^{M}b_kx[n-k]\) ,若取 z 轉換,則

\(Z\{\sum_{k=0}^{N}a_ky[n-k]\} = Z\{\sum_{k=0}^{M}b_kx[n-k]\}\) ,因此

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

可得到系統的轉換程式 Transfer Function 為

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

這是有理式函數(分子分母都是多項式)的表示方法 (page 10-8)

通常假設 \(a_0=1\) ,也就是 IIR Filter。

ex: IIR Filter 為 \(y[n]=0.8 y[n-1]+x[n]\) ,其中濾波器的係數 \(\{a_k\} = \{1, -0.8\}\),\(\{b_k\}=\{1\}\) ,稱為一階 First-Order IIR Filter。求此系統的轉換函式

因為 IIR 可表示為 \(y[n]-0.8y[n-1]=x[n]\) 取z轉換

\(Z\{y[n]-0.8y[n-1]\} = Z\{x[n]\}\)

可得到 \((1-0.8z^{-1})Y(z) = X(z)\),因此系統的轉換函式為:

\(H(z) = \frac{Y(z)}{X(z)} = \frac{1}{1-0.8z^{-1}}\)


IIR Filter 的系統方塊圖為


假設有個輸入的數位訊號:

\(x=\{1,2,1,-1,-2,-1\}, n=0,1,...,5\) ,其中 \(x[0]=1, x[1]=2\)

增加兩個初始靜止條件 Initial Rest Conditions

  1. 在某個時間點 \(n_0\) 以前,輸入訊號都為 0。也就是輸入訊號是在 \(n_0\) 開始發生
  2. 在某個時間點 \(n_0\) 以前,輸出訊號都為 0。也就是系統一開始為靜止狀態

計算輸出訊號:

\(y[0]=0.8y[-1]+x[0] = 1\)

\(y[1]=0.8y[0]+x[1] = 2.8\)

\(y[2]=0.8y[1]+x[2] = 0.8 \cdot 2.8 +1 = 3.24\)

\(y[3]=0.8y[2]+x[3] = 0.8 \cdot 3.24 + (-1) = 1.592\)

\(y[4]=0.8y[3]+x[4] = 0.8 \cdot 1.592 + (-2) = -0.7264\)

接下來雖然輸入訊號為0 但還是一直有產生輸出訊號,因此才稱為 "Infinite" Impulse Response

剛剛已經算出此 IIR 系統的轉換函式為:

\(H(z) = \frac{Y(z)}{X(z)} = \frac{1}{1-0.8z^{-1}}\)

可用轉換函式透過 SciPy 的 lfilter 進行計算

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

# n = np.array( [ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 ] )
n = np.arange(20)
# 輸入訊號,後面補上幾個 0
x = np.array( [ 1, 2, 1, -1, -2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ] )

# 轉換函式
b = np.array( [ 1 ] )
a = np.array( [ 1, -0.8 ] )
y = signal.lfilter( b, a, x )

print( "x =", x )
print( "y =", y )

# plt.figure( 1 )
plt.subplot(121)
plt.stem( n, x, use_line_collection=True)
plt.xlabel( 'n' )
plt.ylabel( 'x[n]' )

plt.subplot(122)
plt.stem( n, y, use_line_collection=True )
plt.xlabel( 'n' )
plt.ylabel( 'y[n]' )

plt.show( )

脈衝響應

若輸入訊號為單位脈衝 Unit Impulse \(x[n] = δ[n]\) ,代入 IIR Filter 求脈衝響應(輸出訊號 \(h[n]\))

因為 IIR Filter 定義為 \(y[n]=a_1y[n-1]+b_0x[n]\)

\(h[n]=a_1h[n-1]+b_0δ[n]\)

依序計算可得

\(h[0]=a_1h[-1]+b_0δ[0] = b_0\)

\(h[1]=a_1h[0]+b_0δ[1]=a_1b_0\)

\(h[2]=a_1h[1]+b_0δ[2]=a_1(a_1b_0)=a_1^2b_0\)

一般式:

\( h[n] = \left\{ \begin{array}{ll} b_0(a_1)^n, & \mbox{if n≥0} \\ 0, & \mbox{if n<0} \end{array} \right.\)


若輸入訊號為單位步階訊號

\( 𝜇[n] = \left\{ \begin{array}{ll} 1, & \mbox{if n≥0} \\ 0, & \mbox{if n<0} \end{array} \right.\)

可進一步表示成

\(h[n]=b_0(a_1)^n𝜇[n]\)


前面的範例中 IIR 定義為 \(y[n]=0.8 y[n-1]+x[n]\) ,其中濾波器的係數 \(\{a_k\} = \{1, -0.8\}\),\(\{b_k\}=\{1\}\)

因此脈衝響應為

\(h[n]=(0.8)^n𝜇[n]\)


也可以透過系統轉換函式的反x轉換求得 IIR 的脈衝響應

ex: IIR 定義為 \(y[n]=0.8 y[n-1]+x[n]\)

IIR 系統的轉換函式為:

\(H(z) = \frac{Y(z)}{X(z)} = \frac{1}{1-0.8z^{-1}}\)

根據反 z 轉換,則查表可得

\(h[n]=Z^{-1}\{\frac{1}{1-0.8z^{-1}}\} = (0.8)^n𝜇[n]\)

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

# n = np.array( [ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 ] )
# x = np.array( [ 1, 0, 0, 0, 0, 0, 0, 0, 0, 0 ] )

n = np.arange(20)
# 輸入訊號,後面補上幾個 0
x = np.zeros(20)
x[0] = 1

b = np.array( [ 1 ] )
a = np.array( [ 1, -0.8 ] )
y = signal.lfilter( b, a, x )

print( "x =", x )
print( "y =", y )

# plt.figure( 1 )
plt.subplot(121)
plt.stem( n, x, use_line_collection=True )
plt.xlabel( 'n' )
plt.ylabel( 'x[n]' )

# plt.figure( 2 )
plt.subplot(122)
plt.stem( n, y, use_line_collection=True)
plt.xlabel( 'n' )
plt.ylabel( 'y[n]' )

plt.show( )

步階響應

若輸入訊號為單位步階 Unit Step \(x[n] = 𝜇[n]\) ,代入 IIR Filter 求步階響應(輸出訊號 \(h[n]\))

因為 IIR Filter 定義為 \(y[n]=a_1y[n-1]+b_0x[n]\)

單位步階 Unit Step定義為

\( 𝜇[n] = \left\{ \begin{array}{ll} 1, & \mbox{if n≥0} \\ 0, & \mbox{if n<0} \end{array} \right.\)

可依序求得

\(y[0]=a_1y[-1]+b_0𝜇[0] = b_0\)

\(y[1]=a_1y[0]+b_0𝜇[1] = a_1b_0+b_0 = b_0(1+a_1)\)

\(y[2]=a_1y[1]+b_0𝜇[2] = a_1(b_0(1+a_1))+b_0 = b_0(1+a_1+a_1^2)\)

一般式:

\(y[n]=b_0(1+a_1+a_1^2+...++a_1^n) = b_0\sum_{k=0}^{n}a_1^k\)

根據等比級數公式可得

\(y[n]=b_0\sum_{k=0}^{n}a_1^k = b_0 \frac{1-a_1^{n+1}}{1-a_1}, n≥0, a_1≠1\)

  1. 當 \(|a_1|>1\),會使得 y[n] 數值指數成長,這是不穩定系統

  2. 當 \(|a_1|<1\),則 \(a_1^{n+1}\) 會逐漸衰減為 0,形成穩定系統

  3. 當 \(|a_1| = 1\)

    1. 當 \(a_1=1\),輸出訊號為 \(y[n]=(n+1)b_0\),y[n]的數值呈現線性成長,是不穩定系統

    2. 當 \(a_1=-1\),輸出訊號為

      \( y[n] = \left\{ \begin{array}{ll} b_0, & \mbox{if n 是偶數} \\ 0, & \mbox{if n 是奇數} \end{array} \right.\)


ex: 計算 IIR 定義為 \(y[n]=0.8 y[n-1]+x[n]\) 的步階響應

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

# n = np.array( [ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 ] )
# x = np.ones( 10 )
n = np.arange(20)
x = np.ones(20)

b = np.array( [ 1 ] )
a = np.array( [ 1, -0.8 ] )
y = signal.lfilter( b, a, x )

print( "x =", x )
print( "y =", y )

# plt.figure( 1 )
plt.subplot(121)
plt.stem( n, x, use_line_collection=True )
plt.xlabel( 'n' )
plt.ylabel( 'x[n]' )

# plt.figure( 2 )
plt.subplot(122)
plt.stem( n, y, use_line_collection=True )
plt.xlabel( 'n' )
plt.ylabel( 'y[n]' )

plt.show( )

IIR Filter 的應用

IIR 是一種回饋系統,在訊號處理系統、自動控制系統中很常見

迴音系統

Echo System

講話者所發出的聲音,在傳遞過程中遇到反射,會傳時會產生時間延遲現象。時間延遲通常跟距離成正比,可利用 IIR 模擬出迴音效果

假設輸入訊號為指數衰減的淡出弦波

\(x(t)=A e^{-1}sin(2 \pi ft)\)

若振幅 A = 1,頻率 f = 5Hz,產生時間長度 1s 的淡出弦波,取樣頻率定為 200Hz,預計產生時間長度 5s 的輸出訊號。因此在輸入訊號後面補上 0 ,讓樣本數為 500

套用 IIR filter,目標是每隔 1s 產生一次回饋訊號,振幅逐漸變小,並與原始輸入訊號進行疊加,得到迴音輸出訊號。設定產生迴音次數為 5。

IIR 可設計為:

\(y[n]=-0.8y[n-100]-0.6y[n-200]-0.4y[n-300]-0.2y[n-400]+x[n]\)

轉換函式為

\(H(z)=\frac{1}{1+0.8z^{-100}+0.6z^{-200}+...+0.2z^{-400}}\)

所以 \(a_0=1, a_{100}=0.8, ..., b_0=1\)

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

fs = 200                                        # 取樣率
t = np.linspace( 0, 1, fs, endpoint = False )   # 定義時間陣列
x = np.exp( -t ) * np.sin( 2 * np.pi * 5 * t )  # 產生淡出弦波
x = np.pad( x, ( 0, fs * 4 ), 'constant' )      # 補零

b = np.array( [ 1 ] )                           # 定義b陣列

num_echos = 5                                   # 迴音次數
a = np.zeros( fs * num_echos )                  # 定義a陣列
for i in range( num_echos ):
    a[i * fs] = 1 - i / num_echos

y = signal.lfilter( x, b, a )                   # IIR濾波器

# plt.figure( 1 )                                   # 繪圖
plt.subplot(121)
plt.plot( x )
plt.xlabel( 'n' )
plt.ylabel( 'Amplitude' )

# plt.figure( 2 )
plt.subplot(122)
plt.plot( y )
plt.xlabel( 'n' )
plt.ylabel( 'Amplitude' )

plt.show( )

將剛剛的迴音寫入 wav

import numpy as np
import wave
import struct
import scipy.signal as signal

file = "sinusoid_echo.wav"  # 檔案名稱

amplitude = 20000           # 振幅
frequency = 200             # 頻率(Hz)
duration = 5                # 時間長度(秒)
fs = 44100                  # 取樣頻率(Hz)
num_samples = duration * fs # 樣本數

num_channels = 1            # 通道數
sampwidth = 2               # 樣本寬度
num_frames = num_samples    # 音框數 = 樣本數
comptype = "NONE"           # 壓縮型態
compname = "not compressed" # 無壓縮

t = np.linspace( 0, 1, fs, endpoint = False )                           # 定義時間陣列
x = np.exp( -t ) * amplitude * np.sin( 2 * np.pi * frequency * t )      # 產生淡出弦波
x = np.pad( x, ( 0, 4 * fs ), 'constant' )                              # 補零

b = np.array( [ 1 ] )           # 定義b陣列

a = np.zeros( duration * fs )   # 定義a陣列
num_echos = 5
for i in range( num_echos ):
    a[ int( i * fs * 5 / num_echos ) ] = 1 - i / num_echos

y = signal.lfilter( x, b, a )
# 避免資料溢位
y = np.clip( y, -30000, 30000 )

# 寫入到 wav
wav_file = wave.open( file, 'w' )
wav_file.setparams(( num_channels, sampwidth, fs, num_frames, comptype, compname ))

for s in y :
   wav_file.writeframes( struct.pack( 'h', int ( s ) ) )

wav_file.close( )

References

數位訊號處理:Python程式實作(附範例光碟)(第二版)

沒有留言:

張貼留言