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

2020/12/21

DSP 相關 Correlation

包含交互相關 cross-correlation 以及 自相關 autocorrelation 及其應用。

correlation 相關就是兩個隨機變數之間,是否具有關聯性。相關係數可用來測量兩個隨機變數集合的相關性,通常介於 -1 (負相關) ~ +1 (正相關) 之間。

交互相關 cross correlation

給定兩個函數 \(x(t), h(t)\),連續時間的 cross correlation 定義為:

\(y(t) = \int_{-∞}^{∞}x^*(τ)h(t+τ)dτ\) 或

\(y(t) = x(t) ⊗ h(t)\)

其中 \(x^*(τ)\) 是共軛複數 complex conjugate 。但數位訊號的輸入資料都是實數,所以 cross correlation 跟卷積運算類似,只是把 \(h(t-τ)\) 換成 \(h(t+τ)\)

數位訊號 x[n], h[n] 的交互相關定義為:

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

\( y[n]=x[n] ⊗ h[n] \)

h[n] 可視為第二個數位訊號,用來進行訊號比對

因輸入訊號的數量 N=7,脈衝響應數量為 M = 5,full cross-correlation 結果長度為 \(M+N-1 = 11\)

計算前先將 \(x[n]\) 兩邊補上 \(M-1=4\) 個 0,稱為 Zero-Padding,但不像卷積運算一樣,不將 h(n) 旋轉 180 度

x h
0 0 0 0 1 2 4 3 2 1 1 0 0 0 0 1 2 3 1 1
1 2 3 1 1
\(y[0] = 0 \cdot 1 + 0 \cdot 2 + 0 \cdot 3 + 0 \cdot 1 + 1 \cdot 1 = 1 \)
x h
0 0 0 0 1 2 4 3 2 1 1 0 0 0 0 1 2 3 1 1
1 2 3 1 1
\(y[1] = 0 \cdot 1 + 0 \cdot 2 + 0 \cdot 3 + 1 \cdot 1 + 2 \cdot 1 = 3 \)

最後卷積運算結果為

n 0 1 2 3 4 5 6 7 8 9 10
y[n] 1 3 9 15 22 22 18 11 7 3 1

實際運算時,會擷取卷積運算的部分結果

n 0 1 2 3 4 5 6
x[n] 1 2 4 3 2 1 1
y[n] 9 15 22 22 18 11 7

結果在 n=2, 3 的位置,數值最大,也就是 x[n] 與 h[n] 在這個位置的相關性最高,波形最相似。

cross correlation 經典的應用是訊號比對 Signal Matching,例如:音樂搜尋系統,聲紋比對系統

numpy 有提供 full, same 兩種 cross correlation 運算的結果

import numpy as np

x = np.array( [ 1, 2, 4, 3, 2, 1, 1 ] )
h = np.array( [ 1, 2, 3, 1, 1 ] )
y = np.correlate( x, h, 'full' )
y1 = np.correlate( x, h, 'same' )

print( "x =", x )
print( "h =", h )
print( "Full Correlation y =", y )
print( "Correlation y =", y1 )

執行結果

$ python correlation.py
x = [1 2 4 3 2 1 1]
h = [1 2 3 1 1]
Full Correlation y = [ 1  3  9 15 22 22 18 11  7  3  1]
Correlation y = [ 9 15 22 22 18 11  7]

如果 h[n] 左右對稱且輸入的數位訊號是實數,則 cross correlation 的結果跟卷積運算一樣。另外 cross correlation 跟卷積一樣,滿足交換率、結合率、分配率。

自相關 autocorrelation

給定兩個函數 \(x(t), h(t)\),連續時間的 autocorrelation 定義為:

\(y(t) = \int_{-∞}^{∞}x^*(τ)x(t+τ)dτ\) 其中 τ 稱為延遲 lag

其中 \(x^*(τ)\) 是共軛複數 complex conjugate ,自相關 autocorrelation 是延遲 τ 的函數。但數位訊號的輸入資料都是實數,所以 autocorrelation 可解釋為:訊號與訊號本身的延遲訊號,進行積分運算的結果。

數位訊號的卷積定義為:

\( R[l] = \sum_{k=-∞}^{∞}x^*[k] \cdot x[n+l] \) 其中 l 是延遲 lag

autocorrelation 是延遲 l 的函數,且 l 為正整數。可用來偵測數位訊號是否具有週期性。通常週期性的數位訊號,自關性較強。 noise 因為每個 sample 在時間軸上都是獨立的,所以不具有 autocorrelation

因輸入訊號的數量 N=5,脈衝響應數量為 M = 5,autocorrelation 結果長度為 5

計算前先將 \(x[n]\) 後面補上 \(N-1=4\) 個 0

