無限脈衝響應 (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 中，Ｍ與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
- 在某個時間點 \(n_0\) 以前，輸入訊號都為 0。也就是輸入訊號是在 \(n_0\) 開始發生
- 在某個時間點 \(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\)
當 \(|a_1|>1\)，會使得 y[n] 數值指數成長，這是不穩定系統
當 \(|a_1|<1\)，則 \(a_1^{n+1}\) 會逐漸衰減為 0，形成穩定系統
當 \(|a_1| = 1\)
當 \(a_1=1\)，輸出訊號為 \(y[n]=(n+1)b_0\)，y[n]的數值呈現線性成長，是不穩定系統
當 \(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( )
