快速傅里叶变换（FFT）原理与实现

什么是FFT？#

快速傅里叶变换（Fast Fourier Transform，FFT）是一种高效计算离散傅里叶变换（DFT）的算法，由Cooley和Tukey于1965年提出。它将DFT的时间复杂度从O(N²)降低到O(N log N)，使傅里叶分析在实时信号处理中变得可行。

核心思想#

FFT利用分治策略和旋转因子的周期性，将大点数的DFT分解为小点数的DFT组合计算。最常见的基2-FFT算法要求信号长度N为2的整数次幂。

FFT应用场景#

信号处理：音频/图像滤波、频谱分析
通信系统：OFDM调制解调、信道估计
医学成像：MRI、CT扫描重建
雷达系统：目标检测与测距
数据压缩：MP3/JPEG等压缩算法

算法原理#

离散傅里叶变换（DFT）#

对于N点离散信号x[n]，其DFT定义为： $X[k] = \sum_{n=0}^{N-1} x[n] \cdot e^{-j2\pi kn/N} \quad k=0,1,...,N-1$

FFT分治过程#

按奇偶分解： $X[k] = \sum_{m=0}^{N/2-1} x[2m] \cdot W_N^{2mk} + W_N^k \sum_{m=0}^{N/2-1} x[2m+1] \cdot W_N^{2mk}$
旋转因子性质： $W_N^{2mk} = W_{N/2}^{mk}$
递归计算： $X[k] = G[k] + W_N^k H[k]$ 其中G[k]是偶数点DFT，H[k]是奇数点DFT

代码实现分析#

数据结构定义#

1
typedef struct {
2
    float real;   // 实部
3
    float imag;   // 虚部
4
} compx;

FFT核心函数#

1
/**
2
 * @brief 快速傅立叶变换-FFT
3
 *
4
 * @param xin  待变换数组
5
 * @param Num  数组元素个数
6
 * @param xout 存储FFT结果的数组
7
 */
8
void FFT(const compx *xin, int Num, compx *xout)
9
{
10
    int f, m, LH, nm, i, k, j, L;
11
    double p, ps; // 双精度 p: 指数，ps: 指数
12
    int le, B, ip; // le: 级数间隔，B: 级数间隔的点数，ip: 级数间隔的点序号
13
    float pi; // 圆周率
14
    compx w, t; // w: 旋转因子，t: 临时变量
15

16
    // 初始化输出数组，复制输入数组的内容
17
    for (i = 0; i < Num; ++i) {
18
        xout[i].real = xin[i].real;
19
        xout[i].imag = xin[i].imag;
20
    }
21

22
    LH = Num / 2;
23
    f = Num;
24
    for (m = 1; (f = f / 2) != 1; m++);
25

26
    for (L = m; L >= 1; L--) {
27
        le = pow(2, L);
28
        B = le / 2;
29
        pi = 3.1415;
30
        for (j = 0; j <= B - 1; j++) {
31
            p = pow(2, m - L) * j;
32
            ps = 2 * pi / Num * p;
33
            w.real = cos(ps);
34
            w.imag = -sin(ps);
35
            for (i = j; i <= Num - 1; i += le) {
36
                ip = i + B;
37
                t = xout[i]; // 临时变量保存xout[i]
38
                xout[i].real = t.real + xout[ip].real; // 蝶形运算
39
                xout[i].imag = t.imag + xout[ip].imag;
40
                xout[ip].real = t.real - xout[ip].real;
41
                xout[ip].imag = t.imag - xout[ip].imag;
42
                xout[ip] = EE(xout[ip], w); // 乘以旋转因子
43
            }
44
        }
45
    }
46

47
    // 变址运算
48
    nm = Num - 2;
49
    j = Num / 2;
50
    for (i = 1; i <= nm; i++) {
51
        if (i < j) {
52
            t = xout[j];
53
            xout[j] = xout[i];
54
            xout[i] = t;
55
        }
56
        k = LH;
57
        while (j >= k) {
58
            j = j - k;
59
            k = k / 2;
60
        }
61
        j = j + k;
62
    }
63
}

关键运算函数#

1
/**
2
 * @brief 复数乘法
3
 *
4
 * @param b1 : 复数b1
5
 * @param b2 : 复数b2
6
 * @return struct compx
7
 */
