如何在Python中绘制功率谱图？

Question

我在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)

绘制时，它看起来像图中左图：

我需要以某种方式将图表放在右侧。如果我没猜错的话，这可能是一个功率谱图。它与快速傅里叶变换有关，但我不太确定如何绘制它，因为网络上的所有内容都会给出不同的结果，不幸的是也与所需的图表不同。

我尝试过 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))
)

但是，我仍然无法实现这一目标：

Answer 1

标准频谱图看起来很像您的上一张图，尽管您的函数在通带中产生的值远多于 10**-4：

import matplotlib.pyplot as plt
import numpy as np


def a(t: np.ndarray) -> np.ndarray:
    return (
        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)
    )


fsample = 1e5
tsample = 1/fsample
t = np.arange(0, 1, tsample)
a_t = a(t)
a_f = np.abs(np.fft.rfft(a_t))
f = np.fft.rfftfreq(n=t.size, d=tsample)

fig, ax = plt.subplots()
ax.loglog(f, a_f)
plt.show()