4. 傅里叶变换与频谱分析：离散傅里叶变换（DFT）、快速傅里叶变换（FFT）、频谱泄漏与窗函数
各位同学，咱们今天聊点硬核的——傅里叶变换。说实话，我刚入行那会儿，看到这一堆公式就头疼。但干我们微机保护这行的，这东西绕不开。你想想看，电力系统里的故障信号，说白了就是一堆正弦波的叠加。你要从波形里把基波、谐波都揪出来，靠的就是傅里叶变换这把“手术刀”。
4.1 从连续到离散：DFT的工程落地
教科书上的傅里叶变换是连续的，但咱们的微机保护装置是数字的，只能处理离散的采样点。所以就有了离散傅里叶变换，简称DFT。
DFT的公式长这样：
X(k) = Σ[n=0 to N-1] x(n) * e^(-j*2π*k*n/N)别被这个公式吓到。我给大家翻译一下：
- x(n)：你采集到的N个采样点，比如一个周波采64点
- X(k)：第k次谐波的幅值和相位
- e^(-j...)：说白了就是一组旋转因子，帮你把时域信号投影到频域
举个例子。假设你采集了一个工频50Hz的信号，每周波采64点。那么：
- k=1 算出来的是基波（50Hz）
- k=2 算出来的是2次谐波（100Hz）
- k=3 算出来的是3次谐波（150Hz）
- 以此类推，直到k=32（奈奎斯特频率，1600Hz）
关键点：DFT的计算量是O(N²)。N=64时还好，如果N=1024，计算量就上去了。我早期做的一个项目，CPU主频才几十兆，算一次DFT要好几毫秒，根本跑不动。所以就有了FFT。
4.2 快速傅里叶变换：把计算量降下来
FFT不是新算法，它是DFT的“快速实现版”。核心思想就四个字：分而治之。
怎么分？把N点DFT拆成两个N/2点的DFT，再拆成四个N/4点的DFT……一直拆到2点DFT。计算量从O(N²)降到了O(N log₂N)。
我给大家算笔账：
| N（采样点数） | DFT计算量 | FFT计算量 | 提速倍数 |
|---|---|---|---|
| 64 | 4096 | 384 | 约10倍 |
| 256 | 65536 | 2048 | 约32倍 |
| 1024 | 1048576 | 10240 | 约102倍 |
看到没？点数越多，FFT的优势越明显。我在做故障录波分析时，经常用1024点FFT，算一次也就几十微秒，实时性完全没问题。
常用的FFT实现是基2时间抽取法。代码实现也不复杂：
// 基2 FFT实现（简化版）
void fft(complex* x, int N) {
// 1. 位反转排序
for (int i = 0; i < N; i++) {
int j = bit_reverse(i, log2(N));
if (i < j) swap(x[i], x[j]);
}
// 2. 蝶形运算
for (int len = 2; len <= N; len <<= 1) {
complex w = exp(-j * 2 * PI / len);
for (int i = 0; i < N; i += len) {
complex wn = 1;
for (int k = 0; k < len/2; k++) {
complex u = x[i + k];
complex v = x[i + k + len/2] * wn;
x[i + k] = u + v;
x[i + k + len/2] = u - v;
wn *= w;
}
}
}
}我的经验：实际工程中，别自己手写FFT。用现成的库，比如CMSIS-DSP或者ARM的DSP库，它们针对特定CPU做了优化，比你自己写的快得多。我曾经在一个项目里手写FFT，结果算出来结果总是不对，折腾了两天才发现是位反转排序写错了。
4.3 频谱泄漏：一个让人头疼的问题
好，现在你学会了用FFT算频谱。但你会发现一个奇怪的现象：明明输入的是50Hz的纯净正弦波，FFT算出来却在49Hz和51Hz也有能量。这就是频谱泄漏。
为什么会这样？说白了，FFT要求你采样的信号是周期性的，且采样长度正好是信号周期的整数倍。但实际中，你很难做到这一点。比如你采了64点，但信号周期是64.3个采样点，那截断后的信号在边界处就不连续了。
这种不连续性，在频域里就表现为能量“泄漏”到了旁边的频率上。
注意：频谱泄漏对微机保护影响很大。比如你算基波幅值，泄漏会导致计算结果波动。我遇到过最严重的情况，泄漏误差达到了5%，直接导致保护误动。嗯，那次之后我就再也不敢忽视窗函数了。
4.4 窗函数：给信号“加个窗户”
解决频谱泄漏的办法，就是给信号加窗。窗函数的作用，是让信号在边界处平滑地衰减到0，减少截断带来的不连续性。
常用的窗函数有几种：
| 窗函数 | 主瓣宽度 | 旁瓣衰减 | 适用场景 |
|---|---|---|---|
| 矩形窗 | 最窄 | -13dB | 频率分辨率要求高，泄漏不敏感 |
| 汉宁窗 | 较宽 | -31dB | 通用场景，兼顾分辨率和泄漏 |
| 海明窗 | 较宽 | -41dB | 旁瓣衰减要求高 |
| 布莱克曼窗 | 最宽 | -57dB | 对泄漏要求极严 |
加窗的数学操作很简单：
// 加汉宁窗
for (int n = 0; n < N; n++) {
window[n] = 0.5 * (1 - cos(2*PI*n/(N-1)));
x_windowed[n] = x[n] * window[n];
}然后对加窗后的信号做FFT就行了。
我的建议：做微机保护，我一般用汉宁窗。它主瓣宽度适中，旁瓣衰减也够用。如果你对谐波分析要求特别高，比如要做5次、7次谐波的精确测量，可以考虑布莱克曼窗。但要注意，窗函数会降低频率分辨率，主瓣越宽，两个相近频率的信号就越难区分。
4.5 实战中的取舍
最后，我给大家总结一下实际工程中的几个要点：
- 采样率选择：至少是最高分析频率的2倍。比如你要分析到20次谐波（1000Hz），采样率至少2000Hz。我一般留3-5倍裕量。
- FFT点数：选2的整数次幂，比如64、128、256。点数越多，频率分辨率越高，但计算量也大。保护装置里常用64点或128点。
- 窗函数选择：保护算法里，汉宁窗是首选。如果做故障测距，对频率精度要求高，可以用矩形窗。
- 实时性：FFT计算要放在中断里还是主循环里？我建议放在主循环里，因为FFT计算时间不确定，放中断里容易出问题。
避坑指南：我曾经在一个项目里，FFT算出来的基波幅值总是偏小。查了半天，发现是采样时钟有抖动，导致采样间隔不均匀。后来加了硬件锁相环，问题才解决。所以，FFT的结果好不好，不光看算法，还得看前端采样电路稳不稳。
好了，这一章的内容就到这儿。傅里叶变换这东西，光看理论是不够的，你得动手写代码、跑数据、看波形。我建议你找个实际的故障录波数据，用Python或者C语言实现一遍FFT加窗，看看频谱长什么样。相信我，亲手做一遍，比看十遍书都管用。