同步相量算法对比：FFT、窗函数、HHT与小波变换的Matlab实践

做电力系统同步相量计算研究，最容易踩的坑就是——把FFT、窗函数法、希尔伯特-黄变换、小波变换四条路线各自跑一遍，得到几张漂亮的对比图，然后发现不知道该信谁。FFT快、窗函数法稳、HHT自适应、小波能抓暂态，这个“常识”谁都会背，但真正到了Matlab里拿同一组测试波形逐一实现时，结果往往跟你背的常识不完全一致。这篇文章既是一次对比研究，也是一份可以直接参照的操作笔记：从IEEE C37.118对同步相量的定义和测试要求出发，用Matlab构造包含谐波、噪声、频率偏移和幅值突变的测试信号，把四种算法全部实现了，并详细拆解每一步的数学逻辑、参数选择和真正决定精度的细节。适合正在做PMU算法仿真、毕设选题或者电力信号处理研究的朋友参考。

1. 同步相量计算的问题本质：四条路线为什么不该被“公平”对比

1.1 同步相量到底在算什么

先说清楚对象。同步相量（Synchrophasor）不是简单地把电压波形幅值和初相角算出来，而是要以UTC绝对时标为参考，给出基波分量在某一时刻的幅值和相角。电网频率在正常运行时会围绕50Hz缓慢波动，动态过程中还会出现谐波、间谐波、次同步振荡甚至幅值突变，所以“算一个相量”这件事天然带有两个互相矛盾的要求——精度和速度。IEEE C37.118标准给出了包括总矢量误差（TVE）、频率误差、频率变化率误差以及响应时间在内的一整套指标体系，这比“看波形对不对”严苛得多。

TVE的定义可以这样理解：把测量得到的相量（幅值和相角）与真实相量做矢量差，再除以真实相量的幅值。哪怕幅值误差1%，或者相角误差大约0.57°，TVE就已经接近1%。这意味着任何算法都不能只盯着“幅值准确”，相位基准错一点，整套数据都没法用。我在做算法对比时，第一步不是跑仿真，而是把TVE的这个几何含义写在纸上——后面所有调参行为，最终都是在跟这个指标较劲。

1.2 四种子方法的技术出身

四种方法出身不同，这是理解后面所有对比的关键。FFT是经典的离散频域分析工具，它假设信号是周期性的，在整周期采样条件下能精确提取各次谐波；窗函数法是在FFT基础上针对非整周期采样和频谱泄漏做的修正；希尔伯特-黄变换（HHT）由Huang等人提出，核心是先通过经验模态分解（EMD）把非平稳信号拆成本征模态函数（IMF），再做Hilbert变换得到瞬时幅值和瞬时频率，它不预设基函数，是自适应的；小波变换则通过尺度伸缩和平移，在时频平面上提供多分辨率分析，对暂态突变尤其敏感。

我把四个方法的技术定位整理成了一张表：

这里要特别提醒：正因为出身不同，四种方法并不能在“谁更准”上做简单排名。准确的问题是——在C37.118的哪个测试场景下，哪种方法能同时满足精度和响应时间要求。我后面第5节会给出同一组测试数据下的实际表现。

2. 经典路线：FFT与窗函数法的组合拳怎么打

2.1 FFT提取相量的两个固有缺陷

先看最简单的实现。对一段N点采样序列x(n)，做N点FFT后，基波对应的谱线就是我们要的那个分量。幅值A=2|X(k0)|/N，相角φ=angle(X(k0))，Matlab代码很短：

fs = 12800; f0 = 50; N = fs/f0*10; % 10个工频周期 x = ...; % 电压采样序列 X = fft(x)/N; k0 = round(f0/fs*N); % 基波对应的谱线索引 A = 2*abs(X(k0+1)); % 直流在X(1)，基波要加1 phi = angle(X(k0+1));

