4. FFT原理：离散傅里叶变换推导、FFT算法核心思想、频谱泄露与窗函数

好，咱们进入第四讲。这一讲，我打算把FFT的底裤扒干净。

很多人用FFT，就是调个库，出个图。但实际项目中，你如果不懂原理，出了问题根本不知道从哪下手。我记得刚入行那会儿，调一个振动信号的频谱，怎么调都有毛刺，折腾了两天，最后发现是窗函数没选对。嗯，从那以后，我再也不敢跳过原理直接调库了。

4.1 从连续到离散：DFT的诞生

咱们先聊聊离散傅里叶变换，也就是DFT。

你想想看，现实世界里的信号，比如声音、振动，都是连续的。但计算机不认识连续信号，它只认离散的点。所以，我们必须把连续信号切成一段一段的，再采样成离散的点。

DFT干的事，说白了就是：给一串离散的点，算出它包含哪些频率成分。

公式长这样：

X[k] = Σ_{n=0}^{N-1} x[n] * e^{-j(2π/N)kn}

其中：

• x[n]：时域的第n个采样点
• X[k]：频域的第k个频率分量
• N：采样点数
• e^{-j(2π/N)kn}：旋转因子，也叫基函数

这个公式怎么理解？我个人的习惯是把它看作一个「匹配度检测器」。你把信号x[n]和不同频率的基函数做内积，内积越大，说明这个频率成分越强。

4.2 FFT：DFT的加速器

DFT虽然好，但计算量太大了。你算一个N点的DFT，需要N²次复数乘法。如果N=1024，那就是一百多万次计算。这在嵌入式系统里根本跑不动。

FFT就是来解决这个问题的。它的核心思想就四个字：分而治之。

怎么分？我举个例子。假设你有8个点的信号：

1. 先把8个点分成两组：偶数点一组，奇数点一组
2. 每组4个点，各自做DFT
3. 再把两组结果合并起来

合并的时候有个技巧，叫「蝶形运算」。你看这个结构像不像蝴蝶的两只翅膀？

// 蝶形运算核心代码（伪代码）
for (int len = 2; len <= N; len *= 2) {
    for (int i = 0; i < N; i += len) {
        for (int j = 0; j < len/2; j++) {
            // 计算旋转因子
            complex w = exp(-2*PI*j/len);
            // 蝶形运算
            complex t = w * x[i + j + len/2];
            complex u = x[i + j];
            x[i + j] = u + t;
            x[i + j + len/2] = u - t;
        }
    }
}

这样一拆分，计算量就从N²降到了N*log₂(N)。N=1024时，FFT只需要约一万次计算，比DFT快了一百多倍。

4.3 频谱泄露：一个躲不开的坑

好，现在你学会了FFT。但别高兴太早，有个大坑等着你——频谱泄露。

为什么会泄露？因为FFT假设你的信号是周期性的。但实际中，你截取的一段信号，首尾往往不连续。这种不连续，在频域里就会表现为能量扩散到旁边的频率上。

举个例子：你采集了一个50Hz的正弦波，采样了1秒钟。按理说频谱上应该只有50Hz一根线。但实际做FFT后你会发现，50Hz旁边也出现了能量。这就是泄露。

4.4 窗函数：给信号做「整形」

解决频谱泄露的办法，就是加窗。

窗函数的作用，是把截取的那段信号的首尾慢慢衰减到0，这样首尾就连续了。说白了，就是给信号做个「整形手术」。

常用的窗函数有这些：

窗函数	主瓣宽度	旁瓣衰减	适用场景
矩形窗	最窄	最差（-13dB）	瞬态信号、频率分辨率要求高
汉宁窗	较宽	较好（-31dB）	一般信号、振动分析
海明窗	较宽	更好（-41dB）	语音信号、窄带信号
布莱克曼窗	最宽	最好（-58dB）	需要极低旁瓣的场景

怎么选？我个人的习惯是：

• 如果信号是瞬态的（比如冲击信号），用矩形窗，因为要保留瞬态信息
• 如果是稳态信号（比如旋转机械的振动），用汉宁窗，兼顾分辨率和泄露抑制
• 如果两个频率靠得很近，需要区分它们，用矩形窗或汉宁窗
• 如果信号很弱，容易被旁瓣淹没，用布莱克曼窗

4.5 实战中的FFT流程

好了，理论讲完了。我分享一下我在项目中常用的FFT处理流程：

1. 数据预处理：去直流分量（减去均值），去除趋势项
2. 选择窗函数：根据信号类型和需求选择
3. 加窗：信号乘以窗函数
4. 补零：如果点数不是2的幂次，补零到最近的2的幂次
5. 做FFT：调用库函数
6. 取幅值：计算|X[k]|，注意要乘以2/N才能得到真实幅值
7. 频率轴：频率分辨率 = 采样率 / FFT点数

import numpy as np

def my_fft_analysis(signal, fs, window='hann'):
    N = len(signal)
    # 1. 去直流
    signal = signal - np.mean(signal)
    
    # 2. 加窗
    if window == 'hann':
        w = np.hanning(N)
    elif window == 'hamming':
        w = np.hamming(N)
    else:
        w = np.ones(N)  # 矩形窗
    
    signal_windowed = signal * w
    
    # 3. 补零到2的幂次
    N_fft = 2 ** int(np.ceil(np.log2(N)))
    signal_padded = np.pad(signal_windowed, (0, N_fft - N), 'constant')
    
    # 4. FFT
    spectrum = np.fft.fft(signal_padded)
    spectrum = spectrum[:N_fft//2]  # 取单边谱
    
    # 5. 幅值修正
    amplitude = np.abs(spectrum) * 2 / N_fft
    
    # 6. 频率轴
    freq = np.linspace(0, fs/2, N_fft//2)
    
    return freq, amplitude

这段代码我用了很多年了，基本没出过问题。你拿去用的时候，注意根据实际情况调整窗函数和补零策略。