电力系统同步相量计算：FFT、窗函数、HHT与小波变换的Matlab仿真对比

做电力系统同步相量计算算法选型时，很多人第一反应就是把FFT搬过来，跑通一次仿真就认为完事了。我在实际对比FFT、窗函数法、希尔伯特-黄变换、小波变换这四类方法时发现，真正决定结果质量的往往不是算法本身的数学公式，而是你如何理解同步相量的定义、如何构造测试信号、又如何评估误差。这篇内容就是从这些角度出发，把我做Matlab仿真时的完整思路、关键代码片段和踩过的坑一次性讲清楚，适合正在做PMU算法研究、课程设计或者准备毕业课题的同学参考。

1. 同步相量的定义与四条技术路线的分工

1.1 同步相量到底在算什么

同步相量（Synchronized Phasor）的核心思想，是把电力系统某一节点上的交流电压或电流，表示成相对于全球统一时间基准的复数形式。可以简单理解成一个旋转向量在某个参考时刻的“快照”：幅值对应电压或电流的RMS值，相位对应它相对于参考余弦波的偏移角。电力系统广域测量系统（WAMS）和同步相量测量单元（PMU）的基础就是这一步计算，IEEE C37.118标准里用总向量误差（TVE）来衡量算法输出的相量与真实相量之间的偏差。

在Matlab里做同步相量计算，通常输入是一段离散采样序列x(n)，采样频率固定为Fs。你要输出的是基波分量的幅值A、相位φ和频率f。如果信号是理想的正弦波，做一次FFT后找到最大谱线就能算出结果。但真实电网里，频率会在50Hz附近波动，还有谐波、噪声、间谐波和暂态分量，这导致直接FFT的结果并不可靠。于是窗函数法、希尔伯特-黄变换、小波变换这些方法才被引入进来，各自解决不同场景下的测量问题。

1.2 为什么FFT在电力系统动态条件下会“失灵”

FFT本身是稳态频谱分析工具，它假设被分析的信号在窗内是周期信号，并且频率严格等于谱线频率。当实际系统频率偏离50Hz时，信号周期与采样窗口不再对齐，频谱就会发生泄漏：原来只应该集中在50Hz附近的能量扩散到相邻谱线上去，幅值测量变低，相位也会出现偏移。

我用一个简单例子验证过：采样率Fs=12800Hz，取0.2秒数据窗，也就是2560点。当信号频率正好是50Hz时，FFT峰值谱线位置正好落在第200条谱线上（50×0.2=10个完整周期对应的谱线），计算出的幅值和相位非常精确。把频率改成50.5Hz后，同样窗长下只有9.8个左右完整周期，峰值能量分散到第201和第202条谱线附近，直接取峰值谱线的幅值会损失大约几个百分点，相位误差也随窗口起始点变化而变化。这就是典型的栅栏效应加频谱泄漏的组合问题。

因此，在同步相量计算里，FFT不能只做一次“峰值查找”就结束，必须搭配窗函数和插值修正，或者改用具备时频分析能力的工具，这就是接下来几个章节要展开的内容。

2. 窗函数法改善FFT泄漏的关键细节

2.1 频谱泄漏与窗函数选择的权衡

加窗是为了压低FFT的旁瓣，让频谱泄漏不污染基波附近的谱线。窗函数种类很多，每一种都对应主瓣宽度和旁瓣衰减之间的权衡。下面的表格是我在电力系统基波相量计算中常对比的几类窗：

为什么同步相量计算通常推荐汉宁窗？因为它主瓣宽度适中，旁瓣衰减足够，而且插值修正公式简单。矩形窗虽然主瓣最窄，频率分辨率最高，但旁瓣太大，谐波很容易把基波附近的小信号淹没。布莱克曼窗旁瓣衰减大，但主瓣也宽，两个相近的间谐波可能会被混在一起。所以做同步相量研究，我一般先把汉宁窗作为基准，等测试信号和误差要求明确之后再换窗型优化。

窗长选择也很有讲究。窗越长，频率分辨率越高，但动态响应越慢。PMU标准里通常要求在不同报告速率下都满足TVE限制，比如每秒50帧报告时，窗长不能太长，否则跟不上系统频率变化。我做实验时常用0.1秒到0.2秒的窗，既保证基波谱线清晰，又不会让时延过大。

2.2 双谱线插值修正的Matlab实现思路

纯加窗虽然压低了旁瓣，但频率偏移导致的栅栏效应并没有消失，峰值谱线未必落在真实频率对应的那条谱线上。双谱线插值就是利用峰值谱线及其相邻谱线的幅值比例，估计出真实频率和幅值。这个思路在Matlab里实现起来并不复杂，核心代码骨架如下：

