## 傅立葉級數

$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_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)。

$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}$

$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.$

$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( )

## 傅立葉轉換

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

$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)]$

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π})$

$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( )

## 平移定理

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(ω) = 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(e^{jω}) = \sum_{n=-∞}^{∞}x[n]e^{-jωn}$

​ $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

​ $|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

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

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ω}}$

## 離散傅立葉轉換

Discrete Fourier Transform (DFT)

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

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

​ $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]

$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\}$

### 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}$

$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_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}$

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

ex: DFT 具有可逆性

$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]

## 自相關 autocorrelation

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

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

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

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$

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]

## 串聯、並聯

$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]$ (卷積運算符合結合率與交換率)

$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

### 平均濾波器 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$

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( )

### 高斯濾波器 Gaussian Filter

$g[n]=e^{-n^2/2𝜎^2}$ 高斯函數的平均值 $𝜇=0$，通常介於 -3𝜎 ~ 3𝜎 之間，因此設定高斯濾波器的大小為 6𝜎+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 濾波

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())


## 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

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) )

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


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

$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())


## 雜訊生成

$y[n]=x[n]+η[n]$ 其中 η[n] 稱為雜訊 Noise

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

## 均勻雜訊 Uniform Noise

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

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

## 高斯雜訊 Gaussian Noise

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

## 布朗尼雜訊 Brownian Noise

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

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

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

## 以 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())


## 訊號雜訊比

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$

## 週期性訊號

$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))$

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.$

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)$

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$是高頻訊號的頻率

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( )

### 振幅調變

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

## 非週期性訊號

### 淡入與淡出

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

Chirp 經常稱為 Sweep Singal

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

• sinusoid() 弦波
• square() 方波
• sawtooth() 鋸齒波
• triangle() 三角波
• harmonic() 諧波
• beat() 節拍波
• 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);

wav_params = WavFileParams(

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()
chirp()

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

## 取樣與量化

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為整數

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

Nyquist-Shannon 取樣定理

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

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

## 數學表示法

$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$

### 單位脈衝

$δ[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.$

ex: $δ[n-2]$ 的圖形為

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

$δ[n]=μ[n]-μ[n-1]$

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

## 數位語音檔

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

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( )

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:
# 字串轉換為 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( )