2021/02/22

頻率響應

DSP (ex: 音響、通訊設備) 的目的都是希望達到無失真的傳遞,在輸出端重現原始訊號,因此,理想的訊號處理系統,在不同頻率範圍,需表現均勻的強度響應。頻率響應 Frequency Response 是用來衡量系統在特定頻率範圍的操作特性。

音響設備可用頻率響應來衡量品質,ex: 理想的喇叭,在人類聽力範圍中 20Hz~20kHz,應該具備等值的頻率響應,藉以達到原音重現。

頻率響應的定義:

系統對於輸入訊號在特定頻率範圍的響應情形,可定義為頻率的函式,響應則包含強度 Magnitude 與相位角 Phase Angle 等等量化數據,通常是以頻譜的方式呈現,用來表示系統的操作特性。


典型的 DSP (LTI) 系統, x[n] 為輸入訊號,經過 h[n] 脈衝響應(impulse response, 也稱為濾波器 filter),得到 y[n] 輸出訊號

如果 DSP 牽涉卷積運算 \(y[n]=h[n]*x[n]\)

根據卷積定理 Covolution Theorem 則

\(F\{y[n]\} = F\{h[n]*x[n]\} = F\{h[n]\} \cdot F\{x[n]\}\)

或是以離散時間傅立葉轉換可表示成:

\(Y(e^{jω}) = H(e^{jω}) \cdot X(e^{jω})\)

因此系統的頻率響應為:

\(H(e^{jω}) = \frac{Y(e^{jω})}{X(e^{jω})}\)

根據 z 轉換,系統的頻率響應跟轉換函式的關係,可表示為:

\(H(e^{jω}) = H(z)|_{z=e^{jω}}\)

因此系統的頻率響應,可表示為

\(H(z) = \frac{Y(z)}{X(z)}\)

其中 X(z), Y(z) 分別為輸入訊號與輸出訊號的 z 轉換結果


DSP 系統的頻率響應 定義為:

\(H(e^{jω}) = \sum_{n=-∞}^{∞} h[n]e^{-jωn}\)

在 LTI 系統,x[n] 是輸入訊號,y[n] 是輸出訊號,卷積運算為:

\(y[n]=\sum_{k=-∞}^{∞}h[k] \cdot x[n-k]\)

假設輸入訊號為複數指數訊號 \(x[n]=e^{jωn}\)

則輸出訊號為 \(y[n]=\sum_{k=-∞}^{∞}h[k] \cdot e^{jω(n-k)} = (\sum_{k=-∞}^{∞}h[k] \cdot e^{jωk})e^{jωn}\)

可改寫為 \(y[n]=H(e^{jω})e^{jωn}\)

因此,可得到以下公式:

\(H(e^{jω}) = \sum_{n=-∞}^{∞}h[n]e^{-jωn}\)

這個公式稱為 DSP 的頻率響應 frequency response

頻率響應就是對脈衝響應 h[n] 求離散時間傅立葉轉換 DTFT 的結果

因為 DTFT 的結果通常是複數,因此經常取其強度 magnitude,與相位角 phase angle

\(|H(e^{jω})|\) 與 \(ang\{H(e^{jω})\}\)

結果都是實數。以圖形表示就是強度頻譜 magnitude spectrum 與相位頻譜 phase spectrum,這是在頻率域的操作特性

頻率響應也常以分貝 decibels, dB 為單位表示,稱為系統的增益 Gain