但这个写法暗藏两个缺陷。第一是栅栏效应：基波频率f0不一定正好落在整数谱线上。当系统频率偏离50Hz时，基波对应谱线会落在两根谱线之间，直接取临近谱线会带来幅值和相角的系统误差。第二是频谱泄漏：采样数据窗长度如果不是信号周期的整数倍，DFT会把基波能量泄漏到旁边的谱线上去。

具体数量级可以算一下。采样率12800Hz，256点/工频周期，10周波窗长正好2560点，频率分辨率Δf = fs/N = 5Hz，基波谱线在k=10这根线上。系统频率偏移0.5Hz时，峰值实际偏移0.1根谱线，看起来不大，但泄漏造成的TVE可以达到百分之几，直接超标。C37.118标准里专门设了“频率偏移”测试场景，就是这个原因——同步相量测量的前提是系统频率未知，必须边测边估。直接做FFT等于默认基波谱线位置已知，这个前提在实际电网中根本不成立。

2.2 加窗、插值、相位补偿：三步修正

解决泄漏问题最成熟的手段就是窗函数法。思路是：先用一个旁瓣衰减较大的窗函数把时域数据截断，抑制非整周期采样造成的频谱拖尾；然后在频谱峰值附近取两根或三根相邻谱线，利用它们的幅值比反推真正的谱峰位置，得到频率偏移量δ；最后用δ修正幅值和相位。

以汉宁窗为例，加窗后的频谱中，峰值谱线k1与相邻谱线k1+1的幅值比β，可以换算成偏移量δ。汉宁窗下有一个近似关系δ≈(2β−1)/(β+1)，更精确的做法是事先用多项式拟合出β到δ的查表曲线。得到δ之后，实际频率f0' = (k1+δ)·fs/N，幅值修正系数由窗函数主瓣形状决定，相角修正则把DFT参考点从数据窗起点移到窗中心。核心代码：

w = hann(N, 'periodic'); xw = x .* w; Xw = fft(xw) / sum(w) * 2; % 加窗后的幅度归一化 [~, k1] = max(abs(Xw(2:N/2))); beta = abs(Xw(k1)) / abs(Xw(k1+1)); % 峰值与右邻谱线比 delta = (2*beta-1)/(beta+1); % 汉宁窗近似插值公式 f_est = (k1+delta) * fs/N; % 估计实际频率

相位修正这一步最容易被忽略。FFT输出的相位对应数据窗第一个采样点，而同步相量要求的是UTC整秒（或PPS边沿）时刻的相位。因此需要做相位补偿：把参考点从窗起点移到窗中心，再对齐到同步时标。忘了这一步，即便幅值算得很准，相角也会整体差一个固定角度，在动态测试里表现为明显的TVE。

2.3 窗长选择和数据窗设计

窗函数法里，窗长和窗型是两个互相制约的参数。长窗旁瓣抑制能力好，但响应时间变长；短窗响应快，但频率分辨率和谐波抑制能力下降。真实PMU常用双窗设计：稳态用10周波窗保证精度，动态事件用2~3周波窗保证响应速度，两个通道并行切换。

窗型选择上，汉宁窗适合一般测量，布莱克曼窗旁瓣衰减更猛但主瓣变宽，平顶窗（flattopwin）的幅值误差对频率不敏感，适合先估频再精确测幅值的场景。我自己的习惯是：仿真阶段先用汉宁窗加双谱线插值跑通整个链路，再去试其他窗型。汉宁窗的解析公式简单、随机误差小，平顶窗虽然幅值误差小，但相位插值处理稍麻烦。如果从零开始做，优先汉宁窗没毛病。

3. 自适应路线：HHT处理非平稳信号的逻辑

3.1 EMD分解：把混合信号拆成“干净成分”

HHT的第一步是EMD。它假设任何复杂信号都能分解成有限个本征模态函数（IMF），每个IMF的上下包络关于时间轴对称，且极值点数和过零点数至多差一个。分解过程直观地说就是“剥洋葱”：对原始信号x(t)，找出所有局部极大值和极小值，用三次样条分别拟合成上包络和下包络，取均值m1，得到第一个分量h1=x−m1。如果h1不满足IMF条件，就重复这个过程，直到满足为止，得到第一个IMF。然后从原信号中减掉这个IMF，对剩余部分重复同样的操作，直到余量单调或足够小。