x h
1 2 1 2 1 0 0 0 0 1 2 1 2 1
1 2 1 2 1
\(R[0] = 1 \cdot 1 + 2 \cdot 2 + 1 \cdot 1 + 2 \cdot 2 + 1 \cdot 1 = 11 \)
x h
1 2 1 2 1 0 0 0 0 1 2 1 2 1
1 2 1 2 1
\(y[1] = 2 \cdot 1 + 1 \cdot 2 + 2 \cdot 1 + 1 \cdot 2 + 0 \cdot 1 = 8 \)

最後 autocorrelation 運算結果為

n 0 1 2 3 4
R[n] 11 8 6 4 1

cross correlation 運算結果的後半部分,就是 autocorrelation 的結果

import numpy as np

def autocorr( x ):
    # cross correlation 運算結果的後半部分,就是 autocorrelation
    R = np.correlate( x, x, 'full' )
    return R[ int( R.size / 2 ) : ]

def main( ):
    x = np.array( [ 1, 2, 1, 2, 1 ] )
    R = autocorr( x )
    print( "x =", x )
    print( "Autocorrelation =", R )

main( )

執行結果

$ python autocorrelation.py
x = [1 2 1 2 1]
Autocorrelation = [11  8  6  4  1]

autocorrelation 的應用

週期性訊號的自相關性較強,雜訊不具有自相關性。自相關適合用來偵測訊號的週期性,即使週期性的訊號混入了雜訊,一樣能找出自相關性。

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

def autocorr( x ):
    R = np.correlate( x, x, 'full' )
    return R[ int( R.size / 2 ) : ]

def main( ):
    t = np.linspace( 0, 1,  100, endpoint = False )     # 定義時間陣列
    x = 10 * np.cos( 2 * np.pi * 5 * t )                # 原始訊號
    noise = random.uniform( -2, 2, 100 )                # 雜訊(均勻分佈)  
    y = x + noise

    auto_corr1 = autocorr( x )
    auto_corr2 = autocorr( noise )
    auto_corr3 = autocorr( y )

    plt.figure( 1 )                                     
    plt.subplot( 121 )
    plt.plot( t, x )
    plt.xlabel( 't (second)' )
    plt.ylabel( 'Amplitude' )

    plt.subplot( 122 )                                  
    plt.plot( auto_corr1 )
    plt.xlabel( 'Lag' )
    plt.ylabel( 'Auto Correlation' )

    plt.figure( 2 )
    plt.subplot( 121 )
    plt.plot( t, noise )
    plt.xlabel( 't (second)' )
    plt.ylabel( 'Amplitude' )

    plt.subplot( 122 )
    plt.plot( auto_corr2 )
    plt.xlabel( 'Lag' )
    plt.ylabel( 'Auto Correlation' )

    plt.figure( 3 ) 
    plt.subplot( 121 )
    plt.plot( t, y )
    plt.xlabel( 't (second)' )
    plt.ylabel( 'Amplitude' )

    plt.subplot( 122 )
    plt.plot( auto_corr3 )
    plt.xlabel( 'Lag' )
    plt.ylabel( 'Auto Correlation' )

    plt.show( )
    
main( )

左圖 x[n] 是原始訊號,右邊是 autocorrelation 的結果

雜訊,及其 autocorrelation

弦波加上雜訊的原始訊號,右圖是 autocorrelation

References

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

2020/12/14

DSP 卷積 Convolution

Convolution 最基本的應用是訊號的 filtering,以下會討論兩種最常見的 filtering:average filters 及 Gaussian filters

卷積運算

Convolution 可解釋為:針對兩個時間函數進行數學運算,產生另一個時間函數的過程。

給定兩個函數 \(x(t), h(t)\),連續時間的卷積定義為:

\(y(t) = \int_{-∞}^{∞}x(τ)h(t-τ)dτ = \int_{-∞}^{∞}h(τ)x(t-τ)dτ\) 或

\(y(t) = x(t) * h(t)\)

卷積運算符合交換率,線性與時間不變性原則,是 LTI 系統。


數位訊號的卷積定義為:

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

\( y[n]=x[n]*h[n] \)

離散的卷積運算,也符合交換率

因為任意數位訊號都可以用單位脈衝 Unit Impulse表示

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

代入 DSP 系統,可得到

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

跟卷積的定義比較

\(h[n-k] = T\{δ[n-k]\}\) ,假設 DSP 系統具有時間不變性,則

\(h[n] = T\{δ[n]\}\)


DSP 的輸入訊號為單位脈衝 Unit Impulse δ[n],輸出訊號為 h[n],因此 h[n] 經常稱為脈衝響應 Impulse Response。如果輸入 \( \{x(n)\}, n=1,2,...N\) 長度為 N,脈衝響應 \( \{h(n)\}, n=1,2,...M\) 長度為 M,則卷積運算結果,長度為 \(M+N-1\)

ex: 數位訊號 \(x=\{1,2,4,3,2,1,1\}\) ,脈衝響應為 \(h=\{1,2,3,1,1\}\),卷積運算

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

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