\(G(ω) = 20 \cdot log_{10}|H(e^{jω}| (dB)\)

濾波器分類

根據系統(濾波器)的頻率響應,可分成以下幾種:

  • 低通濾波器 Lowpass Filter

    使得低頻範圍的訊號通過,抑制高頻範圍的訊號,頻率的閥值稱為截止頻率(Cutoff Frequency) ,定義為 \(f_c\),對應的截止角頻率 \(ω_c\)

  • 高通濾波器 Highpass Filter

    低通濾波器的相反,讓高頻範圍訊號通過,抑制低頻訊號

  • 帶通濾波器 Bandpass Filter

    讓頻率介於 \(f_1 ~ f_2\) 或角頻率 \(ω_1~ω_2\) 之間的訊號通過,抑制其他頻率範圍的訊號

  • 帶阻濾波器 Bandstop Filter

    帶通濾波器的相反,抑制頻率範圍落在 \(f_1 ~ f_2\) 或角頻率 \(ω_1~ω_2\) 之間的訊號,使得其他頻率範圍的訊號通過

  • 陷波濾波器 Notch Filter

    是帶阻濾波器的一種,僅抑制某個特定頻率,ex: 60Hz 的干擾雜訊。陷波濾波器主要是使得某特定頻率的訊號變為 0,因此也稱為歸零濾波器 Nulling Filter

  • 梳狀濾波器 Comb Filter

    Comb Filter 其頻率響應的形狀跟梳子一樣,主要是同時抑制多個頻率範圍,頻率與頻率之間,通常是設定為整數倍數

  • 全通濾波器 All-Pass Filter

    讓所有頻率範圍都通過,通常會調整訊號在某些頻率範圍的相位移

濾波器也可根據訊號的型態分成 類比濾波器 Analog Filter,與數位濾波器 Digital Filter

頻率響應範例

介紹幾種經典的頻率響應

理想濾波器

濾波器的用途是讓某些特定範圍的頻率分量通過,抑制其他範圍的頻率分量,為了不造成訊號失真,頻率通過的頻率分量的頻率響應值是設定為 1,抑制的頻率分量,頻率響應值設定為 0。

使得頻率分量通過的頻率範圍稱為通過頻帶 (Passband),抑制的頻率範圍稱為截止頻帶 (Stopband)

  • 理想低通濾波器 Ideal Lowpass Filter

    轉換函式:\( H_{LP}(e^{jω}) = \left\{ \begin{array}{ll} 1 & \mbox{if 0 ≤ ω ≤ \)ω_c\(} \\ 0 & \mbox{otherwise} \end{array} \right.\)

    在此只定義單邊 Single-Sided 的頻率響應,ω 介於 0 ~ π 之間,且 \(ω_c=2 \pi f_c\),其中 \(f_c\) 稱為截止頻率 Cutoff Frequency

  • 理想高通濾波器 Ideal Highpass Filter

    轉換函式:\( H_{HP}(e^{jω}) = \left\{ \begin{array}{ll} 1 & \mbox{if \)ω_c\( < ω ≤ π} \\ 0 & \mbox{otherwise} \end{array} \right.\)

    低通與高通濾波器的關係為: \(H_{HP}(e^{jω})=1- H_{LP}(e^{jω})\)

  • 理想帶通濾波器 Ideal Bandpass Filter

    轉換函式:\( H_{BP}(e^{jω}) = \left\{ \begin{array}{ll} 1 & \mbox{if \)ω1\( ≤ ω ≤ \)ω2\(} \\ 0 & \mbox{otherwise} \end{array} \right.\)

    允許頻率範圍 \(ω_1 ≤ ω ≤ ω_2\) 或頻率 \(f_1 ≤ f ≤ f_2\) 的訊號通過,\(f_1, f_2\) 稱為截止頻率

  • 理想帶阻濾波器 Ideal Bandstop Filter

    轉換函式:\( H_{BS}(e^{jω}) = \left\{ \begin{array}{ll} 1 & \mbox{if ω < \)ω1\( or ω > \)ω2\(} \\ 0 & \mbox{otherwise} \end{array} \right.\)

    帶通與帶阻濾波器的關係為: \(H_{BS}(e^{jω})=1- H_{BP}(e^{jω})\)


ex: Ideal Lowpass Filter定義如下,求其在離散時間域的脈衝響應 Impulse Response

​ \( H_{LP}(e^{jω}) = \left\{ \begin{array}{ll} ​ 1 & \mbox{if 0 ≤ ω ≤ \)ω_c\(} \\ ​ 0 & \mbox{otherwise} ​ \end{array} \right.\)

根據 Inverse DTFT 公式

\(x[n] = \frac{1}{2π} \int_{-∞}^{∞}X(e^{jωn})e^{jωn}dω\)

Ideal Lowpass Filter 在離散時間域的脈衝響應為:

\(h_{LP}[n] = \frac{1}{2π} \int_{-∞}^{∞}H_{LP}(e^{jωn})e^{jωn}dω \\ = \frac{1}{2π} \int_{-ω_c}^{ω_c}e^{jωn}dω \\ = \frac{1}{2π} [\frac{1}{jn}e^{jωn}]_{-ω_c}^{ω_c} \\ = \frac{1}{2π} [\frac{e^{jω_cn}}{jn}- \frac{e^{-jω_cn}}{jn} ] \\ = \frac{sin(ω_cn)}{πn} \quad\quad \mbox{-∞<n<∞}\)

當 n=0,使用羅必達規則 L'Hopital's Rule

\(h_{LP}[0] = \lim_{n→0} \frac{sin(ω_cn)}{πn} \\ = \lim_{n→0} \frac{ \frac{d}{dn} sin(ω_cn)}{ \frac{d}{dn} (πn) } \\ = \lim_{n→0} \frac{ cos(ω_cn) \cdot ω_c}{ π } \\ = \frac{ ω_c }{ π }\)

因為 \(-∞ < n < ∞\),無法在離散時間域中,定義有限長度的脈衝響應,藉以實現頻率域中的 ideal lowpass filter。脈衝響應在計算時有用到過去、現在、未來的數位訊號,不是因果系統。脈衝響應不符合絕對可加總 Absolutely Summable 原則,系統也不具 BIBO 穩定性。

用 python 處理離散時間域的脈衝響應時,因為離散序列是有限長度,以近似方式接近 ideal lowpass filter。原則上採用的濾波器大小,是使用奇數

當 filter 大小為 5時,-2<n<2

當 filter 大小為 11 時,-5<n<5

ideal lowpass filter 的脈衝響應,具有 sinc 函數的特性,向左右兩邊無限延伸

以下範例選擇 \(ω_c=\frac{\pi}{2}\),因此當 n=0,\(\frac{ω_c}{\pi}=0.5\) ,可調整 \(ω_c\) 觀察脈衝響應的結果

import numpy as np
import matplotlib.pyplot as plt

def plot(filter_size, plotpos):
    filter_half = int( filter_size / 2 )
    wc = np.pi / 2

    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] =  wc/np.pi
        else:
            h[n+filter_half] = np.sin( wc * n ) / ( np.pi * n )

    plt.subplot(plotpos)
    plt.stem( na, h, use_line_collection=True )
    plt.xlabel( 'n (filter_size='+ str(filter_size) +')')
    plt.ylabel( 'h[n]' )

plot(5, 221)
plot(11, 222)
plot(21, 223)
plot(31, 224)

plt.show( )

以下選擇 \(ω_c=\frac{\pi}{10}\),因此當 n=0,\(\frac{ω_c}{\pi}=1/10\)


由於有限長度的脈衝響應,無法提供 ideal filter 的操作特性,在 DSP 實際應用時,通常是讓頻率響應在通過頻帶與截止頻帶之間,由 1 逐漸降為 0,且在通過頻帶或截止頻帶內,頻率響應允許有些微變動。

剛剛是求濾波器大小 5, 11,21, 31 的脈衝響應,現在改求頻率響應。

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


def plot(filter_size, plotpos):
    filter_half = int( filter_size / 2 )
    wc = np.pi / 2

    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] =  wc/np.pi
        else:
            h[n+filter_half] = np.sin( wc * n ) / ( np.pi * n )

    w, H = signal.freqz( h )
    mag = abs( H )

    plt.subplot(plotpos)
    plt.plot( w, mag )
    plt.xlabel( r'$\omega$'+' (filter_size='+ str(filter_size) +')' )
    plt.ylabel( 'Magnitude' )


plot(5, 221)
plot(11, 222)
plot(21, 223)
plot(31, 224)

plt.show( )

結果發現,若 filter size 越大,或脈衝響應樣本越多,越接近理想的 lowpass filter。此外可注意到頻率響應在通過頻帶與截止頻帶之間,出現了 Gibbs 現象。

平均濾波器

average filter 定義為

\(h[n] = \frac{1}{M}{1,1,1,...,1}, n=0,1,...,M-1\) ,M為濾波器大小 filter size

平均濾波器的頻率響應,可根據其 DTFT 得來

\(H(e^{jω}) = \sum_{n=-∞}^{∞}h[n]e^{-jωn} \\ = \frac{1}{M} \sum_{n=0}^{M-1}e^{-jωn} \\ = \frac{1}{M} ( \sum_{n=0}^{∞}e^{-jωn} - \sum_{n=M}^{∞}e^{-jωn}) \\ = \frac{1}{M} ( \sum_{n=0}^{∞}e^{-jωn} ) (1 - e^{-jωM}) \\ = \frac{1}{M} \frac{1 - e^{-jωM}}{1 - e^{-jω}} \\ = \frac{1}{M} \frac{sin(Mω/2)}{sin(ω/2)} e^{-j(M-1)ω/2} \)

取其強度 magnitude

\(|H(e^{jω})| = |\frac{sin(Mω/2)}{sin(ω/2)}|\)

以圖形表示,就是平均濾波器的頻率響應


ex: 若平均濾波器定義 \(h[n]=\frac{1}{M}{1,1,...,1}, n=0,1,...,M-1\)

當 M=5, 15,求其頻率響應

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

def average_filter(filter_size):
    h = np.ones( filter_size ) / filter_size
    return h

def plot(filter_size, plotpos):
    h = average_filter(filter_size)

    w, H = signal.freqz( h )
    mag = abs( H )

    plt.subplot(plotpos)
    plt.plot( w, mag )
    plt.xlabel( r'$\omega$'+' (filter_size='+ str(filter_size) +')' )
    plt.ylabel( 'Magnitude' )

plot(5, 121)
plot(15, 122)

plt.show( )

平均濾波器有 lowpass filter 的特性,當 size 越小,通過頻帶越寬。但平均濾波器有允許少量高頻分量通過

高斯濾波器

Gaussian Filter 可定義為

\(g[n]=e^{-n^2/2σ^2}\) ,其中 σ 為標準差

chap9: 高斯函數的傅立葉轉換,形成另一個高斯函數。因此高斯函數的頻率響應,可用另一個高斯函數表示

ex: Gaussian Filter 定義為 \(g[n]=e^{-n^2/2σ^2}\) ,當 σ=1, 3,求其頻率響應

note: 程式中有對係數進行正規化,讓 \(\sum_{n}g[n]=1\)

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

def plot(sigma, plotpos):
    filter_size = int( 6 * sigma + 1 )              # 濾波器大小
    gauss = signal.gaussian( filter_size, sigma )   # 濾波器係數
    sum = np.sum( gauss )                           # 正規化
    gauss = gauss / sum

    w, H = signal.freqz( gauss )
    mag = abs( H )

    plt.subplot(plotpos)
    plt.plot( w, mag )
    plt.xlabel( r'$\omega$'+' (sigma='+ str(sigma) +')' )
    plt.ylabel( 'Magnitude' )

plot(sigma=1, plotpos=121)
plot(sigma=3, plotpos=122)

plt.show( )

高斯濾波器有 lowpass filter 的特性,標準差越小,通過頻帶較寬。在高頻範圍的抑制效果,比平均濾波器好。

FIR 濾波器

一階 FIR 濾波器有 lowpass, highpass filter 兩種

FIR filters 裡面最簡單的是移動平均濾波器 Moving Averge Filter,濾波器大小為 2

最簡單的 FIR Lowpass Filter 的一階轉換函式定義為:

\(H_0(z)=\frac{1}{2}(1+z^{-1})\)

最簡單的 FIR Highpass Filter 的一階轉換函式定義為:

\(H_1(z)=\frac{1}{2}(1-z^{-1})\)

因為 \(H_1(z) = 1 - H_0(z) = 1-\frac{1}{2}(1+z^{-1}) = \frac{1}{2}(1-z^{-1})\)

ex: 分別求上面兩個 FIR filters 的頻率響應

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

b = np.array( [ 0.5, 0.5 ] )    # 低通濾波器
w, H0 = signal.freqz( b )
H0 = abs( H0 )

b = np.array( [ 0.5, -0.5 ] )   # 高通濾波器
w, H1 = signal.freqz( b )
H1 = abs( H1 )

plt.subplot(121)
plt.plot( w, H0 )
plt.xlabel( r'$\omega$' +' (lowpass)' )
plt.ylabel( 'Magnitude' )

plt.subplot(122)
plt.plot( w, H1 )
plt.xlabel( r'$\omega$' +' (highpass)' )
plt.ylabel( 'Magnitude' )

plt.show( )

根據 Lowpass Filter 的頻率響應,強度的最大/最小值分別是 1 與 0。頻率響應也是從 1 逐漸降到 0。

在討論 Lowpass Filter 時,當截止頻率 \(ω=ω_0\),強度降為

\(|H_0(e^{jω_c})| = \frac{1}{\sqrt{2}}|H_0(e^{j0})|\)

如果以分貝表示則為

\(20 log_{10}|H_0(e^{jω_c})| = 20 log_{10}( \frac{1}{\sqrt{2}}|H_0(e^{j0})| ) \\ = 20 log_{10}|H_0(e^{j0})| - 20 log_{10}\sqrt{2} \\ = 0 - 3.0103 ≈ -3dB\)

\(ω_c\)或 \(f_c\) 稱為 3dB 截止頻率 (3dB Cutoff Frequency),通常被視為是通過頻帶 (Passband) 的邊緣。因此,\(ω_c\)或 \(f_c\) 也常被稱為 通過頻帶邊緣頻率(Passband Edge Frequency)


串接多個 FIR Filter,可用來近似 ideal lowpass filter

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

# order = eval( input( "Enter number of cascade FIR filters: " ) )

def plot(order, plotpos):
    b = np.array( [ 0.5, 0.5 ] )
    w, H = signal.freqz( b )

    for i in range( order - 1 ):
        H *= H

    H0 = abs( H )

    plt.subplot(plotpos)
    plt.plot( w, H0 )
    plt.xlabel( r'$\omega$' +'(order='+str(order)+')' )
    plt.ylabel( 'Magnitude' )

plot(1, 221)
plot(2, 222)
plot(3, 223)
plot(4, 224)

plt.show( )

IIR 濾波器

  • IIR Lowpass Filter

    一階 IIR Lowpass Filter 轉換函式定義為 \(H_0(z)= \frac{1-α}{2} \frac{1+z^{-1}}{1-αz^{-1}}\) ,其中 \(|α|<1\) 為穩定的 IIR Filter

  • IIR Highpass Filter

    一階 IIR Highpass Filter 轉換函式定義為 \(H_1(z)= \frac{1+α}{2} \frac{1-z^{-1}}{1-αz^{-1}}\) ,其中 \(|α|<1\) 為穩定的 IIR Filter

ex: 一階 IIR Lowpass/Highpass Filter,當 α = 0.2, 0.4, 0.6, 0.8 時,求其頻率響應

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

# alpha = eval( input( "Please enter alpha (< 1): " ) )

def iir_lowpass(alpha, plotpos):
    b = np.array( [ 1 - alpha, 1 - alpha ] )
    a = np.array( [ 2, -2 * alpha ] )
    w, H = signal.freqz( b, a )
    H0 = abs( H )

    plt.subplot(plotpos)
    plt.plot( w, H0 )
    plt.xlabel( r'$\omega$' +' (alpha='+str(alpha)+')' )
    plt.ylabel( 'Magnitude' )

def iir_highpass(alpha, plotpos):
    b = np.array( [ 1 + alpha, - ( 1 + alpha ) ] )
    a = np.array( [ 2, -2 * alpha ] )
    w, H = signal.freqz( b, a )
    H1 = abs( H )

    plt.subplot(plotpos)
    plt.plot( w, H1 )
    plt.xlabel( r'$\omega$' +' (alpha='+str(alpha)+')' )
    plt.ylabel( 'Magnitude' )

plt.figure( 1 )
iir_lowpass(0.2, 221)
iir_lowpass(0.4, 222)
iir_lowpass(0.6, 223)
iir_lowpass(0.8, 224)

plt.figure( 2 )
iir_highpass(0.2, 221)
iir_highpass(0.4, 222)
iir_highpass(0.6, 223)
iir_highpass(0.8, 224)

plt.show( )

IIR Lowpass Filters

IIR Highpass Filters

  • IIR 帶通濾波器

    二階 IIR Bandpass Filter 的轉換函式定義為

    \(H_{BP}(z) = \frac{1-α}{2} \frac{1-z^{-2}}{1-β(1+α)z^{-1}+αz^{-2}}\) ,其中 α, β 都小於 1

ex1: 當 α = 0.2, 0.5, 0.8 且 β=0.5 時,求其頻率響應

ex2: 當 α = 0.5 且 β=0.2, 0.5, 0.8 時,求其頻率響應

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

def IIR_bandpass_frequency_response(alpha, beta, plotpos):
    b = np.array( [ 1 - alpha, 0, -( 1 - alpha ) ] )
    a = np.array( [ 2, -2 * beta * ( 1 + alpha ), 2 * alpha ] )
    w, H = signal.freqz( b, a )
    H1 = abs( H )

    plt.subplot(plotpos)

    plt.plot( w, H1, '-', label = r'$\alpha$ = '+str(alpha)+','+r'$\beta$ = '+str(beta) )
    plt.legend( loc = 'upper right' )
    plt.xlabel( r'$\omega$' )
    plt.ylabel( 'Magnitude' )


IIR_bandpass_frequency_response(0.2, 0.5, 121)
IIR_bandpass_frequency_response(0.5, 0.5, 121)
IIR_bandpass_frequency_response(0.8, 0.5, 121)

IIR_bandpass_frequency_response(0.5, 0.2, 122)
IIR_bandpass_frequency_response(0.5, 0.5, 122)
IIR_bandpass_frequency_response(0.5, 0.8, 122)

plt.show( )

α 參數用來調整頻寬,β 參數用來調整頻率中心

梳狀濾波器

  • Comb Lowpass Filter

    簡易梳狀低通濾波器的轉換函式,可定義為

    \(H_{CLP}(z) = \frac{1}{2}(1+z^{-M})\),M 為正整數

  • Comb Highpass Filter

    簡易梳狀高通濾波器的轉換函式,可定義為

    \(H_{CHP}(z) = \frac{1}{2}(1-z^{-M})\) ,M 為正整數

ex: 當 M=10,求頻率響應

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

# M = eval( input( "Please enter M for the Comb Filter: " ) )

def comb_lowpass_filter_frequency_response(M, plotpos):

    b = np.zeros( M + 1 )           # 梳狀低通濾波器
    b[0] = 0.5
    b[M] = 0.5
    w, H = signal.freqz( b, 1 )
    H0 = abs( H )

    plt.subplot(plotpos)
    plt.plot( w, H0 )
    plt.xlabel( r'$\omega$'+' (comb lowpass)' )
    plt.ylabel( 'Magnitude' )

def comb_highpass_filter_frequency_response(M, plotpos):

    b = np.zeros( M + 1 )
    b[0] = 0.5
    b[M] = -0.5                     # 梳狀高通濾波器
    w, H = signal.freqz( b, 1 )
    H1 = abs( H )

    plt.subplot(plotpos)
    plt.plot( w, H1 )
    plt.xlabel( r'$\omega$'+' (comb highpass)' )
    plt.ylabel( 'Magnitude' )

M = 10
comb_lowpass_filter_frequency_response(M, 121)
comb_highpass_filter_frequency_response(M, 122)

plt.show( )

References

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

2021/02/01

頻譜分析

頻譜分析 spectrum analysis 目的是分析訊號在不同頻率下的分量,以圖形表示。頻譜分析的方法可分成兩大類型:傅立葉頻譜 Fourier Spectrum 與 功率頻密度 Power Spectral Density (PSD)

根據傅立葉基本理論:「任意週期性函數,均可表示成不同頻率、不同振幅的正弦函數或餘弦函數,加總而得的無窮級數。」更進一步,傅立葉轉換也可適用於非週期性函數。

頻譜分析 spectrum analysis 的目的是將訊號分解成不同頻率的分量,並以圖形表示,稱為頻譜(spectrum)。頻譜可用來分析訊號中的頻率分量,例如男生的語音聲音低沈,屬於低頻,女生或兒童是高頻,可藉此辨識說話者是男生、女生或兒童。

弦波的定義為

\(x(t)=A cos(ω_0+ϕ)\) 或是 \(x(t)=A cos(2 \pi f_0t+ϕ)\) ,其中振幅 A ,角頻率 \(ω_0\),頻率 \(f_0\)

根據反尤拉公式 \(cosθ = \frac{e^{jθ}+e^{-jθ}}{2}\)

代入可得

\(x(t) = A cos(ω_0t+ϕ) = A [e^{j(ω_0t+ϕ)}+e^{-j(ω_0t+ϕ)}]/2 \)

因此弦波包含了兩個振幅 Amplitude 的分量,當頻率為 \(ω_0, -ω_0\),振幅均為 A/2,如果以圖形表示,就稱為振幅頻譜 Amplitude Spectrum。

同時弦波 x(t) 也包含了兩個相位移 Phase Shift 的分相,當頻率為 \(ω_0, -ω_0\),相位移分別為 \(ϕ, -ϕ\),如果以圖形表示,稱為相位頻譜 Phase Spectrum。

振幅頻譜及相位頻譜,通稱為傅立葉頻譜 Fourier Spectrum 或頻率頻譜 Frequency Spectrum,或簡稱 Spectrum

由於頻譜是以原點為中心分成兩邊,經常稱為雙邊頻譜 Doubl-Sided Spectrum。頻譜的橫軸以角頻率 ω為主,通常介於 -π ~ π 之間,單位為 弧度/秒 (Raidans/Second),此外,頻譜的橫軸也經常以頻率 f 為主,單位為 Hz

若諧波由多個弦波組成,頻率為基礎頻率的整數倍數,定義為

\(x(t)=\sum_{k=1}^{N}A_k cos(2 \pi f_kt)\),其中 \(f_k=k \cdot f_1\),\(f_1\) 稱為基礎頻率

用反尤拉公式展開後得到振幅頻譜。對應不同頻率的分量,其頻率分量是弦波振幅的一半。

傅立葉頻譜

Fourier Spectrum 是指時間域訊號在頻率域的圖形表示法,可使用傅立葉轉換得來

先對輸入的數位訊號 \(x[n], n=0,1,...,N-1\) 進行離散傅立葉轉換

\(X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2 \pi kn/N}, k=0,1,...,N-1\)

輸出結果為 X[k],是複數構成的離散數列。在 DSP 通常採用 Fast Fourier Transform (FFT) 縮短離散副業轉換的計算所需的時間

FFT 結果為複數,強度為 \(|X[k|, k=0,1,...,N-1\),以圖形表示,稱為強度頻譜 Magnitude Spectrum,也就是前述的振幅頻譜 Amplitude Spectrum。

若求複數的幅角 Argument 或相位角 Phase Angle,並以圖形表示,則稱為相位頻譜 Phase Spectrum

\(arg(X[k]), k=0,1,...,N-1\)


ex: 數位訊號 \(x={1,2,4,3}, n=0,1,2,3\) ,求強度頻譜的圖形

離散傅立葉轉換結果為複數

\(X=\{10, -3+j, 0, -3-j\}\)

強度為 \(|X[0]|=10, |X[1]|=\sqrt{(-3)^2+1^2}=\sqrt{10}\)

Magnitude Spectrum:

弦波

定義 \(x(t)=cos(2 \pi 100 \cdot t)\),振幅 A=1,頻率 f=100,相位移 ϕ=0,在 t=0~1s 之間取 1000 個樣本,取樣頻率 1000Hz,求弦波的傅立葉頻譜?

這邊先以強度頻譜為主,先不討論相位頻譜。由圖形結果可發現 Fourier Spectrum 在 k=100, 900 兩處出現峰值 peaks,形成 Double-Sided Spectrum 的對稱結構,但不是以原點為中心,其中 k=100 的峰值就是弦波的頻率,振幅為 \((A/2) \cdot N = 500\)

為方便傅立葉頻譜的觀察與分析,通常會對 FFT 的結果進行平移 shift,讓頻譜的中心為原點。numpy 的 fftshift 有實現平移的功能。因為 N=1000,根據 Nyquist Theorem,可分析的最大頻率為 500 Hz

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

# time,1000 個 sample 點
t = np.linspace( 0, 1, 1000, endpoint = False )
# 弦波 原始訊號
x = np.cos( 2 * np.pi * 100 * t )

## fft 計算
X = fft( ( x ) )
# 強度
Xm = abs( X )

# 強度頻譜
plt.subplot(121)
plt.plot( Xm )
plt.xlabel( 'k' )
plt.ylabel( 'Magnitude' )

# t = np.linspace( 0, 1, 1000, endpoint = False )
# x = np.cos( 2 * np.pi * 100 * t )

# numpy.fft.fftfreq(n, d=1.0) 定義 Fourier Spectrum 的取樣頻率
# n 為樣本數, d 為取樣週期 (1 second)
# fftshift 用來平移,以原點為中心
f = fftshift( fftfreq( 1000, 0.001 ) )
X = fftshift( fft( x ) )
Xm = abs( X )

plt.subplot(122)
plt.plot( f, Xm )
plt.xlabel( 'f' )
plt.ylabel( 'Magnitude' )

plt.show( )

方波

若方波的參數為 A=1, f=10Hz,在一秒內震盪 10 次,若取樣頻率為 1000Hz,求方波的 Fourier Spectrum ?

方波由無限多個弦波加總而成,可視為一種諧波,基礎頻率為 10Hz,也就是方波的頻率

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

# time,1000 個 sample 點
t = np.linspace( 0, 1, 1000, endpoint = False )
# 方波 原始訊號
x = signal.square( 2 * np.pi * 10 * t )

# numpy.fft.fftfreq(n, d=1.0) 定義 Fourier Spectrum 的取樣頻率
# n 為樣本數, d 為取樣週期 (1 second)
# fftshift 用來平移,以原點為中心
f = fftshift( fftfreq( 1000, 0.001 ) )
X = fftshift( fft( x ) )
Xm = abs( X )

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

plt.show( )

節拍波

beat wave是由兩個不同頻率的弦波乘法運算得來

\(x(t)= A \cdot cos(2 \pi f_1 \cdot t) \cdot cos(2 \pi f_2 \cdot t)\)

通常 \(f_1\) 為低頻訊號的頻率。\(f_2\) 是高頻訊號的頻率,稱為載波頻率 carrier frequency,也常表示為 \(f_c\)。乘法運算是調變 modulation 的主軸,常被用在 DSP 系統中

ex: 節拍波定義為

\(x(t)= \cdot cos(2 \pi 20 \cdot t) \cdot cos(2 \pi 200 \cdot t)\),其中 A=1,\(f_1=20Hz, f_2=200Hz\),取樣頻率為 1000Hz,求 Fourier Spectrum?

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

f1 = 20
f2 = 200
t = np.linspace( 0, 1, 1000, endpoint = False )
# 節拍波 原始訊號
x = np.cos( 2 * np.pi * f1 * t ) * np.cos( 2 * np.pi * f2 * t )

f = fftshift( fftfreq( 1000, 0.001 ) )
X = fftshift( fft( x ) )
Xm = abs( X )

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

plt.show( )

啁啾訊號

這是非週期性訊號

ex: 啁啾訊號在時間 0~1s 之間,頻率範圍是 0Hz 線性遞增到 100Hz,取樣頻率 1000Hz,求 Fourier Spectrum?

頻譜呈現 ripple 現象

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

t = np.linspace( 0, 1, 1000, endpoint = False )
x = signal.chirp( t, 0, 1, 100, 'linear' )

f = fftshift( fftfreq( 1000, 0.001 ) )
X = fftshift( fft( x ) ) 
Xm = abs( X )  

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

plt.show( )

功率頻密度 Power Spectral Density (PSD)

數位訊號的能量 Energy 可定義為

\(E=\int_{-∞}^{∞}|x(t)|^2 {\rm d} t\)

但只適合用來描述有限的訊號。當訊號能量集中在有限時間區間時,總能量有限,因此可透過積分計算而得。

以離散時間域來說,數位訊號的能量可用以下公式計算:

\(E=\sum_{-∞}^{∞}x^2[n]\)


帕塞瓦爾定理

若訊號 \(x(t)\) 的傅立葉轉換為 \(X(f)\),也就是 \(X(f)=\int_{-∞}^{∞}x(t)e^{-j(2 \pi ft)} {\rm d}t\)

則訊號的能量滿足以下公式:

\(E=\int_{-∞}^{∞}|x(t)|^2 {\rm d}t = \int_{-∞}^{∞}|X(f)|^2 {\rm d}f\)

「時間域訊號,若取平方總和(或積分),則結果與傅立葉轉換後頻率域訊號的平方總和(或積分)相等」

換句話說,如果想分析訊號的能量,除了在時間域進行積分外,也可以在頻率域進行積分,結果一樣。

證明:

根據反傅立葉轉換

\(x(t)=\frac{1}{2\pi}\int_{-∞}^{∞}X(ω)e^{jωt} {\rm d}ω\)

則 \(E = \int_{-∞}^{∞} x(t)[\frac{1}{2\pi} \int_{-∞}^{∞}X(ω) e^{jωt} {\rm d}ω ] {\rm d}t \\ = \frac{1}{2 \pi} \int_{-∞}^{∞}X(ω) [\int_{-∞}^{∞} x(t) e^{jωt} {\rm d}t]{\rm d}ω \\ = \frac{1}{2 \pi} \int_{-∞}^{∞}X(ω) X(-ω){\rm d}ω \\ = \frac{1}{2 \pi} \int_{-∞}^{∞}X(ω) X^*(ω){\rm d}ω \quad\quad X^*(ω) 是 X(ω) 的共軛複數 \\ = \frac{1}{2 \pi} \int_{-∞}^{∞} |X(ω)|^2 {\rm d}ω \\ = \frac{1}{2 \pi} \int_{-∞}^{∞} |X(f)|^2 {\rm d}f\)


能量譜密度

若類比訊號 x(t) 的傅立葉轉換為 X(f),則能量譜密度 (Energy Spectral Density) 定義為:

\(X_{xx}(f) = |X(f)|^2\)

能量譜密度的目的是用來描述訊號在不同頻率下的能量分佈狀況,以圖形表示

能量譜密度可根據 自相關 Autocorrelation 取其傅立葉轉換得來

假設訊號 x(t) 的自相關定義為

\(R(τ)=\int_{-∞}^{∞}x^*(t) \cdot x(t+τ) {\rm d}t\)

則其傅立葉轉換為:

\(F\{R(τ)\} = \int_{-∞}^{∞}R(τ)e^{-j2 \pi ft} {\rm d}τ \\ = \int_{-∞}^{∞} [\int_{-∞}^{∞} x^*(t) \cdot x(t+τ) {\rm d}t ] e^{-j2 \pi ft} {\rm d}τ \\ = \int_{-∞}^{∞} x^*(t) [\int_{-∞}^{∞} x(t+τ) e^{-j2 \pi ft} {\rm d}τ ] {\rm d}t \quad\quad 假設 \hatτ = t+τ, {\rm d}\hatτ = {\rm d}τ \\ = \int_{-∞}^{∞} x^*(t) [\int_{-∞}^{∞} x(\hatτ) e^{-j2 \pi f \hatτ} \cdot e^{j2 \pi f t} {\rm d}τ ] {\rm d}t \\ = \int_{-∞}^{∞} x^*(t) [\int_{-∞}^{∞} x(\hatτ) e^{-j2 \pi f \hatτ} {\rm d}τ ] e^{j2 \pi f t} {\rm d}t \\ = \int_{-∞}^{∞} x^*(t) X(f) e^{j2 \pi f t} {\rm d}t \\ = X(f) \int_{-∞}^{∞} x^*(t) e^{j2 \pi f t} {\rm d}t \\ = X(f) X^*(f) \\ = |X(f)|^2 \\ = S_{xx}(f)\)


功率頻密度

功率頻密度 Power Spectral Density 就是訊號在不同頻率下的功率分佈圖

上面的 能量譜密度 定義,只適用於能量集中在某個時間區間的訊號,例如:脈衝訊號。對於持續存在的連續訊號,或是平穩過程 Stationary Process 的訊號,要改用 功率頻密度 Power Spectral Density (PSD),用來描述訊號在不同頻率下,功率 Power 的分佈狀況

評估 功率頻密度 有兩個具代表性的方法:

  • 週期圖 Periodogram

    週期圖是評估功率頻密度最基本的方法,主要根據傅立葉轉換的結果,取平方得來

  • Welch Method

    由 P. D. Welch 提出,也是基於週期圖的方法,但同時參考 Barlett's Method,牽涉功率的平均值運算,又稱為平均週期圖法 Method of Averaged Periodograms


ex: 若產生時間介於 0~1s 的均勻雜訊,值介於 -1~1 之間,取樣頻率 1000Hz,使用週期圖語 Welch Method 求功率頻密度 PSD

這兩種方法得到的功率頻密度,頻率範圍是 0~500Hz (因為 Nyquist Theorem,上限是 500 Hz),在不同頻率下,均勻雜訊都含有少量功率

若以比較嚴謹的定義,理想的 white noise,其功率頻密度在不同頻率下,應該要為固定常數。因此均勻雜訊,並不完全符合 white noise 的條件

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

fs = 1000
t = np.linspace( 0, 1, fs, endpoint = False )
# 均勻雜訊
noise = random.uniform( -1, 1, fs )

# 用兩種方法求得 PSD
f1, pxx1 = signal.periodogram( noise, fs )
f2, pxx2 = signal.welch( noise, fs )

plt.subplot(121)
plt.plot( f1, pxx1 )
plt.xlabel( 'frequency (Hz)' )
plt.ylabel( 'PSD' )

plt.subplot(122)
plt.plot( f2, pxx2 )
plt.xlabel( 'frequency (Hz)' )
plt.ylabel( 'PSD' )

plt.show( )


ex: 訊號定義為

\(x(t) = 10 cos(2 \pi 100 t) + 5 cos(2 \pi 200 t) + η(t)\)

其中包含兩個弦波,η(t) 是介於 -1~1 之間的均勻雜訊,取樣頻率 1000 Hz,用 週期圖與 Welch Method 求 PSD

兩種方法呈現的 PSD 都是單邊頻譜 Single-Sided Spectrum,範圍是 0~500 Hz。

雖然有雜訊干擾, PSD 在 100Hz 與 200Hz 呈現功率的峰值

概括來說,Welch Method 在頻率域上的解析度較差,但具有較好的抗雜訊能力

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

fs = 1000
t = np.linspace( 0, 1, fs, endpoint = False )
# 原始訊號
x = 10 * np.cos( 2 * np.pi * 100 * t ) + 5 * np.cos( 2 * np.pi * 200 * t ) + random.uniform( -1, 1, fs )

f1, pxx1 = signal.periodogram( x, fs )
f2, pxx2 = signal.welch( x, fs )

plt.subplot(121)
plt.plot( f1, pxx1 )
plt.xlabel( 'frequency (Hz)' )
plt.ylabel( 'PSD' )

plt.subplot(122)
plt.plot( f2, pxx2 )
plt.xlabel( 'frequency (Hz)' )
plt.ylabel( 'PSD' )

plt.show( )

也可再對 PSD 取對數,查看主要頻率外,其他頻率的分佈狀況

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

fs = 1000
t = np.linspace( 0, 1, fs, endpoint = False )
x = 10 * np.cos( 2 * np.pi * 100 * t ) + 5 * np.cos( 2 * np.pi * 200 * t )
noise = random.uniform( -1, 1, fs )
y = x + noise

f1, pxx1 = signal.periodogram( y, fs )
f2, pxx2 = signal.welch( y, fs )

# 取對數
lp1 = 10 * np.log10( pxx1 )
lp2 = 10 * np.log10( pxx2 )

plt.subplot(121)
plt.plot( f1, lp1 )
plt.xlabel( 'frequency (Hz)' )
plt.ylabel( 'PSD' )
plt.axis( [ 0, 500, -50, 10 ] )

plt.subplot(122)
plt.plot( f2, lp2 )
plt.xlabel( 'frequency (Hz)' )
plt.ylabel( 'PSD' )
plt.axis( [ 0, 500, -40, 10 ] )

plt.show( )

References

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

2021/01/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程式實作(附範例光碟)(第二版)