2021/01/18

FIR 濾波器

Finite Impulse Response Filter: FIR 濾波器是線性時間不變性 Linear Time-Invariant LTI 系統,而且是因果系統 (causal system) (輸出的訊號只根據目前與過去的訊號運算得來)。

FIR Filter

前面討論過移動平均濾波器 (Moving Average Filter),有兩種定義方式:

  1. \(y[n]=\frac{1}{3}(x[n]+x[n+1]+x[n+2])\) ,包含未來的訊號,屬於非因果系統
  2. \(y[n]=\frac{1}{3}(x[n]+x[n-1]+x[n-2])\) ,計算結果跟上面一樣,差異是輸出訊號的時間會比較慢,這個是因果系統,比較適合即時性 DSP應用

FIR Filter 定義:

假設輸入訊號 x[n],輸出訊號 y[n],則 FIR Filter 為

\(y[n]=\sum_{k=0}^{M}b_kx[n-k] = b_0x[n]+b_1x[n-1]+...+b_Mx[n-M]\)

其中 \(b_k, k=0,1,...M\) 稱為 FIR 的係數,M 稱為濾波器的階數 Order

輸出訊號是目前訊號跟過去M個訊號的加權平均(weighted average) 運算而得。FIR Filter 是 causal system 適合即時性的 DSP 應用。

因為 FIR 必須等到 M+1 個輸入訊號後,才能開始輸出訊號,故須具備足夠的記憶能力。


假設 \(x[n]=δ[n]\) 單位脈衝 unit impulse

FIR filter 為 \(y[n]=\sum_{k=0}^{M}b_k \cdot x[n-k] = \sum_{k=0}^{M}b_k \cdot δ[n-k] = b_0δ[n]+b_1δ[n-1]+...+b_Mδ[n-M]\)

若計算 z 轉換

\(Y(z) = Z\{y[n]\} = Z\{\sum_{k=0}^{M}b_k \cdot x[n-k]\} \\ = Z\{b_0δ[n]+b_1δ[n-1]+...+b_Mδ[n-M]\} \\ = Z\{b_0δ[n]\}+Z\{b_1δ[n-1]\}+...+Z\{b_Mδ[n-M]\} \\ = b_0Z\{δ[n]\}+b_1Z\{δ[n-1]\}+...+b_MZ\{δ[n-M]\} \\ = b_0X(z)+b_1z^{-1}X(z)+...+b_Mz^{-M}X(z) \\ = (b_0+b_1z^{-1}+...+b_{M}z^{-M})X(z)\)

FIR Filter 的轉換函式為

\(H(z)=(b_0+b_1z^{-1}+...+b_{M}z^{-M})\)

在複數平面上,FIR filter 只有零點 zeros 沒有極點 poles,因此 FIR filter 為穩定系統。


DSP 系統常會以方塊圖的型態表示,基本建構元件為:加法器 adder、乘法器 multiplier、單位延遲 unit delay

根據 FIR 的定義 \(y[n]=\sum_{k=0}^{M}b_kx[n-k]\)

z轉換 \(Y(z)=(b_0+b_1z^{-1}+...+b_{M}z^{-M})X(z)\)

FIR 濾波器的系統方塊圖為

FIR Filter 的應用

移動平均濾波、股價趨勢分析、歸零濾波器

移動平均濾波

moving average filter 是最簡單的 FIR filter

\(b_k={\frac{1}{3}, \frac{1}{3}, \frac{1}{3}}, k=0,1,2\) 階數為 2

使用 SciPy 的 lfilter 實作 moving average filter

import numpy as np
import scipy.signal as signal

x = np.array( [ 1, 2, 4, 3, 2, 1, 1 ] )
b = np.ones( 3 ) / 3
y = signal.lfilter( b, 1, x )

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

執行結果

$ python FIR_example.py
x = [1 2 4 3 2 1 1]
y = [0.33333333 1.         2.33333333 3.         3.         2. 1.33333333]

股價趨勢分析

將股價趨勢當作數位訊號,套用 FIR filter 可產生均線,觀察股價走勢,如果平均天數為 5,可產生 5 日均線(也就是週線),如果平均天數為20,可產生月線

股價趨勢分析常用 5日均線(週線)、10日均線、20日均線(月線)、60日均線(季線)

到 Yahoo Finance 下載 TSM 的股價走勢 csv 資料,裡面有七個欄位

日期Date,開盤價Open,最高價High,最低價Low,收盤價Close,調整收盤價Adj Close,交易量Volume

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

csvDataFile = open( 'TSM.csv' )
reader = csv.reader( csvDataFile )

data = []                           # 讀取收盤價資料
for row in reader:
    data.append( row[4] )

price = []                          # 去掉第一行的資料,將收盤價轉換為股價數值
for i in range( 1, len( data ) ):
    price.append( eval( data[i] ) )

day = np.arange( len( price ) )
x = np.array( price )               # 轉換成陣列

b1 = np.ones( 5 ) / 5               # 週線
y1 = signal.lfilter ( b1, 1, x )

b2 = np.ones( 20 ) / 20             # 月線
y2 = signal.lfilter ( b2, 1, x )

# 替換中文字型,處理 matplotlib 中文 label 問題
# https://blog.csdn.net/gmr2453929471/article/details/78655834
plt.rcParams['font.sans-serif'] = ['Heiti TC'] # 步驟一(替換sans-serif字型)
plt.rcParams['axes.unicode_minus'] = False  # 步驟二(解決座標軸負數的負號顯示問題)

plt.figure( 1 )                     # 繪圖
plt.subplot(131)
plt.plot( day, x, '-', fillstyle = 'bottom' )
plt.xlabel( 'Day' )
plt.ylabel( 'Price' )
plt.axis( [ 0, len( price), 28, 45 ] )

# plt.figure( 2 )
plt.subplot(132)
plt.plot( day, x, '--', day, y1, '-' )
plt.xlabel( 'Day (週線)' )
plt.ylabel( 'Price' )
plt.axis( [ 0, len( price), 28, 45 ] )

# plt.figure( 3 )
plt.subplot(133)
plt.plot( day, x, '--', day, y2, '-' )
plt.xlabel( 'Day (月線)' )
plt.ylabel( 'Price' )
plt.axis( [ 0, len( price), 28, 45 ] )

plt.show()

歸零濾波器

Nulling Filter,用來消除某個特定頻率的濾波器,也稱為陷波濾波器 Notch Filter

因為 FIR Filter 在 z 轉換後,會產生 zeros,可以使得某個特定頻率歸零,適合用來消除特定干擾訊號 Jamming Signals。例如來自電源線的干擾訊號,通常是 60Hz 的弦波訊號,可在 DSP 中,加入 Nulling Filter 消除特定頻率的干擾訊號。

假設想要消除的干擾訊號定義為

\(\hat{x} = cos(\hat{ω}t)\) ,其中 \(\hat{ω}\) 是干擾訊號的角頻率,且 \(\hat{ω} = 2 \pi \hat{f}\)

假設取樣頻率為 \(f_s\),則干擾訊號的數位訊號可定義為

\(\hat{x}[n]=cos(2\pi\hat{f}n/f_s)\)

根據反尤拉公式

\(cos(2\pi\hat{f}n/f_s) = \frac{1}{2}(e^{j2\pi\hat{f}n/f_s}+e^{-j2\pi\hat{f}n/f_s})\)

可用兩個一階的 FIR filter 串聯成二階 FIR filter,藉以消除這兩個複數指數訊號

(

因為 FIR Filter 的轉換函式為 \(H(z)=(b_0+b_1z^{-1}+...+b_{M}z^{-M})\)

取M=1,就是一階濾波器 \(H(z)=b_0+b_1z^{-1}\)

)

因此二階 FIR filter 要設計為包含兩個零點:

\(z_1= e^{j2\pi\hat{f}n/f_s}\) , \(z_2= e^{-j2\pi\hat{f}n/f_s}\)

