1. 从“看波形”到“看频谱”：为什么我们需要功率谱分析

在信号处理的日常工作中，我们最常打交道的就是信号的时域波形图。一个电压随时间变化的曲线，能告诉我们信号的幅度、频率（如果周期明显）和大概的形态。但很多时候，光看时域波形就像雾里看花。比如，你拿到一段录音，里面混杂着人声、背景音乐和持续的嗡嗡噪声，在时域上它们全部叠加在一起，波形复杂得难以分辨。又或者，在分析机械振动数据时，你想知道设备在哪些频率点上振动能量最强，从而判断是否存在特定部件的故障特征频率。这时，时域分析就显得力不从心了。

我们需要一个工具，能把信号从“时间”这个维度，转换到“频率”维度来观察。这就是频谱分析。而功率谱和功率谱密度，正是频谱分析中两个最核心、最实用的指标。简单来说，它们回答了一个关键问题：信号的能量在不同频率上是如何分布的？

• 功率谱：它告诉我们信号在各个离散频率点上的功率大小。你可以把它想象成一份“能量账单”，清晰地列出了在50Hz、100Hz、150Hz……这些具体频率上，信号分别贡献了多少能量。这对于分析周期性信号、谐波成分特别有用。
• 功率谱密度：当信号中包含连续分布的频率成分（如白噪声）时，用离散点的功率来描述就不够准确了。PSD描述的是单位频率带宽内的功率，单位通常是W/Hz。它反映的是功率的“密度”，更适合分析随机信号或宽带噪声。

那么，如何从一段时域信号计算出它的功率谱或PSD呢？最强大、最通用的工具就是快速傅里叶变换。FFT是离散傅里叶变换的高效算法，它能够将时域信号转换为频域表示。但直接从FFT结果得到功率谱，中间有许多细节和陷阱，比如频谱泄露、栅栏效应、加窗的影响、以及如何从幅度谱换算成功率等。这些细节直接决定了分析结果的准确性和可信度。

接下来，我将结合多年的工程实践，抛开教科书式的理论推导，直接切入如何正确使用FFT来估计这两个关键指标，并分享那些在标准文档里不会写，但能让你少走弯路的实操经验。

2. 理解核心概念：幅度谱、功率谱与功率谱密度到底什么关系？

在动手写代码之前，我们必须理清几个基本但极易混淆的概念。很多初学者直接调用库函数，得到一幅频谱图，却说不清纵坐标到底是什么，这会导致严重的误读。

假设我们有一段离散的时域信号x[n]，长度为N，采样频率为Fs。我们对其做FFT，得到一个复数数组X[k]，其中k=0, 1, ..., N-1对应着从0Hz到Fs（实际上到奈奎斯特频率Fs/2）的频率点。

2.1 从FFT结果到幅度谱

FFT的直接结果X[k]是复数，包含了每个频率成分的幅度和相位信息。我们首先关心幅度：幅度谱[k] = |X[k]|这里的|·|表示取复数的模（绝对值）。但这个幅度谱的物理意义是什么？它对应的是原始信号中该频率成分的振幅。例如，一个纯正弦波A*sin(2πft)，做FFT后在其频率f对应的谱线上，幅度谱值大约为A * (N/2)（对于实数信号，能量会分布在正负频率上，需注意缩放因子）。

2.2 从幅度谱到功率谱（Periodogram）

信号在电阻R上的瞬时功率是v(t)² / R。为简化，通常假设R=1，那么功率就与电压的平方成正比。对于离散信号，一个点的功率近似为x[n]²。整个信号的总能量是Σ x[n]²。

根据帕塞瓦尔定理，信号在时域的总能量等于在频域的总能量。对于FFT，有：Σ |x[n]|² = (1/N) * Σ |X[k]|²因此，|X[k]|² / N可以解释为在第k个频率点上的能量。如果我们关心的是功率（能量除以时间），而信号的总时间为N / Fs，那么第k个频率点上的功率可以表示为：P[k] = (|X[k]|² / N) / (N / Fs) = |X[k]|² / (N * Fs)但更常见的、也是许多教科书和初期谱估计方法（直接法或周期图法）采用的公式是：功率谱[k] = |X[k]|² / N（对于双边谱）功率谱[k] = |X[k]|² / (N/2)（对于单边谱，k=1...N/2-1，且DC和Nyquist成分特殊处理）

