如何用SciPy FFT生成平滑精准的功率谱？

问题分析 #

当前FFT实部谱图出现锯齿状波动与峰前小谷，核心诱因有三点：

• 栅栏效应：FFT频率分辨率由总采样时间决定（$\Delta f=1/T$），有限分辨率导致真实频率落在FFT离散频率点之间，催生虚假波动；
• 信号截断泄漏：未加窗的信号截断会引入频谱泄漏，放大高频杂散波动；
• 实部特性：FFT实部包含正负相位分量，本身比绝对值更易出现波动，直接绘制会放大该现象。

解决方案步骤 #

1. 信号预处理：加窗抑制泄漏 #

在FFT前对信号施加汉宁窗，减少截断带来的频谱泄漏，这是抑制虚假波动最直接的手段：

window = np.hanning(len(x))
x_windowed = x * window
fft_vals = fft(x_windowed)

2. 提升真实频率分辨率 #

放弃线性插值（会引入虚假特征），通过增加总采样时间提升分辨率：

• 当前总采样时间stoptime_main=204.8，频率分辨率约为0.00488Hz；若提升至409.6，分辨率翻倍至0.00244Hz，可显著缓解栅栏效应。
• 若无法增加采样时间，可在加窗后补零（仅用于绘图平滑，不提升真实分辨率）：

# 补零到下一个2的幂次长度，提升FFT点数
n_padded = 2**np.ceil(np.log2(len(x))).astype(int)
fft_vals = fft(x_windowed, n_padded)
frequencies = fftfreq(n_padded, d=timestep)[:n_padded//2]

3. 实部谱图的平滑处理 #

针对实部波动，使用Savitzky-Golay滤波器（在保留峰位的同时平滑噪声）：

from scipy.signal import savgol_filter
# 窗口长度需为奇数，多项式阶数按需调整
smoothed_real = savgol_filter(np.real(fft_vals[:n//2][mask]), window_length=11, polyorder=2)

4. 批量处理优化 #

针对数千条序列，优化代码提升效率：

• 预计算窗函数，避免重复生成；
• 向量化FFT操作，利用NumPy广播特性；
• 预计算频率轴，无需每条序列重复计算。

修改后的完整代码示例 #

import scipy
from scipy.integrate import solve_ivp
from scipy.fft import fft, fftfreq
from scipy.signal import savgol_filter
import numpy as np
import matplotlib.pyplot as plt

def dtr(degrees):
    return (degrees/360)*(2*np.pi)

# 提升总采样时间以优化频率分辨率
stoptime_main = 409.6
interval_main = 0.01

def eqs(t, w):
    z, y, m, n = w
    f = [
         (m - (n * np.cos(z - y))) / (1 + (np.sin(z - y)) ** 2), 
         ((2*n) - (m * np.cos(z - y))) / (1 + (np.sin(z - y)) ** 2), 
         (-2 * np.sin(z)) - ((m * n * np.sin(z - y)) / (1 + (np.sin(z - y)) ** 2)) + ((((m**2) + (2*(n**2)) - (2*m*n*np.cos(z - y))) * np.sin(2*(z - y))) / (2 * ((1 + (np.sin(z - y)) ** 2) ** 2))),
         (-1 * np.sin(y)) + ((m * n * np.sin(z - y)) / (1 + (np.sin(z - y)) ** 2)) - ((((m**2) + (2*(n**2)) - (2*m*n*np.cos(z - y))) * np.sin(2*(z - y))) / (2 * ((1 + (np.sin(z - y)) ** 2) ** 2)))
    ]
    return f

def rungekutta8(theta, phi):
    z1 = theta
    y1 = phi
    m1 = 0
    n1 = 0
    abserr = 1.0e-8
    relerr = 1.0e-6
    stoptime = stoptime_main
    interval = interval_main
    t = np.arange(0.0, stoptime, interval)
    t_span = (0.0, stoptime)
    w0 = [z1, y1, m1, n1]
    wsol = solve_ivp(eqs, t_span, w0, t_eval=t, method='DOP853')
    return wsol, t 

wsol, t = rungekutta8(dtr(32), dtr(0))
z_list = wsol.y[0, :]

def decompose_signal(x, time):
    n = len(x)
    timestep = time[1] - time[0]
    
    # 加窗处理抑制泄漏
    window = np.hanning(n)
    x_windowed = x * window
    
    # 补零到2的幂次提升FFT点数
    n_padded = 2**np.ceil(np.log2(n)).astype(int)
    fft_vals = fft(x_windowed, n_padded)
    
    # 计算并过滤频率轴
    frequencies = fftfreq(n_padded, d=timestep)[:n_padded//2]
    mask = (frequencies >= 0) & (frequencies <= 3)
    frequencies = frequencies[mask]
    fft_real = np.real(fft_vals[:n_padded//2][mask])
    
    # 平滑实部谱图
    smoothed_real = savgol_filter(fft_real, window_length=15, polyorder=2)
    
    return frequencies, smoothed_real

frequencies, smoothed_real = decompose_signal(z_list, t)
plt.plot(frequencies, smoothed_real)
plt.xlim([0, 0.8])
plt.xlabel('频率 (Hz)')
plt.ylabel('FFT实部')
plt.show()

关键注意事项 #

• 加窗会略微降低谱峰幅度，若需保留幅度信息，可对窗函数进行能量归一化；
• 平滑参数（如Savitzky-Golay窗口长度）需根据实际频谱峰宽调整，避免过度平滑丢失细节；
• 批量处理时统一采样时间与FFT参数，确保所有序列谱图具有可比性。

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