Python中FFT的PSD缩放异常问题（汉明窗与矩形窗结果不符）

测试过程 #

我分别对比矩形窗（无窗）和汉明窗的PSD计算，测试步骤如下：

1. 生成测试信号 #

import numpy as np
from numpy.fft import fft, fftfreq
import matplotlib.pyplot as plt

f_aq = 2**18  # 采样率
Np = 2**22    # 数据点数

t = np.arange(0,Np) / f_aq

f1 = 70       # 信号频率（Hz）
A1 = 6
y = A1 * np.sin(2*np.pi*f1 * t)

2. 应用汉明窗并执行FFT #

window = np.hamming(Np)
y_window = y * window

Y = fft(y, n=Np, norm='backward')
Y_window = fft(y_window, n=Np, norm='backward')

freq = fftfreq(Np, 1/f_aq)

3. PSD计算与异常现象 #

• 单边FFT幅值：用1/np.mean(window)缩放汉明窗结果后，与矩形窗幅值匹配。
• 功率谱：通过窗函数求和S1缩放后，两者结果一致。
• PSD异常：使用ENBW（等效噪声带宽）缩放后，矩形窗PSD与手动计算的水平线一致，但汉明窗PSD幅值偏低，且scipy.welch验证结果相同。

PSD计算代码：

# 提取单边FFT结果（假设已完成双边转单边处理）
freq_one_sided = freq[:Np//2]
Y_one_sided = Y[:Np//2]
Y_one_sided_window = Y_window[:Np//2]
# 双边转单边的幅值加倍（直流和Nyquist分量除外）
Y_one_sided[1:-1] *= 2
Y_one_sided_window[1:-1] *= 2

# 计算功率谱
Y_one_sided_fft_power_spectrum = np.abs(Y_one_sided)**2 / (Np**2)
Y_one_sided_fft_power_spectrum_window = np.abs(Y_one_sided_window)**2 / (Np**2)

# 功率谱密度（PSD）
Y_one_sided_fft_psd = Y_one_sided_fft_power_spectrum / (f_aq/Np)

# 窗参数与ENBW计算
S1 = np.sum(window)
S2 = np.sum(window**2)
ENBW = f_aq * S2 / S1**2

Y_one_sided_fft_psd_window = Y_one_sided_fft_power_spectrum_window / ENBW

# 绘图
plt.plot(freq_one_sided, Y_one_sided_fft_psd)
plt.plot(freq_one_sided+2, Y_one_sided_fft_psd_window)

plt.title("one-sided PSD")
plt.xlim(-1, 100)
plt.xlabel("Frequency in Hz")
plt.ylabel("PSD")

# 手动计算PSD参考线
energy = np.dot(y,y) / len(y)
psd = energy / (f_aq/Np)
plt.axhline(psd, color='k', alpha=0.3)
plt.show()

问题根源与解决方法 #

你的核心问题是混淆了加窗信号的PSD和原始信号的PSD，以及缩放逻辑的遗漏：

1. 矩形窗的特殊性质 #

矩形窗的S1=Np、S2=Np，计算得ENBW = f_aq * Np / Np² = f_aq/Np，刚好和你用的频率分辨率(f_aq/Np)一致，所以矩形窗的PSD计算自动匹配原始信号的参考值。

2. 汉明窗的缩放遗漏 #

加窗会导致信号能量损失，你需要先补偿这部分损失，再计算PSD：

• 加窗后信号的FFT幅值是原始信号的S1/Np倍（S1为窗函数求和），因此功率谱需要乘以(Np/S1)²来还原原始信号的功率水平。
• 之后再除以频率分辨率(f_aq/Np)得到PSD，或者等价地，用补偿后的功率谱除以ENBW。

修正后的PSD计算代码 #

# 汉明窗功率谱补偿：还原原始信号功率水平
power_spectrum_hamming_corrected = Y_one_sided_fft_power_spectrum_window * (Np/S1)**2
# 计算修正后的PSD
Y_one_sided_fft_psd_window_corrected = power_spectrum_hamming_corrected / (f_aq/Np)

# 重新绘图对比
plt.plot(freq_one_sided, Y_one_sided_fft_psd, label='矩形窗')
plt.plot(freq_one_sided+2, Y_one_sided_fft_psd_window_corrected, label='汉明窗（修正后）')
plt.axhline(psd, color='k', alpha=0.3, label='手动计算参考线')
plt.title("修正后的one-sided PSD")
plt.xlim(-1, 100)
plt.xlabel("Frequency in Hz")
plt.ylabel("PSD")
plt.legend()
plt.show()

关于scipy.welch的结果 #

scipy.welch默认输出的是加窗信号的PSD，而非原始信号的PSD，所以它的结果会和你未修正的汉明窗PSD一致。如果要让welch输出匹配原始信号的PSD，需要手动乘以补偿因子(Np/S1)²。

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