DSP筆記6_功率頻譜密度 Power Spectral Density(PSD)_Periodogram(週期圖)_用scipy內建的periodogram跟最原始簡單的作法比較
事實上，在上一篇我們完全不依賴scipy手刻的PSD計算的python程式
它的演算法流程就是所謂的Periodogram(週期圖)最原始簡單的作法。
periodogram（週期圖）是一種用來估計訊號功率頻譜密度
- 使用單一區段全部資料來計算 PSD
- 雖然速度快，但具有對雜訊敏感且變異數較大的限制。
- 顯示訊號的功率如何分布在不同頻率上。
- 常用於signal processing、vibration analysis、PPG/EEG/ECG 分析
在數學上，對於一個長度為 的離散時間訊號 ：
會用以下數學式子表示週期圖功率頻譜密度
Step1.計算離散傅立葉轉換(實際應用場景中，可用快速傅立葉轉換更有效率做出來。)
將時域訊號(比方交響樂)轉換到頻域(各不同頻率的純音)
Step2.計算功率譜，也就是計算「振幅(幅度大小)平方」，這邊會去計算每個頻率成分的模數平方。也就是Magnitude Squared，會得到各頻率純音的能量。
Step2.計算功率譜，也就是計算「振幅(幅度大小)平方」，這邊會去計算每個頻率成分的模數平方。也就是Magnitude Squared，會得到各頻率純音的能量。
Step3.除以樣本數N做正規化，確保功率頻譜的總能量和原始訊號的功率一致。
Step4.針對特定頻帶進行積分計算，以取得該頻帶總功率密度特性
原始periodogram功率頻譜密度特性
- 變異數 : 高；即使增加 \(N\)，變異數也不會下降，因此結果通常會比較多雜訊干擾。
- 偏差(Bias)、頻譜洩漏（Spectral Leakage）: 在有限長度訊號中的不連續現象，會使功率擴散到鄰近的頻率。
- 解析度(Resolution)
| ；較長的訊號可以提升解析度。 |
- 一致性（Consistency）: 原始 periodogram是不一致的估計量
原始periodogram功率頻譜密度限制
- High variance : 峰值會產生波動
- Spectral leakage : 有限長度資料會造成頻譜能量向周圍頻率擴散。
- Poor consistency : N→∞時，原始 periodogram 不會收斂。
- Trade-off between bias and resolution : 加窗（Windowing）可以降低頻譜洩漏，但會使頻譜峰值變寬。
常見改善方法
練習程式碼ver1.這邊拿上一篇正弦波ver4案例改寫來觀測
import numpy as np import matplotlib.pyplot as plt from scipy.signal import periodogram # Sinusoid parameters fs = 100 # Sampling freq (Hz)，每秒取樣 100 次，取樣時間間隔為：0.01秒。 t = np.arange(0, 1, 1/fs) #建立時間序列：t=[0,0.01,0.02,...,0.99]一共產生100個取樣點。 f0 = 10 # signal frequency 頻率為 10 Hz，正弦波每秒完成 10 次週期。
# x(t) = A*sin(2*pi*f*t+theta)# Sinusoidal signal，振幅A為1，頻率f為10，初始相位:0 x = np.sin(2 * np.pi * f0 * t) # Periodogram #SciPy 對實數訊號預設使用return_onesided=True也就是只回傳正頻率，並把負頻率那一半的功率合併到正頻率 freq, PSD = periodogram(x,fs=fs,scaling='density')#預設只取前半段正頻率。 # Plot plt.figure(figsize=(8,4)) plt.plot(freq, PSD) plt.title("PSD of a 10 Hz Sinusoid") plt.xlabel("Frequency (Hz)") plt.ylabel("Power") plt.grid(True) plt.show()
運行效果
可以觀察到原先訊號的取樣頻率100Hz會自動頻率範圍轉換為0~50Hz，滿足Nyquist頻率。
根據 Nyquist 頻率，f=fs/2=100/2=50Hz
[0, 1, 2, 3, ..., 49, 50]
此外，periodogram() 在本程式使用的幾個預設參數為：
- window='boxcar' 使用矩形窗，也就是沒有額外加權的窗函數。
- detrend='constant' 先扣除訊號平均值。本範例正弦波平均值約為 0。
- return_onesided=True 因為輸入為實數訊號，所以只傳回單邊頻譜。 會將對應的負頻率貢獻合併到正頻率。
週期圖可以由有限長度自相關函數的 DTFT 求得。反過來，自相關函數則是週期圖的反 DTFT。這也是 Wiener–Khinchin 定理在譜估計上的重要應用。
將程式碼中的變數對應套用式子數學符號如下
核心公式：
\[
P_{x,M}(\omega)
=
\frac{1}{M}
\left|
\sum_{n=0}^{M-1}x_w(n)e^{-j\omega n}
\right|^2
\]
先代入：
\[
M=100
\]
\[
x_w(n)=\sin(0.2\pi n)
\]
得到：
\[
\boxed{
P_{x,100}(\omega)
=
\frac{1}{100}
\left|
\sum_{n=0}^{99}
\sin(0.2\pi n)e^{-j\omega n}
\right|^2
}
\]
接著我們需要找出 10 Hz 對應的數位角頻率。
步驟 1：計算數位角頻率
數位角頻率與 Hz 的轉換關係：
\[
\omega=2\pi\frac{f}{f_s}
\]
代入 \(f=10\)、\(f_s=100\)：
\[
\omega=2\pi\frac{10}{100}
\]
得到：
\[
\boxed{\omega=0.2\pi\ \text{rad/sample}}
\]
步驟 2：將頻率代入 DTFT 公式
\[
P_{x,100}(0.2\pi)
=
\frac{1}{100}
\left|
\sum_{n=0}^{99}
\sin(0.2\pi n)e^{-j0.2\pi n}
\right|^2
\]
這個式子可以繼續進行紙筆推導。
利用歐拉公式：
\[
\sin(\theta)=\frac{e^{j\theta}-e^{-j\theta}}{2j}
\]
將正弦函數展開：
\[
\begin{aligned}
&\sum_{n=0}^{99}\sin(0.2\pi n)e^{-j0.2\pi n}\\
&=\frac{1}{2j}\sum_{n=0}^{99}
\left(1-e^{-j0.4\pi n}\right)
\end{aligned}
\]
拆成兩個求和：
\[
=\frac{1}{2j}
\left[
\sum_{n=0}^{99}1
-
\sum_{n=0}^{99}e^{-j0.4\pi n}
\right]
\]
第一項：
\[
\sum_{n=0}^{99}1=100
\]
第二項是一個完整週期整數倍的複數指數序列之和：
\[
\sum_{n=0}^{99}e^{-j0.4\pi n}=0
\]
因此：
\[
X(0.2\pi)=\frac{100}{2j}
\]
因為 \(1/j=-j\)：
\[
\boxed{X(0.2\pi)=-50j}
\]
步驟 3：計算平方幅值
\[
|X(0.2\pi)|^2=|-50j|^2
\]
因此：
\[
|X(0.2\pi)|^2=2500
\]
帶回 Periodogram 公式：
\[
P_{x,100}(0.2\pi)
=\frac{1}{100}\times2500
\]
得到：
\[
\boxed{P_{x,100}(0.2\pi)=25}
\]
這裡算出來是 25，而不是 Python 圖表的 0.5。
原因在於公式是以數位角頻率 \(\omega\) 表示的雙邊週期圖，SciPy 的輸出則是以 Hz 為單位的單邊 PSD。
從式子計算出來的25轉為Scipy產出的0.5
(從單位角頻率轉換到Hz的單邊PSD，正負加總起來轉為單邊PSD。)
第一步：轉換為每 Hz 的雙邊 PSD
對於本範例的矩形窗：
\[
PSD_{\text{two-sided}}(f)
=
\frac{P_{x,M}(\omega)}{f_s}
\]
代入：
\[
PSD_{\text{two-sided}}(10)
=
\frac{25}{100}
\]
得到：
\[
\boxed{PSD_{\text{two-sided}}(10)=0.25}
\]
而在 \(-10\) Hz 也有相同結果：
\[
PSD_{\text{two-sided}}(-10)=0.25
\]
第二步：轉換為單邊 PSD
須把 10 Hz 與 -10 Hz 的功率譜密度加總：
\[
PSD_{\text{one-sided}}(10)
=0.25+0.25
\]
最後得到：
\[
\boxed{PSD_{\text{one-sided}}(10)=0.5}
\]
Ref:
The Periodogram
Welch's Method
3.2 Non-parametric spectral estimation
What is Periodogram method? | Power Spectral Estimation | Signal processing | Signal Processing
10. The Periodogram
离散时间信号处理DSP-吴炳飞
留言
張貼留言