\(y[1] = \sum_{k=-∞}^{∞}x[k] \cdot h[1-k] = x[0] \cdot h[1] + x[1] \cdot h[0] = 1 \cdot 2 + 2 \cdot 1 = 4\)

\(y[2] = \sum_{k=-∞}^{∞}x[k] \cdot h[2-k] = x[0] \cdot h[2] + x[1] \cdot h[1] + x[2] \cdot h[0] = 1 \cdot 3 + 2 \cdot 2 + 4 \cdot 1 = 11\)

另一種計算方法:

因輸入訊號的數量 N=7,脈衝響應數量為 M = 5,卷積結果長度為 \(M+N-1 = 11\)

計算前先將 \(x[n]\) 兩邊補上 \(M-1=4\) 個 0,稱為 Zero-Padding,再將 h(n) 旋轉 180 度

x h
0 0 0 0 1 2 4 3 2 1 1 0 0 0 0 1 2 3 1 1
1 1 3 2 1
\(y[0] = 0 \cdot 1 + 0 \cdot 1 + 0 \cdot 3 + 0 \cdot 2 + 1 \cdot 1 = 1 \)
x h
0 0 0 0 1 2 4 3 2 1 1 0 0 0 0 1 2 3 1 1
1 1 3 2 1
\(y[1] = 0 \cdot 1 + 0 \cdot 1 + 0 \cdot 3 + 1 \cdot 2 + 2 \cdot 1 = 4 \)

最後卷積運算結果為

n 0 1 2 3 4 5 6 7 8 9 10
y[n] 1 4 11 18 23 20 16 10 6 2 1

實際運算時,會擷取卷積運算的部分結果

n 0 1 2 3 4 5 6
x[n] 1 2 4 3 2 1 1
y[n] 11 18 23 20 16 10 6

numpy 有提供 full, same 兩種卷積運算的結果

import numpy as np

x = np.array( [ 1, 2, 4, 3, 2, 1, 1 ] )
h = np.array( [ 1, 2, 3, 1, 1 ] )
y = np.convolve( x, h, 'full' )
y1 = np.convolve( x, h, 'same' )

print( "x =", x )
print( "h =", h )
print( "Full Convolution y =", y )
print( "Convolution y =", y1 )

執行結果

$ python convolution.py
x = [1 2 4 3 2 1 1]
h = [1 2 3 1 1]
Full Convolution y = [ 1  4 11 18 23 20 16 10  6  2  1]
Convolution y = [11 18 23 20 16 10  6]

交換率、結合率

import numpy as np

x = np.array( [ 1, 2, 4, 3, 2, 1, 1 ] )
h = np.array( [ 1, 1, 1 ] ) / 3

y = np.convolve( x, h, 'full' )
y1 = np.convolve( h, x, 'full' )
z = np.convolve( x, h, 'same' )
z1 = np.convolve( h, x, 'same' )

print( "x =", x )
print( "h =", h )

# 卷積運算滿足 交換率
print( "交換率")
print( "Full Convolution y =", y )
print( "Full Convolution y1 =", y1 )

print( "Convolution z =", z )
print( "Convolution z1 =", z1 )

print( "結合率")
h2 = np.array( [ 1, 2, 1 ] ) / 4
k1 = np.convolve( np.convolve( x, h, 'full' ), h2, 'full' )
k2 = np.convolve( x, np.convolve( h, h2, 'full' ), 'full' )
print( "Full Convolution k1 =", k1 )
print( "Full Convolution k2 =", k2 )

執行結果

$ python convolution.py
x = [1 2 4 3 2 1 1]
h = [0.33333333 0.33333333 0.33333333]
交換率
Full Convolution y = [0.33333333 1.         2.33333333 3.         3.         2.
 1.33333333 0.66666667 0.33333333]
Full Convolution y1 = [0.33333333 1.         2.33333333 3.         3.         2.
 1.33333333 0.66666667 0.33333333]
Convolution z = [1.         2.33333333 3.         3.         2.         1.33333333
 0.66666667]
Convolution z1 = [1.         2.33333333 3.         3.         2.         1.33333333
 0.66666667]
結合率
Full Convolution k1 = [0.08333333 0.41666667 1.16666667 2.16666667 2.83333333 2.75
 2.08333333 1.33333333 0.75       0.33333333 0.08333333]
Full Convolution k2 = [0.08333333 0.41666667 1.16666667 2.16666667 2.83333333 2.75
 2.08333333 1.33333333 0.75       0.33333333 0.08333333]

串聯、並聯

卷積運算的 DSP 系統符合線性時間不變性,是 LTI (Linear Time-Invariant) 系統。LTI 系統可以透過串聯或並聯,組合成另一個 LTI 系統。

串聯可表示為