對應的一階 FIR filter 分別為

\(H_1(z) = 1-z_1z^{-1}\), \(H_2(z)=1-z_2z^{-1}\)

串聯後得到的二階濾波器為

\(H(z) = H_1(z)H_2(z) = (1-z_1z^{-1})(1-z_2z^{-1}) \\ = 1-(z_1+z_2)z^{-1} + (z_1z_2)z^{-2} \\ = 1-(e^{j2\pi\hat{f}n/f_s}+e^{-j2\pi\hat{f}n/f_s})z^{-1} + (e^{j2\pi\hat{f}n/f_s}e^{-j2\pi\hat{f}n/f_s}) z^{-2} \\ = 1-2cos(2\pi\hat{f}n/f_s)z^{-1}+z^{-2} \)


ex: 輸入訊號 \(x(t)=cos(2\pi\cdot10\cdot t)+cos(2\pi\cdot20\cdot t)\)

假設 \(cos(2\pi\cdot10\cdot t)\) 是原始訊號,另一個是干擾訊號,取樣頻率為 100 Hz,設計一個 Nulling Filter

因為干擾訊號的頻率 \(\hat{f}=20Hz\) 取樣頻率 100,剛剛知道Nulling Filter 二階濾波器轉換函式為

\(H(z)=1-2cos(2\pi \hat{f} n/f_s)z^{-1}+z^{-2}\)

分別代入 \(\hat{f}, f_s\)

\(H(z)=1-2cos(2\pi \cdot 20 \cdot n/100)z^{-1}+z^{-2}\)

套用 Nulling Filter 後,輸出訊號會發生時間延遲的現象。

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

fs = 100
t = np.linspace( 0, 1, fs, endpoint = False )   # 定義時間陣列
x1 = np.cos( 2 * np.pi * 10 * t )               # 原始訊號
x2 = np.cos( 2 * np.pi * 20 * t )               # 干擾訊號
x = x1 + x2                                     # 輸入訊號

b = np.array( [ 1, -2 * np.cos( 2 * np.pi * 20 / fs ), 1 ] )
y = signal.lfilter( b, 1, x )                   # FIR濾波, Nulling Filter

# plt.figure( 1 )                                   # 繪圖
plt.subplot(121)
plt.plot( t, x )
plt.xlabel( 't (second)' )
plt.ylabel( 'Amplitude' )
plt.axis( [ 0, 1, -2, 2 ] )

# plt.figure( 2 )
plt.subplot(122)
plt.plot( t, x1, '--', t, y, '-' )
plt.xlabel( 't (second)' )
plt.ylabel( 'Amplitude' )
plt.axis( [ 0, 1, -2, 2 ] )

plt.show( )

References

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

2021/01/11

z 轉換

z Transform 用來將離散的數位訊號表示成複數指數函數的數學工具。以下討論 z transform,並用 z 轉換求 LTI 系統的轉換函式,並求零點與極點等參數,分析 LTI 系統的特性,最後說明反 z 轉換。

z轉換

z轉換 源自 Laplace Transform (適用於分析連續時間域的函數),適用於離散時間域 discrete time domain 的數位訊號分析。跟離散時間傅立葉轉換 DTFT 比較,z轉換提供更廣義的訊號表示法,可用來分析 DSP 系統的操作特性。

DTFT 是針對絕對可加總序列(Absolute Summable Sequence) 進行轉換,僅適合用來分析穩定的 LTI 系統, z 轉換應用範圍較廣,也可分析不穩定的 LTI 系統。


給定離散序列 x[n],z轉換定義為:

​ \(X(z)=Z\{x[n]\}=\sum_{n=-∞}^{∞}x[n]z^{-n}\)

反z轉換定義為:

​ \(x[n]=Z^{-1}\{X(z)\} = \frac{1}{2πj} \oint_CX(z)z^{n-1}dz\) 其中 C 是包含原點的逆時針封閉路徑,落在收斂區域中。


z轉換是將離散時間域 (discrete-time domain) 的序列 x[n],轉換成 z domain 的函數 X(z)。\(Z\{\}\) 代表 z 轉換。z 轉換符合可逆性,可用反 z 轉換將 X(z) 還原為 x[n]。

跟 discrete-time fourier transform (DTFT) 的公式比較

​ \( X(e^{jω}) = \sum_{n=-∞}^{∞}x[n]e^{-jωn} \)

可發現到 \(z=e^{jw}\)

z轉換就是 DTFT 的另一種表示法。


收斂區域 Region of Convergence (ROC)

就是使得 z 轉換收斂的複數平面的 z 點集合,定義為

\( ROC = \{z: |\sum_{n=-∞}^{∞}x[n]z^{-n}|<∞\} \)

其中 C 是包含原點的逆時針封閉路徑,落在收斂區域中。


ex: 求單位脈衝 δ[n] 的 z 轉換,並決定其收斂區域 ROC