这里的关键是缩放因子。不同的软件、库（如MATLAB的periodogram函数，Python的scipy.signal.periodogram）默认的缩放可能不同，但目标都是让频谱图的积分（求和）等于信号的总功率。

实操心得1：归一化因子的混乱这是最大的坑之一。我见过无数团队因为缩放因子不统一，导致对比不同系统或不同参数算出的频谱时，数值差了几个数量级。我的建议是：永远从物理定义出发进行校准。生成一个已知幅度A、频率f的正弦波测试信号，用你的流程计算功率谱，看对应f的谱线值是否等于A²/2（正弦波的平均功率）。如果不是，调整你的缩放因子。建立自己团队内部的“标准流程”比盲目相信默认配置更重要。

2.3 从功率谱到功率谱密度

功率谱P[k]给出了离散频率点k * Fs / N上的功率。它的单位是W（假设信号是电压且电阻为1欧姆）。而功率谱密度描述的是功率的“密度”，单位是W/Hz。

如何从离散的功率谱得到PSD呢？关键在于分辨率带宽。FFT的频率分辨率是Δf = Fs / N。我们可以认为，功率谱P[k]所代表的功率，是均匀分布在以频率f_k为中心、宽度为Δf的这个小小频率区间内的。因此，这个区间内的功率密度（即PSD）为：PSD[k] = P[k] / Δf = P[k] * (N / Fs)

将周期图法的功率谱P_periodogram[k] = |X[k]|² / N代入，得到经典的周期图法PSD估计：PSD_periodogram[k] = (|X[k]|² / N) * (N / Fs) = |X[k]|² / Fs

这个公式非常直观：用FFT结果幅值的平方，除以采样频率，就得到了（双边）功率谱密度的初始估计。同样，转换为单边谱时需要加倍（除DC和Nyquist点外），并注意单位。

注意：这只是最基础的估计方法（周期图法），它的方差性能很差，估计结果波动剧烈，并不是一个“优质”的PSD估计。但它是一切高级方法（如Welch法）的基础。下文我们会深入探讨如何改进它。

3. 周期图法的局限与Welch方法的实战改进

如果我们直接对一段长度为N的信号做FFT，然后套用上述公式计算PSD，这个方法就叫周期图法。它简单直接，但存在两个致命缺点：

1. 高方差：估计结果不稳定。即使你用同一过程生成两段不同的随机信号，它们的周期图也会差异很大。这不利于我们观察信号稳定的频谱特性。
2. 频率分辨率与数据长度的矛盾：高分辨率需要大的N（因为Δf=Fs/N），但大的N会导致计算出的周期图方差更大，曲线更“毛躁”。

工程上，几乎不会直接用原始的周期图法作为最终结果。取而代之的是Welch方法，它是对周期图法的极大改进，也是目前使用最广泛的PSD估计方法。

3.1 Welch方法的核心思想：平均与加窗

Welch方法的精髓在于两点：分段平均和加窗处理。

1. 分段：将长度为N的原始数据分成L段，每段长度为M。允许相邻段之间有部分重叠（通常为50%）。重叠可以减少因分段导致的信息损失。
2. 加窗：对每一段数据，先乘以一个窗函数（如汉宁窗、汉明窗），然后再做FFT。加窗的主要目的是抑制频谱泄露。什么是频谱泄露？简单说，因为FFT默认假设信号是周期性的，且我们截取的长度正好是它的整数个周期。如果不是，截断就会在频域产生额外的、虚假的频率分量，看起来就像能量从主频“泄露”到了旁边。加窗可以使截断处的信号平滑过渡到零，极大减弱这种效应。
3. 计算与平均：对每一段加窗后的数据，计算其修正后的周期图（即PSD估计），然后将这L段的结果平均起来。平均操作能有效降低估计的方差，得到一条更平滑、更稳定的PSD曲线。

3.2 关键参数的选择与权衡

使用Welch方法时，你需要面对几个核心参数的选择，它们共同决定了最终PSD图的质量：

