Python中偏态正态分布信号FFT的振幅校正问题

核心错误定位 #

你的振幅校正问题根源在于FFT计算和信号重构过程中的两处关键错误：

• 无意义的π乘法：fft_own函数里data_fft = fft(data) * np.pi完全多余，直接导致所有频谱振幅被错误放大π倍，后续振幅计算全部失真。
• 重构时错误的振幅缩放：手动重构时amplitude = amp[i] * df是错误操作——你已经通过2/N * np.abs(data_fft)算出了单频率分量的振幅，无需再乘频率间隔df，这会再次缩小振幅，导致重构信号整体偏小。

修正后的完整代码 #

from numpy.fft import fft, fftfreq
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import skewnorm
np.random.seed(1234)


def data():
    height = 1
    # a=0时为标准正态分布PDF，loc设置为时间中点
    data = height * skewnorm.pdf(t, a=0, loc=t[int(N/2)])
    return data


def fft_own(data):
    freq = fftfreq(N, dt)
    data_fft = fft(data)  # 移除错误的*np.pi
    # 直流分量振幅：abs(data_fft[0])/N；正/负频率振幅：2*abs(data_fft[i])/N
    amp = np.abs(data_fft) / N
    amp[1:-1] *= 2  # 对非直流的所有频率应用2倍缩放
    phase = np.angle(data_fft)
    # 保留所有频率分量用于重构
    peaks = np.arange(len(freq))
    return freq, amp, phase, peaks


def rebuild(fft_own):
    freq, amp, phase, peaks = fft_own
    data_rebuild = 0
    for i in peaks:
        # 直接使用计算好的amp[i]作为振幅，无需额外缩放
        data_rebuild += amp[i] * np.exp(1j * (2*np.pi * freq[i] * t + phase[i]))

    f, ax = plt.subplots(1, 1)
    ax.plot(t, data_init, label="initial signal")
    ax.plot(t, np.real(data_rebuild), label="rebuild")
    ax.plot(t, data_init - np.real(data_rebuild), label="diff", alpha=0.5)
    ax.set_xlim(0, t1-1)
    ax.legend()
    plt.show()


t0 = 0
t1 = 10
N = 4096
t = np.linspace(t0, t1, int(N))
dt = (t1 - t0) / N

data_init = data()
fft_init = fft_own(data_init)
rebuild_init = rebuild(fft_init)

高斯信号FFT处理的优化建议 #

你提到的“直接用裸FFT峰值”的建议并不适合高斯信号：

• 高斯信号的频谱本身是连续的高斯分布，没有尖锐单频峰值，FFT仅能得到离散采样点，取峰值会丢失大量频谱信息，导致振幅计算误差。
• 更可靠的方法是使用帕塞瓦尔定理做积分计算：信号能量等于频谱能量，即sum(data**2)*dt = sum(np.abs(data_fft)**2)/(N*dt)，通过能量反推等效振幅，对截断高斯信号的稳定性远高于取峰值。
• 周期信号和矩形脉冲的不稳定问题，窗函数是正确选择——这类信号时域不连续会引发严重频谱泄漏，窗函数可有效抑制；高斯信号时域平滑，泄漏本身较小，若需更高精度，可尝试汉宁窗，但效果提升有限。

效果说明 #

修正后，重构信号与原始信号的差值会趋近于机器精度级别的噪声，不会再出现高斯形态的差值信号。调整skewnorm的a参数（偏度）时，重构精度依然能保持一致。

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