Fs = 12800; T = 0.1; % 窗长 0.1秒 N = round(Fs * T); % 1280点 sig = your_signal(1:N); % 取一帧数据 w = hanning(N, 'periodic'); % 周期汉宁窗 xw = sig(:) .* w(:); X = fft(xw); mag = abs(X); [~, k0] = max(mag(1:floor(N/2))); % 峰值谱线索引，注意下标偏移 % 取相邻谱线 k1 = k0 - 1; k2 = k0 + 1; if k1 < 1 k1 = 1; end % 双谱线比值参数 beta = (mag(k2) - mag(k1)) / (mag(k2) + mag(k1)); % 查表或多项式计算频率偏移量 d，汉宁窗可近似为 % d = 1.5 * beta; % 这是简化形式，实际应用建议用多项式拟合 d = 1.5 * beta; % 频率估计 f_est = (k0 - 1 + d) * Fs / N; % 幅值修正：汉宁窗下系数近似 % A_est = 2 * (mag(k1) + mag(k2)) * (b0 + b2*d^2 + ...) / N % 这里给出一阶近似 A_est = 2 * (mag(k1) + mag(k2)) / N * (0.5 + 0.25*d);

这段代码里的d就是真实谱线与峰值谱线之间的频偏量，取值范围通常在-0.5到0.5之间。得到d之后，真实频率就等于(k0-1+d)*Fs/N。幅值修正系数与窗函数有关，汉宁窗的修正系数可以展开成d的偶次幂多项式。实际工程里更推荐提前用离线方式把修正系数表算好，运行时查表，这样既快又稳定。

相位计算也要注意。FFT结果的相位是窗函数起点处信号的相位，但加窗会引入一个与d相关的相位偏移。对汉宁窗，这个偏移大约等于-pi*d，所以最终相位估计要补偿回来。这个细节很多人第一次做都会漏掉，结果幅值对了，TVE里相位部分始终超标，最终整体误差卡在临界线上。

3. 希尔伯特-黄变换和小波变换处理非平稳信号的方法论

3.1 HHT：自适应分解与瞬时相量提取

希尔伯特-黄变换（HHT）由EMD和Hilbert变换两步组成。EMD把原始信号自适应地分解成若干个本征模态函数（IMF），每个IMF都可以看作一个幅值和频率随时间变化的振荡分量。对含基波的电力信号做EMD之后，基波成分通常会被分解到某一个或两个IMF中，然后对IMF做Hilbert变换，构造解析信号，就能得到瞬时幅值A(t)和瞬时频率f(t)。

在Matlab中，从R2021a开始，官方已经内置了emd函数，用法如下：

imf = emd(signal, 'SiftRelativeTolerance', 0.01); % imf的每一行是一个IMF，最后一列常常是残余项

需要注意的是，EMD对噪声和采样率比较敏感。采样率太低时，基波频率成分没法被有效分离；采样率太高但数据长度不足，又会出现端点飞翼和模式混叠。我测试时发现，信号里如果含有较大的5次谐波，EMD有可能把基波和谐波混在同一个IMF里，这时需要结合信号的能量占比或者与参考基波的相关系数来筛选IMF。

一个实用的筛选策略是：对每个IMF计算其Hilbert瞬时频率的平均值，如果这个平均值接近系统额定频率（比如50Hz），并且该IMF与原始信号的相关性较高，就把它当作基波分量。提取基波IMF后，瞬时幅值就是解析信号的模，瞬时相位就是解析信号相角的连续展开，再去掉2π跳变就能得到相位序列。

HHT最大的优势是自适应，不需要提前选定基函数，能够追踪幅值调制和相位跳变。但它的理论背景不如FFT或小波严谨，分解结果可能随EMD参数变化而变化，所以在同步相量计算中更适合作为动态相量分析的辅助手段，而不是唯一依赖的工具。

3.2 小波变换：多分辨率下的相位追踪

小波变换通过平移和伸缩一个小波基函数，把信号映射到时间-尺度平面，可以获得不同频带上的时变信息。对同步相量计算，工程上常用连续小波变换（CWT）配合复Morlet小波。复小波能同时提供幅值和相位信息，而且可以在特定频率处提取随时间和相位变化的曲线。

Matlab的Wavelet Toolbox里，最简单的方式是直接调用cwt：

[wt, f] = cwt(signal, Fs, 'amor'); % 'amor'是复Morlet小波（Amorphous? 实际是'amor'代表复数Morlet）