\( x[n] → h_1[n] → h_2[n] → y[n] = x[n] → h_1[n] * h_2[n] → y[n]\)

\(y[n]=h_2[n]*(h_1[n]*x[n])=(h_2[n]*h_1[n])*x[n]=(h_1[n]*h_2[n])*x[n]\) (卷積運算符合結合率與交換率)

兩個 LTI 系統串聯後,修改先後順序,不影響結果。

並聯可表示為

\(x[n] → h_1[n]+h_2[n] → y[n] \)

\(y[n]=h_1[n]*x[n]+h_2[n]*x[n]=(h_1[n]+h_2[n])*x[n]\) (滿足分配率)

濾波器 filter

卷積最基本的應用是訊號的濾波 filtering,因此脈衝響應 Impulse Response 也經常稱為濾波器 filter。

平均濾波器 Average Filter

\(h[n]=\frac{1}{M} \{1,1,....1 \}, n=0,1,....,M-1\)

ex: \(x={1,2,4,3,2,1,1}, n=0,1,...6\) ,套用大小為 3 的平均濾波器,\(h[n]=\frac{1}{3} \{1,1,1\}, n=0,1,2\)

輸出訊號為 \(y= \{1,\frac{7}{3}, 3, 3, 2, \frac{4}{3}, \frac{2}{3} \}, n=0,1,...6\)

平均濾波器典型用來去除雜訊 (Noise Removal),也稱為 Smoothing。通常設計 Smoothing 的濾波器,其脈衝響應的總和為 1 ( \(\sum_{n}h[n]=1\))

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

# 產生 200 個點
t = np.linspace( 0, 1, 200, endpoint = False )
# 作 sinusoid() 弦波,並加上 均勻雜訊 Uniform Noise
x = 10 * np.cos( 2 * np.pi * 5 * t ) + random.uniform ( -5, 5, 200 )
#  average filter
h = np.ones( 7 ) / 7
# convolution
y = np.convolve( x, h, 'same' )

plt.figure( 1 )
plt.plot( t, x )
plt.xlabel( 't (second)' )
plt.ylabel( 'Amplitude' )

plt.figure( 2 )
plt.plot( t, y )
plt.xlabel( 't (second)' )
plt.ylabel( 'Amplitude' )

plt.show( )

左圖是加上雜訊的原始訊號,右圖是套用 average filter 的結果

高斯濾波器 Gaussian Filter

\(g[n]=e^{-n^2/2𝜎^2} \) 高斯函數的平均值 \(𝜇=0\),通常介於 -3𝜎 ~ 3𝜎 之間,因此設定高斯濾波器的大小為 6𝜎+1

如果高斯濾波器的 𝜎=1.0 ,則獲取高斯濾波器的數值,介於 -3 ~ 3 之間,濾波器大小為 7,且 \(\sum_{n}g[n]=1\),為了讓輸出訊號的數值範圍更接近輸入訊號,因此在取高斯函數作為濾波器的係數後,再進一步正規化。

高斯濾波器的典型應用也是去雜訊

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

# 產生 200 個點
t = np.linspace( 0, 1, 200, endpoint = False )
# 作 sinusoid() 弦波,並加上 均勻雜訊 Uniform Noise
x = 10 * np.cos( 2 * np.pi * 5 * t ) + random.uniform( -5, 5, 200 )

sigma = 3                                      # 標準差
filter_size = 6 * sigma + 1                    # 濾波器大小  19
gauss = signal.gaussian( filter_size, sigma )  # 濾波器係數

sum = np.sum( gauss )                          # 正規化
gauss = gauss / sum

y = np.convolve( x, gauss, 'same' )

plt.figure( 1 )
plt.plot( t, x )
plt.xlabel( 't (second)' )
plt.ylabel( 'Amplitude' )

plt.figure( 2 )
plt.plot( t, y )
plt.xlabel( 't (second)' )
plt.ylabel( 'Amplitude' )

plt.show( )

wav 濾波

一般來說,平均、高斯濾波的結果較平滑,造成輸出訊號的振幅會變小,使得音量變小。加上正規化 Normalization的處理,可調整輸出音訊檔的訊號範圍,介於 -30000~30000 之間。

import numpy as np
import wave
from scipy.io import wavfile
import struct
import scipy.signal as signal
import sys

def average_filtering( x, filter_size ):
    h = np.ones( filter_size ) / filter_size
    y = np.convolve( x, h, 'same' )
    return y

def gaussian_filtering( x, sigma):
    filter_size = 6 * sigma + 1
    gauss = signal.gaussian( filter_size, sigma )

    sum = np.sum( gauss )
    gauss = gauss / sum

    y = np.convolve( x, gauss, 'same' )
    return y

def normalization( x, maximum ):
    x_abs = abs( x )
    max_value = max( x_abs )
    y = x / max_value * maximum
    return y

