傅立葉級數 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
沒有留言:
張貼留言