使用Matplotlib绘制功率谱密度图：Python数据可视化指南

使用Matplotlib绘制功率谱密度图：Python数据可视化指南

参考：Plot the power spectral density using Matplotlib – Python

功率谱密度（Power Spectral Density，PSD）是信号处理和数据分析中的重要工具，它可以帮助我们了解信号在频域中的能量分布。Matplotlib是Python中最流行的数据可视化库之一，它提供了强大的功能来绘制各种类型的图表，包括功率谱密度图。本文将详细介绍如何使用Matplotlib绘制功率谱密度图，并提供多个实用示例。

1. 功率谱密度简介

功率谱密度是描述信号功率如何随频率分布的一种方法。它在许多领域都有广泛应用，如通信、音频处理、地震学等。PSD可以帮助我们识别信号中的主要频率成分，检测周期性模式，以及分析噪声特性。

在开始绘制PSD图之前，我们需要了解一些基本概念：

1. 频率：信号在单位时间内完成周期性变化的次数。
2. 功率：信号能量在单位时间内的变化率。
3. 谱密度：描述信号功率或能量如何分布在不同频率上。

2. Matplotlib简介

Matplotlib是一个综合性的Python绘图库，它提供了一个类似MATLAB的绘图框架。使用Matplotlib，我们可以创建各种静态、动态和交互式图表。它的主要特点包括：

1. 高度可定制化
2. 支持多种输出格式（PNG, PDF, SVG, EPS等）
3. 良好的跨平台支持
4. 与NumPy和SciPy等科学计算库无缝集成

3. 准备工作

在开始绘制PSD图之前，我们需要安装必要的库并导入它们：

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

print("Welcome to how2matplotlib.com")

这段代码导入了NumPy用于数值计算，Matplotlib.pyplot用于绘图，以及SciPy的signal模块用于信号处理。

4. 生成示例信号

为了演示PSD绘图，我们首先需要创建一个示例信号。让我们生成一个包含多个频率成分的信号：

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

# 生成示例信号
np.random.seed(0)
t = np.linspace(0, 1, 1000, endpoint=False)
signal = np.sin(2 * np.pi * 10 * t) + 0.5 * np.sin(2 * np.pi * 20 * t) + 0.2 * np.random.randn(len(t))

plt.figure(figsize=(10, 4))
plt.plot(t, signal)
plt.title("Example Signal - how2matplotlib.com")
plt.xlabel("Time (s)")
plt.ylabel("Amplitude")
plt.show()

Output:

这段代码生成了一个包含10Hz和20Hz正弦波以及随机噪声的信号。我们使用Matplotlib绘制了时域信号，以便直观地了解信号的特征。

5. 使用Matplotlib绘制基本PSD图

现在，让我们使用Matplotlib绘制这个信号的PSD图：

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

# 计算并绘制PSD
f, Pxx = signal.periodogram(signal, fs=1000)

plt.figure(figsize=(10, 4))
plt.semilogy(f, Pxx)
plt.title("Power Spectral Density - how2matplotlib.com")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Power/Frequency (dB/Hz)")
plt.grid(True)
plt.show()

在这个例子中，我们使用SciPy的signal.periodogram函数计算PSD，然后使用Matplotlib的semilogy函数绘制对数刻度的PSD图。这样可以更清楚地显示不同幅度的频率成分。

6. 自定义PSD图的外观

Matplotlib提供了丰富的选项来自定义图表的外观。让我们来改进我们的PSD图：

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

# 自定义PSD图的外观
f, Pxx = signal.periodogram(signal, fs=1000)

plt.figure(figsize=(12, 6))
plt.semilogy(f, Pxx, color='blue', linewidth=2)
plt.title("Customized Power Spectral Density - how2matplotlib.com", fontsize=16)
plt.xlabel("Frequency (Hz)", fontsize=14)
plt.ylabel("Power/Frequency (dB/Hz)", fontsize=14)
plt.grid(True, which="both", ls="-", alpha=0.5)
plt.xlim(0, 50)
plt.ylim(1e-7, 1e0)
plt.tick_params(axis='both', which='major', labelsize=12)
plt.show()

在这个例子中，我们调整了图表大小、线条颜色和宽度、标题和轴标签的字体大小、网格线样式、x轴和y轴的范围，以及刻度标签的大小。这些自定义选项可以使PSD图更加清晰和专业。

7. 添加多个PSD曲线进行比较