得到的wt是复数矩阵，每一行对应一个频率点。如果想提取接近50Hz处的基波幅值和相位，可以先从f数组里找离50Hz最近的索引，再对该行小波系数取abs和angle。不过需要注意，这样直接取的是“某个尺度”附近的时频系数，它相当于一个带通滤波结果，幅值会受小波基选择影响。要做定量相量测量，需要先对基波频率处的小波系数作归一化标定。

离散小波变换（DWT）也可以用来做谐波分析，但它的频带划分是二进制的，很难正好卡在50Hz，所以不太适合精确相量估计。CWT更灵活，代价是计算量比较大。小波变换对突变和暂态信号非常敏感，适合检测电压暂降、相位跳变过程中的相量变化轨迹。

3.3 两种方法在实际应用中的适用边界

在我自己的对比实验里，HHT和小波的表现呈现出明显的互补性：

• HHT在信号非平稳、幅值缓慢波动时表现更好，因为它能自适应提取时变的幅值和频率。但对含噪信号，EMD容易产生模态混叠，需要额外处理，比如先用阈值滤波或者集合经验模态分解（EEMD/CEEMDAN）。
• 小波变换在暂态突变检测上更直接，频率分辨率也更可控。但它的结果受小波基和尺度范围影响较大，比如用db4和用复Morlet算出来的相位轨迹在突变点附近会有区别，做横向对比时需要固定小波基参数。

从计算时间看，HHT比小波变换慢。一个约0.5秒长、6400个采样点的信号，emd函数在我的电脑上要跑一到两秒，而cwt只需要几百毫秒。如果研究目标是实时相量计算，HHT更适合离线分析或事件后分析，小波还可以通过选择有限尺度和优化实现勉强贴近实时。

因此选型建议是：稳态测量优先用窗函数法，暂态或动态事件分析用小波，复杂非平稳信号的离线精细分析用HHT。

4. 基于Matlab的多算法仿真对比框架

4.1 测试信号与工况设计

做同步相量算法研究，一个规范的测试信号生成模块是跑不掉的。我参考IEEE C37.118里的类型测试和常见文献，设计了四个基础工况：

• 工况A：稳态正弦波，50Hz，幅值1，初相30度
• 工况B：频率斜坡，从49.5Hz线性增加到50.5Hz，斜率0.5Hz/s
• 工况C：基波叠加5次谐波（幅值0.1），并加上白噪声（信噪比40dB）
• 工况D：基波相位在某个时刻阶跃20度，模拟开关操作

生成信号的Matlab代码片段：

Fs = 12800; t = (0: N-1) / Fs; % 工况B f_ramp = 49.5 + 0.5 * t; % 线频率变化 phase = 2 * pi * cumsum(f_ramp) / Fs; % 积分得到相位 sig = cos(phase + pi/6); % 工况C sig = cos(2*pi*50*t + pi/6) + 0.1*cos(2*pi*250*t + pi/3) + 0.01*randn(size(t));

相位初值选择不要总是选0度，因为相位对窗口起点非常敏感。动态工况里用相位累积的方式生成信号，比直接用2pif*t更真实，否则频率变化时相位是断开的，和实际电网连续旋转的相量不同。

4.2 评估指标与实测结果对比

评估指标主要看总向量误差（TVE），定义为：

TVE = sqrt((Xr - Xr_est)^2 + (Xi - Xi_est)^2) / sqrt(Xr^2 + Xi^2)

也就是真实相量和估计相量在复数平面上的距离，除以真实相量的幅值。TVE综合反映了幅值和相位误差，是PMU标准中的核心指标。

我在同样参数下跑过四种算法的仿真，得到一组典型结果代表相对趋势（实际数值会随窗长、小波基和EMD参数变化，只当参照）：

从表格能看出：FFT在纯稳态下精度很高，但频率一变就崩了；加了窗和插值之后，频率斜坡工况下TVE下降了一个数量级。小波在稳态下反而不如窗函数法精确，这是因为小波系数受边界效应影响，需要舍弃数据段两端的部分结果。HHT在含噪工况里TVE优势不明显，因为EMD会把噪声解析成若干虚假IMF，基波IMF的提取不稳定。

4.3 结果背后的物理解释

为什么同一个信号，四种方法结果差这么多？关键在“频率偏差”和“时变信息”这两个维度上。