• 窗函数：最常用的是汉宁窗和汉明窗。汉宁窗在抑制旁瓣（减少泄露）方面表现更优，是通用性最好的选择。汉明窗的主瓣稍窄，频率分辨率理论上略好一点，但旁瓣抑制不如汉宁窗。除非有特殊理由，首选汉宁窗。
• 段长度：这直接决定了频率分辨率Δf = Fs / M。M越大，分辨率越高（能区分更近的两个频率），但段数L会变少，平均效果变差，方差可能增大。你需要根据实际需求权衡。如果关注精细的频谱结构，就选大M；如果追求平滑稳定的谱形，可以选小M。
• 重叠率：通常设置为50%。这能在不减少段数L的前提下，使用更长的窗M，从而在保证一定平均次数的同时，获得更好的频率分辨率。重叠超过50%收益递减，且计算量增加。
• FFT点数：通常直接取为段长度M。但也可以通过补零来增加FFT点数NFFT（NFFT >= M）。补零不会增加真实的频率分辨率（因为信息量没变），但可以让频谱图在频率轴上看起来更“连续”，是一种插值效果，有助于更精确地定位谱峰的位置。

3.3 一个完整的Python (scipy.signal.welch) 实操示例

让我们抛开理论，直接看代码。假设我们有一个混合信号：一个50Hz的正弦波，一个120Hz的正弦波，以及一些高斯白噪声。

import numpy as np import matplotlib.pyplot as plt from scipy import signal # 1. 生成模拟信号 Fs = 1000 # 采样频率 1000 Hz T = 2.0 # 信号时长 2秒 N = int(Fs * T) # 总采样点数 2000 t = np.linspace(0, T, N, endpoint=False) # 信号成分：50Hz和120Hz的正弦波，加噪声 x = 1.0 * np.sin(2 * np.pi * 50 * t) + 0.5 * np.sin(2 * np.pi * 120 * t) x += 0.2 * np.random.randn(N) # 加入高斯白噪声 # 2. 使用Welch方法计算PSD # 关键参数设置 nperseg = 256 # 每段长度，决定频率分辨率 Δf = Fs/256 ≈ 3.9 Hz noverlap = 128 # 重叠点数，50%重叠 nfft = 512 # FFT点数，这里选择补零到512点，使频谱线更密 window = 'hann' # 汉宁窗 # 调用welch函数 frequencies, psd = signal.welch(x, Fs, window=window, nperseg=nperseg, noverlap=noverlap, nfft=nfft, scaling='density') # 返回单边PSD # 3. 绘制结果 plt.figure(figsize=(10, 6)) plt.semilogy(frequencies, psd) # 纵坐标用对数坐标，便于观察不同量级的成分 plt.title('Power Spectral Density (Welch‘s Method)') plt.xlabel('Frequency [Hz]') plt.ylabel('Power/Frequency [V**2/Hz]') plt.grid(True, which='both', linestyle='--', alpha=0.6) plt.xlim([0, Fs/2]) # 只显示0到奈奎斯特频率的部分 plt.tight_layout() plt.show() # 4. 验证：找到峰值对应的频率和功率 # 找到PSD中前两个最大峰值的位置 peaks, properties = signal.find_peaks(psd, height=0.01) # 设置一个最小高度阈值 print(f"Detected peak frequencies: {frequencies[peaks]} Hz") print(f"Peak PSD values: {psd[peaks]} V**2/Hz")

在这段代码中：

• scaling='density'告诉函数我们想要计算的是功率谱密度。
• 对数坐标（semilogy）对于观察同时存在强信号和弱噪声的频谱非常有用。
• 通过signal.find_peaks可以自动识别谱峰，这在自动化分析中很实用。

实操心得2：参数选择的“试凑”与准则面对新数据，如何快速确定nperseg（段长）？我的经验是：先根据你关心的最小频率间隔来定。比如，你想分辨间隔10Hz的两个峰，那么分辨率Δf至少要小于10Hz，最好小于5Hz。根据Δf = Fs / nperseg，可以反推出nperseg至少需要Fs / 5。然后，在保证有足够段数（比如8段以上）进行平均以平滑曲线的条件下，可以适当增加nperseg来提高分辨率。这是一个迭代过程，通常需要根据频谱图的“清晰度”和“平滑度”微调几次。

