STM32F4-Base/docs/算法/傅氏算法详解.md

从《项目详细分析-交流采样与有效值计算》中提取扩展，配图讲解傅氏算法（全周波傅里叶算法）的原理与工程实现。 工程背景：CLS100 DTU 保护测控装置，每周波 20 点采样（1000Hz @ 50Hz），保护有效值/相量由 GlobalUnit/CalFFTVecFull.c 全周傅氏计算。 公式为 LaTeX 排版，需 Markdown 预览（VS Code / Typora 等）渲染。

---

1. 一句话思想

把一周波的采样序列分别和"标准正弦、标准余弦"做内积（投影）：正弦分量投影到正弦上得实部，余弦分量投影到余弦上得虚部——用两根已知相位的"标杆"去量未知信号里基波的幅值和相位。数学上这就是单频 DFT（只计算 k=1 那一根谱线）。

图中含 20% 五次谐波和 1V 直流的畸变波形，经两次投影后仍精确还原出基波：幅值 57.70V、相角 30.0°——谐波和直流在投影中"自己抵消了"，这是傅氏算法的全部魔力，第 3 节用面积图解释为什么。

2. 数学推导

2.1 连续形式

设周期信号 $x(t)$，基波角频率 $\omega = 2\pi f$，其基波分量：

x_1(t) = A_1\sin(\omega t + \varphi_1) = a_1\sin(\omega t) + b_1\cos(\omega t)

在一个整周波 T 上：

a_1 = \frac{2}{T}\int_0^T x(t)\sin(\omega t)\,\mathrm{d}t
\qquad
b_1 = \frac{2}{T}\int_0^T x(t)\cos(\omega t)\,\mathrm{d}t

2.2 离散形式（工程实现）

一个周波采 N 个点 $x(0),\dots,x(N-1)$，积分 → 求和：

a_1 \approx \frac{2}{N}\sum_{k=0}^{N-1} x(k)\,\sin\!\left(\frac{2\pi k}{N}\right)
\qquad
b_1 \approx \frac{2}{N}\sum_{k=0}^{N-1} x(k)\,\cos\!\left(\frac{2\pi k}{N}\right)

幅值、相位、有效值：

A_1 = \sqrt{a_1^2 + b_1^2}
\qquad
\varphi_1 = \arctan\!\frac{b_1}{a_1}
\qquad
U = \frac{A_1}{\sqrt{2}}

2.3 本项目的写法：把 \sqrt{2} 藏进系数表

GIUnit/GLobalFunc.c L337 的系数表预乘了 $\sqrt{2}/N$：

g\_fFFTSin20[k] = \sin\!\left(\frac{2\pi k}{N}\right)\cdot\frac{\sqrt{2}}{N}

于是累加结果直接是有效值相量分量：

f_{Re} = \sum_{k=0}^{N-1} x(k)\,\sin\!\left(\frac{2\pi k}{N}\right)\frac{\sqrt{2}}{N} = \frac{a_1}{\sqrt{2}}
\qquad
f_{Im} = \frac{b_1}{\sqrt{2}}

f_{AmSQ} = \sqrt{f_{Re}^2 + f_{Im}^2} = \frac{A_1}{\sqrt{2}} = U
\quad (\text{直接就是有效值，无需再除}\ \sqrt{2})

3. 为什么能滤除直流和谐波 —— 正交性

整周波上三角函数系正交：

\int_0^T \sin(\omega t)\,\mathrm{d}t = 0,
\qquad
\int_0^T \sin(m\omega t)\sin(\omega t)\,\mathrm{d}t = 0 \quad (m \in \mathbb{Z},\ m \neq 1)

看图最直观：蓝=正面积，红=负面积。

信号成分 × sinθ	面积结果	结论
直流 1V	红蓝完全抵消	完全滤除
3、5、…次整数谐波	红蓝完全抵消	完全滤除
基波	正面积占优，净值 ≠ 0	被提取
衰减直流（非严格周期）	近似抵消，有残余	大部分滤除
非整数次谐波（如 2.7 次）	不严格抵消	按"泄漏"比例部分反映

前提是数据窗恰好一个整周波。窗长偏离整周波（系统频率漂移），正交性被破坏，直流和 2 次谐波就会漏进测量——这就是本工程用"软件测频 + 片内 EHRPWM 调频"严格锁定 20 点/周波的原因。

4. 频域视角：一个梳状带通滤波器

把整周数据窗对基波做相关，等效为滤波器，其频率响应：

|H(f)| = \left|\frac{\sin\!\big(\pi (f/f_1 - 1)\big)}{N\,\sin\!\big(\pi (f/f_1 - 1)/N\big)}\right|

• 中心 50Hz 通带，增益 1；
• 零点恰好落在 0Hz（直流）和 100、150、200…Hz（全部整数谐波）；
• 窗越长（周波数越多）主瓣越窄、选择性越好，但时延越大。

"一窗 = 一个整周波"是保护领域在速度与滤波之间的经典折中，IEC 标准对保护动作时间的考核也按全周波算法进行。

