4. 时域与频域分析：FFT算法实现、频谱泄露与窗函数

说实话，做通信算法这些年，我最大的感触就是——时域和频域就像一枚硬币的两面。你在时域里看半天搞不明白的问题，切到频域里一眼就清楚了。今天咱们就来聊聊这个核心工具：FFT，以及它背后那些让人头疼的坑。

4.1 从DFT到FFT：一个加速的故事

先说说DFT（离散傅里叶变换）。它的数学公式长这样：

X[k] = Σ x[n] · e^(-j·2π·k·n/N)   (n=0到N-1)

你算一下复杂度——O(N²)。N=1024时，要算一百多万次乘法。我早年用DSP做实时处理时，这简直要命。

后来有了FFT（快速傅里叶变换）。它利用旋转因子的对称性和周期性，把复杂度降到了O(N·log₂N)。N=1024时，只需要约一万次运算。快了上百倍。

嗯，这里要注意：FFT不是一种新的变换，它只是DFT的一种高效算法实现。说白了，结果一模一样，只是算得更快。

4.2 基2时间抽取FFT算法

我个人最常用的是基2时间抽取法（DIT-FFT）。它的思路很简单：

1. 把序列按奇偶分成两组
2. 分别计算两组N/2点的DFT
3. 用蝶形运算合并结果

下面是我写的一个Python实现，你直接拿去用：

import numpy as np

def fft_dit(x):
    """基2时间抽取FFT"""
    N = len(x)
    if N <= 1:
        return x
    
    # 按奇偶拆分
    even = fft_dit(x[0::2])
    odd = fft_dit(x[1::2])
    
    # 蝶形运算
    result = np.zeros(N, dtype=complex)
    for k in range(N // 2):
        # 旋转因子
        W = np.exp(-2j * np.pi * k / N)
        result[k] = even[k] + W * odd[k]
        result[k + N//2] = even[k] - W * odd[k]
    
    return result

# 测试一下
x = np.array([1, 2, 3, 4, 4, 3, 2, 1])
X = fft_dit(x)
print(np.round(X, 3))

我在项目中遇到过一个问题：输入长度必须是2的幂次。如果不是，就补零到下一个2的幂。这个操作叫"补零FFT"。

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

你想想看，FFT处理的是有限长序列。但实际信号往往是无限长的。这就好比从窗户里看风景——你只能看到窗框框住的那部分。

频谱泄露，说白了就是：信号的能量从主瓣"漏"到了旁瓣。为什么会这样？

因为截断操作在时域里相当于乘了一个矩形窗。矩形窗的频域响应是sinc函数，它有很高的旁瓣。这些旁瓣会把能量"洒"到不该去的地方。

我举个例子：

import numpy as np
import matplotlib.pyplot as plt

# 生成一个50Hz的正弦波
fs = 1000  # 采样率
t = np.arange(0, 1, 1/fs)
f0 = 50
x = np.sin(2 * np.pi * f0 * t)

# 取前0.1秒（非整数周期）
N = 100
x_truncated = x[:N]

# 做FFT
X = np.fft.fft(x_truncated)
freq = np.fft.fftfreq(N, 1/fs)

# 你会发现：50Hz附近有能量扩散
plt.plot(freq[:N//2], np.abs(X[:N//2]))
plt.xlabel('频率(Hz)')
plt.ylabel('幅度')
plt.show()

运行这段代码，你会看到50Hz附近出现了一个"小山坡"，而不是一根干净的谱线。这就是频谱泄露。

4.4 窗函数：给数据"加个框"

既然矩形窗不好，那就换个窗呗。窗函数的作用就是：让截断变得平滑，减少旁瓣能量。

常用的窗函数有这些：

窗函数	主瓣宽度	旁瓣衰减	适用场景
矩形窗	最窄	-13dB	频率分辨率要求高，无泄露问题
汉宁窗	较宽	-31dB	通用场景，平衡性好
海明窗	较宽	-41dB	旁瓣要求高，可接受分辨率损失
布莱克曼窗	最宽	-57dB	强旁瓣抑制，分辨率要求低
凯泽窗	可调	可调	灵活，通过β参数控制

我个人习惯用汉宁窗。它在主瓣宽度和旁瓣衰减之间取得了很好的平衡。代码实现也很简单：

def hanning_window(N):
    """生成汉宁窗"""
    n = np.arange(N)
    return 0.5 * (1 - np.cos(2 * np.pi * n / (N - 1)))

# 使用窗函数
x_windowed = x_truncated * hanning_window(N)
X_windowed = np.fft.fft(x_windowed)

加了窗之后，你会发现频谱干净多了。但代价是主瓣变宽了——频率分辨率下降了。这就是一个trade-off。

4.5 知识体系总览

下面这张图总结了本章的核心逻辑：

4.6 实战建议

最后，给你几个实战中的建议：

• 先判断是否需要加窗：如果信号是周期性的且截断长度正好是整数周期，可以不用窗。否则，加窗是必须的。
• 窗函数的选择：先试汉宁窗。如果旁瓣还不够低，换海明窗或布莱克曼窗。
• 补零不是万能的：补零只能让频谱看起来更平滑，不能提高真实的分辨率。分辨率由信号长度决定。
• 注意幅度校正：加窗后信号能量会变化，记得做幅度校正。校正系数 = N / sum(window)。

---

公众号:蓝海资料掘金营，微信deep3321