Matlab本身没有官方EMD函数，大多数研究者用的是MathWorks File Exchange上Rilling等人发布的开源emd包，调用方式很直接：

imfs = emd(x); % 每一行是一个IMF分量，自上而下频率由高到低 % 实际使用时建议限制IMF个数，避免分解出过多伪分量 imfs = emd(x, 'MaxNumIMF', 6);

这个“剥洋葱”过程不依赖预设基函数，理论上对频率缓慢变化、幅值调制的信号有很好的自适应性。电力系统里电压幅值波动、功角摇摆等非平稳过程，正是它的目标场景。

3.2 Hilbert变换求出瞬时幅值和瞬时频率

EMD分解出IMF后，对每个IMF做Hilbert变换构造解析信号：z(t)=c(t)+j·H[c(t)]。解析信号的幅值就是瞬时幅值a(t)，相位θ(t)对时间求导就得到瞬时频率f(t)。对电力系统基波相量来说，基波所在IMF的a(t)就是相量幅值的动态变化轨迹，θ(t)经过解卷绕后，减去2πf0·t得到的剩余相位就是相角偏差。

ht = hilbert(imf_50Hz); % imf_50Hz是基波所在IMF A_hht = abs(ht); theta = unwrap(angle(ht)); f_hht = diff(theta)/(2*pi*Ts); % 瞬时频率序列

有个容易犯迷糊的点：瞬时频率的定义虽然数学上简洁，但要求信号是窄带的，否则Hilbert变换得到的相位没有明确物理意义。这恰恰是EMD要先做分离的原因。在相量计算里，如果EMD没能把基波和邻近频率成分干净地分开，瞬时频率会出现明显的抖动甚至跳到负值，这通常是模态混叠的征兆，后面会细说。

3.3 模态混叠、端点效应和计算代价

HHT在实际应用中三个坑最多。一是模态混叠：当信号中有频率接近的分量时，EMD可能把它们分到同一个IMF，或者同一个频率被拆到两个IMF。工程经验是用EEMD或CEEMDAN（加入辅助白噪声再总体平均）来缓解，代价是成倍增加计算量。二是端点效应：三次样条包络在数据两端容易发散，导致IMF两端出现明显摆动，瞬时幅值在边界处会翘起或凹陷。解决办法是数据延拓，比如镜像延拓、极值延拓。做离线分析时，我会故意多采一段数据，分析时只取中间一段，两端直接扔掉，这是最简单粗暴也最有效的方法。三是计算量：EMD是迭代过程，数据一长就慢得让人焦虑。我在12800Hz采样率下分析几秒数据，纯Matlab循环实现的emd包要跑相当久，和FFT完全不在一个量级。

这三个坑决定了HHT更适合离线诊断分析，而不是实时PMU的核心测算法。不过对于研究课题来说，它的价值在于：在频率缓变、幅值调制的场景下，能给出FFT给不出的时变相量轨迹。

4. 多分辨率路线：小波变换对相量计算的独特价值

4.1 为什么电力暂态信号需要“变焦镜头”

FFT和加窗FFT的问题在于：窗长一旦确定，整段数据的频率分辨率就固定了。要抑制谐波就要长窗，要捕捉暂态就要短窗，鱼和熊掌不可兼得。小波变换通过一组可伸缩平移的基函数把信号投影到不同尺度上，低频段用宽窗提高频率分辨率，高频段用窄窗提高时间分辨率。对于电压跌落、相位跳变、短路冲击这类包含突变分量的信号，小波能同时给出突变发生时刻和基波参数的动态变化。

这一点在相量计算里非常实用。PMU不只是稳态仪表，它还要捕捉动态事件。基于短窗FFT的算法在电压突变后需要重新积累一个窗长的数据才能给出稳定读数，而小波方法的时频“变焦”能力可以让突变定位和相量估计同时进行。

