無限脈衝響應 (Infinite Impulse Response) IIR filter 是線性時間不變性 LTI 系統，而且是因果系統，適合用在即時性 DSP 應用
IIR Filter
通常因果的 LTI 系統可用常係數差異方程式(Constant Coefficient Difference Equation)定義
常係數差異方程式(Constant Coefficient Difference Equation) 定義：
∑Nk=0aky[n−k]=∑Mk=0bkx[n−k] ，其中 {ak},k=1,2,...,N {bk},k=0,1,...M 都是常係數
IIR Filter 定義：
如果輸入訊號 x[n], 輸出訊號 y[n]，IIR 定義為：
y[n]=−∑Nk=1aky[n−k]+∑Mk=0bkx[n−k] ，其中 {ak},k=0,1,...,N {bk},k=0,1,...M 都是IIR 的常係數
根據定義，IIR 也可以表示為
y[n]=−a1y[n−1]−a2y[n−2]−...−aNy[n−N]+b0x[n]+b1x[n−1]+...+bMx[n−M]
因此 IIR 在產生出輸出訊號時，同時牽涉 x[n] 及 yn，所以 IIR 是一種回饋系統 Feedback System。IIR 也常被稱為 Recursive Filter
IIR 的係數包含 回饋Feedback 係數 {ak} 及前饋 Feed-Forward 係數 {bk}，共需要 M+N+1 個係數。如果將 回饋Feedback 係數 {ak} 設定為 0，則 IIR 就變成一個 FIR Filter
在 FIR filter 定義中，M 稱為階數 Order，但在 IIR filter 中，Ｍ與N都是計算量，為避免混淆，改用 N 代表 IIR filter 的階數Order
根據差異方程式的定義：∑Nk=0aky[n−k]=∑Mk=0bkx[n−k] ，若取 z 轉換，則
Z{∑Nk=0aky[n−k]}=Z{∑Mk=0bkx[n−k]} ，因此
(a0+a1z−1+...+aNz−N)Y(z)=(b0+b1z−1+...+bMz−M)X(z)
可得到系統的轉換程式 Transfer Function 為
H(z)=Y(z)X(z)=b0+b1z−1+...+bMz−Ma0+z1z−1+...+aNz−N
這是有理式函數(分子分母都是多項式)的表示方法 (page 10-8)
通常假設 a0=1 ，也就是 IIR Filter。
ex: IIR Filter 為 y[n]=0.8y[n−1]+x[n] ，其中濾波器的係數 {ak}={1,−0.8}，{bk}={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)=Y(z)X(z)=11−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
- 在某個時間點 n0 以前，輸入訊號都為 0。也就是輸入訊號是在 n0 開始發生
- 在某個時間點 n0 以前，輸出訊號都為 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⋅2.8+1=3.24
y[3]=0.8y[2]+x[3]=0.8⋅3.24+(−1)=1.592
y[4]=0.8y[3]+x[4]=0.8⋅1.592+(−2)=−0.7264
接下來雖然輸入訊號為0 但還是一直有產生輸出訊號，因此才稱為 "Infinite" Impulse Response
剛剛已經算出此 IIR 系統的轉換函式為：
H(z)=Y(z)X(z)=11−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( )