直接FFT把整窗当作平稳周期信号，频率一旦偏离采样栅格，能量泄漏导致误差。窗函数插值相当于在频域做了更精细的重构，能恢复真实频率和幅值。小波变换相当于一组带通滤波器组，跟踪的是中心频率附近的时变能量，但它的中心频率跟实际频率未必完全重合，而且时间分辨率有限，所以稳态精度不如专门为单频估计设计的插值算法。HHT的EMD没有“基函数频率”概念，它完全依靠信号本身极值点的时间尺度分离模态，在噪声和高次谐波干扰下，IMF的纯度下降，相位计算自然不稳定。

这个对比结果给研究者的启示是：不要只看算法名头，要把算法与信号模型匹配起来。如果只做传统稳态PMU性能测试，窗函数法是性价比最高的选择；如果要做扰动事件分析，小波变换能提供更丰富的时频演化信息；如果研究对象是次同步振荡或非平稳波动，HHT有独有的优势，但必须处理好分解质量和端点问题。

5. Matlab工程实现中的高频坑与优化建议

5.1 采样同步、频率估计与坐标系转换

在Matlab仿真中，大多数人习惯直接设定一个固定采样率，然后认为采样序列和真实时间轴完美对齐。但在工程PMU里，采样必须和GPS/北斗秒脉冲同步，否则时间参考漂移会直接变成相位误差。我做仿真时会把采样率设置成跟50Hz不成整数倍关系的值，比如12800Hz，这样更接近真实非同步采样的情况，算法性能测试也更严格。

频率估计是另一个容易被忽略的环节。窗函数插值法能够同时估计频率，但如果信号频率偏移比较大，双谱线插值的线性近似误差会增加，这时候可以考虑加一个迭代步骤：先粗估计频率，然后以估计频率为参考重新设计采样窗口或修正相位累积，迭代一两次能让TVE进一步下降。不过迭代会增加计算负担，实时系统里要权衡。

还有坐标系转换的问题。三相系统里经常把abc三相变换到dq旋转坐标系，然后再计算同步相量。这个时候要明确参考角度是A相余弦过零还是d轴角度。Matlab里用park函数或者自己写变换矩阵都行，但务必保证测试信号里的相位初值和坐标系参考一致，不然算出来的相位和理论值对不上，容易浪费大量时间找Bug。

5.2 工具箱选型、计算速度与实时化改造

Matlab版本会影响你能否直接用官方函数。官方emd函数从R2021a开始提供，但早期版本只能用第三方工具箱，比如常用的“HHT”开源包。cwt函数在Wavelet Toolbox里，需要注意不同版本中cwt的输入输出格式变化很大，老版本是cwt(signal, scales, 'db4')，新版本是cwt(signal, Fs, 'amor')，我就在版本迁移时被坑过一次。

计算速度方面，我的经验是：FFT加窗函数法在12800Hz采样率、0.1秒窗长下，单帧处理不到几毫秒，完全可以满足每秒50帧的实时PMU计算。小波变换如果要处理所有尺度，计算量会大很多，但可以只计算基波附近几条尺度线，速度会明显提升。HHT则很难做到实时，除非对数据做分段滑窗并且限制迭代次数，否则最好离线使用。

如果后续要把Matlab算法往嵌入式或实时平台迁移，建议先用MATLAB Coder把核心的FFT加窗函数法转成C代码，窗和插值系数预先算好，避免运行时实时生成。滤波器组的系数也可以离线计算并固化，这样移植后代码更稳定。

5.3 代码组织与可复现性经验

我习惯把一个完整的同步相量算法对比工程分成四个模块：信号生成、相量计算、误差评估、结果绘图。信号生成模块单独放，方便切换不同工况；相量计算模块里每个算法做一个函数，输入是采样序列、采样率和算法参数，输出是相量序列；误差评估模块负责计算TVE和频率误差；结果绘图模块统一做图，不然每次手动看数据太低效。

这个过程中一个容易忽略的点是随机数种子。只要信号里加了噪声，就要在生成信号前设置固定的rng(0)之类的种子，否则两次运行结果对不上，算法对比的结论也不可复现。另外，保存结果时我会把Matlab版本、工具箱名称和关键参数写进一个mat格式的结果结构体里，方便后续回溯。

最后再分享一个小技巧：不要只盯着TVE，还要观察相量序列在动态工况下的“相量轨迹”。把实部、虚部画在复数平面上，能直观看到算法是不是发生了振荡或延迟。这个图比堆一堆误差数字更能帮你判断算法的动态行为，尤其是在做小波和HHT对比时，轨迹形态差异非常明显。

做同步相量计算研究，算法本身是工具，对问题的理解和误差评估的设计才决定工作质量。希望这篇内容能帮你少走一些弯路。