cctg: 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 中，Ｍ與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$
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程式實作（附範例光碟）（第二版）

cctg

2021/1/25

IIR 濾波器

IIR Filter

脈衝響應

步階響應

IIR Filter 的應用

迴音系統

References

沒有留言:

張貼留言

analytics

Creative Commons License