在实际应用中，我们可能需要比较多个信号的PSD。Matplotlib可以轻松地在一个图表中绘制多条PSD曲线：

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

# 生成两个不同的信号
t = np.linspace(0, 1, 1000, endpoint=False)
signal1 = np.sin(2 * np.pi * 10 * t) + 0.5 * np.sin(2 * np.pi * 20 * t) + 0.2 * np.random.randn(len(t))
signal2 = np.sin(2 * np.pi * 15 * t) + 0.3 * np.sin(2 * np.pi * 25 * t) + 0.1 * np.random.randn(len(t))

# 计算PSD
f1, Pxx1 = signal.periodogram(signal1, fs=1000)
f2, Pxx2 = signal.periodogram(signal2, fs=1000)

# 绘制比较图
plt.figure(figsize=(12, 6))
plt.semilogy(f1, Pxx1, label='Signal 1', color='blue')
plt.semilogy(f2, Pxx2, label='Signal 2', color='red')
plt.title("Comparison of Power Spectral Densities - how2matplotlib.com", fontsize=16)
plt.xlabel("Frequency (Hz)", fontsize=14)
plt.ylabel("Power/Frequency (dB/Hz)", fontsize=14)
plt.grid(True, which="both", ls="-", alpha=0.5)
plt.xlim(0, 50)
plt.ylim(1e-7, 1e0)
plt.legend(fontsize=12)
plt.tick_params(axis='both', which='major', labelsize=12)
plt.show()

Output:

这个例子展示了如何在同一图表中比较两个不同信号的PSD。我们使用不同的颜色和图例来区分两条曲线，使比较更加直观。

8. 使用Welch方法计算PSD

Welch方法是一种改进的PSD估计技术，它通过将信号分成多个重叠的段，然后平均这些段的周期图来减少噪声的影响。让我们使用Welch方法来计算和绘制PSD：

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

# 使用Welch方法计算PSD
f, Pxx_welch = signal.welch(signal, fs=1000, nperseg=256)

plt.figure(figsize=(12, 6))
plt.semilogy(f, Pxx_welch)
plt.title("PSD using Welch's Method - how2matplotlib.com", fontsize=16)
plt.xlabel("Frequency (Hz)", fontsize=14)
plt.ylabel("Power/Frequency (dB/Hz)", fontsize=14)
plt.grid(True, which="both", ls="-", alpha=0.5)
plt.xlim(0, 50)
plt.ylim(1e-7, 1e0)
plt.tick_params(axis='both', which='major', labelsize=12)
plt.show()

在这个例子中，我们使用signal.welch函数来计算PSD。Welch方法通常会产生更平滑的PSD估计，特别是对于含有噪声的信号。

9. 添加频率标记

在某些情况下，我们可能想要突出显示某些特定的频率。Matplotlib允许我们轻松地添加这些标记：

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

# 计算PSD
f, Pxx = signal.periodogram(signal, fs=1000)

# 绘制PSD并添加频率标记
plt.figure(figsize=(12, 6))
plt.semilogy(f, Pxx)
plt.title("PSD with Frequency Markers - how2matplotlib.com", fontsize=16)
plt.xlabel("Frequency (Hz)", fontsize=14)
plt.ylabel("Power/Frequency (dB/Hz)", fontsize=14)
plt.grid(True, which="both", ls="-", alpha=0.5)
plt.xlim(0, 50)
plt.ylim(1e-7, 1e0)

# 添加频率标记
for freq in [10, 20]:
    plt.axvline(x=freq, color='r', linestyle='--')
    plt.text(freq, 1e-1, f'{freq} Hz', rotation=90, verticalalignment='bottom')

plt.tick_params(axis='both', which='major', labelsize=12)
plt.show()

这个例子展示了如何在PSD图上添加垂直线和文本标签来标记特定的频率。这对于强调信号中的重要频率成分非常有用。

10. 创建子图比较不同的PSD估计方法

有时我们可能想要比较不同的PSD估计方法。Matplotlib的子图功能使这变得简单：

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

# 生成示例信号
np.random.seed(0)
t = np.linspace(0, 1, 1000, endpoint=False)
signal = np.sin(2 * np.pi * 10 * t) + 0.5 * np.sin(2 * np.pi * 20 * t) + 0.2 * np.random.randn(len(t))