單位脈衝的定義:\( δ[n] = \left\{ \begin{array}{ll} 1, & \mbox{if n=0} \\ 0, & \mbox{if n≠0} \end{array} \right.\)

\(Z\{δ[n]\}=\sum_{n=-∞}^{∞}δ[n]z^{-n} = δ[0]z^{-0} = 1\)

ROC 為複數平面上所有 z 點集合


ex: 求單位步階 𝜇[n] 的轉換及其收斂區域

單位步階的定義:\( 𝜇[n] = \left\{ \begin{array}{ll} 1, & \mbox{if t≥0} \\ 0, & \mbox{if t<0} \end{array} \right.\)

\(Z\{𝜇[n]\}=\sum_{n=-∞}^{∞}𝜇[n]z^{-n} = \sum_{n=0}^{∞}z^{-n} = \frac{1}{1-z^{-1}}\)

ROC 為 \( |z^{-1}|<1 \) 或 \( |z|>1\)


ex: 給定單位脈衝時間延遲 δ[n-k] ,其中 k>0,求 z 轉換及其 ROC

\(Z\{δ[n-k]\}=\sum_{n=-∞}^{∞}δ[n-k]z^{-n} = z^{-k}\)

ROC 為複數平面上除了原點以外的所有 z 點集合


ex: 數位訊號 \(x={x[n]}, n=1,2,3,4,5\) 或 \(x={1,2,4,3,2,1}\) ,求z轉換及 ROC

\(X(z)=Z\{x[n]\}=\sum_{n=-∞}^{∞}x[n]z^{-n} \\ = 1 + 2z^{-1}+ 4z^{-2}+ 3z^{-3}+ 2z^{-4}+ z^{-5}\)

ROC 為複數平面上除了原點以外的所有 z 點集合


ex: 數位訊號\(x[n] = (0.5)^n 𝜇[n]\) ,求z轉換及其 ROC

\(X(z)=Z\{x[n]\}=\sum_{n=-∞}^{∞}x[n]z^{-n} \\ = \sum_{n=-∞}^{∞}(0.5)^n𝜇[n]z^{-n} \\ = \sum_{n=0}^{∞}(0.5)^nz^{-n} \\ = \sum_{n=0}^{∞}(0.5z^{-1})^{n} \\ = \frac{1}{1-0.5z^{-1}}\)

要讓無窮等比級數收斂的條件是 \(|r|<1\),因此

ROC 為 \(|0.5z^{-1}|<1\) 或 \(|z|>0.5\)


常見的 z 轉換對應表

離散序列 z轉換 ROC
δ[n] 1 all z
δ[n-k] \(z^{-k}\) z≠0, k>0
𝜇[n] \(\frac{1}{1-z^{-1}}\) 或 \(\frac{z}{z-1}\) $
-𝜇[-n-1] \(\frac{1}{1-z^{-1}}\) 或 \(\frac{z}{z-1}\) $
\(a^n𝜇[n]\) \(\frac{1}{1-az^{-1}}\) 或 \(\frac{z}{z-a}\) $
\(-a^n𝜇[-n-1]\) \(\frac{1}{1-az^{-1}}\) 或 \(\frac{z}{z-a}\) $
\(n𝜇[n]\) \(\frac{z^{-1}}{(1-az^{-1})^2}\) 或 \(\frac{z}{(z-1)^2}\) $
\(n^2𝜇[n]\) \(z^{-1}\frac{(1+z^{-1})}{(1-z^{-1})^3}\) 或 \(\frac{z(z+1)}{(z-1)^3}\) $
\(e^{-an}\) \(\frac{1}{1-e^{-a}z^{-1}}\) 或 \(\frac{z}{z-e^{-a}}\) $
\(sinω_0n\) \(\frac{z sinω_0}{z^2-2zcosω_0+1}\) $
\(cosω_0n\) \(\frac{z(z-cosω_0)}{z^2-2zcosω_0+1}\) $

z轉換的性質

線性運算原則

\(Z\{αx_1[n]+βx_2[n]\} = αX_1(z) + βX_2(z)\)

證明:

\(Z\{αx_1[n]+βx_2[n]\} = \sum_{n=-∞}^{∞}(αx_1[n]+βx_2[n])z^{-n} \\ = α \sum_{n=-∞}^{∞}x_1[n]z^{-n}+β\sum_{n=-∞}^{∞}x_2[n]z^{-n} \\ = αZ\{x_1[n]\} + βZ\{ x_2[n] \} \\ = αX_1(z) + βX_2(z)\)

時間延遲

若 x[n] 的 z 轉換為 X(z),時間延遲 x[n-k] 的 z 轉換為

\(Z\{x[n-k]\}=z^{-k}X(z)\)

證明:

\(Z\{x[n-k]\} = \sum_{-∞}^{∞}x[n-k]z^{-n} \\ = \sum_{j=-∞}^{∞}x[j]z^{-(j+k)} \quad\quad 假設 j=n-k, n=j+k \\ = \sum_{j=-∞}^{∞}x[j]z^{-j}z^{-k)} \\ = z^{-k} \sum_{j=-∞}^{∞}x[j]z^{-j} \\ = z^{-k} X(z)\)

轉換函式

LTI 系統方塊圖

\(y[n]=h[n]*x[n]\) 其中 \(h[n]\) 稱為脈衝響應 (Impulse Response)

卷積定理也成立:

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

或 \(Y(z) = H(z) * X(z)\)

因此轉換函式(Transform Function) 或 系統函式 (System Function)可表示為

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

零點與極點

LTI 的轉換函式,可以用有理式函數 (Rational Function) 的型態表示,讓分子與分母都是 \(z^{-1}\) 的多項式函數

\(H(z) = \frac{b_0+b_1z^{-1}+...+b_Mz^{-M}}{a_0+a_1z^{-1}+...+a_Nz^{-N}}\) 其中的係數包含 \(\{a_k\}, k=0,1,2,...N\) 與 \(\{b_k\}, k=0,1,2,...M\)

系統的轉換函式可因式分解為:

\(H(z) = \frac{b_0}{a_0}z^{N-M}\frac{\prod_{l=1}^{M}(z-z_l)}{\prod_{l=1}^{N}(z-p_l)}\)

讓分子多項式為 0 的所有根,也就是 \(z=z_l, l=1,2,...M\) ,稱為零點 Zeros

讓分母多項式為 0 的所有根,也就是\(z=p_l, l=1,2,...N\),稱為極點 Poles

\(\frac{b_0}{a_0}\) 稱為系統的增益 gain


ex: 如 LTI 系統的轉換函式定義為 \(H(z)=\frac{1}{1-z^{-1}}\) ,求收斂區域、零點、極點

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

ROC 為 \(|z|>1 \),零點為 z=0 (讓分子為0), 極點為 z=1 (讓分母為0)


定理:穩定系統與極點

若 LTI 系統的轉換函式 H(z) 可表示成有理式函數,則 LTI 系統為穩定系統,若且惟若 H(z) 的所有極點都落在複數平面的單位圓內。

換句話說:LTI 要是穩定系統的充要條件為:LTI 的轉換函式,所有極點都要落在複數平面的單位圓內。若條件不成立,則構成不穩定的 LTI 系統。上面範例中,極點落在單位圓上,因此是不穩定系統。


ex: 轉換函式 \(H(z) = \frac{0.8-0.16z^{-1}-0.64z^{-2}}{1-0.2z^{-1}-0.2z^{-2}+z^3}\) ,求系統的零點、極點、增益

可使用 SciPy Signal 的 tf2zpk (Transfer Funciton to Zeros, Poles, Gain)

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
from matplotlib import patches
from matplotlib.markers import MarkerStyle

def zplane(z, p):
    fig = plt.figure( )
    ax = plt.subplot( 1, 1, 1 )

    # 單位圓
    unit_circle = patches.Circle( ( 0,0 ), radius = 1, fill = False, color = 'black', ls = 'dashed' )
    ax.add_patch( unit_circle )

    # 軸線
    plt.axvline( 0, color = 'black' )
    plt.axhline( 0, color = 'black' )
    # 橫軸 -2~2
    plt.xlim( ( -2, 2 ) )
    # 縱軸 -1.5 ~ 1.5
    plt.ylim( ( -1.5, 1.5 ) )
    plt.grid( )

    # 畫上 zeros, ko 為圓圈
    plt.plot( z.real, z.imag, 'ko', fillstyle = 'none', ms = 12 )
    # 畫上 poles, ko 為叉叉
    plt.plot( p.real, p.imag, 'kx', fillstyle = 'none', ms = 12 )
    return fig

def main( ):
    # 分子的係數
    b = np.array( [ 0.8, -0.16, -0.64 ] )
    # 分母的係數
    a = np.array( [ 1, -0.2, -0.2, 1 ] )

    # 用 tf2zpk 找到 zeros, poles, gain
    z, p, k = signal.tf2zpk( b, a )

    print( "Zeros =", z )
    print( "Poles =", p )
    print( "Gain =", k )

    zplane( z, p )
    plt.show( )

main( )

執行結果

$ python tf2zpk.py
Zeros = [ 1.  -0.8]
Poles = [ 0.6+0.8j  0.6-0.8j -1. +0.j ]
Gain = 0.8


反過來可用 zp2tf,根據零點及極點,找到轉換函式的係數

import numpy as np
import scipy.signal as signal

z = np.array( [ -0.8, 1 ] )
p = np.array( [ 0.6 + 0.8j, 0.6 - 0.8j, -1 ] )
k = 0.8

b, a = signal.zpk2tf( z, p, k )

print( "Numerator Polynomial Coefficients =", b )
print( "Denominator Polynomial Coefficients =", a )

執行結果

$ python zp2tf.py
Numerator Polynomial Coefficients = [ 0.8  -0.16 -0.64]
Denominator Polynomial Coefficients = [ 1.  -0.2 -0.2  1. ]

反z轉換

反z轉換定義為:

​ \(x[n]=Z^{-1}\{X(z)\} = \frac{1}{2πj} \oint_CX(z)z^{n-1}dz\) 其中 C 是封閉曲線,落在收斂區域中。

雖然公視是使用封閉曲線積分,實際的反z轉換運算,會採用以下三種方法:

  1. 長除法 long division method
  2. 部分分式展開法 partial fraction expansion method
  3. 餘數法 residue method

長除法

訊號或系統的 z 轉換,通常表示成兩個多項式相除的型態,因此可表示成幂級數 Power Series:

\(X(Z) = \frac{N(z)}{D(z)} = \sum_{n=0}^{∞}a_nz^{-n} = a+a_1z^{-1}+a_2z^{-2}+....\)

其中 N(z) 是 Numerator 分子多項式,D(z) 是 Denominator 分母多項式。

求反 z 轉換時,可以使用長除法求幂級數的係數


ex: 訊號的 z 轉換函式如下,求 反 z 轉換

\(X(z)=\frac{1+z^{-1}+2z^{-2}-z^{-3}+3z{-4}}{1-z^{-1}+z^{-2}}\)

\(\begin{split} \\ &\underline {\ 1 \ +2^z{-1} \ +3^{z-2} \ } \\ 1 \ -z^{-1} \ +z^{-2}\big)& 1 \ +z^{-1} \ +2z^{-2} \ -z^{-3} \ +3z^{-4}\ \\ &\underline{1 \ -z^{-1} \ +z^{-2}\ \ \ \ \ \ }\ \\ & \ \ \ \ 2z^{-1} \ +z^{-2} \ -z^{-3} \ \ \\ &\underline{\ \ \ \ 2z^{-1} \ -2z^{-2} \ +2z^{-3}\ \ \ \ \ \ }\ \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ 3z^{-2} \ -3z^{-3} \ +3z^{-4} \ \ \\ &\underline{\ \ \ \ \ \ \ \ \ \ \ \ \ \ 3z^{-2} \ -3z^{-3} \ +3z^{-4}\ \ \ \ \ \ }\ \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 0 \\ \end{split}\)

