DSP筆記6_功率頻譜密度 Power Spectral Density(PSD)_Periodogram(週期圖)_用scipy內建的periodogram跟最原始簡單的作法比較

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

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

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

• periodogram（週期圖）是一種非參數方法（nonparametric method）
• 最早由 Schuster 於 1898 年提出，用於研究太陽黑子數量（sunspot numbers）中的週期性。
• 使用單一區段全部資料來計算 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頻率。

會發覺跟原始手刻方法不一樣

根據原文書教材Digital Signal Processing,Fourth Edition(2006) John G. Proakis, Dimitris K Manolakis

在第14章的說明14.1.3

如果我們有N個資料點，則至少需要計算N-point DFT。

其對應的頻率為 

$f_k=k/N$。

然而，在實務上，對 spectrum 進行如此稀疏的取樣，無法對連續的 spectrum estimate $P_{xx}(f)$ 提供非常良好的表示或圖像。這個問題可以很容易地藉由在額外的頻率上計算 $P_{xx}(f)$ 來改善。

我們可以透過 zero padding，有效地增加 sequence 的長度，然後在一組更密集的頻率上計算 $P_{xx}(f)$。因此，如果我們透過 zero padding 將 data sequence 的長度增加到 $L$ 點，並計算 $L$-point DFT。

要再次強調 zero padding，以及在 $L>N$ 個點上計算 DFT，並不會改善 spectral estimate 中的 frequency resolution。它只是提供一種方法，讓我們可以在更多頻率位置上對量測到的 spectrum 數值進行 interpolation。spectral estimate $P_{xx}(f)$ 中的 frequency resolution，是由 data record 的長度 $N$ 所決定。

Ref:

Digital Signal Processing,Fourth Edition(2006) John G. Proakis, Dimitris K Manolakis-Ch14

(主要學習教材)

https://www.tenlong.com.tw/products/9780131873742

https://www.scribd.com/document/949970176/Digital-Signal-Processing-4th-Edition-Proakis-and-Manolakis

10. The Periodogram

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

What is Periodogram method? | Power Spectral Estimation | Signal processing | Signal Processing

https://www.youtube.com/watch?v=Wtudah2Nr_M

麻省理工學院2008年教材-2.161 Signal Processing: Continuous and Discrete-Lecture 23

https://ocw.mit.edu/courses/2-161-signal-processing-continuous-and-discrete-fall-2008/c7d1951d44fe78ea97d0e14ea5158aee_lecture_23.pdf