简介本资源系统讲解全周傅氏算法的数学原理与工程实现面向电子信息工程、计算机及数学专业本科生适用于课程设计、期末大作业及毕业设计中对电力系统谐波分析、信号频谱提取等典型场景的需求。压缩包共3个文件862KB含1份详尽PDF理论文档、1个Simulink仿真模型.mdl用于算法动态验证、1个结构清晰的MATLAB主程序.m支持matlab2014a至2024a多版本直接运行。代码采用参数化设计关键变量如采样频率、信号周期、谐波阶数等均集中定义并配有中文注释便于理解算法逻辑与快速调整实验条件。已有123人学习下载配套案例数据开箱即用无需额外配置可帮助学习者从理论推导、公式编程到仿真验证全流程掌握全周傅氏算法的核心思想与实践要点。1. 全周傅氏算法到底解决什么问题打开这个压缩包的时候我猜你大概率是在做微机继电保护相关的课设、毕设或者真的需要在一台保护装置上算基波。全周傅氏算法在电力系统里算得上是元老级算法了线路保护、变压器差动保护、母差保护、发电机保护只要涉及从故障电流电压里提取基波幅值和相位基本都绕不开它。这篇博文不打算铺陈太多背景直接把算法本身拆开理论基础怎么来的、Matlab代码怎么落地、实际用起来会踩哪些坑一次讲清楚。先说结论全周傅氏算法本质上就是一个基于傅里叶级数的基波提取器它用一个完整工频周期20毫秒的数据窗把故障信号里的基波分量“挑”出来同时把恒定直流和各次整数谐波统统挡在外面。对继电保护来说基波幅值决定过流保护动不动基波相位决定距离保护测的阻抗方向所以“提取基波”是整个保护逻辑的第一步也是最重要的一步。适合看这篇的人有两类一类是刚接触保护算法、需要快速理解原理并交作业的学生另一类是已经在做保护算法验证或者电能质量分析的工程师。前者可以重点看理论推导和代码后者可以直接跳到第四节以后的工程坑那里面的东西常规教材不会写。1.1 故障电流里的“成分”有多复杂电力系统正常运行时电流电压是标准的50Hz正弦波。但一旦发生短路故障波形就变脏了。最常见的故障暂态信号里至少有四类成分基波50Hz、整数次谐波主要是3次、5次变压器励磁涌流里2次谐波也很常见、非周期分量也就是衰减直流短路瞬间往往很大按指数衰减时间常数在几十毫秒量级以及各种高频噪声。保护装置要做的第一件事就是从这一锅粥里把基波分量单独盛出来。为什么偏偏是基波因为正常运行、故障时的稳态短路、系统振荡等工况下基波分量的大小和相位携带了最能反映系统状态的信息。比如短路故障后基波电流会突然增大而振荡时基波电流是缓慢波动。谐波和直流分量在这些判断里不仅没用还会干扰保护逻辑。比如变压器励磁涌流里2次谐波含量高微机差动保护恰恰是靠识别2次谐波来制动的所以谐波提取本身也是个核心技术但那是另一个话题。这就引出一个关键问题用什么办法能从混杂信号里干净地取出基波傅里叶变换当然是标准答案但保护装置里不能按教科书那样做全局变换它需要快速、递归、可嵌入微处理器计算的离散实现。于是全周傅氏算法就成了最合适的方案。1.2 为什么偏偏是“全周”工程上有几种从故障信号里提取基波的路线。最传统的是模拟窄带滤波器用RLC网络把基波附近的频率放行、其他频率衰减。这种方案在早期模拟式保护里用得很多缺点是元器件温漂、老化后中心频率会偏而且模拟滤波器对衰减直流几乎没有抑制作用相位特性也不好控制。数字化以后大家几乎都转向了基于离散傅里叶变换的算法这就是全周傅氏算法登场的背景。全周傅氏算法全称是“全周波傅里叶算法”数据窗长度是一个完整的工频周期。与之对应的是半周波傅里叶算法只取半个周期10毫秒。两者对比很有意思半周算法速度快故障后10毫秒就能出结果适合对速动性要求极高的场合但半个周期内的数据无法区分偶次谐波对恒定直流也完全无能为力抗干扰能力明显不足。全周算法虽然要等满20毫秒才能给出第一个计算结果但它在原理上能完全滤除所有整数次谐波和恒定直流抗衰减直流的能力也更强。对大多数保护而言20毫秒的等待完全可接受所以全周算法成了微机保护的主流选择。还有一个优势经常被忽略全周傅氏算法不仅能算出基波幅值还能同步算出基波相位。相位这个信息太重要了距离保护要靠电压电流的相位差算测量阻抗方向保护要靠相位判断故障方向。如果用普通带通滤波器相位延迟是频率相关的需要额外校准而全周傅氏算法直接输出正交分解后的实部和虚部幅值相位一次拿到手干净利落。2. 理论基础五分钟从傅里叶推到这里2.1 核心武器正交性要理解全周傅氏算法不需要把傅里叶变换的所有性质都背下来只需要抓住一个核心机制三角函数的正交性。简单说在一个周期内不同频率的cos和cos相乘再积分结果是零只有同频率的cos和cos相乘积分才不为零。这句话往深了想其实就是一个筛子我用基波频率的cos去和信号做内积信号里所有非基波频率的成分都会被筛掉只剩基波分量里与cos同相的那部分再用基波频率的sin去筛一遍就得到基波分量里与sin同相的那部分。把这两个值组合起来基波幅值和相位就都出来了。打个比方这就好比你在一群身高各异的人里面找特定身高的人你用一个正好卡在那个高度的板子去比划能通过的就是你要找的其他都被挡在外面。三角函数的正交性就是这个“正好卡住”的板子。用数学表达设信号x(t)是一个周期为T的周期信号角频率ω 2π/T。把它展开成傅里叶级数x(t) X0 Σ[Xk·cos(kωt) Yk·sin(kωt)]k从1到无穷其中X0是直流分量Xk和Yk是第k次谐波的余弦项和正弦项系数。根据正交性这些系数可以这么求Xk (2/T)·∫x(t)·cos(kωt)dt积分区间是一个周期 Yk (2/T)·∫x(t)·sin(kωt)dt这里的“2/T”是归一化系数因为同频cos和cos在一个周期内的积分是T/2乘以2/T后正好变成1这样提取出来的就是真实的幅值。2.2 离散化全过程保护装置里没有连续的积分器只有经过采样保持和A/D转换得到的离散序列。所以要把上面的连续积分改成离散求和。设每个工频周期采样N个点采样周期Ts T/N采样时刻t n·Tsn 0,1,...,N-1。用矩形积分代替连续积分就得到全周傅氏算法的基本公式Xr (2/N)·Σ x(n)·cos(2πn/N) Xi -(2/N)·Σ x(n)·sin(2πn/N)这里n从0到N-1一共N个采样点。注意第二个式子里我加了个负号这是很多教材里容易绕晕的地方。负号的作用是让最后计算相位时数学表达更简洁如果输入信号是x(t) A·cos(ωt φ)那么推导出来Xr A·cosφXi A·sinφ这样相位φ atan2(Xi, Xr)一步到位不用再手动调象限。如果不用负号公式变成Xi A·sinφ那边就会差个负号相位算出来要自己再补。有了Xr和Xi基波幅值和相位就是A sqrt(Xr² Xi²) φ atan2(Xi, Xr)这套算法之所以能把恒定直流完全滤掉是因为直流分量在求和时cos(2πn/N)和sin(2πn/N)在一个完整周期内的和正好是零。同理对于任意整数次谐波k它与基频cos、sin的内积在整周期内也恒为零所以全周算法天然免疫整数次谐波。这个性质在原理上就被保证了不需要任何附加滤波器。2.3 固定窗长和下标的细节实现时有个容易被忽视的细节离散化公式里的角度是2πn/N而不是2πf0·n/fs。两者在理论上等价因为N fs/f0。但在代码实现层面直接用一个N点长度的系数表最省事也最不容易出错。先把cos(2πn/N)和sin(2πn/N)的N个值算好存成表然后每次计算就是查表乘累加速度快到能在低端DSP上跑实时计算。另一个细节是起点的选择。如果实际系统有一定相位偏移公式依然成立因为正交性是平移不变的。也就是说不管信号是cos(ωt)还是cos(ωt φ)只要数据窗恰好是一个整周期算法输出的Xr和Xi就分别是A·cosφ和A·sinφ误差为零。但前提是数据窗长度必须严格等于一个工频周期。这句“前提”在工程里非常关键后面我会专门讲。3. Matlab代码直接能跑的两个版本3.1 测试信号构造写算法之前先造一个测试信号。真实故障电流没有教科书那么干净所以我构造的信号里包含了基波、3次谐波、5次谐波、恒定直流和衰减直流模拟一个比较典型的短路暂态电流%% 参数设置 fs 2000; % 采样率 2000 Hz f0 50; % 基波频率 50 Hz N fs / f0; % 一个周期采样点数 40 t (0:4*N-1) / fs; % 仿真4个工频周期 %% 故障电流基波50A(30度) 3次谐波5A 5次谐波2A 恒定直流1A 衰减直流10A(τ50ms) x 50*cos(2*pi*f0*t pi/6) ... 5*cos(2*pi*3*f0*t) ... 2*cos(2*pi*5*f0*t pi/4) ... 1 ... 10*exp(-t/0.05);这里N 2000/50 40正好是整数所以“一个周期取40个点”数据窗长度就是40个采样点。实际装置里N不一定取40常见的有12、20、24、32、48、64取决于采样率和保护采样规格。N越大高频抗扰能力越强但数据窗对应的时间不变都是20毫秒。3.2 基本版算法实现从公式到代码非常直接。我把算法封装成一个函数fdfull输入采样序列、采样率和基波频率输出基波幅值和相位function [A, phi] fdfull(x, fs, f0) % 全周傅氏算法提取基波分量幅值和相位 % 输入 % x - 采样序列长度至少为一个工频周期 % fs - 采样率 % f0 - 基波频率 % 输出 % A - 基波幅值 % phi - 基波相位弧度 N round(fs / f0); x x(:); if length(x) N error(输入数据长度不足一个工频周期); end n (0:N-1); Xr 2/N * sum(x(1:N) .* cos(2*pi*n/N)); Xi -2/N * sum(x(1:N) .* sin(2*pi*n/N)); A sqrt(Xr^2 Xi^2); phi atan2(Xi, Xr); end调用方式很简单[A, phi] fdfull(x, fs, f0); fprintf(基波幅值: %.4f A\n, A); fprintf(基波相位: %.4f rad (理论值 %.4f rad)\n, phi, pi/6);我在仿真里用了4个周期的数据但基本版只取前40个点也就是第一个周期的数据就算出了结果。基波理论幅值是50A相位是π/6≈0.5236弧度。没有衰减直流干扰的理想情况下输出会精确到浮点数精度。3.3 递推版算法实现上面的基本版有个局限它必须等到一个完整数据窗收集完毕才能计算一次而且每次计算都要重新累加40个点。在保护装置里我们希望每来一个新采样点就能更新一次基波计算结果这样故障后20毫秒窗口滑动过程中我们能连续观察幅值的变化轨迹。这就需要递推形式也叫滑动窗全周傅氏算法。递推的核心思想很简单新窗口比旧窗口多了一个新采样点、少了一个旧采样点所以实部和虚部都可以在旧值基础上做增量更新。关键是对应关系旧窗口的数据点是x(1)到x(N)新窗口是x(2)到x(N1)。计算时新窗口 旧窗口 x(N1)的贡献 - x(1)的贡献。由于cos表是周期性的x(N1)对应的角度和x(1)对应的角度相差2π在三角函数表里的索引是同一个。这就是为什么能用一个固定的N点系数表反复循环使用。代码实现如下function [A, phi] fdfull_recursive(x, fs, f0) % 递推全周傅氏算法每来一个新采样点输出一次基波 % 输入 % x - 采样序列长度至少为 N1 % fs - 采样率 % f0 - 基波频率 % 输出 % A - 每个滑窗位置的基波幅值 % phi - 每个滑窗位置的基波相位弧度 N round(fs / f0); x x(:); L length(x); if L N 1 error(数据长度需要至少 N1 个点); end c cos(2*pi*(0:N-1)/N); s sin(2*pi*(0:N-1)/N); % 用第一个窗口初始化 Xr 2/N * sum(x(1:N) .* c); Xi -2/N * sum(x(1:N) .* s); A zeros(L-N1, 1); phi zeros(L-N1, 1); A(1) sqrt(Xr^2 Xi^2); phi(1) atan2(Xi, Xr); for k 2 : L-N1 % 旧点 x(k-1) 和新点 x(kN-1) 对应的系数表索引 old_idx mod(k-2, N) 1; new_idx mod(kN-2, N) 1; Xr Xr 2/N * (x(kN-1)*c(new_idx) - x(k-1)*c(old_idx)); Xi Xi - 2/N * (x(kN-1)*s(new_idx) - x(k-1)*s(old_idx)); A(k) sqrt(Xr^2 Xi^2); phi(k) atan2(Xi, Xr); end end这段代码里有几个细节要强调第一old_idx和new_idx的下标计算。Matlab的索引从1开始而公式里的n从0开始所以取模之后要加1。old_idx mod(k-2, N) 1new_idx mod(kN-2, N) 1这两个式子不好背但逻辑很清晰信号点x(m)在Matlab里的位置是m对应的采样序号是m-1角度是2π(m-1)/N在系数表里的下标就是mod(m-1, N)1。把m分别换成k-1和kN-1就是上面两个式子。写代码的时候实在不确定就先用基本版比对结果确保递推版和基本版输出一致。第二递推版输出的第一个值A(1)和phi(1)和基本版完全一样因为都是用第一个窗口的数据算的。后续每一个新采样点都会推一个新的输出。这就是保护装置里常用的“滑窗更新”思路。4. 仿真验证算法到底有多能打4.1 谐波和恒定直流完全滤除先用不含衰减直流的信号做一组验证。把测试信号改成基波50A30度加3次谐波5A加5次谐波2A加恒定直流1A理论上全周算法应该精确输出50A、0.5236弧度。运行基本版代码输出结果在浮点精度内与理论值一致误差在10的负15次方量级。这不是巧合而是正交性的必然结果整数次谐波和恒定直流在整周期窗内完全被cos和sin的累积抵消。这个特性在继电保护里非常值钱。比如区外故障时CT饱和会产生大量谐波如果算法不能干净地滤除谐波保护就可能误判。全周算法在这一点上的表现几乎无可挑剔。但注意它只能滤除整数次谐波和恒定直流如果信号里有间谐波比如变频器输出的非整数倍频谱波全周算法也无能为力因为间谐波与基频在整周期内不正交。这个局限性在某些电能质量分析场景下需要注意。4.2 衰减直流测试把前面构造的完整故障信号含衰减直流送入算法结果就不是完美的50A了。用起始幅值10A、时间常数50ms的衰减直流量做测试算法输出的基波幅值大约在49.6A上下相位也有约1度左右的偏差。这个偏差看起来不大但如果衰减直流起始幅值更大比如短路发生在电压过零点附近时衰减直流分量可能达到上百安培此时对基波提取的影响就不能忽略了。衰减直流为什么治不了从频域看衰减直流不是一条在零频的线而是一个连续频谱它在50Hz处有泄漏分量。一个整周期矩形窗的正交性只能保证对恒定直流完全滤除对非周期的衰减直流残余误差取决于衰减时间常数、初相和幅值。时间常数越短衰减越快误差越大数据窗越靠近故障发生时刻误差也越大。工程上解决这个问题的常见手段是前置差分滤波器。在做全周傅氏之前先对采样序列做一次一阶差分差分会大大衰减低频直流分量然后再做傅氏提取。代价是差分会改变信号的幅值相位需要做系数补偿。这个方案简单有效在很多保护装置里都有应用。4.3 频率偏移测试电力系统的额定频率是50Hz但实际运行中会有偏差国标允许49.5到50.5Hz。如果系统频率漂到49.5Hz而算法仍然按50Hz的周期取40个点数据窗就不是完整的信号周期了正交性被破坏输出会出现误差。从仿真看频率偏差0.5Hz时幅值误差一般还在几个百分点以内但相位误差会随着滑窗位置变化呈周期性波动这一点在距离保护里可能引发测距误差。要缓解频率偏移的影响工程上有两个思路一是加频率跟踪先用测频算法估计实时频率动态调整采样率或重采样保证数据窗始终对应一个完整周期二是用更长数据窗、加窗函数来处理频谱泄漏但保护实时性要求下用长窗不现实。多数保护装置采用第一种思路配合锁相环或软件测频来锁定真实频率。5. 工程实战这些坑我帮你踩过了5.1 采样率必须和工频严格匹配这是最常踩的坑没有之一。全周傅氏算法理论上要求数据窗长度严格等于一个工频周期也就是N fs/f0必须是整数。如果N不是整数比如fs 1000Hz、f0 50Hz倒好说N20但有些装置采样率是1200Hz除以50也是24也没问题。真正出问题的是系统频率偏移比如实际49.8HzN算出来不是整数代码里用了round取整结果数据窗和信号周期相差一点点长期运行下来累积误差会很明显。我自己的经验是做MATLAB验证时可以在函数入口加一个检查如果fs/f0不是整数就提示用户确认是否要重采样或者改用频率跟踪。虽然加了这行代码会让函数显得啰嗦但它能拦住不少低级错误。5.2 递推版本的边界处理递推版代码里最容易出现索引越界或结果突跳。我在调试早期版本时就遇到过一个问题滑窗滑到数据末尾时old_idx和new_idx的取模运算一旦写错结果会突然跳到错误值。排查方法很简单把递推版和基本版在相同输入下的输出画在一起对比如果两条曲线重合说明递推逻辑没问题中间某个点突然分叉就去查那个点对应的索引对不对。另外注意递推版不能一上来就出结果。它的第一个有效输出需要对前N个采样点做完整初始化计算所以在实时采样系统里程序启动后要先攒够一个周期的数据再进入递推循环。这个“启动暂态”需要配套的状态机逻辑别在代码里省掉判断。5.3 算法响应速度和灵敏度的取舍全周算法的最大短板是固有延迟故障发生后要满20ms才能拿到第一个准确的基波值。这在某些超高速保护场景下是不够的。我做过一次对比半周算法故障后10ms就能出结果但谐波抑制能力差而且对直流更敏感。工程上常见的折中方案是双算法并行启动元件用半周算法或瞬时值快速动作确认故障后再切换到全周算法做精确测量。这个套路在数字式线路保护里很常见理解全周算法的性能边界才能更好地理解整个保护方案为什么这么设计。5.4 常见问题速查表现象可能原因解决办法输出幅值偏小/偏大采样率与信号频率不匹配N取整引入周期误差检查fs/f0是否为整数必要时重采样或加频率跟踪相位结果反复跳动频率偏移导致数据窗非整周期用锁相环测频、动态调整采样本文还有配套的精品资源点击获取
全周傅氏算法：从原理到Matlab实现的基波提取指南
简介本资源系统讲解全周傅氏算法的数学原理与工程实现面向电子信息工程、计算机及数学专业本科生适用于课程设计、期末大作业及毕业设计中对电力系统谐波分析、信号频谱提取等典型场景的需求。压缩包共3个文件862KB含1份详尽PDF理论文档、1个Simulink仿真模型.mdl用于算法动态验证、1个结构清晰的MATLAB主程序.m支持matlab2014a至2024a多版本直接运行。代码采用参数化设计关键变量如采样频率、信号周期、谐波阶数等均集中定义并配有中文注释便于理解算法逻辑与快速调整实验条件。已有123人学习下载配套案例数据开箱即用无需额外配置可帮助学习者从理论推导、公式编程到仿真验证全流程掌握全周傅氏算法的核心思想与实践要点。1. 全周傅氏算法到底解决什么问题打开这个压缩包的时候我猜你大概率是在做微机继电保护相关的课设、毕设或者真的需要在一台保护装置上算基波。全周傅氏算法在电力系统里算得上是元老级算法了线路保护、变压器差动保护、母差保护、发电机保护只要涉及从故障电流电压里提取基波幅值和相位基本都绕不开它。这篇博文不打算铺陈太多背景直接把算法本身拆开理论基础怎么来的、Matlab代码怎么落地、实际用起来会踩哪些坑一次讲清楚。先说结论全周傅氏算法本质上就是一个基于傅里叶级数的基波提取器它用一个完整工频周期20毫秒的数据窗把故障信号里的基波分量“挑”出来同时把恒定直流和各次整数谐波统统挡在外面。对继电保护来说基波幅值决定过流保护动不动基波相位决定距离保护测的阻抗方向所以“提取基波”是整个保护逻辑的第一步也是最重要的一步。适合看这篇的人有两类一类是刚接触保护算法、需要快速理解原理并交作业的学生另一类是已经在做保护算法验证或者电能质量分析的工程师。前者可以重点看理论推导和代码后者可以直接跳到第四节以后的工程坑那里面的东西常规教材不会写。1.1 故障电流里的“成分”有多复杂电力系统正常运行时电流电压是标准的50Hz正弦波。但一旦发生短路故障波形就变脏了。最常见的故障暂态信号里至少有四类成分基波50Hz、整数次谐波主要是3次、5次变压器励磁涌流里2次谐波也很常见、非周期分量也就是衰减直流短路瞬间往往很大按指数衰减时间常数在几十毫秒量级以及各种高频噪声。保护装置要做的第一件事就是从这一锅粥里把基波分量单独盛出来。为什么偏偏是基波因为正常运行、故障时的稳态短路、系统振荡等工况下基波分量的大小和相位携带了最能反映系统状态的信息。比如短路故障后基波电流会突然增大而振荡时基波电流是缓慢波动。谐波和直流分量在这些判断里不仅没用还会干扰保护逻辑。比如变压器励磁涌流里2次谐波含量高微机差动保护恰恰是靠识别2次谐波来制动的所以谐波提取本身也是个核心技术但那是另一个话题。这就引出一个关键问题用什么办法能从混杂信号里干净地取出基波傅里叶变换当然是标准答案但保护装置里不能按教科书那样做全局变换它需要快速、递归、可嵌入微处理器计算的离散实现。于是全周傅氏算法就成了最合适的方案。1.2 为什么偏偏是“全周”工程上有几种从故障信号里提取基波的路线。最传统的是模拟窄带滤波器用RLC网络把基波附近的频率放行、其他频率衰减。这种方案在早期模拟式保护里用得很多缺点是元器件温漂、老化后中心频率会偏而且模拟滤波器对衰减直流几乎没有抑制作用相位特性也不好控制。数字化以后大家几乎都转向了基于离散傅里叶变换的算法这就是全周傅氏算法登场的背景。全周傅氏算法全称是“全周波傅里叶算法”数据窗长度是一个完整的工频周期。与之对应的是半周波傅里叶算法只取半个周期10毫秒。两者对比很有意思半周算法速度快故障后10毫秒就能出结果适合对速动性要求极高的场合但半个周期内的数据无法区分偶次谐波对恒定直流也完全无能为力抗干扰能力明显不足。全周算法虽然要等满20毫秒才能给出第一个计算结果但它在原理上能完全滤除所有整数次谐波和恒定直流抗衰减直流的能力也更强。对大多数保护而言20毫秒的等待完全可接受所以全周算法成了微机保护的主流选择。还有一个优势经常被忽略全周傅氏算法不仅能算出基波幅值还能同步算出基波相位。相位这个信息太重要了距离保护要靠电压电流的相位差算测量阻抗方向保护要靠相位判断故障方向。如果用普通带通滤波器相位延迟是频率相关的需要额外校准而全周傅氏算法直接输出正交分解后的实部和虚部幅值相位一次拿到手干净利落。2. 理论基础五分钟从傅里叶推到这里2.1 核心武器正交性要理解全周傅氏算法不需要把傅里叶变换的所有性质都背下来只需要抓住一个核心机制三角函数的正交性。简单说在一个周期内不同频率的cos和cos相乘再积分结果是零只有同频率的cos和cos相乘积分才不为零。这句话往深了想其实就是一个筛子我用基波频率的cos去和信号做内积信号里所有非基波频率的成分都会被筛掉只剩基波分量里与cos同相的那部分再用基波频率的sin去筛一遍就得到基波分量里与sin同相的那部分。把这两个值组合起来基波幅值和相位就都出来了。打个比方这就好比你在一群身高各异的人里面找特定身高的人你用一个正好卡在那个高度的板子去比划能通过的就是你要找的其他都被挡在外面。三角函数的正交性就是这个“正好卡住”的板子。用数学表达设信号x(t)是一个周期为T的周期信号角频率ω 2π/T。把它展开成傅里叶级数x(t) X0 Σ[Xk·cos(kωt) Yk·sin(kωt)]k从1到无穷其中X0是直流分量Xk和Yk是第k次谐波的余弦项和正弦项系数。根据正交性这些系数可以这么求Xk (2/T)·∫x(t)·cos(kωt)dt积分区间是一个周期 Yk (2/T)·∫x(t)·sin(kωt)dt这里的“2/T”是归一化系数因为同频cos和cos在一个周期内的积分是T/2乘以2/T后正好变成1这样提取出来的就是真实的幅值。2.2 离散化全过程保护装置里没有连续的积分器只有经过采样保持和A/D转换得到的离散序列。所以要把上面的连续积分改成离散求和。设每个工频周期采样N个点采样周期Ts T/N采样时刻t n·Tsn 0,1,...,N-1。用矩形积分代替连续积分就得到全周傅氏算法的基本公式Xr (2/N)·Σ x(n)·cos(2πn/N) Xi -(2/N)·Σ x(n)·sin(2πn/N)这里n从0到N-1一共N个采样点。注意第二个式子里我加了个负号这是很多教材里容易绕晕的地方。负号的作用是让最后计算相位时数学表达更简洁如果输入信号是x(t) A·cos(ωt φ)那么推导出来Xr A·cosφXi A·sinφ这样相位φ atan2(Xi, Xr)一步到位不用再手动调象限。如果不用负号公式变成Xi A·sinφ那边就会差个负号相位算出来要自己再补。有了Xr和Xi基波幅值和相位就是A sqrt(Xr² Xi²) φ atan2(Xi, Xr)这套算法之所以能把恒定直流完全滤掉是因为直流分量在求和时cos(2πn/N)和sin(2πn/N)在一个完整周期内的和正好是零。同理对于任意整数次谐波k它与基频cos、sin的内积在整周期内也恒为零所以全周算法天然免疫整数次谐波。这个性质在原理上就被保证了不需要任何附加滤波器。2.3 固定窗长和下标的细节实现时有个容易被忽视的细节离散化公式里的角度是2πn/N而不是2πf0·n/fs。两者在理论上等价因为N fs/f0。但在代码实现层面直接用一个N点长度的系数表最省事也最不容易出错。先把cos(2πn/N)和sin(2πn/N)的N个值算好存成表然后每次计算就是查表乘累加速度快到能在低端DSP上跑实时计算。另一个细节是起点的选择。如果实际系统有一定相位偏移公式依然成立因为正交性是平移不变的。也就是说不管信号是cos(ωt)还是cos(ωt φ)只要数据窗恰好是一个整周期算法输出的Xr和Xi就分别是A·cosφ和A·sinφ误差为零。但前提是数据窗长度必须严格等于一个工频周期。这句“前提”在工程里非常关键后面我会专门讲。3. Matlab代码直接能跑的两个版本3.1 测试信号构造写算法之前先造一个测试信号。真实故障电流没有教科书那么干净所以我构造的信号里包含了基波、3次谐波、5次谐波、恒定直流和衰减直流模拟一个比较典型的短路暂态电流%% 参数设置 fs 2000; % 采样率 2000 Hz f0 50; % 基波频率 50 Hz N fs / f0; % 一个周期采样点数 40 t (0:4*N-1) / fs; % 仿真4个工频周期 %% 故障电流基波50A(30度) 3次谐波5A 5次谐波2A 恒定直流1A 衰减直流10A(τ50ms) x 50*cos(2*pi*f0*t pi/6) ... 5*cos(2*pi*3*f0*t) ... 2*cos(2*pi*5*f0*t pi/4) ... 1 ... 10*exp(-t/0.05);这里N 2000/50 40正好是整数所以“一个周期取40个点”数据窗长度就是40个采样点。实际装置里N不一定取40常见的有12、20、24、32、48、64取决于采样率和保护采样规格。N越大高频抗扰能力越强但数据窗对应的时间不变都是20毫秒。3.2 基本版算法实现从公式到代码非常直接。我把算法封装成一个函数fdfull输入采样序列、采样率和基波频率输出基波幅值和相位function [A, phi] fdfull(x, fs, f0) % 全周傅氏算法提取基波分量幅值和相位 % 输入 % x - 采样序列长度至少为一个工频周期 % fs - 采样率 % f0 - 基波频率 % 输出 % A - 基波幅值 % phi - 基波相位弧度 N round(fs / f0); x x(:); if length(x) N error(输入数据长度不足一个工频周期); end n (0:N-1); Xr 2/N * sum(x(1:N) .* cos(2*pi*n/N)); Xi -2/N * sum(x(1:N) .* sin(2*pi*n/N)); A sqrt(Xr^2 Xi^2); phi atan2(Xi, Xr); end调用方式很简单[A, phi] fdfull(x, fs, f0); fprintf(基波幅值: %.4f A\n, A); fprintf(基波相位: %.4f rad (理论值 %.4f rad)\n, phi, pi/6);我在仿真里用了4个周期的数据但基本版只取前40个点也就是第一个周期的数据就算出了结果。基波理论幅值是50A相位是π/6≈0.5236弧度。没有衰减直流干扰的理想情况下输出会精确到浮点数精度。3.3 递推版算法实现上面的基本版有个局限它必须等到一个完整数据窗收集完毕才能计算一次而且每次计算都要重新累加40个点。在保护装置里我们希望每来一个新采样点就能更新一次基波计算结果这样故障后20毫秒窗口滑动过程中我们能连续观察幅值的变化轨迹。这就需要递推形式也叫滑动窗全周傅氏算法。递推的核心思想很简单新窗口比旧窗口多了一个新采样点、少了一个旧采样点所以实部和虚部都可以在旧值基础上做增量更新。关键是对应关系旧窗口的数据点是x(1)到x(N)新窗口是x(2)到x(N1)。计算时新窗口 旧窗口 x(N1)的贡献 - x(1)的贡献。由于cos表是周期性的x(N1)对应的角度和x(1)对应的角度相差2π在三角函数表里的索引是同一个。这就是为什么能用一个固定的N点系数表反复循环使用。代码实现如下function [A, phi] fdfull_recursive(x, fs, f0) % 递推全周傅氏算法每来一个新采样点输出一次基波 % 输入 % x - 采样序列长度至少为 N1 % fs - 采样率 % f0 - 基波频率 % 输出 % A - 每个滑窗位置的基波幅值 % phi - 每个滑窗位置的基波相位弧度 N round(fs / f0); x x(:); L length(x); if L N 1 error(数据长度需要至少 N1 个点); end c cos(2*pi*(0:N-1)/N); s sin(2*pi*(0:N-1)/N); % 用第一个窗口初始化 Xr 2/N * sum(x(1:N) .* c); Xi -2/N * sum(x(1:N) .* s); A zeros(L-N1, 1); phi zeros(L-N1, 1); A(1) sqrt(Xr^2 Xi^2); phi(1) atan2(Xi, Xr); for k 2 : L-N1 % 旧点 x(k-1) 和新点 x(kN-1) 对应的系数表索引 old_idx mod(k-2, N) 1; new_idx mod(kN-2, N) 1; Xr Xr 2/N * (x(kN-1)*c(new_idx) - x(k-1)*c(old_idx)); Xi Xi - 2/N * (x(kN-1)*s(new_idx) - x(k-1)*s(old_idx)); A(k) sqrt(Xr^2 Xi^2); phi(k) atan2(Xi, Xr); end end这段代码里有几个细节要强调第一old_idx和new_idx的下标计算。Matlab的索引从1开始而公式里的n从0开始所以取模之后要加1。old_idx mod(k-2, N) 1new_idx mod(kN-2, N) 1这两个式子不好背但逻辑很清晰信号点x(m)在Matlab里的位置是m对应的采样序号是m-1角度是2π(m-1)/N在系数表里的下标就是mod(m-1, N)1。把m分别换成k-1和kN-1就是上面两个式子。写代码的时候实在不确定就先用基本版比对结果确保递推版和基本版输出一致。第二递推版输出的第一个值A(1)和phi(1)和基本版完全一样因为都是用第一个窗口的数据算的。后续每一个新采样点都会推一个新的输出。这就是保护装置里常用的“滑窗更新”思路。4. 仿真验证算法到底有多能打4.1 谐波和恒定直流完全滤除先用不含衰减直流的信号做一组验证。把测试信号改成基波50A30度加3次谐波5A加5次谐波2A加恒定直流1A理论上全周算法应该精确输出50A、0.5236弧度。运行基本版代码输出结果在浮点精度内与理论值一致误差在10的负15次方量级。这不是巧合而是正交性的必然结果整数次谐波和恒定直流在整周期窗内完全被cos和sin的累积抵消。这个特性在继电保护里非常值钱。比如区外故障时CT饱和会产生大量谐波如果算法不能干净地滤除谐波保护就可能误判。全周算法在这一点上的表现几乎无可挑剔。但注意它只能滤除整数次谐波和恒定直流如果信号里有间谐波比如变频器输出的非整数倍频谱波全周算法也无能为力因为间谐波与基频在整周期内不正交。这个局限性在某些电能质量分析场景下需要注意。4.2 衰减直流测试把前面构造的完整故障信号含衰减直流送入算法结果就不是完美的50A了。用起始幅值10A、时间常数50ms的衰减直流量做测试算法输出的基波幅值大约在49.6A上下相位也有约1度左右的偏差。这个偏差看起来不大但如果衰减直流起始幅值更大比如短路发生在电压过零点附近时衰减直流分量可能达到上百安培此时对基波提取的影响就不能忽略了。衰减直流为什么治不了从频域看衰减直流不是一条在零频的线而是一个连续频谱它在50Hz处有泄漏分量。一个整周期矩形窗的正交性只能保证对恒定直流完全滤除对非周期的衰减直流残余误差取决于衰减时间常数、初相和幅值。时间常数越短衰减越快误差越大数据窗越靠近故障发生时刻误差也越大。工程上解决这个问题的常见手段是前置差分滤波器。在做全周傅氏之前先对采样序列做一次一阶差分差分会大大衰减低频直流分量然后再做傅氏提取。代价是差分会改变信号的幅值相位需要做系数补偿。这个方案简单有效在很多保护装置里都有应用。4.3 频率偏移测试电力系统的额定频率是50Hz但实际运行中会有偏差国标允许49.5到50.5Hz。如果系统频率漂到49.5Hz而算法仍然按50Hz的周期取40个点数据窗就不是完整的信号周期了正交性被破坏输出会出现误差。从仿真看频率偏差0.5Hz时幅值误差一般还在几个百分点以内但相位误差会随着滑窗位置变化呈周期性波动这一点在距离保护里可能引发测距误差。要缓解频率偏移的影响工程上有两个思路一是加频率跟踪先用测频算法估计实时频率动态调整采样率或重采样保证数据窗始终对应一个完整周期二是用更长数据窗、加窗函数来处理频谱泄漏但保护实时性要求下用长窗不现实。多数保护装置采用第一种思路配合锁相环或软件测频来锁定真实频率。5. 工程实战这些坑我帮你踩过了5.1 采样率必须和工频严格匹配这是最常踩的坑没有之一。全周傅氏算法理论上要求数据窗长度严格等于一个工频周期也就是N fs/f0必须是整数。如果N不是整数比如fs 1000Hz、f0 50Hz倒好说N20但有些装置采样率是1200Hz除以50也是24也没问题。真正出问题的是系统频率偏移比如实际49.8HzN算出来不是整数代码里用了round取整结果数据窗和信号周期相差一点点长期运行下来累积误差会很明显。我自己的经验是做MATLAB验证时可以在函数入口加一个检查如果fs/f0不是整数就提示用户确认是否要重采样或者改用频率跟踪。虽然加了这行代码会让函数显得啰嗦但它能拦住不少低级错误。5.2 递推版本的边界处理递推版代码里最容易出现索引越界或结果突跳。我在调试早期版本时就遇到过一个问题滑窗滑到数据末尾时old_idx和new_idx的取模运算一旦写错结果会突然跳到错误值。排查方法很简单把递推版和基本版在相同输入下的输出画在一起对比如果两条曲线重合说明递推逻辑没问题中间某个点突然分叉就去查那个点对应的索引对不对。另外注意递推版不能一上来就出结果。它的第一个有效输出需要对前N个采样点做完整初始化计算所以在实时采样系统里程序启动后要先攒够一个周期的数据再进入递推循环。这个“启动暂态”需要配套的状态机逻辑别在代码里省掉判断。5.3 算法响应速度和灵敏度的取舍全周算法的最大短板是固有延迟故障发生后要满20ms才能拿到第一个准确的基波值。这在某些超高速保护场景下是不够的。我做过一次对比半周算法故障后10ms就能出结果但谐波抑制能力差而且对直流更敏感。工程上常见的折中方案是双算法并行启动元件用半周算法或瞬时值快速动作确认故障后再切换到全周算法做精确测量。这个套路在数字式线路保护里很常见理解全周算法的性能边界才能更好地理解整个保护方案为什么这么设计。5.4 常见问题速查表现象可能原因解决办法输出幅值偏小/偏大采样率与信号频率不匹配N取整引入周期误差检查fs/f0是否为整数必要时重采样或加频率跟踪相位结果反复跳动频率偏移导致数据窗非整周期用锁相环测频、动态调整采样本文还有配套的精品资源点击获取