4.2 用小波重构提取基波分量

具体到同步相量计算，最常用的套路不是直接拿小波系数当相量，而是分两步：先用离散小波分解，把基波所在的频带单独重构出来，得到一个“干净的基波时域波形”，再对这个波形用Hilbert变换或最小二乘正弦拟合法求瞬时幅值和相位。

Matlab里用wavedec做多级分解：

[C, L] = wavedec(x, 6, 'db4'); base = wrcoef('a', C, L, 'db4', 6); % 第6层近似分量重构 h = hilbert(base); A_wave = abs(h); theta_wave = unwrap(angle(h));

分解层数的选择必须和采样率对应。以fs=12800Hz为例，第6层近似分量对应的频带大约是0~100Hz，基波50Hz正好落在频带中部，而3次谐波150Hz、5次谐波250Hz都落在更高频带的细节分量里。这样重构出来的时域波形基本就是纯净的基波分量。实际使用前，可以用freqz看一下所选小波滤波器的实际通带，做到心里有数。

如果信号里存在100Hz以下的间谐波或次同步分量，第6层近似就不够干净，这时可以用小波包分解做更细的频带划分，让50Hz单独落在一个子带里。Matlab里用wpdec和wprcoef，节点选择取决于采样率和分解层数，需要根据频率分辨率逐节点核对。

4.3 小波基选择、边界效应与稳态纹波

小波基的选择会影响重构质量。db4是电力暂态分析的经典选择，紧支撑、与突变信号匹配较好；sym8对称性好、相位失真相对小；如果用连续小波变换（CWT）做时频脊提取，常用复Morlet小波，相位信息更完整。

和HHT一样，小波重构也有边界效应，因为滤波器在数据两端拿不到完整上下文。处理办法不外乎延拓或丢弃两端数据。另外还有个容易被忽略的问题：小波滤波器通带内的波纹会让重构基波幅值出现小幅振荡，如果要做高精度TVE评估，这个纹波会直接变成误差。我的做法是在重构信号之后加一个窄带平滑滤波器，或者用最小二乘拟合正弦参数，把纹波对幅值和相位的影响压下去。

5. 同一组测试数据，四个方法的高下之分

5.1 测试信号与评价指标设计

我在Matlab里构造了一组与C37.118测试思路对齐的信号，统一采样率12800Hz，时长1s。场景包括：

• 场景A：纯50Hz稳态。
• 场景B：频率偏移到50.5Hz。
• 场景C：50Hz基波叠加5次、7次、11次谐波，幅度分别为基波的5%、3%、2%。
• 场景D：幅值在0.4s时从1.0 pu跌落到0.8 pu，0.6s恢复。
• 场景E：基波幅值按1Hz正弦调制（±5%），模拟动态摆动。

每个场景先算出理论相量（幅值、相角随时间变化），再计算各算法的TVE、频率误差和响应时间。所有算法都统一使用相同的输入数据，不额外做预处理，这样对比才有意义。

5.2 各方法的实测表现对比

我把实测结果汇总成表（结果基于我的仿真条件，参数不同会有合理浮动）：

这个结果符合理论预期，但也有几个反直觉的地方。比如HHT在稳态场景下并不比加窗FFT更准，甚至由于端点效应和EMD分解的不确定性，TVE反而偏高；又比如小波在幅值突变時定位很准，但重构出来的过渡段波形会有振铃，不能直接拿来做相量输出。仿真跑完之后，我对“哪种方法最好”这个问题基本没兴趣了，更想搞清楚的是——每个场景下哪个环节在拖后腿。

5.3 选型建议

根据实测结果，我给出的选型建议很简单：

• 目标是工程标准PMU：直接选窗函数FFT，配合GPS授时和多窗协同，这是主流路线。
• 做动态振荡、次同步振荡分析：HHT最能给出模态包络变化，但要做好边界延拓和EEMD去模态混叠。
• 做暂态扰动检测、行波故障测距：小波优势明显，DWT或CWT能精确定位突变时刻。
• 做课程设计或科研演示：FFT加窗函数法作基础对照组，HHT作为特色，小波作为补充，三组都做，论文结构会很完整。