4. 从单边谱到物理单位：工程中的校准与解读

我们计算出了PSD数组，也画出了图，但纵坐标的数值到底意味着什么？如何把它和真实的物理世界联系起来？这是工程应用的最后一步，也是至关重要的一步。

4.1 单边谱与双边谱

对于实数信号，其频谱具有共轭对称性，负频率部分是正频率部分的镜像，不包含新信息。因此，我们通常只显示和关心从0Hz到奈奎斯特频率（Fs/2）的单边谱。

• 双边PSD：计算时使用了从 -Fs/2 到 Fs/2 的所有频率点。总功率等于对所有频率点PSD求和。
• 单边PSD：只显示0到Fs/2的部分。为了保持总功率不变，除了直流（0Hz）和奈奎斯特频率（Fs/2，如果存在）分量外，其他频率点的PSD值需要乘以2。scipy.signal.welch函数在scaling='density'时，默认返回的就是单边PSD。

4.2 单位换算与物理意义

我们的原始信号x通常是从数据采集卡（DAQ）或传感器读出的电压值，单位是伏特（V）。那么计算出的PSD单位就是V²/Hz。这代表了在1Hz带宽内，信号功率的期望值。

如何理解这个值？假设我们在100Hz处读到的PSD值是S_xx(100) = 0.05 V²/Hz。 这意味着，如果我们用一个中心频率为100Hz、带宽为1Hz的理想带通滤波器去过滤这个信号，那么滤波器输出信号的平均功率大约是0.05瓦特（假设负载电阻为1欧姆）。

如果我们的信号代表其他物理量，比如加速度（m/s²）、速度（m/s）、位移（m）呢？这时，PSD的单位就变成了(m/s²)²/Hz等。它直接反映了振动能量在不同频率上的分布密度，是故障诊断、模态分析中判断异常频带的核心依据。

4.3 校准：从ADC读数到真实物理量

很多时候，我们拿到的是数据采集系统的模数转换（ADC）后的数字码值（LSB），而不是直接的电压。这就需要校准。

1. 灵敏度校准：传感器或测量链通常有一个灵敏度系数，单位可能是 mV/g（加速度计）、mV/(m/s)（速度传感器）等。假设灵敏度是Sens = 100 mV/g，采集卡量程是 ±5V，对应16位ADC的 ±32768 LSB。
2. 换算关系：那么，一个数字码值LSB对应的真实物理量（如加速度a）为：a = (LSB / 32768) * (5V) / (0.1 V/g) = LSB * (5 / (32768 * 0.1)) g即存在一个缩放因子K = 5 / (32768 * 0.1)。
3. 应用到PSD：PSD是功率（幅值的平方）的密度。因此，如果时域信号乘以了因子K，那么其PSD需要乘以 K²。PSD_physical(f) = PSD_digital(f) * K²单位也从 (LSB)²/Hz 转换为了 (g)²/Hz。

实操心得3：建立校准流程文档这是团队协作中最容易出错的地方。强烈建议为每一个数据采集系统建立一份《PSD计算校准手册》，明确记录：传感器灵敏度、采集卡量程、ADC位数、软件中设置的工程单位换算关系，以及最终PSD图的纵坐标单位。在分享频谱图时，务必在图的纵坐标轴标签或标题中注明单位，例如 “Acceleration PSD [(m/s²)²/Hz]”。缺少单位的频谱图，其价值大打折扣。

5. 高级话题：窗函数的影响、频谱泄露与参数化方法初探

掌握了Welch方法，你已经能解决90%的工程频谱估计问题。但要成为专家，还需要理解一些更深层次的影响和知道其他工具的存在。

5.1 窗函数选择的深层影响

我们之前推荐了汉宁窗，但窗函数的选择本质上是主瓣宽度与旁瓣衰减之间的权衡。

