如何在Python中绘制功率谱图？附信号函数及尝试记录

我定义了如下时域信号：

omega = (
    5.4547990171221e-20 * np.sin(4.35893728452007e-5 * t) ** 2
    - 6.80677962622999e-40 * np.sin(4.35893728452007e-5 * t)
                           * np.cos(4.35893728452007e-5 * t)
    + 3.84259827520928e-19 * np.sin(2831.25294891184 * t) ** 2
    - 2.77555756156289e-17 * np.sin(2831.25294891184 * t)
                           * np.cos(2831.25294891184 * t)
    + 5.45479901712209e-20 * np.cos(4.35893728452007e-5 * t) ** 2
    - 3.87725551857472e-22 * np.cos(2831.25294891184 * t) ** 2
    + 1391.74933087011 * np.exp(-40.0880015429485 * t)
                       * np.sin(2831.25294891184 * t)
    + 9.5003864918335e-13 * np.exp(-40.0880015429485 * t)
                          * np.cos(2831.25294891184 * t)
    - 0.000158807497479067 * np.exp(-9.500167e-15 * t)
                           * np.sin(4.35893728452007e-5 * t)
    + 2094.39510239319 * np.exp(-9.500167e-15 * t)
    * np.cos(4.35893728452007e-5 * t)
)

绘制时域波形后，我需要生成对应的功率谱图，推测需通过FFT实现，但尝试了numpy的fft及scipy的welch方法均未得到预期图形，尝试的代码如下：

freqs, psd = sig.welch(signal, 10000)

fig, ax = plt.subplots(figsize=(8, 6))
ax.semilogx(
    freqs, np.log10(np.abs(psd))
)

问题分析与解决思路 #

1. 信号分量的核心特征

你的信号包含两个差异极大的频段：

• 低频分量：角频率4.3589e-5 rad/s，对应频率≈6.937e-6 Hz（周期约144150秒/40小时），幅值约2094（主导项）
• 高频分量：角频率2831.25 rad/s，对应频率≈450.6 Hz，幅值初始上千但随指数衰减（衰减系数40，10秒后幅值几乎为0）
• 其余小量级项（e-17~e-20）幅值可忽略，对功率谱无显著贡献

2. 现有方法失效的原因

• 采样时长不足：低频分量周期极长，若你的t仅取了几秒/几十秒，FFT/Welch无法捕捉到这个低频信号的周期特征
• 幅值差异过大：高低频分量幅值差达6个数量级以上，直接合并绘制会导致低频小项（或高频衰减后的信号）被完全淹没
• Welch参数不合理：默认窗口大小/重叠率不适合这种跨多个数量级频段的信号

3. 修正方案

将高低频分量分开处理，针对各自特征调整采样与分析参数：

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

# ---------------------- 处理高频分量 ----------------------
fs_high = 10000  # 高频采样率
t_high = np.arange(0, 10, 1/fs_high)  # 取前10秒，覆盖高频衰减过程
omega_high = (
    3.84259827520928e-19 * np.sin(2831.25294891184 * t_high) ** 2
    - 2.77555756156289e-17 * np.sin(2831.25294891184 * t_high) * np.cos(2831.25294891184 * t_high)
    - 3.87725551857472e-22 * np.cos(2831.25294891184 * t_high) ** 2
    + 1391.74933087011 * np.exp(-40.0880015429485 * t_high) * np.sin(2831.25294891184 * t_high)
    + 9.5003864918335e-13 * np.exp(-40.0880015429485 * t_high) * np.cos(2831.25294891184 * t_high)
)

# 用Welch计算高频功率谱，调整窗口大小提升分辨率
freqs_high, psd_high = sig.welch(omega_high, fs_high, nperseg=4096, noverlap=2048)

# ---------------------- 处理低频分量 ----------------------
fs_low = 1e-3  # 降采样到0.001Hz，每1000秒采1个点
t_low = np.arange(0, 1.5e6, 1/fs_low)  # 取约17天，覆盖10个低频周期
omega_low = (
    5.4547990171221e-20 * np.sin(4.35893728452007e-5 * t_low) ** 2
    - 6.80677962622999e-40 * np.sin(4.35893728452007e-5 * t_low) * np.cos(4.35893728452007e-5 * t_low)
    + 5.45479901712209e-20 * np.cos(4.35893728452007e-5 * t_low) ** 2
    - 0.000158807497479067 * np.exp(-9.500167e-15 * t_low) * np.sin(4.35893728452007e-5 * t_low)
    + 2094.39510239319 * np.exp(-9.500167e-15 * t_low) * np.cos(4.35893728452007e-5 * t_low)
)

# 用Welch计算低频功率谱，窗口取半长保证频率分辨率
freqs_low, psd_low = sig.welch(omega_low, fs_low, nperseg=len(omega_low)//2, noverlap=0)

# ---------------------- 分开可视化 ----------------------
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8))

# 高频功率谱（用10log10标准dB表示）
ax1.semilogx(freqs_high, 10*np.log10(psd_high))
ax1.set_title('高频分量功率谱')
ax1.set_xlabel('频率 (Hz)')
ax1.set_ylabel('功率谱密度 (dB/Hz)')
ax1.grid(True)

# 低频功率谱
ax2.semilogx(freqs_low, 10*np.log10(psd_low))
ax2.set_title('低频分量功率谱')
ax2.set_xlabel('频率 (Hz)')
ax2.set_ylabel('功率谱密度 (dB/Hz)')
ax2.grid(True)

plt.tight_layout()
plt.show()

关键调整点

• 高频部分取前10秒数据，覆盖衰减过程，用较大窗口提升频率分辨率
• 低频部分降采样并延长采样时长，确保能捕捉到完整周期
• 用10*np.log10(psd)替代np.log10(np.abs(psd))，符合功率谱的标准dB可视化规范
• 分开绘制高低频，避免幅值差异导致的信息丢失

内容的提问来源于stack exchange，提问作者ms4rkl