因此 X(z) 可化簡為 \(X(z) = 1+2z^{-1}+3z^{-2}\)

求反 z 轉換可得下列結果:

\(x[n]=Z^{-1}\{X(z)\} = Z^{-1}\{1+2z^{-1}+3z^{-2}\} = δ[n]+2δ[n-1]+3δ[n-2]\)

或 \(x[n]=\{1,2,3\}\)


若X(z) 定義為 \( X(z) = \frac{b_0+b_1z^{-1}+...+b_Mz^{-M}}{a_0+a_1z^{-1}+...+a_Nz^{-N}}\)

長除法可用遞迴方法計算:

\(x[n]=[b_n-\sum_{i=1}^{n}x[n-i]a_i]/a_0, n=1,2,...\) ,其中 \(x[0]=b_0/a_0\)

import numpy as np

b = np.array( [ 1, 1, 2, -1, 3 ] )
a = np.array( [ 1, -1, 1, 0, 0 ] )

M = b.size
N = a.size
x = np.zeros( M )
x[0] = b[0] / a[0]
for n in range( 1, M ):
    sum = 0
    k = n
    if n > N:
        k = N
    for i in range( 1, k + 1 ):
        sum = sum + x[n-i] * a[i]
    x[n] = ( b[n] - sum ) / a[0]

print( x )

執行結果

$ python long_division.py
[1. 2. 3. 0. 0.]

部分分式展開法

如果 z 轉換中的分母多項式 D(z) 可進一步因式分解,則可使用部分分式展開法求反z轉換

ex: 訊號的 z 轉換 \(X(z)=\frac{1}{1-3z^{-1}+2z^{-2}}\) ,求反z轉換

分母可分解為 \((1-z^{-1})(1-2z^{-1})\)

因此可假設 \(X(z)=\frac{1}{1-3z^{-1}+2z^{-2}}= \frac{A}{1-z^{-1}} + \frac{B}{1-2z^{-1}}\)

A, B 為常數。通分後可得

\(1=A(1-2z^{-1})+B(1-z^{-1})\)

假設 \(z=1\),則 \(1=A(1-2)+B(1-1) \Rightarrow A=-1\)

假設 \(z=2\),則 \(1=A(1-1)+B(1-1/2) \Rightarrow B=2\)

因此 \(X(z)= \frac{-1}{1-z^{-1}} + \frac{2}{1-2z^{-1}}\)

求反z轉換

\( x[n]=Z^{-1}\{X(z)\} = Z^{-1} \{ \frac{-1}{1-z^{-1}} + \frac{2}{1-2z^{-1}} \} \quad\quad 查表 \\ = -𝜇[n]+2 \cdot (2)^n𝜇[n] = [-1+2^{n+1}]𝜇[n]\)

餘數法

反z轉換定義為:

​ \(x[n]=Z^{-1}\{X(z)\} = \frac{1}{2πj} \oint_CX(z)z^{n-1}dz\) 其中 C 是包含 X(z) 所有極點的封閉曲線,落在收斂區域中。餘數法根據複變分析的柯西餘數定理 Cauchy's Residue Theorem。

柯西餘數定理 Cauchy's Residue Theorem:

\(x[n]=Z^{-1}\{X(z)\} = \frac{1}{2πj} \oint_CX(z)z^{n-1}dz \\ = 封閉曲線C內所有極點,取 z^{n-1}X(z) 的餘數總和\)

若 X(z) 某個極點 \(P_k\),則該極點的餘數 Residue 為:

\( Residue[F(z), P_k] = (z-P_k)F(z)|_{z=P_k} = (z-P_k)z^{n-1}X(z)|_{z=P_k} \)

其中 \(F(z)=z^{n-1}X(z)\)


若訊號 z 轉換函式為 \(X(z)= \frac{z}{(z-0.75)(z+0.5)}\),求反z轉換

因為 \( z=0.75, z= -0.5 \) 讓分母為 0,是 X(z) 的極點

當 \(z=0.75 (P_1=0.75)\) 時

\(Residue[F(z), P_1] = \\ =(z-P_1)z^{n-1}X(z)|_{z=P_1} \\ = (z-0.75)z^{n-1} \frac{z}{(z-0.75)(z+0.5)}|_{z=0.75} \\ = \frac{z^n}{(z+0.5)}|_{z=0.75} = \frac{(0.75)^n}{(0.75+0.5)} = \frac{4}{5}(0.75)^n \)

當 \(z=-0.5 (P_1=-0.5)\) 時

\(Residue[F(z), P_2] = \\ =(z-P_2)z^{n-1}X(z)|_{z=P_2} \\ = (z+0.5)z^{n-1} \frac{z}{(z-0.75)(z+0.5)}|_{z=-0.5} \\ = \frac{z^n}{(z-0.75)}|_{z=-0.5} = \frac{(-0.5)^n}{(-0.5-0.75)} = -\frac{4}{5}(-0.5)^n \)

根據柯西餘數定理,反z轉換為餘數總和:

\(x[n] = Residue[F(z), P_1] + Residue[F(z), P_2]\) 因此反z轉換為

\(x[n]=[\frac{4}{5}(0.75)^n-\frac{4}{5}(-0.5)^n]𝜇[n]\)

References

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

2020/12/28

傅立葉級數及轉換

傅立葉級數 Fourier Series 及轉換 Fourier Transforms 是重要的數學工具,用來進行訊號的頻率分析 frequency analysis。

傅立葉基本理論是:任何週期性函數,可表示成不同頻率、不同振幅的餘弦函數或正弦函數,所加總而得的無窮級數。這個無窮級數就稱為 Fourier Series。

傅立葉級數

基本定義:

若函數 x(t) 在區間 -L ≤ x ≤ L,則傅立葉級數定義為

\( x(t)=\frac{1}{2}a_0 + \sum_{n=1}^{∞}[a_n cos\frac{nπt}{L}+b_n sin\frac{nπt}{L}] \)

其中 \(a_0 = \frac{1}{L} \int_{-L}^{L}x(t)dt\)

​ \(a_n = \frac{1}{L} \int_{-L}^{L}x(t) cos\frac{nπt}{L}dt\)

​ \(b_n = \frac{1}{L} \int_{-L}^{L}x(t) sin\frac{nπt}{L}dt\)

x(t) 是週期性函數,滿足 \(x(t)=x(t+T)\) 其中週期 \(T=2L\),係數 \(a_0, a_n, b_n\) 稱為 x(t) 的傅立葉係數(Fourier Coefficients)。因為是無窮級數,傅立葉級數必須要收斂,才能用來表示 x(t)。

因為 \( T=2L \) 所以傅立葉級數也可以寫成

\( x(t)=\frac{1}{2}a_0 + \sum_{n=1}^{∞}[a_n cos\frac{n2πt}{T}+b_n sin\frac{n2πt}{T}] \)

