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

2020/11/30

數位訊號雜訊 in python

雜訊生成

\(y[n]=x[n]+η[n]\) 其中 η[n] 稱為雜訊 Noise

DSP 中,如果資料分布情形或統計參數(例如標準差、平均值),不會隨著時間改變,則該訊號為平穩訊號 Stationary Signals,反之稱為非平穩訊號。例如馬達產生的聲音是平穩訊號,一般語音是非平穩訊號

均勻雜訊 Uniform Noise

以 uniform distribution function 產生隨機亂數資料,因為是均勻分佈,平均值 𝜇 為 0,雜訊連續樣本之間互相獨立,沒有相關。在 DSP 這種雜訊稱為 White Noise

均勻分佈函數:

給定隨機變數 random variable z,Uniform Distribution Function 定義為

\(p(z) = \left\{ \begin{array}{ll} 1/2 & \mbox{if -1≤z≤1} \\ 0 & \mbox{otherwise} \end{array} \right.\) ,其中 z 介於 -1~1 之間

均勻分佈的機率密度函數,必須滿足總機率值 (積分) 為 1。曲線下的面積為 1

\(\int_{-∞}^{∞} p(z) {\rm d}z = 1\)

高斯雜訊 Gaussian Noise

高斯分佈函數定義

\(p(z) = \frac{1}{\sqrt{2 \pi σ^2}} e^{(z-μ)^2/(2σ^2)}\)

其中 μ 為平均,σ 是標準差

噪音的隨機亂數資料以高斯分佈函數也就是常態分佈函數 Normal Distribution Function 產生,其平均值 𝜇 為 0,標準差 𝜎 為 1,且滿足 \(\int_{-∞}^{∞}p(z)dz = 1\),就是函數曲線下的面積總和為 1(所有機率的總和為 1)

雜訊連續樣本之間互相獨立,沒有自相關 autocorrelation。在 DSP 這種雜訊稱為 Gaussian White Noise

布朗尼雜訊 Brownian Noise

根據布朗尼運動 Brownian Motion 產生的亂數資料,也稱為 Random Walk Noise

Brownian Noise 也稱為 Brown Noise,因為 Brownian Motion 的發現者 Robert Brown 而這樣命名。

Brownian Noise 的每個 sample,取自先前所有 samples 的總和(積分)。

程式中產生 Brownian Noise 後,會再針對數值做正規化 Normalization,讓結果介於 -1~1 之間

Brownian Noise 的變動幅度比均勻雜訊緩慢,連續的 samples 之間有相關性,結果會造成數位訊號的振幅微幅改變的狀況,也就是發生訊號飄移的狀況

脈衝雜訊 Impulse Noise

瞬間的脈衝產生的雜訊,當訊號受到電磁干擾、或是唱片(CD)有刮痕,會產生 Impulse Noise。可用振幅 Amplitude 與發生機率 Probability 來模擬產生 Impulse Noise

以 python 程式碼產生 noise

需要的套件

pip install numpy
pip install scipy

程式裡面有以下這些雜訊的產生方式

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

class WavFileParams:
    def __init__(self, filename, amplitude, frequency, duration, fs, num_channels, sampwidth, comptype, compname):
        self.filename = filename
        self.amplitude = amplitude                      # 振幅
        self.frequency = frequency                      # 頻率(Hz)
        self.duration = duration                        # 時間長度(秒)
        self.fs = fs                                    # 取樣頻率(Hz)
        self.num_samples = self.duration * self.fs      # 樣本數

        self.num_channels = num_channels                # 通道數
        self.sampwidth = sampwidth                      # 樣本寬度
        self.num_frames = self.num_samples              # 音框數 = 樣本數
        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 uniform_noise():
    wav_params = WavFileParams(
        filename = "uniform_noise.wav",     # 檔案名稱

        amplitude = 30000,          # 振幅
        frequency = 100,            # 頻率(Hz)
        duration = 1,               # 時間長度(秒)
        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, wav_params.duration, wav_params.num_samples, endpoint = False );
    # 原始訊號
    x = wav_params.amplitude * np.cos( 2 * np.pi * wav_params.frequency * t );
    # uniform noise, 在 -1000 ~ 1000 之間
    noise = np.random.uniform( -1000, 1000, wav_params.num_samples )
    y = x + noise

    write_wav_file(wav_params, y);

def gaussian_noise():
    wav_params = WavFileParams(
        filename = "gaussian_noise.wav",        # 檔案名稱

        amplitude = 30000,          # 振幅
        frequency = 100,            # 頻率(Hz)
        duration = 1,               # 時間長度(秒)
        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, wav_params.duration, wav_params.num_samples, endpoint = False );
    # 原始訊號
    x = wav_params.amplitude * np.cos( 2 * np.pi * wav_params.frequency * t );
    # 雜訊(高斯分佈), 在 0 ~ 1000 之間
    noise = np.random.normal( 0, 1000, wav_params.num_samples )
    y = x + noise

    # 為避免相加後超過 16 bits 資料長度的值,選擇使用比較小的振幅,且在相加後,用 np.clip 抑制數值範圍
    y_clip = np.clip(y, -32768, 32767)

    write_wav_file(wav_params, y_clip);

def brownian_noise():
    wav_params = WavFileParams(
        filename = "brownian_noise.wav",        # 檔案名稱

        amplitude = 30000,          # 振幅
        frequency = 100,            # 頻率(Hz)
        duration = 1,               # 時間長度(秒)
        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, wav_params.duration, wav_params.num_samples, endpoint = False );
    # 原始訊號
    x = wav_params.amplitude * np.cos( 2 * np.pi * wav_params.frequency * t );

    # 布朗尼雜訊 brownian_noise, 在 -1 ~ 1 之間
    n1 = np.random.uniform( -1, 1, wav_params.num_samples )
    # 累加
    ns = np.cumsum( n1 )
    # 正規化
    mean = np.mean( ns )
    max = np.max( np.absolute( ns - mean ) )
    noise = (( ns - mean ) / max)
    # 把最後正規化的數值,放大 10000 倍,強化 noise
    noise = noise * 10000.0

    y = x + noise

    # 為避免相加後超過 16 bits 資料長度的值,選擇使用比較小的振幅,且在相加後,用 np.clip 抑制數值範圍
    y_clip = np.clip(y, -32768, 32767)

    write_wav_file(wav_params, y_clip);

def impulse_noise():
    wav_params = WavFileParams(
        filename = "impulse_noise.wav",     # 檔案名稱

        amplitude = 30000,          # 振幅
        frequency = 100,            # 頻率(Hz)
        duration = 1,               # 時間長度(秒)
        fs = 44100,                 # 取樣頻率(Hz)
        # num_samples = duration * fs,  # 樣本數

        num_channels = 1,           # 通道數
        sampwidth = 2,              # 樣本寬度
        # num_frames = num_samples, # 音框數 = 樣本數
        comptype = "NONE",          # 壓縮型態
        compname = "not compressed" # 無壓縮
    )

    # impulse noise 發生機率 (%)
    probability = 5

    t = np.linspace( 0, wav_params.duration, wav_params.num_samples, endpoint = False );
    # 原始訊號
    x = wav_params.amplitude * np.cos( 2 * np.pi * wav_params.frequency * t );

    noise = np.zeros( x.size )                          # 脈衝雜訊
    for i in range( x.size ):
        p1 = np.random.uniform( 0, 1 )
        if p1 < probability / 100:
            p2 = np.random.uniform( 0, 1 )
            if p2 < 0.5:
                noise[i] = wav_params.amplitude
            else:
                noise[i] = -wav_params.amplitude

    y = x + noise

    # 為避免相加後超過 16 bits 資料長度的值,選擇使用比較小的振幅,且在相加後,用 np.clip 抑制數值範圍
    y_clip = np.clip(y, -32768, 32767)

    write_wav_file(wav_params, y_clip);

def main():
    uniform_noise()
    gaussian_noise()
    brownian_noise()
    impulse_noise()

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

以下分別是產生出來的 wav 的波形,第三個 Brownian Noise 可看出有上下飄移的狀況

訊號雜訊比

Signal-to-Noise Ratio 也就是 SNR 定義為

\( SNR = \frac{P_{signal}}{P_{noise}} \) ,\(P_{signal}\) 是訊號功率,\(P_{noise}\) 是雜訊功率

SNR 也可以定義為

\( SNR = (\frac{A_{signal}}{A_{noise}})^2 \) ,\(A_{signal}\) 是訊號振幅,\(A_{noise}\) 是雜訊振幅

SNR 通常用分貝 Decibels dB 為單位,故定義為

\( SNR = 10log_{10}(\frac{P_{signal}}{P_{noise}}) dB\) 或

\( SNR = 10log_{10}(\frac{A_{signal}}{A_{noise}})^2 dB = 20log_{10}(\frac{A_{signal}}{A_{noise}}) dB \)

如果 SNR 值越高,表示訊號的品質越好, SNR 越低則訊號受到雜訊干擾越嚴重

References

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

2020/11/23

數位訊號生成 in python

週期性訊號

\(x(t)=x(t+T)\) ,其中 T 為週期,是固定值

弦波 Sinusoids

\(x[n] = x(nT_s) = A cos(2 \pi f nT_s)= A cos(2 \pi fn/f_s)\)

ex: \(x(t)=cos(2 \pi \cdot 5 t)\) 其中 A=1, f=5Hz, 時間 0~1s

import numpy as np
import matplotlib.pyplot as plt

t = np.linspace( 0, 1, 1000, endpoint = False ) # 定義時間陣列
x = np.cos( 2 * np.pi * 5 * t )                 # 產生弦波

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

plt.show( )

方波 Square

\(x(t) = A \cdot sgn(sin(ωt))\) 或 \(x(t)=A \cdot sgn(sin(2 \pi ft))\)

其中 A 為振幅,ω 是角頻率 Angular Frequency,f 為頻率 frequency

sgn 是符號函數 \(sgn(x) = \left\{ \begin{array}{ll} 1 & \mbox{if x>0} \\ 0 & \mbox{if x = 0} \\ -1 & \mbox{if x<0} \end{array} \right.\)

可用 SciPy 的 signal.square 產生方波

ex: A =1, f=5Hz, t=0~1 的方波

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

t = np.linspace( 0, 1, 1000 )           # 定義時間陣列
x = signal.square( 2 * np.pi * 5 * t )  # 產生方波

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

plt.show( )

鋸齒波 Sawtooh Wave

A = 1, f=5Hz, t=0~1s 的鋸齒波

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

t = np.linspace( 0, 1, 1000 )            # 定義時間陣列
x = signal.sawtooth( 2 * np.pi * 5 * t ) # 產生鋸齒波

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

plt.show( )

三角波 Triangle Wave

A = 1, f=5Hz, t=0~1s

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

t = np.linspace( 0, 1, 1000 )                 # 定義時間陣列
x = signal.sawtooth( 2 * np.pi * 5 * t, 0.5 ) # 產生三角波

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

plt.show( )

諧波 Harmonic

兩個以上的弦波組成

\(x(t) = \sum_{k=1}^{N}A_k cos(2 \pi f_k t)\)

\(A_k\) 是第 k 個弦波的振幅,\(f_1\) 稱為基礎頻率 Fundamental Frequency (基頻), \(f_k\) 是第 k 個弦波的頻率,是 \(f_1\)的整數倍。第一個弦波稱為基礎弦波 Fundamental Sinusoid

ex:

\(x_1(t) = cos(2 \pi 2 t)\)

\(x_2(t) = cos(2 \pi 4 t)\)

諧波 \(x(t) = x_1(t)+x_2(t)\)

import numpy as np
import matplotlib.pyplot as plt

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

f1 = 2                                          # 定義基礎頻率
x1 = np.cos( 2 * np.pi * f1 * t )               # 產生第1個弦波
x2 = np.cos( 2 * np.pi * 2 * f1 * t )           # 產生第2個弦波
x = x1 + x2                                     # 產生諧波

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

plt.subplot(132)
plt.plot( t, x2 )
plt.xlabel( 't (second)' )
plt.ylabel( 'Amplitude' )
plt.axis( [ 0, 1, -2, 2 ] )

plt.subplot(133)
plt.plot( t, x )
plt.xlabel( 't (second)' )
plt.ylabel( 'Amplitude' )
plt.axis( [ 0, 1, -2, 2 ] )

plt.show( )

節拍波 beat

兩個不同頻率的弦波,相乘得來

\(x(t) = A cos(2 \pi f_1 t) \cdot cos(2 \pi f_2 t)\) ,A 為振幅, \(f_1\) 是低頻訊號的頻率,\(f_2\)是高頻訊號的頻率

低頻的訊號,構成訊號的 Envelop,高頻訊號,稱為載波訊號 Carrier Signal

ex: \(x(t) = A cos(2 \pi 20 t) \cdot cos(2 \pi 200 t)\)

import numpy as np
import matplotlib.pyplot as plt

t = np.linspace( 0, 0.1, 1000, endpoint = False )   # 定義時間陣列    

f1 = 20                                             # 低頻頻率
f2 = 200                                            # 高頻頻率
x = np.cos( 2 * np.pi * f1 * t ) * np.cos( 2 * np.pi * f2 * t )
envelop1 =  np.cos( 2 * np.pi * f1 * t )            # 包絡
envelop2 = -np.cos( 2 * np.pi * f1 * t )

plt.plot( t, x, '-' )                               # 繪圖
plt.plot( t, envelop1, '--', color = 'b' )
plt.plot( t, envelop2, '--', color = 'b' )
plt.xlabel( 't (second)' )
plt.ylabel( 'Amplitude' )
plt.axis( [ 0, 0.1, -1, 1 ] )

plt.show( )

振幅調變

調幅 Amplitude Modulation (AM) 可定義為

\(y(t)=x(t) \cdot cos(2 \pi f_ct)\) 其中 x(t) 是輸入訊號, \(f_c\) 稱為載波頻率 Carrier Frequency

跟 beat 乘法運算類似。常應用於 AM 收音機、無線電對講機。通常輸入訊號 x(t) 頻率較低,屬於人類的聽力範圍 20Hz ~ 20kHz,載波頻率比較高,例如低頻無線電波,落在 30k ~ 300k Hz

非週期性訊號

淡入與淡出

隨著時間,讓 amplitude 逐漸增加/減少的訊號

ex: \(x(t) = cos(2 \pi 5 t)\) 套用淡入、淡出效果

import numpy as np
import matplotlib.pyplot as plt

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

x = np.cos( 2 * np.pi * 5 * t )                 # 產生弦波
a = np.linspace( 1, 0, 1000, endpoint = False ) # 產生淡出陣列 
x = x * a                                       # 套用淡出效果 

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

plt.show( )

剛剛是用線性方式產生淡出陣列,也可以用指數衰減方式產生

\(x(t) = A e^{-αt} cos(2 \pi f t)\)

A 為振幅,f 是頻率,α 為衰減參數,可調整作用時間

啁啾訊號 Chirp

隨著時間,頻率逐漸增加/減少的訊號,分別稱為 Up-Chirp/Down-Chirp

Chirp 經常稱為 Sweep Singal

常應用於 radar 或 sonar,在不同時間點掃描不同的頻率區間/頻帶,藉以偵測目標物

ex: 產生 0~5 秒,由 0Hz 線性增加至 5Hz

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

t = np.linspace( 0, 5, 1000, endpoint = False ) # 定義時間陣列
x = signal.chirp( t, 0, 5, 5, 'linear' )        # 產生啁啾訊號

plt.plot( t, x )                                # 繪圖
plt.xlabel( 't (second)' )
plt.ylabel( 'Amplitude' )
plt.axis( [ 0, 5, -1.5, 1.5 ] )

plt.show( )

wave

需要的套件

pip install numpy
pip install scipy

程式裡面有以下這些波形 wav 的產生方式

  • sinusoid() 弦波
  • square() 方波
  • sawtooth() 鋸齒波
  • triangle() 三角波
  • harmonic() 諧波
  • beat() 節拍波
  • fadeout() 淡出
  • chirp() 啁啾
import numpy as np
import scipy.signal as signal
import wave
import struct
import sys

class WavFileParams:
    def __init__(self, filename, amplitude, frequency, duration, fs, num_channels, sampwidth, comptype, compname):
        self.filename = filename
        self.amplitude = amplitude                      # 振幅
        self.frequency = frequency                      # 頻率(Hz)
        self.duration = duration                        # 時間長度(秒)
        self.fs = fs                                    # 取樣頻率(Hz)
        self.num_samples = self.duration * self.fs      # 樣本數

        self.num_channels = num_channels                # 通道數
        self.sampwidth = sampwidth                      # 樣本寬度
        self.num_frames = self.num_samples              # 音框數 = 樣本數
        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( )

# sinusoid 弦波
def sinusoid():
    wav_params = WavFileParams(
        filename = "sinusoid.wav",      # 檔案名稱

        amplitude = 30000,          # 振幅
        frequency = 100,            # 頻率(Hz)
        duration = 3,               # 時間長度(秒)
        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, wav_params.duration, wav_params.num_samples, endpoint = False );
    x = wav_params.amplitude * np.cos( 2 * np.pi * wav_params.frequency * t );

    write_wav_file(wav_params, x);

def square():
    # square 方波
    wav_params = WavFileParams(
        filename = "square.wav",    # 檔案名稱

        amplitude = 30000,          # 振幅
        frequency = 100,            # 頻率(Hz)
        duration = 3,               # 時間長度(秒)
        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, wav_params.duration, wav_params.num_samples, endpoint = False );
    x = wav_params.amplitude * signal.square( 2 * np.pi * wav_params.frequency * t );

    write_wav_file(wav_params, x);

def sawtooth():
    # sawtooth 鋸齒波
    wav_params = WavFileParams(
        filename = "sawtooth.wav",  # 檔案名稱

        amplitude = 30000,          # 振幅
        frequency = 100,            # 頻率(Hz)
        duration = 3,               # 時間長度(秒)
        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, wav_params.duration, wav_params.num_samples, endpoint = False );
    x = wav_params.amplitude * signal.sawtooth( 2 * np.pi * wav_params.frequency * t )

    write_wav_file(wav_params, x);

def triangle():
    # triangle 三角波
    wav_params = WavFileParams(
        filename = "triangle.wav",  # 檔案名稱

        amplitude = 30000,          # 振幅
        frequency = 100,            # 頻率(Hz)
        duration = 3,               # 時間長度(秒)
        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, wav_params.duration, wav_params.num_samples, endpoint = False );
    x = wav_params.amplitude * signal.sawtooth( 2 * np.pi * wav_params.frequency * t, 0.5 )

    write_wav_file(wav_params, x);

def harmonic():
    # harmonic 諧波
    wav_params = WavFileParams(
        filename = "harmonic.wav",  # 檔案名稱

        amplitude = 10000,          # 振幅
        frequency = 100,            # 頻率(Hz)
        duration = 3,               # 時間長度(秒)
        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, wav_params.duration, wav_params.num_samples, endpoint = False );
    x1 = wav_params.amplitude * np.cos( 2 * np.pi * wav_params.frequency * t );
    x2 = wav_params.amplitude * np.cos( 2 * np.pi * (2 * wav_params.frequency) * t );

    x = x1 + x2

    # 為避免相加後超過 16 bits 資料長度的值,選擇使用比較小的振幅,且在相加後,用 np.clip 抑制數值範圍
    x_clip = np.clip(x, -32768, 32767)

    write_wav_file(wav_params, x_clip);

def beat():
    # beat 節拍波
    wav_params = WavFileParams(
        filename = "beat.wav",      # 檔案名稱

        amplitude = 30000,          # 振幅
        frequency = 0,              # 頻率(Hz)
        duration = 10,              # 時間長度(秒)
        fs = 44100,                 # 取樣頻率(Hz)
        # num_samples = duration * fs,  # 樣本數

        num_channels = 1,           # 通道數
        sampwidth = 2,              # 樣本寬度
        # num_frames = num_samples, # 音框數 = 樣本數
        comptype = "NONE",          # 壓縮型態
        compname = "not compressed" # 無壓縮
    )
    f1 = 20                     # 低頻頻率(Hz)
    f2 = 200                    # 高頻頻率(Hz)

    t = np.linspace( 0, wav_params.duration, wav_params.num_samples, endpoint = False )
    x = wav_params.amplitude * np.cos( 2 * np.pi * f1 * t ) * np.cos( 2 * np.pi * f2 * t )

    write_wav_file(wav_params, x);

def fadeout():
    # fadeout 淡出
    wav_params = WavFileParams(
        filename = "fadeout.wav",   # 檔案名稱

        amplitude = 30000,          # 振幅
        frequency = 300,            # 頻率(Hz)
        duration = 3,               # 時間長度(秒)
        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, wav_params.duration, wav_params.num_samples, endpoint = False )
    # 產生淡出陣列, 1 -> 0, num_samples 個值
    a = np.linspace( 1, 0, wav_params.num_samples, endpoint = False )
    x = wav_params.amplitude * a * np.cos( 2 * np.pi * wav_params.frequency * t )

    write_wav_file(wav_params, x);

def chirp():
    # chirp 啁啾, 頻率隨著時間, 逐漸增加或減少
    wav_params = WavFileParams(
        filename = "chirp.wav",     # 檔案名稱

        amplitude = 30000,          # 振幅
        frequency = 0,              # 頻率(Hz)
        duration = 10,              # 時間長度(秒)
        fs = 44100,                 # 取樣頻率(Hz)
        # num_samples = duration * fs,  # 樣本數

        num_channels = 1,           # 通道數
        sampwidth = 2,              # 樣本寬度
        # num_frames = num_samples, # 音框數 = 樣本數
        comptype = "NONE",          # 壓縮型態
        compname = "not compressed" # 無壓縮
    )

    f0 = 0                      # 初始頻率(Hz)
    f1 = 1000                   # 終止頻率(Hz)

    t = np.linspace( 0, wav_params.duration, wav_params.num_samples, endpoint = False )
    # 以線性方式遞增頻率
    x = wav_params.amplitude * signal.chirp( t, f0, wav_params.duration, f1, 'linear' )

    write_wav_file(wav_params, x);

def main():
    sinusoid()
    square()
    sawtooth()
    triangle()
    harmonic()
    beat()
    fadeout()
    chirp()

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

References

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

2020/11/16

數位訊號

類比訊號透過取樣 sampling 與量化 quantization 兩個步驟,擷取數位訊號,並以數學表示法討論數位訊號。

數位訊號定義為:隨著時間改變的離散訊號,也稱為離散時間訊號 Discrete-Time Signals

自然界的訊號通常都是類比訊號,因此必須先透過類比/數位轉換器 Analog/Digital Converter, A/D Converter 將類比訊號轉成數位訊號,以便電腦系統處理。過程是根據類比訊號,在時間軸上擷取離散的樣本 samples,藉以產生數位訊號。x(t) 代表類比訊號,x[n] 代表數位訊號,t 為時間,n 是索引。

取樣與量化

取樣 sampling 可定義為 \(x[n]=x(nT_s)\),\(T_s\) 稱為取樣週期 sampling period 或取樣間隔 sampling interval,單位為 second,另外 \(f_s=\frac{1}{T_s}\) 稱為取樣頻率 sampling frequency,或取樣率 sampling rate

取樣是訊號在時間軸上的數位化過程,每隔 \(T_s\) 時間,擷取一個樣本 sample,產生一個數字序列,因每個樣本來自短暫的時間框,也稱為音框 frame

取樣率 sampling frequency 解釋為:每秒取樣的次數。由於取樣頻率固定不變,又稱為均勻取樣 Uniform Sampling

ex: 弦波 \(x(t)=A cos(ωt+ϕ) = Acos(2 \pi ft+ϕ)\) ,若取樣週期 \(T_s\),取樣頻率 \(f_s\),則弦波的數位訊號可表示為

\(x[n]=x(nT_s)=A cos(2 \pi fnT_s+ϕ)=A cos(2 \pi fn/f_s+ϕ)\),n為整數

通常可假設 \(\hat{ω} = 2 \pi f/f_s\) 稱為正規化角頻率 Normalized Angular Frequency,單位是弧度 Radians,因此數位訊號也可表示為

\(x[n] = A cos(\hat{ω}n+ϕ)\)

Nyquist-Shannon 取樣定理

假設原始訊號為頻帶限制訊號 Band-Limited Signal,最高頻率為 \(f_H\),若取樣頻率為 \(f_s\) ,則 \(f_s > 2 f_H\) 時,才能保證可重建原始訊號

若取樣頻率超過原始訊號的最高頻率的兩倍,可充分表示原始訊號,如取樣頻率不足,會產生混疊 aliasing 現象。

ex: 人類聽力範圍 20Hz ~ 20k Hz,最高為 20kHz,如果 \(f_s>40kHz\) 可充分表示人類可感知的原始訊號。目前 mp3 通常設定為 44kHz or 48kHz,符合基本條件,不會發生 aliasing

根據取樣定理,取樣頻率決定了可取樣的原始最高頻率,又稱為 Nyquist 頻率。例如取樣頻率 44kHz,則 Nyquist 頻率為 22kHz

量化 Quantization

量化是將訊號的振幅 Amplitude 經過數位化轉換為數值,常用的數值範圍:

  • 8 bits: 0~255 or -128~127
  • 16 bits: -32768 ~ 32767

每個樣本可使用的位元數稱為位元解析度 bit resolution,或稱為位元深度 bit depth。

數學表示法

\(x=\{x[n]\}, -∞<n<∞\),n 為整數

DSP 系統中,數位訊號通常為有限 finite 數字序列,可表示為

\(x=\{x[n]\}, n=0,1,2..., N-1\) ,這是從 t=0 開始的樣本,N為樣本數

import numpy as np
import matplotlib.pyplot as plt

n = np.array( [ 0, 1, 2, 3, 4, 5 ] )    # 定義n陣列
x = np.array( [ 1, 2, 4, 3, 2, 1 ] )    # 定義x陣列

plt.stem( n, x, use_line_collection=True )                      # 繪圖
plt.xlabel( 'n' )
plt.ylabel( 'x[n]' )
plt.show( )

基本的數位訊號

單位脈衝函數

Unit Impulse Function 也稱為狄拉克 δ 函數 (Dirac Delta Function),以英國物理學家 Paul Dirac 命名,定義:

\(δ(t)=0, t \neq 0\) ,且 \(\int_{-∞}^{∞}δ(t) {\rm d}t = 1\)

單位脈衝函數是連續時間函數,除了原點 t=0 以外,函數值均為 0,在時間域的積分(總面積)為 1,且集中在原點

單位脈衝

在離散時間域的 Unit Impulse 定義為

\( δ[n] = \left\{ \begin{array}{ll} 1, & \mbox{if n=0} \\ 0, & \mbox{if n≠0} \end{array} \right.\)

單位步階函數

Unit Step Function 也稱為黑維塞步階函數 (Heaviside Step Function),以英國科學家 Oliver Heaviside 命名,定義為

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

單位步階

Unit Step 定義為

\( μ[n] = \left\{ \begin{array}{ll} 1, & \mbox{if n ≥0} \\ 0, & \mbox{if n<0} \end{array} \right.\)


單位脈衝的時間延遲 Time Delay 定義為 \(δ[n-n_0]\)

ex: \(δ[n-2]\) 的圖形為

因此 Unit Step 也可表示為

\( μ[n] = δ[n]+δ[n-1]+δ[n-2]+... = \sum_{k=0}^{n}δ[n-k]\)

單位脈衝 Unit Impulse 也可表示為

\(δ[n]=μ[n]-μ[n-1]\)

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

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

數位語音檔

波型音訊檔案格式 Waveform Audio File Format、wave,是 Microsoft 與 IBM 制定的音訊編碼格式。另外 Apple Mac 定義了 aiff 音訊格式。

wav 是根據 Resource Interchange File Format(RIFF) ,為了存取 CD 數位音樂而設計的。前面是 header 共44 bytes,內容由多個區塊 chunks 組成,每個區塊是 4 bytes

起始位址(Bytes) 位址範圍 區塊名稱 區塊大小 ( Bytes) 內容 說明
0 1 ~ 4 區塊編號 4 RIFF Marks the file as a riff file.
4 5 ~ 8 總區塊大小 4 N + 36 File size (integer) Size of the overall file - 8 bytes
8 9 ~ 12 檔案格式 4 "WAVE" File Type Header. For our purposes, it always equals "WAVE".
12 13~16 子區塊 1 標籤 4 "fmt " Format chunk marker. Includes trailing null
16 17~20 子區塊 1 大小 4 16 Length of format data as listed above
20 21~22 音訊格式 2 1(PCM) Type of format (1 is PCM)
22 23~24 通道數量 2 1: (mono), 2: stereo Number of Channels
24 25~28 取樣頻率 4 Hz Sample Rate, Common values are 44100 (CD), 48000 (DAT). Sample Rate = Number of Samples per second, or Hertz
28 29~32 位元組 ByteRate 4 ByteRate = (Sample Rate * BitsPerSample * Channels) / 8 取樣頻率 * 位元深度 / 8
32 33~34 區塊對齊 BlockAlign 2 4 (BitsPerSample * Channels) / 8 The number of bytes for one sample including all channels ex: 1 -> 8 bit mono, 2 -> 8 bit stereo/16 bit mono, 4 -> 16 bit stereo
34 35~36 位元深度 2 16 Bits per sample
36 37~40 子區塊2標籤 4 "data" "data" chunk header. Marks the beginning of the data section.
40 41~44 子區塊2大小 4 N Size of the data section, i.e. file size - 44 bytes header.
44 45~ 資料 N 音訊資料

wav 沒有採用壓縮技術,檔案大小比其他格式大,但 wav 不會發生失真的狀況。

import wave

filename = input( "Please enter file name: " )
wav = wave.open( filename, 'rb' )

num_channels = wav.getnchannels( )  # 通道數
sampwidth   = wav.getsampwidth( )   # 樣本寬度
frame_rate  = wav.getframerate( )   # 取樣率
num_frames  = wav.getnframes( )     # 音框數
comptype    = wav.getcomptype( )    # 壓縮型態
compname    = wav.getcompname( )    # 壓縮名稱

print( "Number of Channels =", num_channels )
print( "Sample Width =", sampwidth )
print( "Sampling Rate =", frame_rate )
print( "Number of Frames =", num_frames )
print( "Comptype =", comptype )
print( "Compname =", compname )

wav.close( )

要讀取 wav 的資料,可改用 SciPy 的 io.wavfile 套件,讀取後以 numpy 的陣列回傳

from scipy.io import wavfile
import matplotlib.pyplot as plt

filename = input( "Please enter file name: " )
sampling_rate, x = wavfile.read( filename )

plt.plot( x )
plt.xlabel( 'n' )
plt.ylabel( 'Amplitude' )

plt.show( )

可直接取得 microphone 語音的波形

pip install PyAudio
import numpy as np
import pyaudio
import matplotlib.pyplot as plt

# sampling rate
fs = 11000
# 每次擷取樣本數 1024
CHUNK = 1024
pa = pyaudio.PyAudio( )
stream = pa.open( format = pyaudio.paInt16, channels = 1, rate = fs, 
                  input = True, output = False, frames_per_buffer = CHUNK )
            
try:            
    while True:                     
        data = stream.read( CHUNK )
    # 字串轉換為 int16
        x = np.fromstring( data, np.int16 )

        plt.clf( )
        plt.plot( x )
        plt.axis( [ 0, CHUNK, -30000, 30000 ] )

        plt.pause( 0.1 )

except KeyboardInterrupt:
    print( "Quit" )
    pa.close( stream )
    quit( )

References

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