DSP筆記6_功率頻譜密度 Power Spectral Density(PSD)_Periodogram(週期圖)最原始簡單的作法

事實上，在上一篇我們完全不依賴scipy手刻的PSD計算的python程式

它的演算法流程就是所謂的Periodogram(週期圖)最原始簡單的作法。

periodogram（週期圖）是一種用來估計訊號功率頻譜密度

• 使用單一區段全部資料來計算 PSD
• 雖然速度快，但具有對雜訊敏感且變異數較大的限制。
• 顯示訊號的功率如何分布在不同頻率上。
• 常用於signal processing、vibration analysis、PPG/EEG/ECG 分析

在數學上，對於一個長度為 $N$ 的離散時間訊號 $x[n]$：

會用以下數學式子表示週期圖功率頻譜密度

推導過程

Step1.計算離散傅立葉轉換(實際應用場景中，可用快速傅立葉轉換更有效率做出來。)

Step2.計算功率譜，也就是計算「振幅(幅度大小)平方」

Step3.除以樣本數N做正規化

原始periodogram功率頻譜密度

提供每個樣本的功率估計值。

原始 periodogram就是最初版本、最簡單的 PSD 估計方法。

在實務上，通常使用 FFT 來有效率地進行計算。

原始periodogram功率頻譜密度特性

• 變異數 : 高；即使增加 \(N\)，變異數也不會下降，因此結果通常會比較多雜訊干擾。
• 偏差(Bias)、頻譜洩漏（Spectral Leakage）: 在有限長度訊號中的不連續現象，會使功率擴散到鄰近的頻率。
• 解析度(Resolution)

$\mathrm{\Delta }f=\frac{f_s}{N}$；較長的訊號可以提升解析度。

• 一致性（Consistency）: 原始 periodogram是不一致的估計量

原始periodogram功率頻譜密度限制

• High variance : 峰值會產生波動
• Spectral leakage : 有限長度資料會造成頻譜能量向周圍頻率擴散。
• Poor consistency : N→∞時，原始 periodogram 不會收斂。
• Trade-off between bias and resolution : 加窗（Windowing）可以降低頻譜洩漏，但會使頻譜峰值變寬。

常見改善方法

加窗（Windowing）： Hann、Hamming、Blackman 視窗可降低頻譜洩漏（spectral leakage）。

Welch 方法（Welch's Method）： 將訊號分割成彼此重疊的多個區段（overlapping segments），再將各區段的 periodogram 取平均，以降低變異數（variance）。

多窗法（Multitaper Method）： 使用多個彼此正交的 taper，以在偏差與變異數之間取得平衡。

補零（Zero-padding）： 改善 FFT 視覺化時的頻率解析度（frequency resolution）。

練習程式碼ver1.這邊拿上一篇正弦波ver4案例改寫來觀測

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

# Sinusoid parameters
fs = 100  # Sampling freq (Hz)
t = np.arange(0, 1, 1/fs)
f0 = 10   # signal frequency

# Sinusoidal signal
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頻率。

Ref:

10. The Periodogram

https://www.youtube.com/watch?v=G0u3W5fEbg8&list=PLMxwp487ysbNGGLE5rtgfyvRk1u08rbfU&index=11