Python使用matplotlib绘制了包含两个信号的时域图和它们的交叉功率谱密度CSD

使用 matplotlib.pyplot 绘制了包含两个信号的时域图和它们的交叉功率谱密度（Cross Spectral Density，CSD）图。

1. 导入库和设置随机种子：

   import matplotlib.pyplot as plt
   import numpy as np
   
   fig, (ax1, ax2) = plt.subplots(2, 1, layout='constrained')
   dt = 0.01
   t = np.arange(0, 30, dt)
   np.random.seed(19680801)
   

   这里导入了 matplotlib.pyplot 和 numpy 库，并创建了包含两个子图的画布。

2. 生成白噪声和带颜色的噪声：

   nse1 = np.random.randn(len(t))  # 白噪声 1
   nse2 = np.random.randn(len(t))  # 白噪声 2
   r = np.exp(-t / 0.05)
   cnse1 = np.convolve(nse1, r, mode='same') * dt   # 带颜色的噪声 1
   cnse2 = np.convolve(nse2, r, mode='same') * dt   # 带颜色的噪声 2
   

   生成两个白噪声信号 nse1 和 nse2，然后通过卷积操作得到带颜色的噪声信号 cnse1 和 cnse2。

3. 生成两个信号：

   s1 = 0.01 * np.sin(2 * np.pi * 10 * t) + cnse1
   s2 = 0.01 * np.sin(2 * np.pi * 10 * t) + cnse2
   

   生成两个信号 s1 和 s2，其中包含相干部分和带颜色的噪声。

4. 绘制时域图：

   ax1.plot(t, s1, t, s2)
   ax1.set_xlim(0, 5)
   ax1.set_xlabel('Time (s)')
   ax1.set_ylabel('s1 and s2')
   ax1.grid(True)
   

   在第一个子图中绘制了两个信号 s1 和 s2 的时域图，并设置了横轴范围、标签和网格。

5. 计算并绘制交叉功率谱密度（CSD）图：

   cxy, f = ax2.csd(s1, s2, 256, 1. / dt)
   ax2.set_ylabel('CSD (dB)')
   

   使用 csd 函数计算并绘制 s1 和 s2 的交叉功率谱密度图，其中 256 是 FFT 窗口大小，1. / dt 是采样频率。

6. 显示图形：

   plt.show()
   

   • 完整代码

import matplotlib.pyplot as plt
import numpy as np

# 创建包含两个子图的画布
fig, (ax1, ax2) = plt.subplots(2, 1, layout='constrained')

# 定义时间步长和时间向量
dt = 0.01
t = np.arange(0, 30, dt)

# 设置随机种子以确保可重现性
np.random.seed(19680801)

# 生成白噪声
nse1 = np.random.randn(len(t))  # 白噪声 1
nse2 = np.random.randn(len(t))  # 白噪声 2

# 定义指数衰减函数 r
r = np.exp(-t / 0.05)

# 使用卷积得到带颜色的噪声
cnse1 = np.convolve(nse1, r, mode='same') * dt   # 带颜色的噪声 1
cnse2 = np.convolve(nse2, r, mode='same') * dt   # 带颜色的噪声 2

# 生成两个信号，包含相干部分和随机部分
s1 = 0.01 * np.sin(2 * np.pi * 10 * t) + cnse1
s2 = 0.01 * np.sin(2 * np.pi * 10 * t) + cnse2

# 绘制时域图
ax1.plot(t, s1, t, s2)
ax1.set_xlim(0, 5)
ax1.set_xlabel('Time (s)')
ax1.set_ylabel('s1 and s2')
ax1.grid(True)

# 计算并绘制交叉功率谱密度（CSD）图
cxy, f = ax2.csd(s1, s2, 256, 1. / dt)
ax2.set_ylabel('CSD (dB)')

# 显示图形
plt.show()