8
compx EE(compx b1, compx b2)
9
{
10
    compx b3;
11
    b3.real = b1.real*b2.real - b1.imag*b2.imag;
12
    b3.imag = b1.real*b2.imag + b1.imag*b2.real;
13
    return (b3);
14
}

测试信号生成与验证#

测试信号配置#

1
#define SIGNAL_LENGTH 1024 // 信号长度，必须为 2 的幂
2
#define SAMPLING_RATE 1024  // 采样率，单位 Hz
3
#define SIGNAL_FREQ 74.6      // 信号频率，单位 Hz
4

5
// 生成测试信号
6
void generate_test_signal(compx *signal)
7
{
8
  int i;
9
    for (i = 0; i < SIGNAL_LENGTH; ++i) {
10
        signal[i].real = sin(2 * M_PI * SIGNAL_FREQ * i / SAMPLING_RATE);
11
        signal[i].imag = 0.0f;
12
    }
13
}

FFT结果分析#

1
void test_fft() {
2
    int i = 0;
3
    compx test_signal[SIGNAL_LENGTH];
4
    compx fft_result[SIGNAL_LENGTH];
5

6
    generate_test_signal(test_signal);
7
    FFT(test_signal, SIGNAL_LENGTH, fft_result);
8

9
    // 打印结果
10
    printf("FFT Result:\n");
11
    for (i = 0; i < SIGNAL_LENGTH / 2; ++i) {
12
        float magnitude = sqrt(fft_result[i].real * fft_result[i].real
13
                             + fft_result[i].imag * fft_result[i].imag);
14
        float freq = (float)i * SAMPLING_RATE / SIGNAL_LENGTH;
15
        printf("Index: %d, Frequency: %f Hz, Magnitude: %f\n", i, freq, magnitude);
16
    }
17

18
    // 检查最大值
19
    float max_magnitude = 0.0f;
20
    int max_index = 0;
21

22
    for (i = 0; i < SIGNAL_LENGTH / 2; ++i) {
23
        float magnitude = sqrt(fft_result[i].real * fft_result[i].real
24
                             + fft_result[i].imag * fft_result[i].imag);
25
        if (magnitude > max_magnitude) {
26
            max_magnitude = magnitude;
27
            max_index = i;
28
        }
29
    }
30

31
    // 使用max_index计算频率
32
    float max_freq = (float)max_index * SAMPLING_RATE / SIGNAL_LENGTH;
33
    printf("\nPeak at Index: %d, Frequency: %f Hz, Magnitude: %f\n",
34
           max_index, max_freq, max_magnitude);
35
}

常见问题及解决方案#

1. 频率分辨率不足#

问题：无法区分相近频率的信号
解决：增加采样点数

1
#define SIGNAL_LENGTH 4096  // 提高频率分辨率

2. 频谱泄漏#

问题：非整数周期截断导致能量分散
解决：使用窗函数

1
// 汉宁窗应用
2
for (int i = 0; i < SIGNAL_LENGTH; i++) {
3
    float window = 0.5 - 0.5 * cos(2 * M_PI * i / (SIGNAL_LENGTH-1));
4
    signal[i].real *= window;
5
}

3. 幅值校准#

问题：FFT结果需要缩放才能反映真实幅值
解决：对结果进行归一化

1
magnitude = sqrt(real*real + imag*imag) / (SIGNAL_LENGTH / 2);
2
// 直流分量特殊处理
3
if (i == 0) magnitude /= 2;

4. 复数信号处理#

问题：实信号FFT的对称性
现象：频谱关于Nyquist频率对称
利用：只需计算前半部分频谱（0 ~ N/2）

实际运行结果#

1
Peak at Index: 75, Frequency: 75.000000 Hz, Magnitude: 386.443970

FFT优化技巧#

查表法：预计算旋转因子
并行化：利用SIMD指令加速蝶形运算
混合基算法：支持非2的幂长度
迭代实现：避免递归开销

总结#

FFT是数字信号处理的基石算法，理解其原理和实现细节对于开发高性能信号处理系统至关重要。通过合理选择参数（采样率、点数）和应用窗函数等技术，可以获得精确的频谱分析结果。

音乐

音乐