# 计算不同的PSD估计
f_periodogram, Pxx_periodogram = signal.periodogram(signal, fs=1000)
f_welch, Pxx_welch = signal.welch(signal, fs=1000, nperseg=256)

# 创建子图
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 10))

# 绘制周期图方法的PSD
ax1.semilogy(f_periodogram, Pxx_periodogram)
ax1.set_title("Periodogram PSD - how2matplotlib.com", fontsize=14)
ax1.set_xlabel("Frequency (Hz)", fontsize=12)
ax1.set_ylabel("Power/Frequency (dB/Hz)", fontsize=12)
ax1.grid(True)
ax1.set_xlim(0, 50)
ax1.set_ylim(1e-7, 1e0)

# 绘制Welch方法的PSD
ax2.semilogy(f_welch, Pxx_welch)
ax2.set_title("Welch's Method PSD - how2matplotlib.com", fontsize=14)
ax2.set_xlabel("Frequency (Hz)", fontsize=12)
ax2.set_ylabel("Power/Frequency (dB/Hz)", fontsize=12)
ax2.grid(True)
ax2.set_xlim(0, 50)
ax2.set_ylim(1e-7, 1e0)

plt.tight_layout()
plt.show()

这个例子创建了两个子图，分别显示使用周期图方法和Welch方法计算的PSD。这种并排比较可以帮助我们更好地理解不同PSD估计方法的特点。

11. 使用颜色映射绘制时频PSD图

对于非平稳信号，我们可能对信号的时频特性感兴趣。Matplotlib可以用来创建时频PSD图（也称为声谱图）：

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

# 生成一个频率随时间变化的信号
t = np.linspace(0, 10, 1000)
w = 2 * np.pi * (3 + t/2)  # 频率随时间线性增加
signal = np.sin(w * t) + 0.1 * np.random.randn(len(t))

# 计算声谱图
f, t, Sxx = signal.spectrogram(signal, fs=100)

# 绘制声谱图
plt.figure(figsize=(12, 6))
plt.pcolormesh(t, f, 10 * np.log10(Sxx), shading='gouraud', cmap='viridis')
plt.title("Time-Frequency PSD (Spectrogram) - how2matplotlib.com", fontsize=16)
plt.xlabel("Time (s)", fontsize=14)
plt.ylabel("Frequency (Hz)", fontsize=14)
plt.colorbar(label='Power/Frequency (dB/Hz)')
plt.ylim(0, 20)
plt.show()

这个例子展示了如何创建一个时频PSD图（声谱图）。我们使用signal.spectrogram函数计算时频PSD，然后使用pcolormesh函数将其可视化。颜色映射显示了不同时间和频率下的功率强度。

12. 3D PSD图

Matplotlib还支持创建3D图表，这可以用来以不同的方式可视化PSD：

from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# 生成示例信号
t = np.linspace(0, 1, 1000, endpoint=False)
signal = np.sin(2 * np.pi * 10 * t) + 0.5* np.sin(2 * np.pi * 20 * t) + 0.2 * np.random.randn(len(t))

# 计算PSD
f, Pxx = signal.periodogram(signal, fs=1000)

# 创建3D图
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')

# 绘制3D PSD图
ax.plot(f, Pxx, zs=0, zdir='y', color='b')
ax.set_title("3D Power Spectral Density - how2matplotlib.com", fontsize=16)
ax.set_xlabel("Frequency (Hz)", fontsize=14)
ax.set_ylabel("Power/Frequency (dB/Hz)", fontsize=14)
ax.set_zlabel("Amplitude", fontsize=14)
ax.set_xlim(0, 50)
ax.set_ylim(0, 0.1)
ax.view_init(elev=20, azim=45)

plt.show()

这个例子展示了如何创建3D PSD图。我们使用Matplotlib的3D绘图功能，将频率和功率绘制在3D空间中。这种可视化方法可以提供PSD的另一种视角。

13. 交互式PSD图

Matplotlib还支持创建交互式图表，这在探索PSD数据时特别有用：

from matplotlib.widgets import Slider
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# 生成示例信号
t = np.linspace(0, 1, 1000, endpoint=False)
signal = np.sin(2 * np.pi * 10 * t) + 0.5 * np.sin(2 * np.pi * 20 * t) + 0.2 * np.random.randn(len(t))

# 计算初始PSD
f, Pxx = signal.periodogram(signal, fs=1000)

# 创建图形和轴
fig, ax = plt.subplots(figsize=(12, 6))
plt.subplots_adjust(bottom=0.25)