• 矩形窗（即不加窗）：主瓣最窄，频率分辨率最高。但旁瓣衰减很差（仅-13dB），频谱泄露非常严重。除非你确信信号截断正好是周期的整数倍，否则不要用。
• 汉宁窗：旁瓣衰减好（-31dB），主瓣宽度适中。通用性最佳。
• 汉明窗：类似汉宁窗，但设计上优化了第一个旁瓣的抵消，其旁瓣衰减更均匀，但第一个旁瓣后的衰减不如汉宁窗。主瓣宽度与汉宁窗几乎一样。
• 平顶窗：主瓣非常宽，频率分辨率差。但其巨大优势在于幅度精度高。如果你需要非常精确地测量某个频率成分的幅度（而不是分辨两个很近的频率），比如在校准系统中，平顶窗是首选。它用分辨率换取了更小的幅度估计误差。

如何选择？记住这个口诀：分辨频率用汉宁，测量幅值用平顶，周期整数用矩形（慎用）。

5.2 频谱泄露的直观演示与应对

让我们用代码直观感受泄露。生成一个非整数周期的正弦波。

# 生成一个非整数周期截断的正弦波 F0 = 50.5 # 频率不是Fs/N的整数倍 A = 1.0 x_leak = A * np.sin(2 * np.pi * F0 * t) # 计算不加窗和加汉宁窗的频谱 f, Pxx_rect = signal.welch(x_leak, Fs, window='boxcar', nperseg=256, scaling='spectrum') # ‘boxcar‘即矩形窗 f, Pxx_hann = signal.welch(x_leak, Fs, window='hann', nperseg=256, scaling='spectrum') plt.figure(figsize=(12, 4)) plt.subplot(1,2,1) plt.plot(f, Pxx_rect) plt.title('PSD with Rectangular Window (Leakage Severe)') plt.xlabel('Frequency [Hz]'); plt.ylabel('Power'); plt.grid(True); plt.xlim([40, 70]) plt.subplot(1,2,2) plt.plot(f, Pxx_hann) plt.title('PSD with Hann Window (Leakage Suppressed)') plt.xlabel('Frequency [Hz]'); plt.ylabel('Power'); plt.grid(True); plt.xlim([40, 70]) plt.tight_layout() plt.show()

运行这段代码，你会看到不加窗时，50.5Hz的单频信号能量“泄露”到了周围很多频率上，形成了虚假的谱线。而加汉宁窗后，能量基本集中在主瓣内，旁瓣泄露被极大抑制。这就是加窗的核心价值。

5.3 超越Welch：参数化谱估计简介

Welch方法属于非参数化谱估计，它不假设信号的生成模型，直接对数据操作。但对于短数据记录或信噪比很低的情况，它的分辨率可能不够。

参数化谱估计（如AR模型、Music算法）假设信号是由一个参数模型（如自回归模型）产生的，然后通过估计模型参数来间接得到频谱。它的优势在于，在数据量较少时，有可能获得比Welch方法更高的频率分辨率，尤其适合分析由多个正弦波组成的信号。

例如，使用Yule-Walker方法估计AR模型的PSD：

from scipy.signal import ar_model # 估计AR模型参数（阶数p需要选择，这是一个难点） ar_order = 30 ar_coeffs, noise_variance, _ = ar_model.aryule(x, ar_order) # 根据AR参数计算PSD w, ar_psd = ar_model.arma2psd(ar_coeffs, worN=1024, whole=False) freqs_ar = w * Fs / (2 * np.pi) # 将角频率转换为Hz plt.plot(freqs_ar, ar_psd)

参数化方法的关键和难点在于模型阶数p的选择，选低了分辨率不够，选高了会产生虚假峰。这通常需要基于信息准则（如AIC、BIC）或经验来判断。

实操心得4：何时考虑更高级的方法？我的建议是，永远先从Welch方法开始。它鲁棒、直观、易于理解。只有当Welch方法给出的频谱图分辨率明显不足（比如两个已知的、靠近的频率峰在Welch谱中融合成一个胖峰），且你确信信号符合某种参数模型（如主要由几个正弦波组成）时，才考虑尝试AR模型等参数化方法。对于大多数工程噪声和振动信号，Welch方法配合精心选择的参数，已经完全够用。不要为了“高级”而使用高级方法，引入不必要的复杂性。