DSP筆記5_功率頻譜密度 Power Spectral Density(PSD)_從時域訊號到幅度大小的平方再除以資料總長度

功率頻譜密度（Power Spectral Density, PSD）

用來描述訊號在不同頻率下的功率分布情形。

功率頻譜密度（Power Spectral Density, PSD） $S_x(f)$ 適用於功率訊號（power signals）以及廣義平穩程序（wide-sense stationary processes）。

定義，透過自相關函數（Autocorrelation）根據（Wiener–Khinchin theorem，維納–辛欽定理）可得到如下式子

另一種定義(此公式十分重要)：

總功率與性質（Total Power and Properties）

換言之，總功率（Total Power）可以由功率頻譜密度 $S_x(f)$ 對所有頻率積分後得到。

把給定區間的功率譜密度全部加總起來

練習程式ver1.建立一個典型的10Hz正弦波訊號

import numpy as np
import matplotlib.pyplot as plt

# Sinusoid parameters
fs = 100  # Sampling freq (Hz)每秒取樣 100 次。
t = np.arange(0, 1, 1/fs)#[0.   0.01 0.02 ... 0.98 0.99]，每次增加 0.01。
f0 = 10   # signal frequency

# Sinusoidal signal
x = np.sin(2 * np.pi * f0 * t)

# Plot
plt.figure(figsize=(8,4))
plt.plot(t, x)

plt.title("10 Hz Sinusoid")
plt.xlabel("Time (s)")
plt.ylabel("Amplitude")

plt.grid(True)
plt.show()

現在產生的是：

$$x(t)=\sin(2\pi f_0t)$$

其中

$$f_0=10\text{ Hz}$$

意思是一秒鐘會完成 10 個週期。

因此每兩個取樣點之間的時間間隔為：

$$T_s=\frac{1}{f_s} =\frac{1}{100} =0.01\text{ 秒}$$

也就是每隔 0.01 秒取一個資料點。

運行效果

練習程式ver2.做FFT觀察它由哪些頻率成分組成。

import numpy as np
import matplotlib.pyplot as plt

# 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)

# FFT
X = np.fft.fft(x)
# Frequency axis
freq = np.fft.fftfreq(len(t), 1/fs)#手動 FFT 版本會得到：0, 1, 2, ..., 49, -50, -49, ..., -1

# Plot
plt.figure(figsize=(8,4))
plt.plot(freq[:len(freq)//2], np.abs(X[:len(X)//2]))#只取前半段正頻率。

plt.title("FFT of a 10 Hz Sinusoid")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Magnitude")

plt.grid(True)
plt.show()

由於我們產生的是 10 Hz 正弦波，所以應該會看到頻譜主要集中在10Hz

運行效果

練習程式ver3.從 FFT 進一步到 ∣X(f)∣取平方

import numpy as np
import matplotlib.pyplot as plt

# 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)

# FFT
X = np.fft.fft(x)
# Frequency axis
freq = np.fft.fftfreq(len(t), 1/fs)
# Magnitude squared
Power = np.abs(X)**2

# Plot
plt.figure(figsize=(8,4))
plt.plot(freq[:len(freq)//2], Power[:len(freq)//2])

plt.title("Squared Magnitude Spectrum")
plt.xlabel("Frequency (Hz)")
plt.ylabel("|X(f)|^2")

plt.grid(True)
plt.show()

此階段是為了算出PSD之前的前一個動作，也就是取出FFT magnitude 的平方。

可理解為：

$$x(t) \rightarrow X(f) \rightarrow |X(f)| \rightarrow |X(f)|^2$$

但此時結果還不是PSD

運行結果

練習程式ver4.Periodogram estimate

(加入除以資料長度 N，得到最終PSD。)

import numpy as np
import matplotlib.pyplot as plt

# 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)

# FFT
X = np.fft.fft(x)

# Frequency axis
freq = np.fft.fftfreq(len(t), 1/fs)

# Magnitude squared
Power = np.abs(X)**2

# PSD
PSD = Power / len(t)

# Plot
plt.figure(figsize=(8,4))
plt.plot(freq[:len(freq)//2], PSD[:len(freq)//2])

plt.title("PSD of a 10 Hz Sinusoid")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Power")

plt.grid(True)
plt.show()

運行效果

整個算法順序可以理解如下：

$$x(t) \rightarrow X(f) \rightarrow |X(f)| \rightarrow |X(f)|^2 \rightarrow \frac{|X(f)|^2}{N}$$

也就是：

$$\text{時域訊號}\to \text{FFT}\to \text{Magnitude}\to \text{Squared Magnitude}\to \text{PSD}$$

根據取樣頻率，每秒100次取樣可得知，經過這一秒時間區間中資料總長度為100。

根據上述計算公式可得知幅度大小平方為2500除以資料總長100就又回到25了。

那事實上我們第四版的程式碼，還有另一種寫法就是去直接用Scipy內建的periodogram的方法。

會在下一篇做更進一步的介紹。

Ref:

9. Concept of the Power Spectrum

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

10. The Periodogram

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