# 绘制初始PSD
line, = ax.semilogy(f, Pxx)
ax.set_title("Interactive PSD - how2matplotlib.com", fontsize=16)
ax.set_xlabel("Frequency (Hz)", fontsize=14)
ax.set_ylabel("Power/Frequency (dB/Hz)", fontsize=14)
ax.grid(True)
ax.set_xlim(0, 50)
ax.set_ylim(1e-7, 1e0)

# 创建滑块
axfreq = plt.axes([0.25, 0.1, 0.65, 0.03])
freq_slider = Slider(axfreq, 'Frequency', 1, 50, valinit=10, valstep=1)

# 更新函数
def update(val):
    freq = freq_slider.val
    new_signal = np.sin(2 * np.pi * freq * t) + 0.5 * np.sin(2 * np.pi * 20 * t) + 0.2 * np.random.randn(len(t))
    f, Pxx = signal.periodogram(new_signal, fs=1000)
    line.set_ydata(Pxx)
    fig.canvas.draw_idle()

# 连接滑块到更新函数
freq_slider.on_changed(update)

plt.show()

这个例子创建了一个交互式PSD图，其中包含一个频率滑块。用户可以通过调整滑块来改变信号的主要频率，并实时看到PSD的变化。这种交互式图表对于探索不同参数对PSD的影响非常有用。

14. 使用窗函数改善PSD估计

窗函数可以用来改善PSD估计的质量。让我们比较不同窗函数的效果：

from matplotlib.widgets import Slider
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# 生成示例信号
t = np.linspace(0, 1, 1000, endpoint=False)
signal = np.sin(2 * np.pi * 10 * t) + 0.5 * np.sin(2 * np.pi * 20 * t) + 0.2 * np.random.randn(len(t))

# 计算使用不同窗函数的PSD
f, Pxx_rect = signal.periodogram(signal, fs=1000, window='boxcar')
f, Pxx_hann = signal.periodogram(signal, fs=1000, window='hann')
f, Pxx_blackman = signal.periodogram(signal, fs=1000, window='blackman')

# 绘制比较图
plt.figure(figsize=(12, 6))
plt.semilogy(f, Pxx_rect, label='Rectangular')
plt.semilogy(f, Pxx_hann, label='Hann')
plt.semilogy(f, Pxx_blackman, label='Blackman')
plt.title("PSD with Different Window Functions - how2matplotlib.com", fontsize=16)
plt.xlabel("Frequency (Hz)", fontsize=14)
plt.ylabel("Power/Frequency (dB/Hz)", fontsize=14)
plt.legend()
plt.grid(True)
plt.xlim(0, 50)
plt.ylim(1e-7, 1e0)
plt.show()

这个例子比较了矩形窗（无窗）、汉宁窗和布莱克曼窗对PSD估计的影响。不同的窗函数可以在频率分辨率和谱泄漏之间进行权衡。

15. 绘制累积PSD

累积PSD可以帮助我们理解信号功率如何随频率累积：

from matplotlib.widgets import Slider
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# 生成示例信号
t = np.linspace(0, 1, 1000, endpoint=False)
signal = np.sin(2 * np.pi * 10 * t) + 0.5 * np.sin(2 * np.pi * 20 * t) + 0.2 * np.random.randn(len(t))

# 计算PSD
f, Pxx = signal.periodogram(signal, fs=1000)

# 计算累积PSD
cumulative_psd = np.cumsum(Pxx)
cumulative_psd /= cumulative_psd[-1]  # 归一化

# 绘制累积PSD
plt.figure(figsize=(12, 6))
plt.plot(f, cumulative_psd)
plt.title("Cumulative PSD - how2matplotlib.com", fontsize=16)
plt.xlabel("Frequency (Hz)", fontsize=14)
plt.ylabel("Cumulative Power", fontsize=14)
plt.grid(True)
plt.xlim(0, 50)
plt.ylim(0, 1)
plt.show()

这个例子展示了如何计算和绘制累积PSD。累积PSD图显示了信号功率如何随频率累积，这对于识别包含大部分信号能量的频率范围很有用。

16. 使用对数频率轴

在某些情况下，使用对数频率轴可以更好地显示宽频带信号的PSD：