我个人最常用的组合是：加窗FFT负责常规相量输出，小波负责事件诊断，HHT只在需要细致刻画振荡模式时才调用。这样既保证基础精度，又覆盖动态场景。

6. Matlab工程化中的几个实操提醒

6.1 数据预处理：直流分量和频率模糊

无论用哪种算法，预处理都可能比算法本身更影响结果。第一件是去直流：直流分量对FFT的零频附近有泄漏，对EMD分解也会产生额外IMF，先用均值或者detrend把直流去掉，后面会省很多事。第二件是频率粗估：做谱线插值之前，如果完全没有频率先验信息，建议先用过零检测或自相关粗略估计基频，把粗估值作为谱线搜索的引导，否则在多谐波场景下可能把峰值谱线认到高次谐波上去。第三件是采样率偏移：真实采集系统里采样钟和GPS时标不一定严格对齐，要先做重采样对齐，否则相位会产生系统性偏差。

6.2 相位参考点和时标的处理

这里值得单独说。FFT输出的相位对应数据窗第一个采样点，而同步相量要的是同步时标时刻的相位，两者之间差一个固定角度。我在代码里会强行统一约定：所有算法输出的相角，在处理链路的最后一步都修正到以数据窗中心为参考点，再对齐PPS时标。仿真验证时，生成信号的时候就把参考时间0设为PPS时刻，直接对“窗中心=PPS时刻”的数据段计算，这样便于后面对比理论值。这个约定看着简单，但能省掉大量排查“为什么相角对不上”的时间。

6.3 计算性能与实时性的取舍

性能对比也是选型的重要依据。FFT和窗函数法计算量最小，一个N点FFT毫秒级完成，MCU和FPGA都能跑实时。HHT最慢，多轮迭代加三次样条拟合，数据稍长就会让人怀疑人生，离线分析时建议先用粗采样或分段处理。CWT在信号很长时也慢，可以限制voices数量和分析频率范围来加速。实际PMU的实时算法几乎都用加窗DFT或递归DFT，HHT和小波大多用于后台分析，不是前端量测。这个现状短期内不会变。

6.4 一个验证技巧：先构造“完美数据”再逐步加干扰

调试算法时最容易犯的错是直接用复杂信号，出错后不知道是哪个环节的问题。我的习惯是：先给50Hz纯信号、整周期窗，确认相位输出和理论一致；然后加频率偏移，检查插值是否修正；再加谐波，检查泄漏是否被抑制；最后加幅值突变和噪声，评估动态响应。每一步都保存中间结果的曲线，出了问题能直接定位到算法模块。这套流程虽然啰嗦，但比直接跑完整仿真再猜问题高效得多。

6.5 关于emd包的使用建议

最后说下emd包。开源的Rilling版emd，默认参数对短信号很容易分解出一堆伪IMF。建议限制MaxNumIMF，并且对每个IMF都看一眼频谱，只挑频谱峰值在50Hz附近的那个作为基波分量。不要默认第一个IMF一定是基波——第一个IMF往往包含的是噪声或高次谐波。另外，EMD对采样点数比较敏感，点数太少时包络拟合极不稳定，做EMD前尽量保证数据长度覆盖足够多的基波周期。

回头看我这次把四条路线完整实现下来的体会，最重要的不是哪个方法赢了，而是你得先清楚自己手头的问题是稳态测量、动态跟踪还是暂态诊断，再决定用哪把尺子。我自己重做一遍的话，会先花80%精力把加窗FFT的工具链打磨扎实，包括频率粗略估计、谱线插值、相位参考修正这一整套——这套东西在任何方案里都绕不开。然后再把HHT和小波作为动态场景的补充武器逐个加进来。如果你也在做类似的相量计算研究，不妨按这个顺序走，能少走不少弯路。