功率谱密度图示例

在Matplotlib中绘制功率谱密度(PSD)。

PSD是信号处理领域中常见的图形。NumPy有许多用于计算PSD的有用库。下面，我们演示一些如何使用Matplotlib实现和可视化这一点的示例。

  1. import matplotlib.pyplot as plt
  2. import numpy as np
  3. import matplotlib.mlab as mlab
  4. import matplotlib.gridspec as gridspec
  6. # Fixing random state for reproducibility
  7. np.random.seed(19680801)
  9. dt = 0.01
  10. t = np.arange(0, 10, dt)
  11. nse = np.random.randn(len(t))
  12. r = np.exp(-t / 0.05)
  14. cnse = np.convolve(nse, r) * dt
  15. cnse = cnse[:len(t)]
  16. s = 0.1 * np.sin(2 * np.pi * t) + cnse
  18. plt.subplot(211)
  19. plt.plot(t, s)
  20. plt.subplot(212)
  21. plt.psd(s, 512, 1 / dt)
  23. plt.show()

功率谱密度图示例

将其与等效的Matlab代码进行比较，以完成相同的任务：

  1. dt = 0.01;
  2. t = [0:dt:10];
  3. nse = randn(size(t));
  4. r = exp(-t/0.05);
  5. cnse = conv(nse, r)*dt;
  6. cnse = cnse(1:length(t));
  7. s = 0.1*sin(2*pi*t) + cnse;
  9. subplot(211)
  10. plot(t,s)
  11. subplot(212)
  12. psd(s, 512, 1/dt)

下面，我们将展示一个稍微复杂一些的示例，演示填充如何影响产生的PSD。

  1. dt = np.pi / 100.
  2. fs = 1. / dt
  3. t = np.arange(0, 8, dt)
  4. y = 10. * np.sin(2 * np.pi * 4 * t) + 5. * np.sin(2 * np.pi * 4.25 * t)
  5. y = y + np.random.randn(*t.shape)
  7. # Plot the raw time series
  8. fig = plt.figure(constrained_layout=True)
  9. gs = gridspec.GridSpec(2, 3, figure=fig)
  10. ax = fig.add_subplot(gs[0, :])
  11. ax.plot(t, y)
  12. ax.set_xlabel('time [s]')
  13. ax.set_ylabel('signal')
  15. # Plot the PSD with different amounts of zero padding. This uses the entire
  16. # time series at once
  17. ax2 = fig.add_subplot(gs[1, 0])
  18. ax2.psd(y, NFFT=len(t), pad_to=len(t), Fs=fs)
  19. ax2.psd(y, NFFT=len(t), pad_to=len(t) * 2, Fs=fs)
  20. ax2.psd(y, NFFT=len(t), pad_to=len(t) * 4, Fs=fs)
  21. plt.title('zero padding')
  23. # Plot the PSD with different block sizes, Zero pad to the length of the
  24. # original data sequence.
  25. ax3 = fig.add_subplot(gs[1, 1], sharex=ax2, sharey=ax2)
  26. ax3.psd(y, NFFT=len(t), pad_to=len(t), Fs=fs)
  27. ax3.psd(y, NFFT=len(t) // 2, pad_to=len(t), Fs=fs)
  28. ax3.psd(y, NFFT=len(t) // 4, pad_to=len(t), Fs=fs)
  29. ax3.set_ylabel('')
  30. plt.title('block size')
  32. # Plot the PSD with different amounts of overlap between blocks
  33. ax4 = fig.add_subplot(gs[1, 2], sharex=ax2, sharey=ax2)
  34. ax4.psd(y, NFFT=len(t) // 2, pad_to=len(t), noverlap=0, Fs=fs)
  35. ax4.psd(y, NFFT=len(t) // 2, pad_to=len(t),
  36.         noverlap=int(0.05 * len(t) / 2.), Fs=fs)
  37. ax4.psd(y, NFFT=len(t) // 2, pad_to=len(t),
  38.         noverlap=int(0.2 * len(t) / 2.), Fs=fs)
  39. ax4.set_ylabel('')
  40. plt.title('overlap')
  42. plt.show()

功率谱密度图示例2

这是一个来自信号处理工具箱的MATLAB示例的移植版本，它显示了Matplotlib和MATLAB对PSD的缩放之间的一些差异。

  1. fs = 1000
  2. t = np.linspace(0, 0.3, 301)
  3. A = np.array([2, 8]).reshape(-1, 1)
  4. f = np.array([150, 140]).reshape(-1, 1)
  5. xn = (A * np.sin(2 * np.pi * f * t)).sum(axis=0)
  6. xn += 5 * np.random.randn(*t.shape)
  8. fig, (ax0, ax1) = plt.subplots(ncols=2, constrained_layout=True)
  10. yticks = np.arange(-50, 30, 10)
  11. yrange = (yticks[0], yticks[-1])
  12. xticks = np.arange(0, 550, 100)
  14. ax0.psd(xn, NFFT=301, Fs=fs, window=mlab.window_none, pad_to=1024,
  15.         scale_by_freq=True)
  16. ax0.set_title('Periodogram')
  17. ax0.set_yticks(yticks)
  18. ax0.set_xticks(xticks)
  19. ax0.grid(True)
  20. ax0.set_ylim(yrange)
  22. ax1.psd(xn, NFFT=150, Fs=fs, window=mlab.window_none, pad_to=512, noverlap=75,
  23.         scale_by_freq=True)
  24. ax1.set_title('Welch')
  25. ax1.set_xticks(xticks)
  26. ax1.set_yticks(yticks)
  27. ax1.set_ylabel('')  # overwrite the y-label added by `psd`
  28. ax1.grid(True)
  29. ax1.set_ylim(yrange)
  31. plt.show()

功率谱密度图示例3

这是一个来自信号处理工具箱的MATLAB示例的移植版本，它显示了Matplotlib和MATLAB对PSD的缩放之间的一些差异。

它使用了一个复杂的信号，所以我们可以看到，复杂的PSD的工作正常。

  1. prng = np.random.RandomState(19680801)  # to ensure reproducibility
  3. fs = 1000
  4. t = np.linspace(0, 0.3, 301)
  5. A = np.array([2, 8]).reshape(-1, 1)
  6. f = np.array([150, 140]).reshape(-1, 1)
  7. xn = (A * np.exp(2j * np.pi * f * t)).sum(axis=0) + 5 * prng.randn(*t.shape)
  9. fig, (ax0, ax1) = plt.subplots(ncols=2, constrained_layout=True)
  11. yticks = np.arange(-50, 30, 10)
  12. yrange = (yticks[0], yticks[-1])
  13. xticks = np.arange(-500, 550, 200)
  15. ax0.psd(xn, NFFT=301, Fs=fs, window=mlab.window_none, pad_to=1024,
  16.         scale_by_freq=True)
  17. ax0.set_title('Periodogram')
  18. ax0.set_yticks(yticks)
  19. ax0.set_xticks(xticks)
  20. ax0.grid(True)
  21. ax0.set_ylim(yrange)
  23. ax1.psd(xn, NFFT=150, Fs=fs, window=mlab.window_none, pad_to=512, noverlap=75,
  24.         scale_by_freq=True)
  25. ax1.set_title('Welch')
  26. ax1.set_xticks(xticks)
  27. ax1.set_yticks(yticks)
  28. ax1.set_ylabel('')  # overwrite the y-label added by `psd`
  29. ax1.grid(True)
  30. ax1.set_ylim(yrange)
  32. plt.show()

功率谱密度图示例4

下载这个示例