from matplotlib.widgets import Slider
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# 生成宽频带信号
t = np.linspace(0, 1, 10000, endpoint=False)
signal = np.sin(2 * np.pi * 10 * t) + 0.5 * np.sin(2 * np.pi * 100 * t) + 0.1 * np.sin(2 * np.pi * 1000 * t)

# 计算PSD
f, Pxx = signal.periodogram(signal, fs=10000)

# 绘制使用对数频率轴的PSD
plt.figure(figsize=(12, 6))
plt.semilogx(f, 10 * np.log10(Pxx))
plt.title("PSD with Logarithmic Frequency Axis - how2matplotlib.com", fontsize=16)
plt.xlabel("Frequency (Hz)", fontsize=14)
plt.ylabel("Power/Frequency (dB/Hz)", fontsize=14)
plt.grid(True, which="both", ls="-", alpha=0.5)
plt.xlim(1, 5000)
plt.ylim(-80, 0)
plt.show()

这个例子展示了如何使用对数频率轴绘制PSD。对于包含宽范围频率成分的信号，这种方法可以更清晰地显示所有频率范围的细节。

17. 绘制相位谱

除了功率谱，相位谱也是信号分析中的重要工具：

from matplotlib.widgets import Slider
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# 生成示例信号
t = np.linspace(0, 1, 1000, endpoint=False)
signal = np.sin(2 * np.pi * 10 * t) + 0.5 * np.sin(2 * np.pi * 20 * t + np.pi/4)

# 计算FFT
fft = np.fft.fft(signal)
freqs = np.fft.fftfreq(len(t), t[1] - t[0])

# 计算相位谱
phase_spectrum = np.angle(fft)

# 绘制相位谱
plt.figure(figsize=(12, 6))
plt.plot(freqs[:len(freqs)//2], phase_spectrum[:len(freqs)//2])
plt.title("Phase Spectrum - how2matplotlib.com", fontsize=16)
plt.xlabel("Frequency (Hz)", fontsize=14)
plt.ylabel("Phase (radians)", fontsize=14)
plt.grid(True)
plt.xlim(0, 30)
plt.ylim(-np.pi, np.pi)
plt.show()

Output:

这个例子展示了如何计算和绘制信号的相位谱。相位谱可以提供关于信号各频率成分相位关系的信息，这在某些应用中非常重要。

18. 使用Matplotlib绘制瀑布图

瀑布图是另一种可视化时变信号频谱的方法：

from matplotlib.widgets import Slider
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# 生成时变信号
t = np.linspace(0, 10, 1000)
frequencies = [5, 10, 15]
signal = np.zeros_like(t)
for i, freq in enumerate(frequencies):
    signal += np.sin(2 * np.pi * freq * t) * np.exp(-0.1 * i * t)

# 计算短时傅里叶变换
f, t, Sxx = signal.spectrogram(signal, fs=100, nperseg=100)

# 绘制瀑布图
plt.figure(figsize=(12, 8))
plt.pcolormesh(t, f, 10 * np.log10(Sxx), shading='gouraud', cmap='viridis')
plt.title("Waterfall Plot - how2matplotlib.com", fontsize=16)
plt.xlabel("Time (s)", fontsize=14)
plt.ylabel("Frequency (Hz)", fontsize=14)
plt.colorbar(label='Power/Frequency (dB/Hz)')
plt.ylim(0, 20)

# 调整视角以创建瀑布效果
ax = plt.gca()
ax.view_init(elev=60, azim=-45)
ax.set_proj_type('ortho')

plt.show()

这个例子展示了如何创建一个瀑布图。瀑布图可以直观地显示信号频谱随时间的变化，特别适合分析非平稳信号。

总结

本文详细介绍了如何使用Matplotlib绘制功率谱密度图，并提供了多个实用示例。我们探讨了基本的PSD绘图、自定义图表外观、比较多个PSD曲线、使用不同的PSD估计方法、创建时频PSD图、3D PSD图和交互式PSD图等多个方面。此外，我们还讨论了窗函数的使用、累积PSD、对数频率轴的应用以及相位谱的绘制。

通过这些示例，读者应该能够掌握使用Matplotlib绘制各种类型的PSD图的技能，并能够根据具体需求自定义和优化这些图表。这些技能在信号处理、数据分析和科学研究等多个领域都有广泛的应用。

记住，在实际应用中，选择合适的PSD估计方法和可视化技术取决于具体的信号特性和分析目的。通过实践和经验，你将能够更好地理解和应用这些技术，从而更有效地分析和解释信号数据。