5. 滑动窗：每 1ms 一个新有效值

数据窗定义为"以当前时刻为终止的前一周波"。每来一个新采样点（1ms），窗口滑动一格、丢弃最老点、20 点全部重算（"暴力法"）：

• 暴力法：每点重算 20 次乘加 ×2（sin/cos），实现简单、无累积误差；
• 迭代法（CalIterFFTVecFull.c，本工程未启用）：增量更新

S_{\text{new}} = S_{\text{old}} - x_{\text{old}}\cdot\sin\!\left(\frac{2\pi k}{N}\right) + x_{\text{new}}\cdot\sin\!\left(\frac{2\pi k}{N}\right)

每次只 2 次乘加，但浮点舍入误差随时间累积，需定期用暴力法清洗。本工程 CPU 余量足够，选暴力法。

6. 数值算例：三种算法各测到什么

同一组 20 点采样（57.7V 基波 + 20% 五次谐波 + 1V 直流）：

算法	读数	说明
直流平均	1.00 V	只剩直流
真有效值 \sqrt{\overline{x^2}}	≈58.85 V	$\sqrt{57.7^2 + 11.5^2 + 1^2}$，谐波直流全算进去
全周傅氏	57.70 V ∠30°	精确回到基波，且给出相位

保护要的是"基波有效值 + 相位"（方向元件、同期、功率都要相角），且故障时短路电流含衰减直流——傅氏一个整周波积分把它抹掉，真 RMS 则会虚高导致误动。这就是保护偏爱傅氏、计量偏爱真 RMS 的根本原因（方法全景见主文档附录C）。

7. 工程实现：本工程数据流

核心循环只有 6 行（GlobalUnit/CalFFTVecFull.c L392-405）：

for( dwDot=0; dwDot<dwCalRate; dwDot++) {
    fVal  = pfCalUse[dwAdr];              // 第n个采样点(已乘通道系数)
    fRe  += fVal * pfFFTSin[dwFFTCnt];    // Σ x·sin·(√2/N)
    fIm  += fVal * pfFFTCos[dwFFTCnt];    // Σ x·cos·(√2/N)
    dwAdr += dwLen;                       // 滑到下一个采样点
}
Fun_VectorAmp(ptVecSave);                 // fAmSQ = √(fRe²+fIm²) = 有效值

配套工程细节：

• 任意次谐波：系数索引步进量从 1 改为 k 即算 k 次谐波相量 $f_{Re}(k)=\sum x(n)\sin!\big(\tfrac{2\pi k n}{N}\big)\tfrac{\sqrt2}{N}$，同一循环复用于 2 次/5 次谐波（涌流判别、选线）；遥测对 k=1\sim11 全算再方和根 $U_{\text{rms}}=\sqrt{\sum_{k=1}^{11}U_k^2}$；
• 双区切换：结果写入向量区 tVectorRly[2][N] 的计算区，算完原子交换"计算区/使用区"指针，保护元件读使用区，无锁并发；
• 4 个实例（参考项目口径）：C01 基波 I/U、C02 零序基波、C03 二次谐波、C04 五次谐波（APP/CLS100_Dsp/InstCalCell.c）。本项目按 V0.1 有效通道展开为 5 个（C01-IA-F1/C02-UA-F1/C03-UN-F1/C04-IA-F2/C05-IA-F5），见 docs/算法/继电保护功能复现规划.md §3.2。

8. 全周波 vs 半周波

	全周傅氏（本工程用）	半周傅氏（代码有但未启用）
数据窗	1 个周波（20 点）	半个周波（10 点）
滤波能力	直流 + 全部整数谐波	半周相减消直流和偶次谐波，奇次保留
速度	1 周波	快一倍
精度	高	低（受奇次谐波影响）

本工程 CN_ANA_DIF=FALSE，半周波/差分组件全部注释掉——DTU 场景更看重精度和抗衰减直流能力。

9. 常见疑问

Q1：傅氏算的是有效值还是幅值？ 取决于系数表。标准系数 2/N 得幅值 $A_1$；本工程预乘 \sqrt2/N 直接得有效值 $U=A_1/\sqrt2$。看 GLobalFunc.c 系数表注释"折算成有效值"即可确认。

Q2：和 FFT 什么关系？ 傅氏算法 = 只算 $k=1$（或少数几次）谱线的 DFT。FFT 是同时算全部 N 根谱线的快速算法。保护只需要基波（+少数谐波），逐次 DFT 比 FFT 更省算力。

Q3：频率偏离 50Hz 还准吗？ 本工程靠闭环调频保证窗长恒等于实际周波（20 点/周波随频率微调），所以准。固定 1000Hz 采样不调整的方案，频偏 ±2Hz 时幅值误差约 1~2% 且出现相角偏差。

Q4：衰减直流滤不干净怎么办？ 整周波对"严格周期"信号才完全正交，衰减直流有残余。对策：加长数据窗（1.5~2 周波）、微分/余弦滤波改进傅氏、或最小二乘把 A\mathrm{e}^{-t/\tau} 作为未知量一并回归。