def main( ):
    inputfile = "input.wav"
    output1 = "output1.wav"
    output2 = "output2.wav"

    # 讀取 wav
    wav = wave.open( inputfile, 'rb' )
    num_channels = wav.getnchannels( )  # 通道數
    sampwidth    = wav.getsampwidth( )      # 樣本寬度
    fs           = wav.getframerate( )          # 取樣頻率(Hz)
    num_frames   = wav.getnframes( )        # 音框數 = 樣本數
    comptype     = wav.getcomptype( )           # 壓縮型態
    compname     = wav.getcompname( )           # 無壓縮
    wav.close( )

    sampling_rate, x = wavfile.read( inputfile )    # 輸入訊號

    ##########################
    # DSP Average Filtering
    print( "Filtering" )
    print( "(1) Average Filtering" )

    # filter_size = eval( input( "Filter Size = " ) )
    filter_size = 31
    y = average_filtering( x, filter_size )

    # 正規化 調整音量
    y = normalization( x, 30000 )

    # 寫入 wav
    wav_file = wave.open( output1, '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( )

    ##########################
    # DSP Gaussian Filtering
    print( "(2) Gaussian Filtering" )

    # sigma = eval( input( "Sigma = " ) )
    sigma = 31
    y = gaussian_filtering( x, sigma )

    # 正規化 調整音量
    y = normalization( x, 30000 )

    # 寫入 wav
    wav_file = wave.open( output2, '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( )


if __name__ == "__main__":
    sys.exit(main())
    

References

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

2020/12/07

數位訊號 DSP

數位訊號處理系統 DSP,也就是將輸入的數位訊號,轉換成另一種數位訊號,\(y[n]= T(x[n])\)

DSP 系統分類

根據數位訊號處理性質可分為

  • 靜態 Static /動態 Dynamic System

    靜態系統 static system 就是輸出訊號只根據目前輸入訊號計算得來,也稱為 Memoryless System

    例如:Scalar Multiplication, \(y[n] = α*x[n]\),其中 α 為 scaling factor,可調整音量大小

    或是 \(y[n] = (x[n])^2\) 可計算功率 Power

​ 動態系統 dynamic system,就是輸出訊號根據過去、現在、未來的數位訊號計算得來,也稱為 Memory System

​ 例如: Moving Average, \(y[n] = \frac{1}{3}(x[n-1]+x[n]+x[n+1])\) 或是 \(y[n] = \frac{1}{3}(x[n-2]+x[n-1]+x[n])\)

  • 線性 Linear /非線性 Nonlinear System

    線性原則就是數位訊號先進行 Linear Combination,再經過 DSP 處理,其結果跟分別對數位訊號進行 DSP 處理,再進行 Linear Combination 的結果一樣

    線性系統: \(T(α*x_1[n]+𝛽*x_2[n]) = α*T(x_1[n]) + 𝛽*T(x_2[n])\)

    ​ ex: Moving Average \(y[n] = \frac{1}{3}(x[n-1]+x[n]+x[n+1])\)

    非線性系統:\(T(α*x_1[n]+𝛽*x_2[n]) ≠ α*T(x_1[n]) + 𝛽*T(x_2[n])\)

    ​ ex: 功率計算 \(y[n] = (x[n])^2\)

  • 時間不變性 Time-Invariant / 時變性 Time-Varying System

    時間不變性 Time-Invariant System 就是 DSP 運算方式不會隨著時間而改變。如果 \(y[n]=T(x[n])\) 則 \(y[n-n_0] = T(x[n-n_0])\)

    ​ ex: Scalar Multiplication

    時變性系統,DSP 運算方式隨著時間改變。因運算方式根據目前觀察到的數位訊號隨時更新,因此具有適應性

  • 因果 Causal / 非因果 Non-Causal System

    因果 Causal System,輸出訊號只根據目前與過去的數位訊號運算而得來

    ​ ex: Moving Average \(y[n] = \frac{1}{3}(x[n-2]+x[n-1]+x[n])\)

    非因果 Non-Causal System,輸出訊號會根據目前、過去及未來的數位訊號運算而得來

    ​ ex: Moving Average, \(y[n] = \frac{1}{3}(x[n-1]+x[n]+x[n+1])\)

    因果系統符合 Real-Time Application 的需求,因此 DSP 通常使用因果系統 Causal System 設計

  • 穩定 Stable / 不穩定 Non-Stable System

    穩定 Stable System,就是系統具有 Bounded Input - Bounded Output (BIBO) 的特性。當訊號落在限定範圍,則輸出訊號也會落在限定範圍內

    當 \( |x[n]| < B_x 則 |y[n]| < B_x \)

    ex: Scalar Multiplication, \(y[n] = α*x[n]\),如果 \(|x[n]|<1\) 則 \(|y[n|< α\)

    不穩定 Non-Stable System,就是系統不具有 Bounded Input - Bounded Output (BIBO) 的特性。

Linear Time-Invariant System (LTI) 線性時間不變性系統,是最具代表性的 DSP

DSP 基本運算

  • Scalar Multiplication

    \(y[n] = α \cdot x[n]\),其中 α 為縮放因子 scaling factor

  • 加/減法

    \( y[n]=x_1[n] ± x_2[n] \)

  • 乘法

    \( y[n]=x_1[n] \cdot x_2[n] \)

  • 時間延遲 time delay

    \( y[n]=x[n-n_0]\)

    當 \(n_0=1\) ,則 \(y[n]=x[n-1]\),又稱為單位延遲 unit-delay

  • Downsampling 降低取樣率,也稱為 Decimation

    \(y[n] = x[2n]\) 將 sample rate 降為一半

    最簡單的方式,是每兩個 sample,直接選用其中一個,但如果 N 比較大,可以用 N 個 sample 的平均值,作為新的取樣結果

  • Upsampling 增加取樣率

    通常是用內插法 interpolation 計算

    \(y[n] = x[n/2]\) 將 sample rate 兩倍

    最簡單的方法,稱為 Zero-Order Hold 內插法。也就是複製前一個 sample value。如果 N 比較大,可改用 Linear interpolation、Polynomimal Interpolation、Spline Interpolation

  • Resampling

    前面的 upsampling/downsampling 都是整數倍,如果是非整數倍,則稱為 resampling。

    Nyquist-Shannon 取樣定理:如果週期函數 x(t) 不包含高於 B cps(次/秒)的頻率,那麼,一系列小於 1/(2B) 秒的x(t)函數值將會受到前一個週期的x(t)函數值影響。因此 2B 樣本/秒或更高的取樣頻率將能使函數不受干擾。相對的,對於一個給定的取樣頻率 fs,完全重構的頻帶限制為 Bfs/2。

    因為在 downsampling 時,可能會發生無法滿足 Nyquist-Shannon 取樣定理的狀況,而發生取樣後訊號的頻率重疊(即高於取樣頻率一半的頻率成分將被重建成低於取樣頻率一半的訊號),也就是混疊 aliasing。

    解決方式是在 downsampling 前,先對原始訊號進行前處理,例如:low-pass filtering,降低原始數位訊號的最高頻率後,再進行 downsampling

    resampling 後的波形必須跟原始訊號的波形相似


import numpy as np
import scipy.signal as signal
import sys

def scalar_multiplication(x, alpha):
    y = alpha * x
    return y

def add(x1, x2):
    y = x1 + x2
    return y

def multiplication(x1, x2):
    y = x1 * x2
    return y

def time_delay(x, n0):
    y = x
    for i in range( n0 ):
        y = np.insert ( y, 0, 0 )   # 在0的位置插入 1 個 0
    return y

def downsampling( x, method = 1 ):
    N = int( len( x ) / 2 )
    y = np.zeros( N )

    if method == 1:             # Decimation
        for n in range( N ):
            y[n] = x[2*n]
    else:                       # Average
        for n in range( N ):
            y[n] = ( x[2*n] + x[2*n+1] ) / 2

    return y

def upsampling( x, method = 1 ):
    N = len( x ) * 2
    y = np.zeros( N )

    if method == 1:             # Zero-Order Hold
        for n in range( N ):
            y[n] = x[int( n / 2 )]
    else:                       # Linear Interpolation
        for n in range( N ):
            if int( n / 2 ) == n / 2:
                y[n] = x[int( n / 2 )]
            else:
                n1 = int( n / 2 )
                n2 = n1 + 1
                if n2 < len( x ):
                    y[n] = ( x[n1] + x[n2] ) / 2
                else:
                    y[n] = x[n1] / 2

    return y

# 利用 Scipy 的 signal.resample 處理
def resampling( x, sampling_rate ):
    num = int( len(x) * sampling_rate )
    y = signal.resample( x, num )
    return y

def main():
    np.set_printoptions(formatter={'all': lambda x: str(x)+","})

    x = np.array ( [ 1, 2, 4, 3, 2, 1, 1 ] )
    x1 = np.array ( [ 1, 2, 4, 3, 2, 1, 1 ] )
    x2 = np.array ( [ 0, 0, 1, 2, 4, 0, -1 ] )

    alpha = 2
    y = scalar_multiplication(x, alpha)
    print ( str(x) + " * " + str(alpha) + " -> " + str(y) )

    y = add(x1, x2)
    print ( str(x1) + " + " + str(x2) + " -> " + str(y) )

    y = multiplication(x1, x2)
    print ( str(x1) + " * " + str(x2) + " -> " + str(y) )

    n0 = 2
    y = time_delay(x, n0)
    print ( str(x) + " time_delay " + str(n0) + " -> " + str(y) )

    # downsampling
    y1 = downsampling( x, 1 )
    y2 = downsampling( x, 2 )
    # upsampling
    z1 = upsampling( x, 1 )
    z2 = upsampling( x, 2 )

    print ("")
    print ("original                        : " + str(x) )
    print ("downsampling Decimation         : " + str(y1))
    print ("downsampling Average            : " + str(y2))
    print ("upsampling Zero-Order Hold      : " + str(z1))
    print ("upsampling Linear Interpolation : " + str(z2))

    y = resampling( x, 1.5 )
    print ("")
    print ( str(x) + " resampling 1.5 -> " + str(y) )

if __name__ == "__main__":
    sys.exit(main())

執行結果

[1, 2, 4, 3, 2, 1, 1,] * 2 -> [2, 4, 8, 6, 4, 2, 2,]
[1, 2, 4, 3, 2, 1, 1,] + [0, 0, 1, 2, 4, 0, -1,] -> [1, 2, 5, 5, 6, 1, 0,]
[1, 2, 4, 3, 2, 1, 1,] * [0, 0, 1, 2, 4, 0, -1,] -> [0, 0, 4, 6, 8, 0, -1,]
[1, 2, 4, 3, 2, 1, 1,] time_delay 2 -> [0, 0, 1, 2, 4, 3, 2, 1, 1,]

original                        : [1, 2, 4, 3, 2, 1, 1,]
downsampling Decimation         : [1.0, 4.0, 2.0,]
downsampling Average            : [1.5, 3.5, 1.5,]
upsampling Zero-Order Hold      : [1.0, 1.0, 2.0, 2.0, 4.0, 4.0, 3.0, 3.0, 2.0, 2.0, 1.0, 1.0, 1.0, 1.0,]
upsampling Linear Interpolation : [1.0, 1.5, 2.0, 3.0, 4.0, 3.5, 3.0, 2.5, 2.0, 1.5, 1.0, 1.0, 1.0, 0.5,]

[1, 2, 4, 3, 2, 1, 1,] resampling 1.5 -> [1.0000000000000007, 1.384837322270305, 3.027762174254682,
 4.024749771849889, 3.3105260815101465, 2.3971667816643563,
 1.8264208230878693, 1.0876017695626832, 0.8352909211473022,
 1.1056443546527661,]

wav DSP

讀取 wav file 的資料,進行 downsampling, upsampling, resampling 後,刻意以原本的 sample rate 數值,寫入新的 wav file,聽起來就像是將原始 wav 進行加速或別的特效處理。


import numpy as np
import wave
from scipy.io import wavfile
import scipy.signal as signal
import struct
import sys

class WavFileParams:
    def __init__(self, filename, num_channels, sampwidth, fs, num_frames, comptype, compname):
        self.filename = filename

        self.num_channels = num_channels                # 通道數
        self.sampwidth = sampwidth                      # 樣本寬度
        self.fs = fs                                    # 取樣頻率(Hz)
        self.num_frames = num_frames                    # 音框數 = 樣本數
        self.comptype = comptype                        # 壓縮型態
        self.compname = compname                        # 壓縮名稱

def write_wav_file(wav_params, digital_signel_points):
    wav_file = wave.open( wav_params.filename, 'w' )
    wav_file.setparams(( wav_params.num_channels, wav_params.sampwidth, wav_params.fs, wav_params.num_frames, wav_params.comptype, wav_params.compname ))

    for s in digital_signel_points :
        # 以 struct.pack 將 int 資料轉換為 16 bits
        wav_file.writeframes( struct.pack( 'h', int ( s ) ) )

    wav_file.close( )

def downsampling( x, method = 1 ):
    N = int( len( x ) / 2 )
    y = np.zeros( N )

    if method == 1:             # Decimation
        for n in range( N ):
            y[n] = x[2*n]
    else:                       # Average
        for n in range( N ):
            y[n] = ( x[2*n] + x[2*n+1] ) / 2

    return y

def upsampling( x, method = 1 ):
    N = len( x ) * 2
    y = np.zeros( N )

    if method == 1:             # Zero-Order Hold
        for n in range( N ):
            y[n] = x[int( n / 2 )]
    else:                       # Linear Interpolation
        for n in range( N ):
            if int( n / 2 ) == n / 2:
                y[n] = x[int( n / 2 )]
            else:
                n1 = int( n / 2 )
                n2 = n1 + 1
                if n2 < len( x ):
                    y[n] = ( x[n1] + x[n2] ) / 2
                else:
                    y[n] = x[n1] / 2

    return y

# 利用 Scipy 的 signal.resample 處理
def resampling( x, sampling_rate ):
    num = int( len(x) * sampling_rate )
    y = signal.resample( x, num )
    return y

def main():
    input_wav_file = "input.wav"
    output_wav_file_prefix = "output"

    wav = wave.open( input_wav_file, 'rb' )
    num_channels = wav.getnchannels( )  # 通道數
    sampwidth    = wav.getsampwidth( )  # 樣本寬度
    fs           = wav.getframerate( )  # 取樣頻率(Hz)
    num_frames   = wav.getnframes( )    # 音框數 = 樣本數
    comptype     = wav.getcomptype( )   # 壓縮型態
    compname     = wav.getcompname( )   # 無壓縮
    wav.close( )

    sampling_rate, x = wavfile.read( input_wav_file )   # 輸入訊號

    # downsampling Decimation
    y = downsampling( x, 1 )
    wav_params = WavFileParams(
        filename = output_wav_file_prefix+"_downsampling_Decimation.wav",       # 檔案名稱

        num_channels = num_channels,    # 通道數
        sampwidth = sampwidth,          # 樣本寬度
        fs = fs,                        # 取樣頻率(Hz)
        num_frames = num_frames,        # 音框數 = 樣本數
        comptype = comptype,            # 壓縮型態
        compname = compname             # 無壓縮
    )
    write_wav_file(wav_params, y);

    # downsampling Average
    y = downsampling( x, 2 )
    wav_params = WavFileParams(
        filename = output_wav_file_prefix+"_downsampling_Average.wav",      # 檔案名稱

        num_channels = num_channels,    # 通道數
        sampwidth = sampwidth,          # 樣本寬度
        fs = fs,                        # 取樣頻率(Hz)
        num_frames = num_frames,        # 音框數 = 樣本數
        comptype = comptype,            # 壓縮型態
        compname = compname             # 無壓縮
    )
    write_wav_file(wav_params, y);

    # upsampling Zero-Order Hold
    y = upsampling( x, 1 )
    wav_params = WavFileParams(
        filename = output_wav_file_prefix+"_upsampling_Zero-OrderHold.wav",     # 檔案名稱

        num_channels = num_channels,    # 通道數
        sampwidth = sampwidth,          # 樣本寬度
        fs = fs,                        # 取樣頻率(Hz)
        num_frames = num_frames,        # 音框數 = 樣本數
        comptype = comptype,            # 壓縮型態
        compname = compname             # 無壓縮
    )
    write_wav_file(wav_params, y);

    # upsampling Linear Interpolation
    y = upsampling( x, 2 )
    wav_params = WavFileParams(
        filename = output_wav_file_prefix+"_upsampling_LinearInterpolation.wav",        # 檔案名稱

        num_channels = num_channels,    # 通道數
        sampwidth = sampwidth,          # 樣本寬度
        fs = fs,                        # 取樣頻率(Hz)
        num_frames = num_frames,        # 音框數 = 樣本數
        comptype = comptype,            # 壓縮型態
        compname = compname             # 無壓縮
    )
    write_wav_file(wav_params, y);

    # resampling
    sampling_rate = 1.5
    y = resampling( x, sampling_rate )
    wav_params = WavFileParams(
        filename = output_wav_file_prefix+"_resampling.wav",        # 檔案名稱

        num_channels = num_channels,    # 通道數
        sampwidth = sampwidth,          # 樣本寬度
        fs = fs,                        # 取樣頻率(Hz)
        num_frames = num_frames,        # 音框數 = 樣本數
        comptype = comptype,            # 壓縮型態
        compname = compname             # 無壓縮
    )
    write_wav_file(wav_params, y);

if __name__ == "__main__":
    sys.exit(main())

振幅調變 AM

振幅調變 Amplitude Modulation (AM),或簡稱調幅,定義為

\(y(t)=x(t)*cos(2*pi*f_c*t)\) 其中 \(f_c\) 為載波頻率 carrier frequency,x(t) 為輸入訊號, y(t) 為輸出訊號

import numpy as np
import wave
from scipy.io import wavfile
import struct
import sys

def AM( x, f, fs ):
    t = np.zeros( len( x ) )
    for i in range( len( x ) ):
        t[i] = i / fs
    carrier = np.cos( 2 * np.pi * f * t )
    return x * carrier

def main( ):
    # infile  = input( "Input File: " )
    # outfile = input( "Output File: " )
    infile = "input.wav"
    outfile = "output.wav"

    wav = wave.open( infile, 'rb' )
    num_channels = wav.getnchannels( )  # 通道數
    sampwidth    = wav.getsampwidth( )  # 樣本寬度
    fs           = wav.getframerate( )  # 取樣頻率(Hz)
    num_frames   = wav.getnframes( )    # 音框數 = 樣本數
    comptype     = wav.getcomptype( )   # 壓縮型態
    compname     = wav.getcompname( )   # 無壓縮
    wav.close( )

    sampling_rate, x = wavfile.read( infile )   # 輸入訊號

    # AM
    # fc = eval( input( "Enter carrier frequency (Hz): " ) )
    fc = 1000
    y = AM( x, fc, fs )

    # output
    wav_file = wave.open( outfile, '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( )

if __name__ == "__main__":
    sys.exit(main())

References

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