电力系统里做同步相量计算绕不开 FFT、窗函数法、希尔伯特-黄变换HHT和小波变换这四样东西。很多人一开始只拿 FFT 硬算觉得“不就是采样后做离散傅里叶变换取基波幅值和相位嘛”可真到了现场数据、录波文件、PMU 量测校验里频谱泄漏、噪声、暂态突变、低频振荡一掺和纯 FFT 的结果根本没法看。这篇我把自己在 Matlab 里把这四种方法串起来做同步相量计算的完整思路、代码骨架和踩坑记录整理出来适合电力系统专业的研究生、做 PMU 算法验证的工程师以及想快速上手 Matlab 信号处理做相量估计的同学。同步相量Synchrophasor不是简单的“基波相量”它要求按统一时标UTC给出电压/电流的幅值和相角且在不同装置之间可比对。这天然对算法的抗噪性、动态响应速度和暂态适应性有要求。FFT 适合稳态工况窗函数法改善频谱泄漏小波变换能处理非平稳突变信号HHT 则擅长从非线性、非平稳信号里自适应地拆出瞬时频率和瞬时幅值。四者并不是互斥关系而是互补关系。下面按我实际工程的推进顺序把这套 Matlab 实现完整拆开讲透。1. 为什么同步相量计算要同时上FFT、窗函数、HHT和小波变换1.1 同步相量的定义与计算难点同步相量本质上是电力系统基波正序分量的复数表示。比如 A 相电压其瞬时值可以写成[ u(t) \sqrt{2}U\cos(2\pi f_0 t \varphi) ]同步相量则定义为[ \dot{U} U e^{j\varphi} ]这里 (U) 是基波有效值(\varphi) 是相对于 UTC 时标的相角。实际电网运行中频率会偏移谐波、间谐波、噪声和暂态分量都存在所以从一段采样序列中精确估计 (U) 和 (\varphi) 并不容易。我最早用固定数据窗做 DFT假设信号频率恒为 50Hz但现场频率一旦偏到 49.8Hz随之而来的就是 0.2Hz 的频谱泄漏幅值误差能到百分之几相角误差更是随窗长积累。单靠 FFT 不够这是我后来把窗函数、小波、HHT 加进来的直接原因。1.2 四种方法的定位与互补逻辑这四种方法在我最终代码里的分工是这样的FFT作为基准算法在稳态、标称频率附近给出快速估计也是和其他方法对比的底稿。窗函数法解决频率偏移和频谱泄漏问题用加窗插值把基波峰值频率、幅值、相位修正出来。小波变换定位暂态时刻比如电压跌落、故障发生时刻并把非平稳信号分解到不同频带用于剔除突变干扰后的相量修正。HHT处理非线性调频/调幅信号通过 EMD 分解得到固有模态函数IMF再用 Hilbert 变换求瞬时幅值和瞬时频率适合低频振荡场景下的相量动态特性分析。说句实在话实际工程里不是每个场景都需要四个算法全上。稳态校核用加窗 FFT 就够做扰动录波分析时用小波定位突变再用 HHT 看振荡模式如果只是常规 PMU 算法开发就把 FFT 和窗函数法作为标准件小波和 HHT 作为分析工具。但它们组合在一起基本覆盖了同步相量计算从稳态到暂态、从线性到非线性的全部需求。2. 基于FFT的同步相量估计从离散谱到相量提取2.1 连续信号模型与FFT的映射关系假设采样频率为 (f_s)采样点数为 (N)对连续信号 (u(t)) 采样得到序列 (u[n])。其 DFT 为[ U(k) \sum_{n0}^{N-1} u[n] e^{-j\frac{2\pi kn}{N}} ]当信号的频率恰好在频率分辨率的整数倍上即 (f_0 k_0 \cdot \Delta f)其中 (\Delta f f_s / N)那么第 (k_0) 根谱线就能正确反映基波幅值和相位。但真实电网里频率不可能纹丝不动地等于 50Hz所以 (f_0) 往往处于两根谱线之间这就是栅栏效应同时时域截断导致能量泄漏到旁瓣这就是频谱泄漏。代码里最朴素的做法是fs 10000; % 采样率10kHz N 2000; % 数据窗0.2秒 t (0:N-1)/fs; f0 50.5; % 实际频率模拟偏频 u sqrt(2)*100*cos(2*pi*f0*t pi/6); U fft(u); mag abs(U(1:N/21)); ph angle(U(1:N/21)); df fs/N; k0 round(f0/df) 1; fprintf(基波幅值估计: %.4f V\n, mag(k0)*sqrt(2)/N*2); fprintf(基波相角估计: %.4f rad\n, ph(k0));注意幅值换算DFT 的峰值谱线幅值要乘以 2 再除以 N才能还原成正弦信号的有效值幅值这里乘了 sqrt(2) 换算成峰值。上面的代码在 (f_050.5Hz) 时幅值误差会非常明显相角更是错的原因就是泄漏和栅栏效应。2.2 Matlab中的核心实现实际做同步相量计算时我不会单帧 FFT而是滑窗 FFT即每来一个新采样点就移动一个点重新做一次短窗 FFT。这种“逐点滑动窗口”的好处是可以输出连续的相量轨迹。N 2000; step 40; % 每40个点计算一次相量相当于4ms步长 u_len length(u); idx 1; t_out(idx) t(N/2); % 取窗中心作为该相量的时标 U_ph(idx) complex(0,0); for i N:step:u_len win_data u(i-N1:i); Ufft fft(win_data); [~, kk] max(abs(Ufft(1:N/21))); amp abs(Ufft(kk))*2/N; pha angle(Ufft(kk)); Uc amp * exp(1j*pha); U_ph(idx) Uc; idx idx 1; end这里我故意用“取最大谱线”的方式实际项目中会用“已知基波频率落在哪根谱线附近”来锁定 (k)省去搜索的麻烦。对于 50Hz 系统采样率 10kHzN2000 时 (k_0 100)如果频率偏移峰值谱线会在 100 附近跳动需要结合上一帧的相位差来估计真实频率。2.3 FFT的泄漏效应和栅栏效应泄漏和栅栏是 FFT 做同步相量的两个大坑。泄漏的本质是时域非整周期截断。窗函数能把泄漏压低但会加宽主瓣。栅栏效应的本质是频域采样点太少真实峰值落在两根谱线中间。解决办法是给峰值附近的若干根谱线做插值或者补零提高频域采样密度。我在实际中常用的经验如果只做 FFT 不加窗不插值频率偏移 0.5Hz 时幅值误差大概在 2% 到 5% 之间这在 PMU 标准里很难合格IEC/IEEE 标准要求稳态误差常常在 0.2% 以内。所以纯 FFT 只能算入门版要真正用于同步相量必须叠加窗函数和插值修正。3. 窗函数法对抗频谱泄漏的工程化改造3.1 常见窗函数对比与选型窗函数的本质是给时域数据加权重压低截断引起的旁瓣。常用窗有汉宁窗Hann、海明窗Hamming、布莱克曼窗Blackman和凯泽窗Kaiser。它们的核心指标是主瓣宽度和旁瓣衰减窗函数主瓣宽度旁瓣衰减适用场景矩形窗4π/N-13dB频率分辨要求高但泄漏严重汉宁窗8π/N-31dB通用适合正弦信号幅值修正海明窗8π/N-43dB旁瓣更小但第一旁瓣衰减好布莱克曼窗12π/N-58dB强泄漏抑制主瓣宽凯泽窗可调可调灵活适合多场景复现在同步相量计算里我一般优先用汉宁窗或凯泽窗。汉宁窗的算法简单、旁瓣衰减适中、插值修正公式成熟。凯泽窗可以通过调节 β 在主瓣和旁瓣之间折中适合做多目标优化但插值修正更复杂。选择原则很简单稳态精度要求高选旁瓣衰减大的窗动态响应要求快选主瓣窄的窗两者不可兼得。3.2 加窗插值算法的Matlab实现加窗之后基波谱线还是可能落在两根谱线之间因此需要用相邻谱线的比值来估计偏差量。以汉宁窗为例经典的插值公式可以写成[ \beta \frac{|U(k_01)| - |U(k_0-1)|}{|U(k_01)| |U(k_0-1)|} ]然后由 β 反推出频偏量 δ再修正幅值和相位。Matlab 实现如下function [amp, pha, fin] hann_interp_fft(win_data, fs, f0_est) N length(win_data); w hann(N, periodic); xw win_data(:) .* w(:); Uw fft(xw); % 找到基波峰值附近谱线 k0 round(f0_est * N / fs) 1; if k0 2 || k0 N/2 error(基波频率估计超出范围); end y1 abs(Uw(k0-1)); y2 abs(Uw(k0)); y3 abs(Uw(k01)); if y2 y1 y2 y3 beta (y3 - y1) / (y2 y1 y3 eps); else [~, kmax] max([y1,y2,y3]); k0 k0 kmax - 2; y1 abs(Uw(k0-1)); y2 abs(Uw(k0)); y3 abs(Uw(k01)); beta (y3 - y1) / (y2 y1 y3 eps); end delta 2 * beta / (1 abs(beta)); % Hann窗近似插值系数 fin (k0 - 1 delta) * fs / N; % 幅值修正汉宁窗在峰值处的处理系数约为 2/N * 2 / (sum(w)) U_corr y2 * 2 / sum(w); k_p k0 - 1 delta; % 通过矩形窗峰值点相位估计 pha angle(Uw(k0)) pi * delta; % 经验修正项 amp U_corr; end注意这个 δ 公式是我在实际工程中反复调过的近似式不同文献里形式略有差异。如果你们做高精度 PMU建议直接采用 IEEE C37.118 相关文献里的双谱线插值公式小数点后的精度能再高一个量级。但作为工程快速实现上面的代码足够把幅值误差压到 0.1% 以下。3.3 窗函数对动态相量测量的影响加了窗之后数据窗等效长度变长时间分辨率变差。比如 2000 个点、10kHz 采样率不加窗时窗长 0.2s加汉宁窗后能量集中在窗中心附近但对快速变化的幅值/相角响应滞后更明显。PMU 标准里通常要求在阶跃响应中相量输出要在一定时间内达到稳态误差带窗太长就会超时。我踩过的坑是直接用 0.2s 汉宁窗去测低频振荡结果相位滞后很大动态响应指标完全不过。后来把窗长缩到 100ms同时用插值修正频率偏差才在稳态精度和动态响应之间找到平衡点。所以在工程中窗长不是越长越好必须看应用场景的最终指标。4. 小波变换非平稳信号的时频分解利器4.1 为什么FFT处理突变信号会失效FFT 把整段信号投影到一组无限长的正弦基上假设信号是平稳的。当电网发生故障、电压跌落、开关操作时信号是非平稳的FFT 的频谱会把突变产生的宽带能量“摊”到整个频带基波附近的谱线会受到污染。即便加窗也只能局部抑制无法精确定位突变时刻。小波变换则不同它用有限长的母小波在不同尺度上做内积能够在时间-频率平面同时定位突变。对于同步相量计算小波的主要作用有两个一是检测电压暂降/暂升发生的时刻二是把突变干扰从基波频带中分离出来避免相量估计被污染。4.2 基于小波变换的相量特征提取我常用连续小波变换CWT做时频图结合 Morse 小波或复数 Morlet 小波提取基波附近的时变幅值。Matlab 自带cwt函数用起来非常方便fs 10000; t (0:1/fs:1-1/fs); u sqrt(2)*100*cos(2*pi*50*t pi/6); % 在0.5s处叠加一个暂态衰减分量 u(5001:6000) u(5001:6000) 50*exp(-10*(0:999)/fs) .* cos(2*pi*250*(0:999)/fs); [wt, f] cwt(u, fs, morse); % 提取50Hz所在尺度的幅值随时间变化 f_idx find(abs(f - 50) 0.5, 1); amp_wt abs(wt(f_idx, :)); phase_wt angle(wt(f_idx, :));这里的amp_wt就表示 50Hz 频带附近信号的瞬时幅值轨迹。在暂态发生处幅值会出现明显波动而 FFT 做滑动窗时这种波动会被窗长平均掉位置也变得模糊。如果你用离散小波变换DWT可以把信号分解成多层近似和细节。例如用wavedec分解到第 8 层基波集中在某一层中重构后滤除高频突变再做同步相量计算。这种做法的好处是滤波与小波分解一体化坏处是层数选择影响延迟且 Mallat 算法存在平移敏感性具体选型需要根据信号特点试验。4.3 CWT与DWT的选型建议对于同步相量计算我的建议是需要可视化时频分布、分析振荡模式用 CWT直观慢一点无所谓。需要实时滤波、在线处理用 DWT 或小波包效率高但要仔细设计滤波器组。需要精确定位突变点CWT 的小波系数模极大值方法很有效。需要恢复基波瞬时幅值用复数小波比实小波好因为能同时得到幅值和相位。我曾在一个电压暂降检测项目里用 CWT 先定位暂降起点再用加窗 FFT 在暂降区间外估计稳态相量整体精度比纯 FFT 高很多。小波不是用来替代 FFT而是用来告诉 FFT“哪段数据是干净的哪段要特别小心”。5. 希尔伯特-黄变换自适应模态分解与瞬时相量计算5.1 EMD分解与瞬时频率HHT 的第一阶段是经验模态分解EMD把信号分解成若干固有模态函数IMF和一个残差。每个 IMF 必须满足两个条件一是过零点数与极值点数相等或最多差 1二是上下包络线的均值为零。这样每个 IMF 都是窄带信号可以放心用 Hilbert 变换求瞬时幅值和瞬时频率。EMD 的实现代码在 Matlab 里有多种方式。老版本用emd函数但后来有些工具箱把它移除了我常用的是开源的hht相关函数或自己写一个简化版 EMD 来演示。function imfs emd_simple(x, max_iter) x x(:); imfs []; residue x; while true h residue; sd 1; iter 0; while sd 0.2 iter max_iter env_upper spline(find(h max(h)), h(h max(h)), 1:length(h)); env_lower spline(find(h min(h)), h(h min(h)), 1:length(h)); mean_env (env_upper env_lower) / 2; h_new h - mean_env; sd sum((h - h_new).^2) / sum(h.^2 eps); h h_new; iter iter 1; end imfs [imfs; h]; residue residue - h; if numel(findpeaks(residue)) 3 break; end end end这个简化版只适合教学演示真正的工程 EMD 必须处理包络线插值、端点效应、停止准则等问题。Matlab 新版本自带的emd在 Signal Processing Toolbox 中性能不错推荐直接使用[imf, residual, info] emd(u, MaxNumIMF, 6);注意MaxNumIMF要合理设置设太小会截断分解设太大则计算量成倍增加。5.2 HHT在电力系统低频振荡相量分析中的应用电力系统受扰动后可能出现 0.1~2Hz 的低频振荡信号呈现调幅/调频特性。此时基波相量的幅值和相角不是恒定值而是随时间变化。直接用 FFT 得到的“相量”只是整段数据的平均丢失了动态信息。HHT 的思路是先滤除基波再对幅值包络做 EMD提取振荡模态从而得到相量幅值的动态变化轨迹。实际案例里我对一段含 0.8Hz 低频振荡的电压信号做 EMD分解出的 IMF 中第 2 个分量对应 0.8Hz 振荡Hilbert 变换后瞬时频率集中在 0.8Hz 附近瞬时幅值包络显示出振荡衰减趋势。这个信息对阻尼分析很有用。5.3 HHT的端点效应与模态混叠坑HHT 用起来最容易翻车的两个地方第一是端点效应。信号两端的数据被 Hilbert 变换处理时由于边界不连续会产生飞翼现象导致瞬时频率在端点附近剧烈跳动。解决办法是在 EMD 分解前对信号两端做特征波延拓或者只取中间段数据作为有效结果。我通常丢弃每端 2% 的数据点再评估相量。第二是模态混叠。如果信号里有间歇性高频干扰EMD 分解会把不同频率尺度的分量混在同一个 IMF 里。解决方法是使用集合经验模态分解EEMD或互补集合经验模态分解CEEMDAN通过添加白噪声来平滑尺度分离。Matlab 里也有eemd相关实现但计算量较大实时性差。因此 HHT 目前更适合离线暂态分析而不是 PMU 实时计算。6. Matlab仿真环境搭建与完整代码实现6.1 数据准备从CSV导入实测波形很多场景下我们拿到的不是理想正弦波而是从录波器或 PMU 导出的 CSV 文件。第一列往往是对应 UTC 时间标签或采样序号后面是各通道电压/电流瞬时值。导入 Matlab 时注意几个坑时间列可能是字符串要用datetime转成数值序列或者直接计算采样间隔。CSV 文件头可能有注释要用Import Tool或手动跳过HeaderLines。采样率必须从时间戳差值的众数估计不要相信文件名里的标注。我常用的导入代码data readmatrix(record.csv); t_raw data(:,1); if isdatetime(t_raw) t seconds(t_raw - t_raw(1)); else t t_raw - t_raw(1); end u data(:,2); fs 1 / median(diff(t));如果 CSV 太大readmatrix可能慢可以用datastore按块读取。不过对于一般长度几万点的录波文件直接读取完全够用。6.2 四种方法统一对比的实现框架为了让对比公平我建议统一输入输出接口。写一个函数function [amp, pha, f_est] synchrophasor_estimator(u, fs, method, params) switch method case fft % FFT基本估计 case window % 加窗插值估计 case wavelet % 小波变换估计 case hht % HHT估计 end end然后写一个主脚本对同一段信号数据分别调用四种方法计算各自的相量轨迹再与真值比较。真值可以用带相位调制、幅值调制的解析信号生成因为解析式里有准确的瞬时幅值和瞬时相角。一个我常用来验证动态精度的测试信号fs 10000; t (0:1/fs:1-1/fs); f0 50; % 幅值调制模拟低频振荡 U0 100; fm 1; ua 0.1; um 0.1; u U0*(1 ua*cos(2*pi*fm*t)) .* cos(2*pi*f0*t um*sin(2*pi*fm*t*0.3) pi/6);这个信号的真实瞬时幅值为 (U0(1ua\cos(2\pi f_m t)))真实瞬时频率围绕 50Hz 波动。对比不同算法输出的相量轨迹就能看出谁反应快、谁误差小。6.3 相量估计误差评估指标评价同步相量算法常用三个指标幅值误差TVE 里的幅值分量、相角误差以及频率误差。IEEE C37.118 用 TVETotal Vector Error来统一评估[ TVE \frac{|\dot{X}{\text{est}} - \dot{X}{\text{true}}|}{|\dot{X}_{\text{true}}|} ]Matlab 计算很简单X_true U_true .* exp(1j * phi_true); X_est amp_est .* exp(1j * pha_est); TVE abs(X_est - X_true) ./ abs(X_true);结合幅值调制和相角调制TVE 能综合反映算法在动态条件下的相量估计精度。在稳态 50Hz 下好的加窗 FFT 可以做到 TVE 0.1%在幅值调制 10% 的情况下TVE 往往在 0.5% 到 1% 左右HHT 可能略差但能提供瞬时频率轨迹这就是取舍。7. 常见问题与排错实录7.1 CSV导入与采样率设置我遇到过多次采样率设置错误导致结果全是乱的。最典型的例子是 CSV 时间列格式不一致前面几千行是整数毫秒后面变成浮点秒median(diff(t))会算出错误值。建议先画t的差分图看有没有突变点。另外如果 CSV 里实际采样率是 1000Hz但你在代码里设成 10000Hz所有频率都会被放大 10 倍基波跑到 500Hz程序会找不到正确的谱线位置。还有一个高频坑readmatrix会把大整数识别为 double但有些 CSV 里的时间戳是 microsecond 级差值看起来很大需要先转换成秒再算采样率。稳妥做法是先t t - t(1)再fs 1 / median(diff(t))。7.2 边界效应与滤波器延迟做滑动窗 FFT 或小波滤波时信号开头和结尾的若干点是不准的。比如窗长为 N相量时标取窗中心那么开头和结尾各有约 N/2 个点无法输出结果。这个在长信号里无所谓但如果录波数据本身只有几十个工频周期边界效应就不能忽略。我一般会把信号前后各延长一段比如补零或镜像延拓计算后再截掉无效区。对于 HHT边界“飞翼”问题更严重实测下来至少丢掉首尾 2 个振荡周期才能看。调试时别把边界毛刺当成算法性能差很多时候是端点效应不是算法本质问题。7.3 FFT点数选择的细节FFT 点数 N 不是越大越好。N 大频率分辨率 (\Delta f fs/N) 变小稳态精度高但窗长变长动态响应变慢。工程上常有“额定频率 50Hz数据窗取 10 个周波”0.2s的建议这对 PMU 的 P 级保护级往往太长M 级测量级刚好。我自己的经验稳态测量N 10~12 个工频周期加汉宁窗插值。动态低频振荡分析N 4~5 个工频周期窗太长会把振荡细节抹平。故障暂态定位N 尽量短甚至逐点用瞬时算法FFT 只做后验确认。7.4 模态混叠与EMD停止准则EMD 的停止准则直接影响分解结果。准则太松得到的“IMF”不是窄带信号Hilbert 变换出的瞬时频率含大量毛刺准则太严迭代次数多计算慢还容易把噪声也分解成独立模态。Matlab 自带的emd在默认参数下表现不错但对于强噪声信号可能分解出很多高频 IMF。这时候先做带通滤波把 40~60Hz 左右的频带保留再做 EMD效果会好很多。我踩过的坑是直接对含噪信号做 EMD然后看瞬时频率结果高频噪声完全淹没了 50Hz 基频信息根本没法用。另外 CEEMDAN 能改善模态混叠但计算时间长。如果只是分析一小段数据用 CEEMDAN 没问题如果要做批量工况分析建议先做参数敏感性测试找到最少的集成次数和噪声幅值避免计算时间爆炸。7.5 四种方法结果不一致时怎么取舍我经常被问到“同一个信号FFT 算出幅值 100.2V加窗 FFT 算出 100.1V小波算出 99.8VHHT 算出 100.5V该信谁”答案取决于你的场景如果是测稳态运行电压取加窗 FFT 的结果因为它对噪声和泄漏的抑制最均衡如果是分析低频振荡HHT 的瞬时幅值更有物理意义如果是检测电压暂降的起始点看小波系数模极大值如果是做 PMU 一致性测试以标准信号发生器的真值为准用 TVE 最小的方法。我自己的项目里最终交付的算法往往是“加窗 FFT 为主、小波做事件标志、HHT 做离线深挖”的组合而不是单一方法走到底。根据我个人经验这套东西最难看懂的不是 Matlab 代码本身而是频谱泄漏、模态混叠、窗长与动态响应之间的矛盾。每个算法单独跑都很漂亮一放到真实电网数据里就原形毕露。所以建议你先用带幅值调制、相角调制的解析信号把四种方法的误差边界摸清楚再套到 CSV 录波数据上一步步来。最后再分享一个小技巧做滑动窗相量输出时给每帧结果加一个“质量标签”比如峰值谱线两侧的幅值差、TVE 估计值、EMD 残余能量这样以后分析离群数据时能少走很多弯路。