頻譜分析 spectrum analysis 目的是分析訊號在不同頻率下的分量，以圖形表示。頻譜分析的方法可分成兩大類型：傅立葉頻譜 Fourier Spectrum 與 功率頻密度 Power Spectral Density (PSD)
根據傅立葉基本理論：「任意週期性函數，均可表示成不同頻率、不同振幅的正弦函數或餘弦函數，加總而得的無窮級數。」更進一步，傅立葉轉換也可適用於非週期性函數。
頻譜分析 spectrum analysis 的目的是將訊號分解成不同頻率的分量，並以圖形表示，稱為頻譜(spectrum)。頻譜可用來分析訊號中的頻率分量，例如男生的語音聲音低沈，屬於低頻，女生或兒童是高頻，可藉此辨識說話者是男生、女生或兒童。
弦波的定義為
\(x(t)=A cos(ω_0+ϕ)\) 或是 \(x(t)=A cos(2 \pi f_0t+ϕ)\) ，其中振幅 A ，角頻率 \(ω_0\)，頻率 \(f_0\)
根據反尤拉公式 \(cosθ = \frac{e^{jθ}+e^{-jθ}}{2}\)
代入可得
\(x(t) = A cos(ω_0t+ϕ) = A [e^{j(ω_0t+ϕ)}+e^{-j(ω_0t+ϕ)}]/2 \)
因此弦波包含了兩個振幅 Amplitude 的分量，當頻率為 \(ω_0, -ω_0\)，振幅均為 A/2，如果以圖形表示，就稱為振幅頻譜 Amplitude Spectrum。
同時弦波 x(t) 也包含了兩個相位移 Phase Shift 的分相，當頻率為 \(ω_0, -ω_0\)，相位移分別為 \(ϕ, -ϕ\)，如果以圖形表示，稱為相位頻譜 Phase Spectrum。
振幅頻譜及相位頻譜，通稱為傅立葉頻譜 Fourier Spectrum 或頻率頻譜 Frequency Spectrum，或簡稱 Spectrum
由於頻譜是以原點為中心分成兩邊，經常稱為雙邊頻譜 Doubl-Sided Spectrum。頻譜的橫軸以角頻率 ω為主，通常介於 -π ~ π 之間，單位為 弧度/秒 (Raidans/Second)，此外，頻譜的橫軸也經常以頻率 f 為主，單位為 Hz
若諧波由多個弦波組成，頻率為基礎頻率的整數倍數，定義為
\(x(t)=\sum_{k=1}^{N}A_k cos(2 \pi f_kt)\)，其中 \(f_k=k \cdot f_1\)，\(f_1\) 稱為基礎頻率
用反尤拉公式展開後得到振幅頻譜。對應不同頻率的分量，其頻率分量是弦波振幅的一半。
傅立葉頻譜
Fourier Spectrum 是指時間域訊號在頻率域的圖形表示法，可使用傅立葉轉換得來
先對輸入的數位訊號 \(x[n], n=0,1,...,N-1\) 進行離散傅立葉轉換
\(X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2 \pi kn/N}, k=0,1,...,N-1\)
輸出結果為 X[k]，是複數構成的離散數列。在 DSP 通常採用 Fast Fourier Transform (FFT) 縮短離散副業轉換的計算所需的時間
FFT 結果為複數，強度為 \(|X[k|, k=0,1,...,N-1\)，以圖形表示，稱為強度頻譜 Magnitude Spectrum，也就是前述的振幅頻譜 Amplitude Spectrum。
若求複數的幅角 Argument 或相位角 Phase Angle，並以圖形表示，則稱為相位頻譜 Phase Spectrum
\(arg(X[k]), k=0,1,...,N-1\)
ex: 數位訊號 \(x={1,2,4,3}, n=0,1,2,3\) ，求強度頻譜的圖形
離散傅立葉轉換結果為複數
\(X=\{10, -3+j, 0, -3-j\}\)
強度為 \(|X[0]|=10, |X[1]|=\sqrt{(-3)^2+1^2}=\sqrt{10}\)
Magnitude Spectrum:
弦波
定義 \(x(t)=cos(2 \pi 100 \cdot t)\)，振幅 A=1，頻率 f=100，相位移 ϕ=0，在 t=0~1s 之間取 1000 個樣本，取樣頻率 1000Hz，求弦波的傅立葉頻譜？
這邊先以強度頻譜為主，先不討論相位頻譜。由圖形結果可發現 Fourier Spectrum 在 k=100, 900 兩處出現峰值 peaks，形成 Double-Sided Spectrum 的對稱結構，但不是以原點為中心，其中 k=100 的峰值就是弦波的頻率，振幅為 \((A/2) \cdot N = 500\)
為方便傅立葉頻譜的觀察與分析，通常會對 FFT 的結果進行平移 shift，讓頻譜的中心為原點。numpy 的 fftshift 有實現平移的功能。因為 N=1000，根據 Nyquist Theorem，可分析的最大頻率為 500 Hz
import numpy as np
from numpy.fft import fft, fftshift, fftfreq
import matplotlib.pyplot as plt
# time，1000 個 sample 點
t = np.linspace( 0, 1, 1000, endpoint = False )
# 弦波 原始訊號
x = np.cos( 2 * np.pi * 100 * t )
## fft 計算
X = fft( ( x ) )
# 強度
Xm = abs( X )
# 強度頻譜
plt.subplot(121)
plt.plot( Xm )
plt.xlabel( 'k' )
plt.ylabel( 'Magnitude' )
# t = np.linspace( 0, 1, 1000, endpoint = False )
# x = np.cos( 2 * np.pi * 100 * t )
# numpy.fft.fftfreq(n, d=1.0) 定義 Fourier Spectrum 的取樣頻率
# n 為樣本數， d 為取樣週期 (1 second)
# fftshift 用來平移，以原點為中心
f = fftshift( fftfreq( 1000, 0.001 ) )
X = fftshift( fft( x ) )
Xm = abs( X )
plt.subplot(122)
plt.plot( f, Xm )
plt.xlabel( 'f' )
plt.ylabel( 'Magnitude' )
plt.show( )
方波
若方波的參數為 A=1, f=10Hz，在一秒內震盪 10 次，若取樣頻率為 1000Hz，求方波的 Fourier Spectrum ?
方波由無限多個弦波加總而成，可視為一種諧波，基礎頻率為 10Hz，也就是方波的頻率
import numpy as np
import scipy.signal as signal
from numpy.fft import fft, fftshift, fftfreq
import matplotlib.pyplot as plt
# time，1000 個 sample 點
t = np.linspace( 0, 1, 1000, endpoint = False )
# 方波 原始訊號
x = signal.square( 2 * np.pi * 10 * t )
# numpy.fft.fftfreq(n, d=1.0) 定義 Fourier Spectrum 的取樣頻率
# n 為樣本數， d 為取樣週期 (1 second)
# fftshift 用來平移，以原點為中心
f = fftshift( fftfreq( 1000, 0.001 ) )
X = fftshift( fft( x ) )
Xm = abs( X )
plt.plot( f, Xm )
plt.xlabel( 'f' )
plt.ylabel( 'Magnitude' )
plt.show( )
節拍波
beat wave是由兩個不同頻率的弦波乘法運算得來
\(x(t)= A \cdot cos(2 \pi f_1 \cdot t) \cdot cos(2 \pi f_2 \cdot t)\)
通常 \(f_1\) 為低頻訊號的頻率。\(f_2\) 是高頻訊號的頻率，稱為載波頻率 carrier frequency，也常表示為 \(f_c\)。乘法運算是調變 modulation 的主軸，常被用在 DSP 系統中
ex: 節拍波定義為
\(x(t)= \cdot cos(2 \pi 20 \cdot t) \cdot cos(2 \pi 200 \cdot t)\)，其中 A=1，\(f_1=20Hz, f_2=200Hz\)，取樣頻率為 1000Hz，求 Fourier Spectrum?
import numpy as np
import scipy.signal as signal
from numpy.fft import fft, fftshift, fftfreq
import matplotlib.pyplot as plt
f1 = 20
f2 = 200
t = np.linspace( 0, 1, 1000, endpoint = False )
# 節拍波 原始訊號
x = np.cos( 2 * np.pi * f1 * t ) * np.cos( 2 * np.pi * f2 * t )
f = fftshift( fftfreq( 1000, 0.001 ) )
X = fftshift( fft( x ) )
Xm = abs( X )
plt.plot( f, Xm )
plt.xlabel( 'f' )
plt.ylabel( 'Magnitude' )
plt.show( )
啁啾訊號
這是非週期性訊號
ex: 啁啾訊號在時間 0~1s 之間，頻率範圍是 0Hz 線性遞增到 100Hz，取樣頻率 1000Hz，求 Fourier Spectrum?
頻譜呈現 ripple 現象
import numpy as np
import scipy.signal as signal
from numpy.fft import fft, fftshift, fftfreq
import matplotlib.pyplot as plt
t = np.linspace( 0, 1, 1000, endpoint = False )
x = signal.chirp( t, 0, 1, 100, 'linear' )
f = fftshift( fftfreq( 1000, 0.001 ) )
X = fftshift( fft( x ) )
Xm = abs( X )
plt.plot( f, Xm )
plt.xlabel( 'f' )
plt.ylabel( 'Magnitude' )
plt.show( )
功率頻密度 Power Spectral Density (PSD)
數位訊號的能量 Energy 可定義為
\(E=\int_{-∞}^{∞}|x(t)|^2 {\rm d} t\)
但只適合用來描述有限的訊號。當訊號能量集中在有限時間區間時，總能量有限，因此可透過積分計算而得。
以離散時間域來說，數位訊號的能量可用以下公式計算：
\(E=\sum_{-∞}^{∞}x^2[n]\)
帕塞瓦爾定理
若訊號 \(x(t)\) 的傅立葉轉換為 \(X(f)\)，也就是 \(X(f)=\int_{-∞}^{∞}x(t)e^{-j(2 \pi ft)} {\rm d}t\)
則訊號的能量滿足以下公式：
\(E=\int_{-∞}^{∞}|x(t)|^2 {\rm d}t = \int_{-∞}^{∞}|X(f)|^2 {\rm d}f\)
「時間域訊號，若取平方總和（或積分），則結果與傅立葉轉換後頻率域訊號的平方總和（或積分）相等」
換句話說，如果想分析訊號的能量，除了在時間域進行積分外，也可以在頻率域進行積分，結果一樣。
證明：
根據反傅立葉轉換
\(x(t)=\frac{1}{2\pi}\int_{-∞}^{∞}X(ω)e^{jωt} {\rm d}ω\)
則 \(E = \int_{-∞}^{∞} x(t)[\frac{1}{2\pi} \int_{-∞}^{∞}X(ω) e^{jωt} {\rm d}ω ] {\rm d}t \\ = \frac{1}{2 \pi} \int_{-∞}^{∞}X(ω) [\int_{-∞}^{∞} x(t) e^{jωt} {\rm d}t]{\rm d}ω \\ = \frac{1}{2 \pi} \int_{-∞}^{∞}X(ω) X(-ω){\rm d}ω \\ = \frac{1}{2 \pi} \int_{-∞}^{∞}X(ω) X^*(ω){\rm d}ω \quad\quad X^*(ω) 是 X(ω) 的共軛複數 \\ = \frac{1}{2 \pi} \int_{-∞}^{∞} |X(ω)|^2 {\rm d}ω \\ = \frac{1}{2 \pi} \int_{-∞}^{∞} |X(f)|^2 {\rm d}f\)
能量譜密度
若類比訊號 x(t) 的傅立葉轉換為 X(f)，則能量譜密度 (Energy Spectral Density) 定義為：
\(X_{xx}(f) = |X(f)|^2\)
能量譜密度的目的是用來描述訊號在不同頻率下的能量分佈狀況，以圖形表示
能量譜密度可根據 自相關 Autocorrelation 取其傅立葉轉換得來
假設訊號 x(t) 的自相關定義為
\(R(τ)=\int_{-∞}^{∞}x^*(t) \cdot x(t+τ) {\rm d}t\)
則其傅立葉轉換為：
\(F\{R(τ)\} = \int_{-∞}^{∞}R(τ)e^{-j2 \pi ft} {\rm d}τ \\ = \int_{-∞}^{∞} [\int_{-∞}^{∞} x^*(t) \cdot x(t+τ) {\rm d}t ] e^{-j2 \pi ft} {\rm d}τ \\ = \int_{-∞}^{∞} x^*(t) [\int_{-∞}^{∞} x(t+τ) e^{-j2 \pi ft} {\rm d}τ ] {\rm d}t \quad\quad 假設 \hatτ = t+τ, {\rm d}\hatτ = {\rm d}τ \\ = \int_{-∞}^{∞} x^*(t) [\int_{-∞}^{∞} x(\hatτ) e^{-j2 \pi f \hatτ} \cdot e^{j2 \pi f t} {\rm d}τ ] {\rm d}t \\ = \int_{-∞}^{∞} x^*(t) [\int_{-∞}^{∞} x(\hatτ) e^{-j2 \pi f \hatτ} {\rm d}τ ] e^{j2 \pi f t} {\rm d}t \\ = \int_{-∞}^{∞} x^*(t) X(f) e^{j2 \pi f t} {\rm d}t \\ = X(f) \int_{-∞}^{∞} x^*(t) e^{j2 \pi f t} {\rm d}t \\ = X(f) X^*(f) \\ = |X(f)|^2 \\ = S_{xx}(f)\)
功率頻密度
功率頻密度 Power Spectral Density 就是訊號在不同頻率下的功率分佈圖
上面的 能量譜密度 定義，只適用於能量集中在某個時間區間的訊號，例如：脈衝訊號。對於持續存在的連續訊號，或是平穩過程 Stationary Process 的訊號，要改用 功率頻密度 Power Spectral Density (PSD)，用來描述訊號在不同頻率下，功率 Power 的分佈狀況
評估 功率頻密度 有兩個具代表性的方法：
週期圖 Periodogram
週期圖是評估功率頻密度最基本的方法，主要根據傅立葉轉換的結果，取平方得來
Welch Method
由 P. D. Welch 提出，也是基於週期圖的方法，但同時參考 Barlett's Method，牽涉功率的平均值運算，又稱為平均週期圖法 Method of Averaged Periodograms
ex: 若產生時間介於 0~1s 的均勻雜訊，值介於 -1~1 之間，取樣頻率 1000Hz，使用週期圖語 Welch Method 求功率頻密度 PSD
這兩種方法得到的功率頻密度，頻率範圍是 0~500Hz (因為 Nyquist Theorem，上限是 500 Hz)，在不同頻率下，均勻雜訊都含有少量功率
若以比較嚴謹的定義，理想的 white noise，其功率頻密度在不同頻率下，應該要為固定常數。因此均勻雜訊，並不完全符合 white noise 的條件
import numpy as np
import numpy.random as random
import scipy.signal as signal
import matplotlib.pyplot as plt
fs = 1000
t = np.linspace( 0, 1, fs, endpoint = False )
# 均勻雜訊
noise = random.uniform( -1, 1, fs )
# 用兩種方法求得 PSD
f1, pxx1 = signal.periodogram( noise, fs )
f2, pxx2 = signal.welch( noise, fs )
plt.subplot(121)
plt.plot( f1, pxx1 )
plt.xlabel( 'frequency (Hz)' )
plt.ylabel( 'PSD' )
plt.subplot(122)
plt.plot( f2, pxx2 )
plt.xlabel( 'frequency (Hz)' )
plt.ylabel( 'PSD' )
plt.show( )
ex: 訊號定義為
\(x(t) = 10 cos(2 \pi 100 t) + 5 cos(2 \pi 200 t) + η(t)\)
其中包含兩個弦波，η(t) 是介於 -1~1 之間的均勻雜訊，取樣頻率 1000 Hz，用 週期圖與 Welch Method 求 PSD
兩種方法呈現的 PSD 都是單邊頻譜 Single-Sided Spectrum，範圍是 0~500 Hz。
雖然有雜訊干擾， PSD 在 100Hz 與 200Hz 呈現功率的峰值
概括來說，Welch Method 在頻率域上的解析度較差，但具有較好的抗雜訊能力
import numpy as np
import numpy.random as random
import scipy.signal as signal
import matplotlib.pyplot as plt
fs = 1000
t = np.linspace( 0, 1, fs, endpoint = False )
# 原始訊號
x = 10 * np.cos( 2 * np.pi * 100 * t ) + 5 * np.cos( 2 * np.pi * 200 * t ) + random.uniform( -1, 1, fs )
f1, pxx1 = signal.periodogram( x, fs )
f2, pxx2 = signal.welch( x, fs )
plt.subplot(121)
plt.plot( f1, pxx1 )
plt.xlabel( 'frequency (Hz)' )
plt.ylabel( 'PSD' )
plt.subplot(122)
plt.plot( f2, pxx2 )
plt.xlabel( 'frequency (Hz)' )
plt.ylabel( 'PSD' )
plt.show( )
也可再對 PSD 取對數，查看主要頻率外，其他頻率的分佈狀況
import numpy as np
import numpy.random as random
import scipy.signal as signal
import matplotlib.pyplot as plt
fs = 1000
t = np.linspace( 0, 1, fs, endpoint = False )
x = 10 * np.cos( 2 * np.pi * 100 * t ) + 5 * np.cos( 2 * np.pi * 200 * t )
noise = random.uniform( -1, 1, fs )
y = x + noise
f1, pxx1 = signal.periodogram( y, fs )
f2, pxx2 = signal.welch( y, fs )
# 取對數
lp1 = 10 * np.log10( pxx1 )
lp2 = 10 * np.log10( pxx2 )
plt.subplot(121)
plt.plot( f1, lp1 )
plt.xlabel( 'frequency (Hz)' )
plt.ylabel( 'PSD' )
plt.axis( [ 0, 500, -50, 10 ] )
plt.subplot(122)
plt.plot( f2, lp2 )
plt.xlabel( 'frequency (Hz)' )
plt.ylabel( 'PSD' )
plt.axis( [ 0, 500, -40, 10 ] )
plt.show( )
沒有留言:
張貼留言