\( x(t)=\frac{1}{2}a_0 + \sum_{n=1}^{∞}[a_n cos(nωt)+b_n sin(nωt)] \) 其中 \(ω=2πf= \frac{2π}{T}\)


由於 n 是正整數,符合諧波的特性(頻率是正整數倍數),也就是說,傅立葉級數就是無限多個弦波的總和。也就是說,任何週期性函數,都可以表示為無限多個弦波(正弦或餘弦) 組成的諧波。


公式中 \(\frac{1}{2}a_0\) 是固定值常數, \(a_n, b_n\) 可是為第n個弦波的振幅,對應的角頻率為 nω。傅立葉級數也可以用電學的方式表示:

\( x(t)=\frac{1}{2}a_0 + \sum_{n=1}^{∞}[a_n cos(nωt)+b_n sin(nωt)] \) 前面 \(\frac{1}{a_0}\) 是直流分量 DC components,後面另一項是交流分量 AC components


ex: 求函數的傅立葉級數

\( x(t) = \left\{ \begin{array}{ll} -1, & \mbox{if -1<t<0} \\ 1, & \mbox{if 0≤t<1} \end{array} \right.\)

假設 x(t) 是滿足 \(x(t)=x(t+2)\) 的週期函數,會形成 square wave,週期 T=2 (L=1)

\(a_0 = \frac{1}{L} \int_{-L}^{L}x(t)dt = \int_{-1}^{1}x(t)dt = \int_{-1}^{0}(-1)dt + \int_{0}^{1}(1)dt = 0\)

\(a_n = \frac{1}{L} \int_{-L}^{L}x(t) cos\frac{nπt}{L}dt \\ = \int_{-1}^{1}x(t)cos(nωt)dt \\ = \int_{-1}^{0}(-1)cos(nωt)dt + \int_{0}^{1}(1)cos(nωt)dt \\ = [-\frac{1}{nπ}sin(nπt)]_{-1}^{0} + [\frac{1}{nπ}sin(nπt)]_{0}^{1} \\ = 0\)

\(b_n = \frac{1}{L} \int_{-L}^{L}x(t) sin\frac{nπt}{L}dt \\ = \int_{-1}^{1}x(t)sin(nωt)dt \\ = \int_{-1}^{0}(-1)sin(nωt)dt + \int_{0}^{1}(1)sin(nωt)dt \\ = [\frac{1}{nπ}cos(nπt)]_{-1}^{0} + [-\frac{1}{nπ}cos(nπt)]_{0}^{1} \\ = \frac{2}{nπ} - \frac{2}{nπ}cos(nπ) \\ = \frac{2}{nπ}[1-(-1)^n] \)

帶入傅立葉級數

\( x(t)=\frac{1}{2}a_0 + \sum_{n=1}^{∞}[a_n cos\frac{nπt}{L}+b_n sin\frac{nπt}{L}] \\ = \sum_{n=1}^{∞}\frac{2}{nπ}[1-(-1)^n]sin(nπt)\)

展開可得

\( x(t) = = \frac{4}{π}sin(πt) + \frac{4}{3π}sin(3πt)+ \frac{4}{5π}sin(5πt) + ... \) 其中偶數項都是 0

如果取得項數越多,最後加總結果越接近方波

import numpy as np
import matplotlib.pyplot as plt

def series(N, t):
    x = np.zeros( 1000 )            # 方波的傅立葉級數
    for n in range( 1, N + 1 ):
        x += 2 / ( n * np.pi ) * ( 1 - np.power( -1, n ) ) * np.sin( n * np.pi * t )
    return x

def subplot(plotindex, t, x , xlabel, ylabel):
    plt.subplot( plotindex )
    plt.plot( t, x )
    plt.xlabel( xlabel )
    plt.ylabel( ylabel )

t = np.linspace( -1, 1, 1000 )  # 定義時間陣列

N = 1
x = series(N, t)

plt.figure( 1 )
subplot(321, t, x, 't (second)', 'Amplitude '+str(N) )

N = 5
x = series(N, t)
subplot(322, t, x, 't (second)', 'Amplitude '+str(N) )

N = 10
x = series(N, t)
subplot(323, t, x, 't (second)', 'Amplitude '+str(N) )

N = 50
x = series(N, t)
subplot(324, t, x, 't (second)', 'Amplitude '+str(N) )

N = 100
x = series(N, t)
subplot(325, t, x, 't (second)', 'Amplitude '+str(N) )

N = 200
x = series(N, t)
subplot(326, t, x, 't (second)', 'Amplitude '+str(N) )

plt.show( )

在方波的不連續處,可發現明顯地跳動現象,稱為 Gibbs Phenomenon 現象。在分段連續 piecewise-continuous 函數 (ex: 鋸齒波、三角波),如果以傅立葉級數表示時,在不連續處都會發生 Gibbs

傅立葉轉換

傅立葉級數可將任意週期性函數都可分解成不同振幅、不同頻率的弦波,因此是頻率分析 frequency analysis 的重要工具。

傅立葉轉換 Fourier Transforms 延伸了傅立葉級數,可同時適用於週期及非週期性連續函數

基本定義:函數 x(t) 的 Fourier Transform 定義為

\(X(ω) = F\{x(t)\} = \int_{-∞}^{∞}x(t)e^{-jωt}dt\)

反傅立葉轉換 Inverse Fourier Transform 定義為

\(x(t)=F^{-1}\{X(ω)\} = \frac{1}{2π}\int_{-∞}^{∞}X(ω)e^{jωt}dt\) 其中 \( \int_{-∞}^{∞}|x(t)|dt \) 收斂,\(j=\sqrt{-1}\)


Fourier Transform 可將時間函數 x(t) 轉換成頻率函數 X(ω)。反傅立葉轉換 Inverse Fourier Transform 可將 X(ω) 還原為 x(t)

\( x(t) → Fourier Transform → X(ω) \)

\( x(t) ← Inverse Fourier Transform ← X(ω) \)


ex: 弦波 \(x(t)=Acos(ω_0t)\) ,其振幅為 A,角頻率 \(ω_0\),相位移 \(φ=0\),求傅立葉轉換

\(X(ω) = F\{x(t)\} = \int_{-∞}^{∞}x(t)e^{-jωt}dt \\ = \int_{-∞}^{∞}Acos(ω_0t)e^{-jωt}dt \\ = A \int_{-∞}^{∞}(\frac{e^{jω_0t}+e^{-jω_0t}}{2})e^{-jωt}dt \\ = \frac{A}{2} [\int_{-∞}^{∞}e^{jω_0t}e^{-jωt}dt + \int_{-∞}^{∞}e^{-jω_0t}e^{-jωt}dt] \\ = \frac{A}{2} [F(e^{jω_0t}) + F(e^{-jω_0t})] \\ = Aπ[δ(ω-ω_0) + δ(ω+ω_0)] \)

其中 \(F(e^{jω_0t})=2πδ(ω-ω_0)\)


Euler's Formula:

\( e^{j𝜃} = cos𝜃+jsin𝜃, j=\sqrt{-1} \)

反尤拉公式

\(cos𝜃 = \frac{e^{j𝜃}+e^{-j𝜃}}{2} \)

\(sin𝜃 = \frac{e^{j𝜃}-e^{-j𝜃}}{2j} \)


ex: Ideal Pulse Function (非週期性函數) 的傅立葉轉換

\( x(t) = \left\{ \begin{array}{ll} A, & \mbox{if -T/2<t<T/2} \\ 0, & \mbox{otherwise} \end{array} \right.\)

\(X(ω) = F\{x(t)\} = \int_{-∞}^{∞}x(t)e^{-jωt}dt \\ = \int_{-T/2}^{T/2}Ae^{-jωt}dt \\ = A \int_{-T/2}^{T/2}e^{-jωt}dt \\ = A [-\frac{1}{jω}e^{-jωt}]_{-T/2}^{T/2} \\ = A [-\frac{1}{jω}e^{-\frac{jωT}{2}} + \frac{1}{jω}e^{\frac{jωT}{2}} ] \\ = \frac{A}{jω} [e^{\frac{jωT}{2}} - e^{-\frac{jωT}{2}} ] \\ = \frac{A}{jω} [ cos(\frac{ωT}{2}) + j sin(\frac{ωT}{2}) - cos(\frac{ωT}{2}) + j sin(\frac{ωT}{2}) ] \\ = \frac{A}{ω} [2 sin(\frac{ωT}{2})] \\ = AT [\frac{sin(\frac{ωT}{2})}{\frac{ωT}{2}}] \\ = AT sinc(\frac{ωT}{2π}) \)

數學的 sinc 定義為 \(sinc(x)= \frac{sin(x)}{x}\)

工程數學的 sinc 定義為 \(sinc(x)= \frac{sin(πx)}{πx}\)

當 \(ω=0\) 時,需要使用 L'Hopital's Rule 羅必達規則,可以求出特定函數趨近於某數的極限值

\(X(ω=0) = \lim_{w→0}AT \cdot sinc(\frac{ωT}{2}) \\ = AT \cdot \lim_{ω→0} [ \frac{sin(\frac{ωT}{2})}{(\frac{ωT}{2})} ] \\ = AT \cdot \lim_{ω→0} [ \frac{ \frac{d}{dω}sin(\frac{ωT}{2}) }{ \frac{d}{dω}(\frac{ωT}{2}) } ] \\ = AT \cdot \lim_{ω→0} [ \frac{\frac{T}{2} \cdot cos(\frac{ωT}{2}) }{\frac{T}{2}} ] \\ = AT \)

import numpy as np
import matplotlib.pyplot as plt

def ideal_pulse_fourier_transform(w, A, T):
    X = A * T * np.sinc( w * T / ( 2 * np.pi ) )
    return X

A = 1
T = 2
w = np.linspace( -20, 20, 1000 )

X = ideal_pulse_fourier_transform(w, A, T)

plt.subplot( 131 )
plt.plot( w, X )
plt.xlabel( r'$\omega$' + ', T=' + str(T) )
plt.ylabel( r'X($\omega$)')

A = 1
T = 5
X = ideal_pulse_fourier_transform(w, A, T)

plt.subplot( 132 )
plt.plot( w, X )
plt.xlabel( r'$\omega$' + ', T=' + str(T) )
plt.ylabel( r'X($\omega$)')

A = 1
T = 10
X = ideal_pulse_fourier_transform(w, A, T)

plt.subplot( 133 )
plt.plot( w, X )
plt.xlabel( r'$\omega$' + ', T=' + str(T) )
plt.ylabel( r'X($\omega$)')

plt.show( )

當 T 越大,傅立葉轉換後的 sinc 函數週期變小

平移定理

Fourier Transform Shifting Theorems

  • 第一平移定理(時間平移定理 Time-Shifting Theorem)

    f 為時間函數,F{} 為傅立葉轉換,則

    \(F\{f(t-t_0)\} = F(ω) \cdot e^{jωt_0}\),其中 \(t_0\) 為平移的時間 且 \(t_0>0\)

    函數 f 在時間域平移,若取其傅立葉轉換,則結果相當於原始函數的傅立葉轉換,乘上一個複數指數函數

​ 證明:

​ 基本定義:函數 x(t) 的 Fourier Transform 定義為

​ \(X(ω) = F\{x(t)\} = \int_{-∞}^{∞}x(t)e^{-jωt}dt\)

​ 根據定義:

​ \( F\{f(t-t_0)\} = \int_{-∞}^{∞}f(t-t_0)e^{-jωt}dt \)

​ 假設 \( τ = t-t_0 \) 則 \( t = τ-t_0, dτ=dt\),因此

​ \( 原式 = \int_{-∞}^{∞}f(τ)e^{-jω(t-t_0)}dτ \\ ​ = \int_{-∞}^{∞}f(τ)e^{-jωt} \cdot e^{jωt_0} dτ \\ ​ = \{\int_{-∞}^{∞}f(τ)e^{-jωt} dτ\} \cdot e^{jωt_0} \\ ​ = F(ω) \cdot e^{jωt_0} \)

​ 上述的第一平移定理,也可以表示成

​ \(F\{f(t+t_0)\} = F(ω) \cdot e^{-jωt_0}\)

  • 第二平移定理(頻率平移定理 Frequency-Shifting Theorem)

    f 為時間函數,F{} 為傅立葉轉換,則

    \(F\{f(t)\cdot e^{jω_0t}\} = F(ω-ω_0) \),其中 \(ω_0\) 為平移的角頻率 且 \(ω_0>0\)

    函數 f 的傅立葉轉換,其在頻率域的平移,則結果相當於原始的時間函數,乘上一個複數指數函數

​ 證明:

​ 基本定義:函數 x(t) 的 Fourier Transform 定義為

​ \(F(ω) = F\{x(t)\} = \int_{-∞}^{∞}x(t)e^{-jωt}dt\)

​ 根據定義:

​ \( F\{f(t)\cdot e^{jω_0t}\} = \int_{-∞}^{∞}f(t)e^{jω_0t}e^{-jωt}dt = \int_{-∞}^{∞}f(t)e^{j(ω-ω_0)t}dt \)

​ 跟傅立葉轉換的定義比較後,得到

​ \( F\{f(t)\cdot e^{jω_0t}\} = \int_{-∞}^{∞}f(t)e^{-j(ω-ω_0)t}dt = F(ω-ω_0) \)

​ 上述第二平移定理也可以表示為

​ \(F\{f(t)\cdot e^{-jω_0t}\} = F(ω+ω_0) \)

  • 卷積定理 (Convolution Theorem)

    時間域的卷積運算結果,與其在頻率域的點對點乘法運算結果相同

​ 若 f, g 為兩個時間函數, F{} 為傅立葉轉換,則

​ \( F\{ f*g\} = F\{f\} \cdot F\{g\} \) 其中 \(*\) 為卷積運算

​ DSP 領域中,卷積定理是很重要的定理,也是在訊號處理時,可採用時間域及頻率域兩種不同方法的理論依據

​ 證明:

​ \( F\{ f*g\} = \int_{-∞}^{∞}[\int_{-∞}^{∞}f(τ)g(t-τ)dτ] e^{-jωt} dt \\ ​ = \int_{-∞}^{∞} f(τ)[\int_{-∞}^{∞}g(t-τ)e^{-jωt}dt] dτ \\ % 利用Fubini's theorem ​ = \int_{-∞}^{∞} f(τ)[\int_{-∞}^{∞}g(ť)e^{-jω(ť+τ)}dt] dτ \quad\quad ( 假設 ť = t-τ, dť = dt ) \\ ​ = \int_{-∞}^{∞} f(τ)e^{-jωt}[\int_{-∞}^{∞}g(ť)e^{-jωť}dt] dτ \\ ​ = G(ω) \int_{-∞}^{∞}f(τ)e^{-jωτ}dτ \\ ​ = F(ω) \cdot G(ω) \\ ​ = F\{f\} \cdot G\{g\} \)

​ Fubini's theorem: 富比尼定理給出了使用逐次積分的方法計算雙重積分的條件。在這些條件下,不僅能夠用逐次積分計算雙重積分,而且交換逐次積分的順序時,積分結果不變。

  • 高斯函數(常態分佈) 的傅立葉轉換

    高斯函數的傅立葉轉換會形成另一個高斯函數。兩個高斯函數的標準差呈現反比關係。如果高斯函數在時間域的標準差越大,其在頻率域的高斯函數的標準差越小。

    \( F\{e^{-\frac{t^2}{2σ^2}}\} = \sqrt{2πσ^2} e^{-\frac{1}{2}σ^2ω^2} \)

​ 證明:

​ 高斯函數的傅立葉轉換為

​ \( X(ω) = F\{x(t\} = \int_{-∞}^{∞}x(t)e^{-jωt}dt \\ ​ = \int_{-∞}^{∞}e^{-\frac{r^2}{2σ^2}}e^{-jωt}dt \\ ​ = \int_{-∞}^{∞}e^{-(\frac{r^2}{2σ^2}+jωt)}dt \)

​ 將指數的幂次方配成平方型態

​ \( \frac{r^2}{2σ^2}+jωt = (\frac{t}{\sqrt{2σ^2}}+\frac{\sqrt{2σ^2}}{2}jω)^2+\frac{1}{2}σ^2ω^2 \)

​ \( 原式 = \int_{-∞}^{∞}e^{-(\frac{t}{\sqrt{2σ^2}}+\frac{\sqrt{2σ^2}}{2}jω)^2} \cdot e^{-\frac{1}{2}σ^2ω^2} dt \\ ​ = e^{-\frac{1}{2}σ^2ω^2} \cdot \int_{-∞}^{∞}e^{-(\frac{t}{\sqrt{2σ^2}}+\frac{\sqrt{2σ^2}}{2}jω)^2} dt \)

​ 假設 \( 𝜇=\frac{t}{\sqrt{2σ^2}}+\frac{\sqrt{2σ^2}}{2}jω \) ,則 \(d𝜇=\frac{1}{\sqrt{2𝜎^2}}dt, dt= \sqrt{2𝜎^2}d𝜇\)

​ \(原式 = e^{-\frac{1}{2}σ^2ω^2} \cdot \int_{-∞}^{∞}e^{-𝜇^2}\sqrt{2𝜎^2}d𝜇 \\ ​ = \sqrt{2𝜎^2} \cdot e^{-\frac{1}{2}σ^2ω^2} \cdot \int_{-∞}^{∞}e^{-𝜇^2} d𝜇 \\ ​ = \sqrt{2𝜎^2} \cdot e^{-\frac{1}{2}σ^2ω^2} \cdot 2 \int_{0}^{∞}e^{-𝜇^2} d𝜇 \quad\quad (其中已知 \int_{0}^{∞}e^{-𝜇^2} d𝜇 = \frac{\sqrt{π}}{2} ) \\ ​ = \sqrt{2𝜎^2} \cdot e^{-\frac{1}{2}σ^2ω^2} \cdot 2 \cdot \frac{\sqrt{π}}{2} \\ ​ = \sqrt{2πσ^2} e^{-\frac{1}{2}σ^2ω^2} \)


如果高斯函數定義為 \(x(t)= \frac{1}{\sqrt{2πσ^2}}e^{-\frac{t^2}{2σ^2}} \) 並滿足以下積分條件 \(\int_{-∞}^{∞}x(t)dt = 1 \) ,也就是機率總和為 1 ,則其傅立葉轉換為:

\(X(ω) = e^{-\frac{1}{2}σ^2ω^2} \)


ex: 高斯函數 \( x(t) = e^{-\frac{t^2}{2σ^2}} \) ,其標準差為 1, 2, 3,計算其時間域的高斯函數,以及在頻率域的傅立葉轉換 \(X(ω)\)

import numpy as np
import matplotlib.pyplot as plt

def gaussian(t, sigma):
    x = np.exp( - ( t * t ) / ( 2 * sigma * sigma ) )
    return x

def gaussian_fourier_transform(w, sigma ):
    w = np.linspace( -7, 7, 1000 )
    X = np.exp( - ( sigma * sigma * w * w ) / 2 )
    X = np.sqrt( 2 * np.pi * sigma * sigma ) * X
    return X

def plot(x, y, subplot, xlabel, ylabel):
    plt.subplot( subplot )
    plt.plot( x, y )
    plt.xlabel( xlabel )
    plt.ylabel( ylabel )


sigma = 1
# 時間域的高斯函數 x
t = np.linspace( -7, 7, 100 )
x = gaussian(t, sigma)
# 頻率域的傅立葉轉換
w = np.linspace( -7, 7, 1000 )
X = gaussian_fourier_transform( w, sigma )

plot(t, x, 231, 't'+", sigma="+str(sigma), 'x(t)' )
plot(w, X, 234, r'$\omega$'+", sigma="+str(sigma), r'X($\omega$)' )


sigma = 2
# 時間域的高斯函數 x
t = np.linspace( -7, 7, 100 )
x = gaussian(t, sigma)
# 頻率域的傅立葉轉換
w = np.linspace( -7, 7, 1000 )
X = gaussian_fourier_transform( w, sigma )

plot(t, x, 232, 't'+", sigma="+str(sigma), 'x(t)' )
plot(w, X, 235, r'$\omega$'+", sigma="+str(sigma), r'X($\omega$)' )

sigma = 3
# 時間域的高斯函數 x
t = np.linspace( -7, 7, 100 )
x = gaussian(t, sigma)
# 頻率域的傅立葉轉換
w = np.linspace( -7, 7, 1000 )
X = gaussian_fourier_transform( w, sigma )

plot(t, x, 233, 't'+", sigma="+str(sigma), 'x(t)' )
plot(w, X, 236, r'$\omega$'+", sigma="+str(sigma), r'X($\omega$)' )

plt.show( )

離散時間傅立葉轉換

Discrete-Time Fourier Transform (DTFT),目標是將離散的序列 \({x[n]}\) 轉換為複數指數 \(\{e^{-jωn}\}\) 的序列

給定離散的序列 x[n],DTFT 定義為

​ \( X(e^{jω}) = \sum_{n=-∞}^{∞}x[n]e^{-jωn} \)

反離散時間傅立葉轉換 Inverse DTFT 定義為

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


\(X(e^{jω})\) 是 ω 的複數函數,可表示成 \( X(e^{jω}) = |X(e^{jω})| \cdot e^{jθ(ω)} \),其中

​ \(|X(e^{jω})|\) 稱為強度 Magnitude

​ \(θ(ω)\) 稱為幅角 Argument 或相位角 Phase Angle

如果以圖形表示,則分別是 強度頻譜 Magnitude Spectrum,與相位頻譜 Phase Spectrum

如果離散序列是絕對可加總序列 (Absolutely Summable Sequence) \(\sum_{n=-∞}^{∞}|x[n]|<∞\),則 DTFT 的無窮級數才會收斂。

換句話說,必須滿足以下條件,DTFT 才會存在:

​ \(|X(e^{jω})|<∞ \) , for all ω


\(X(e^{jω})\) 是 ω 的複數函數,同時也是週期函數,週期為 2π

證明:

\(X(e^{j(ω+2πk)}) = \sum_{n=-∞}^{∞}x[n] \cdot e^{-j(ω+2πk)n} \\ = \sum_{n=-∞}^{∞}x[n] \cdot e^{-jωn}e^{-j2πkn} \\ = \sum_{n=-∞}^{∞}x[n] \cdot e^{-jωn} = X(e^{jω}) \)


ex: 數位訊號 \(x[n] = δ[n]\) ,單位脈衝,求其 DTFT

單位脈衝的定義:\( δ[n] = \left\{ \begin{array}{ll} 1, & \mbox{if n=0} \\ 0, & \mbox{if n≠0} \end{array} \right.\)

DTFT \(X(e^{jω}) = \sum_{-∞}^{∞}x[n]e^{-jωn} = \sum_{-∞}^{∞}δ[n]e^{-jωn} = δ[0]e^{-jω0}= 1\)


ex: 數位訊號 \(x[n] = (0.5)^n 𝜇[n]\) ,單位步階,求其 DTFT

單位步階的定義:\( 𝜇[n] = \left\{ \begin{array}{ll} 1, & \mbox{if t≥0} \\ 0, & \mbox{if t<0} \end{array} \right.\)

DTFT \(X(e^{jω}) = \sum_{-∞}^{∞}x[n]e^{-jωn} = \sum_{-∞}^{∞}(0.5)^n 𝜇[n]e^{-jωn} = \sum_{0}^{∞}(0.5)^n e^{-jωn} = \sum_{0}^{∞}(0.5e^{-jωn})^n = \frac{1}{1-0.5 e^{-jω}} \)

這個 DTFT 存在的條件為 \( \sum_{n=-∞}^{∞}|x[n]| < ∞ \)


離散傅立葉轉換

Discrete Fourier Transform (DFT)

給定離散序列 \(x[n], n=0,1,2..,N-1\) ,則 DFT 定義為

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

反 DFT 為:

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


跟 DTFT 定義 \( X(e^{jω}) = \sum_{n=-∞}^{∞}x[n]e^{-jωn} \) 比較後可發現,離散傅立葉轉換以有限的離散序列為主,且在 ω 軸進行取樣

​ \( X[k]=X(e^{jω})|_{ω=2πk/N} = \sum_{n=0}^{N-1}x[n]e^{-j2πkn/N} \)


DFT 符合可逆性,N個樣本數的 x[n] 經過 DFT 後,可得到 N 個樣本的 X[k]。X[k] 經過反轉換,可還原為 x[n]


離散傅立葉轉換通常會用這個方式表示:\(W_N=e^{-j(2π/N)}\)

離散傅立葉轉換也可以表示成

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

反轉換為

\(x[n]= \frac{1}{N}\sum_{k=0}^{N-1}X[k]W_N^{kn}, n=0,1,...,N-1 \)


ex: 數位訊號定義為 \(x=\{x[n]\}, n=0,1,2,3\) 或 \(x=\{1,2,4,3\}, n=0,1,2,3\) ,其中 N=4,求 DFT

DFT 定義為

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

\(X[0]= \sum_{n=0}^{3}x[n]=1+2+4+3=10\)

\(X[1]= \sum_{n=0}^{3}x[n]e^{-jπn/2} = 1 \cdot e^0+2 \cdot e^{-j(π/2)} +4 \cdot e^{-jπ}+3 \cdot e^{-j(3π/2)} \\ = 1 \cdot 1 + 2 \cdot (-j) + 4 \cdot (-1) + 3 \cdot j = -3+j\)

\(X[2]= \sum_{n=0}^{3}x[n]e^{-jπn} = 0 \)

\(X[3]= \sum_{n=0}^{3}x[n]e^{-jπ3n/2} = -3-j \)

結果為

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

本範例中,X[0] 稱為直流分量 DC component, X[1]~X[3] 稱為交流分量 (AC Compomenets),另外 X[1] 與 X[3] 為共軛複數


DFT 矩陣

DFT 可用矩陣的方式表示

\(X = D_N \cdot x\)

\(x = [x[0], x[1],...., x[N-1]]^T \) 為輸入向量

\(X = [X[0], X[1],...., X[N-1]]^T \) 為輸出向量

\(D_N\) 為 NxN 的 DFT 轉換矩陣,定義為

\( D_N = \begin{bmatrix} 1 & 1 & 1 & ... & 1 \\ 1 & W_N^1 & W_N^2 & ... & W_N^{N-1} \\ 1 & W_N^2 & W_N^4 & ... & W_N^{2(N-1)} \\ . & . & . & ... & . \\ . & . & . & ... & . \\ 1 & W_N^{N-1} & W_N^{2(N-1)} & ... & W_N^{(N-1) \cdot (N-1)} \\ \end{bmatrix} \)

反離散傅立葉轉換 Inverse Discrete Fourier Transform, Inverse DFT,定義為

\(x = D_N^{-1} \cdot X\)

\(D_N^{-1}\) 為 NxN 的 DFT 反轉換矩陣,定義為

\( D_N^{-1} = \begin{bmatrix} 1 & 1 & 1 & ... & 1 \\ 1 & W_N^{-1} & W_N^{-2} & ... & W_N^{-(N-1)} \\ 1 & W_N^{-2} & W_N^{-4} & ... & W_N^{-2(N-1)} \\ . & . & . & ... & . \\ . & . & . & ... & . \\ 1 & W_N^{-(N-1)} & W_N^{-2(N-1)} & ... & W_N^{-(N-1) \cdot (N-1)} \\ \end{bmatrix} \)

因為 \(W_k=e^{-j(2π/N)}\)

如果 N=4,則 \(W_N^k, k=0,1,2,3\) 分別為

\(W_4^0 = e^0 =1 \) \(W_4^{1} = e^{-j(2π/4)} = -j\)

\(W_4^{2} = e^{-j(4π/4)} = -1\)

\(W_4^{3} = e^{-j(6π/4)} = j\)

\(W_N^k\) 落在單位圓上,依順時針方向旋轉。反傅立葉轉換則是依逆時針方向旋轉。


ex: 若數位訊號定義為 \(x={x[n]}, n=0,1,2,3\) 或 \(x={1,2,4,3}, n=0,1,2,3\),其中 N =4,則 DFT 為?

\( X = \begin{bmatrix} 1 & 1 & 1 & 1 \\ 1 & W_4^1 & W_4^2 & W_4^3 \\ 1 & W_4^2 & W_4^4 & W_4^6 \\ 1 & W_4^3 & W_4^6 & W_4^9 \\ \end{bmatrix}x \\ = \begin{bmatrix} 1 & 1 & 1 & 1 \\ 1 & -j & -1 & j \\ 1 & -1 & 1 & -1 \\ 1 & j & -1 & -j \\ \end{bmatrix} \begin{bmatrix} 1 \\ 2 \\ 4 \\ 3 \\ \end{bmatrix} \\ = \begin{bmatrix} 10 \\ (-3+j) \\ 0 \\ (-3-j) \\ \end{bmatrix}\)


根據 DFT 定義,計算的時間複雜度為 \(O(N^2)\),實際使用時,N 值很大,導致 DFT 計算量龐大。

J. W. Cooley, John Turkey 修改 DFT 演算法,改用 divide-and-conquer 方法提出快速傅立葉轉換 Fast Fourier Transforms (FFT)。FFT 可得到與 DFT 相同的結果,時間複雜度為 \(O(N \cdot log_2N)\)


ex: DFT 具有可逆性

如果 \( X=\{10, -3+j, 0, -3-j\} \) 其中 N=4, 求其 Inverse DFT ?

可根據定義計算,或是用 Inverse DFT Matrix 計算

反 DFT 為:\(x[n]= \frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{-j2πkn/N}, n=0,1,...,N-1\)

\( x = \frac{1}{4}\begin{bmatrix} 1 & 1 & 1 & 1 \\ 1 & W_4^{-1} & W_4^{-2} & W_4^{-3} \\ 1 & W_4^{-2} & W_4^{-4} & W_4^{-6} \\ 1 & W_4^{-3} & W_4^{-6} & W_4^{-9} \\ \end{bmatrix}X \\ = \frac{1}{4}\begin{bmatrix} 1 & 1 & 1 & 1 \\ 1 & j & -1 & -j \\ 1 & -1 & 1 & -1 \\ 1 & -j & -1 & j \\ \end{bmatrix} \begin{bmatrix} 10 \\ (-3+j) \\ 0 \\ (-3-j) \\ \end{bmatrix} \\ = \begin{bmatrix} 1 \\ 2 \\ 4 \\ 3 \\ \end{bmatrix}\)


import numpy as np
from numpy.fft import fft, ifft

x = np.array( [ 1, 2, 4, 3 ] )
X = fft( x )
Xm = abs( X )
xx = ifft ( X )

print( "x =", x )
print( "X =", X )
print( "Magnitude of X =", Xm )
print( "Inverse FFT of X =", xx )

執行結果

$ python FFT_example.py
x = [1 2 4 3]
X = [10.+0.j -3.+1.j  0.+0.j -3.-1.j]
Magnitude of X = [10.          3.16227766  0.          3.16227766]
Inverse FFT of X = [1.+0.j 2.+0.j 4.+0.j 3.+0.j]

SciPy 也有提供 fft, ifft 函數

# from numpy.fft import fft, ifft
from scipy.fftpack import fft, ifft

References

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