C语言实现FFT算法详解从入门到精通掌握信号处理中的快速傅里叶变换编程技巧

引言

快速傅里叶变换（Fast Fourier Transform, FFT）是数字信号处理领域中最重要的算法之一。它是一种高效计算离散傅里叶变换（Discrete Fourier Transform, DFT）及其逆变换的算法，将计算复杂度从O(N²)降低到O(N log N)，极大地提高了信号处理的效率。FFT在音频处理、图像处理、通信系统、雷达信号处理等众多领域有着广泛的应用。

本文将详细介绍FFT算法的原理，并通过C语言实现从入门到精通的过程，帮助读者掌握信号处理中的快速傅里叶变换编程技巧。我们将从基础的傅里叶变换概念开始，逐步深入到FFT算法的实现细节和优化技巧，并提供完整的C语言代码示例。

傅里叶变换基础

离散傅里叶变换(DFT)

离散傅里叶变换(DFT)是傅里叶分析在离散域中的表现形式。对于一个长度为N的复数序列x[n]（n = 0, 1, …, N-1），其DFT定义为：

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

其中，X[k]是频域表示，x[n]是时域表示，j是虚数单位。

DFT的逆变换(IDFT)定义为：

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

DFT将信号从时域转换到频域，使我们能够分析信号的频率成分。然而，直接计算DFT需要O(N²)的运算复杂度，对于大的N值，计算量非常大。

快速傅里叶变换(FFT)的提出

1965年，James Cooley和John Tukey提出了一种高效计算DFT的算法，即快速傅里叶变换(FFT)。FFT利用了DFT计算中的对称性和周期性，将DFT分解为多个较小的DFT，从而将计算复杂度从O(N²)降低到O(N log N)。

FFT算法原理

Cooley-Tukey算法

最常用的FFT算法是Cooley-Tukey算法，它采用分治策略将DFT分解为较小的DFT。该算法要求序列长度N是2的幂次方（即N = 2^m），这种情况下称为基-2 FFT。

基-2 FFT的基本思想是将长度为N的序列x[n]分解为两个长度为N/2的子序列：

偶数索引序列：x[0], x[2], …, x[N-2]
奇数索引序列：x[1], x[3], …, x[N-1]

然后，分别计算这两个子序列的DFT，最后通过组合得到原序列的DFT。

数学上，可以表示为：

X[k] = E[k] + W_N^k * O[k] X[k + N/2] = E[k] - W_N^k * O[k], k = 0, 1, ..., N/2-1

其中，E[k]是偶数索引序列的DFT，O[k]是奇数索引序列的DFT，W_N^k = e^(-j * 2π * k / N)是旋转因子（也称为蝶形因子）。

这个过程可以递归地进行，直到子序列长度为1，此时DFT就是序列本身。

蝶形运算

FFT算法中的基本运算单元是蝶形运算，它描述了如何将两个较小的DFT组合成一个较大的DFT。基-2 FFT的蝶形运算如下图所示：

 ----->(+)----> X[k] / x[k] W_N^k / ----->(-)----> X[k + N/2]

蝶形运算包括一次复数乘法和两次复数加法（或减法）。

位反转

在FFT算法中，输入序列通常需要按照位反转的顺序重新排列。位反转是指将序列索引的二进制表示反转。例如，对于N=8，索引3（二进制011）的位反转是6（二进制110）。

位反转的目的是确保在递归分解过程中，数据能够正确地组合。

C语言实现FFT

基本数据结构

在开始实现FFT之前，我们需要定义一些基本的数据结构和常量：

#include <stdio.h> #include <stdlib.h> #include <math.h> #include <complex.h> #define PI 3.14159265358979323846 #define MAX_SIZE 1024 typedef double complex Complex; // 使用C99的复数类型

复数运算函数

虽然C99标准支持复数运算，但为了清晰起见，我们可以定义一些基本的复数运算函数：

// 复数加法 Complex complex_add(Complex a, Complex b) { return a + b; } // 复数减法 Complex complex_sub(Complex a, Complex b) { return a - b; } // 复数乘法 Complex complex_mul(Complex a, Complex b) { return a * b; } // 计算旋转因子 W_N^k = e^(-j * 2π * k / N) Complex twiddle_factor(int N, int k) { double angle = -2 * PI * k / N; return cos(angle) + I * sin(angle); }

位反转算法

位反转是FFT算法中的重要步骤，下面是实现位反转的函数：

// 计算一个数的位反转 int bit_reverse(int num, int log2n) { int reversed = 0; for (int i = 0; i < log2n; i++) { reversed |= ((num >> i) & 1) << (log2n - 1 - i); } return reversed; } // 对数组进行位反转重排 void bit_reverse_reorder(Complex *x, int N) { int log2n = 0; while ((1 << log2n) < N) { log2n++; } for (int i = 0; i < N; i++) { int reversed = bit_reverse(i, log2n); if (reversed > i) { Complex temp = x[i]; x[i] = x[reversed]; x[reversed] = temp; } } }

递归实现FFT

递归实现FFT是最直观的方法，它直接反映了FFT的分治思想：

// 递归实现FFT void fft_recursive(Complex *x, Complex *X, int N, int stride) { if (N == 1) { X[0] = x[0]; return; } // 分配临时数组 Complex *E = (Complex *)malloc(N/2 * sizeof(Complex)); Complex *O = (Complex *)malloc(N/2 * sizeof(Complex)); Complex *even_X = (Complex *)malloc(N/2 * sizeof(Complex)); Complex *odd_X = (Complex *)malloc(N/2 * sizeof(Complex)); // 分离偶数和奇数索引的元素 for (int i = 0; i < N/2; i++) { E[i] = x[2 * i * stride]; O[i] = x[(2 * i + 1) * stride]; } // 递归计算偶数和奇数序列的FFT fft_recursive(E, even_X, N/2, stride); fft_recursive(O, odd_X, N/2, stride); // 组合结果 for (int k = 0; k < N/2; k++) { Complex W = twiddle_factor(N, k); Complex temp = complex_mul(W, odd_X[k]); X[k] = complex_add(even_X[k], temp); X[k + N/2] = complex_sub(even_X[k], temp); } // 释放临时数组 free(E); free(O); free(even_X); free(odd_X); } // 包装函数，简化调用 void fft(Complex *x, Complex *X, int N) { fft_recursive(x, X, N, 1); }

非递归实现FFT

递归实现虽然直观，但效率不高，因为函数调用开销较大。下面是非递归实现FFT的方法：

// 非递归实现FFT void fft_iterative(Complex *x, Complex *X, int N) { // 首先进行位反转重排 for (int i = 0; i < N; i++) { X[i] = x[i]; } bit_reverse_reorder(X, N); // 进行FFT计算 for (int s = 1; s <= log2(N); s++) { int m = 1 << s; // 2^s Complex Wm = twiddle_factor(m, 1); // W_m^1 for (int k = 0; k < N; k += m) { Complex W = 1.0 + 0.0 * I; // W_m^0 = 1 for (int j = 0; j < m/2; j++) { Complex t = complex_mul(W, X[k + j + m/2]); Complex u = X[k + j]; X[k + j] = complex_add(u, t); X[k + j + m/2] = complex_sub(u, t); W = complex_mul(W, Wm); // W = W * W_m^1 } } } }

逆FFT实现

逆FFT（IFFT）与FFT类似，只有以下几点不同：

旋转因子的符号相反
最后需要除以N

// 逆FFT void ifft(Complex *X, Complex *x, int N) { // 首先进行位反转重排 for (int i = 0; i < N; i++) { x[i] = X[i]; } bit_reverse_reorder(x, N); // 进行IFFT计算 for (int s = 1; s <= log2(N); s++) { int m = 1 << s; // 2^s Complex Wm = twiddle_factor(m, -1); // 注意这里是-1 for (int k = 0; k < N; k += m) { Complex W = 1.0 + 0.0 * I; for (int j = 0; j < m/2; j++) { Complex t = complex_mul(W, x[k + j + m/2]); Complex u = x[k + j]; x[k + j] = complex_add(u, t); x[k + j + m/2] = complex_sub(u, t); W = complex_mul(W, Wm); } } } // 除以N for (int i = 0; i < N; i++) { x[i] /= N; } }

实数序列的FFT优化

在实际应用中，输入信号通常是实数序列，我们可以利用这一特性进行优化。一种常用的方法是使用一个N点的FFT来同时计算两个N点实数序列的FFT：

// 同时计算两个实数序列的FFT void fft2real(double *x1, double *x2, Complex *X1, Complex *X2, int N) { Complex *x = (Complex *)malloc(N * sizeof(Complex)); // 构造复数序列 x = x1 + j*x2 for (int i = 0; i < N; i++) { x[i] = x1[i] + I * x2[i]; } // 计算x的FFT Complex *X = (Complex *)malloc(N * sizeof(Complex)); fft_iterative(x, X, N); // 从X中提取X1和X2 X1[0] = creal(X[0]) + 0.0 * I; // X1[0]是实数 X2[0] = cimag(X[0]) + 0.0 * I; // X2[0]是实数 for (int k = 1; k < N/2; k++) { Complex Xk = X[k]; Complex X_N_minus_k = conj(X[N - k]); // 共轭 X1[k] = 0.5 * (Xk + X_N_minus_k); X1[N - k] = conj(X1[k]); // 共轭对称 X2[k] = -0.5 * I * (Xk - X_N_minus_k); X2[N - k] = conj(X2[k]); // 共轭对称 } X1[N/2] = creal(X[N/2]) + 0.0 * I; // X1[N/2]是实数 X2[N/2] = cimag(X[N/2]) + 0.0 * I; // X2[N/2]是实数 free(x); free(X); }

实际应用案例

信号频谱分析

FFT最常见的应用之一是分析信号的频谱。下面是一个简单的例子，展示如何使用FFT分析一个包含多个频率成分的信号：

#include <stdio.h> #include <stdlib.h> #include <math.h> #include <complex.h> #define PI 3.14159265358979323846 #define N 256 int main() { // 生成测试信号：包含50Hz和120Hz的正弦波 double signal[N]; double Fs = 1000.0; // 采样频率1000Hz double T = 1.0 / Fs; // 采样间隔 for (int i = 0; i < N; i++) { double t = i * T; signal[i] = 0.7 * sin(2 * PI * 50 * t) + sin(2 * PI * 120 * t); } // 准备FFT输入 Complex x[N]; for (int i = 0; i < N; i++) { x[i] = signal[i] + 0.0 * I; } // 计算FFT Complex X[N]; fft_iterative(x, X, N); // 计算幅度谱 double magnitude[N/2]; for (int i = 0; i < N/2; i++) { magnitude[i] = cabs(X[i]) / (N/2); } // 找出峰值频率 double max_magnitude = 0; int peak_index = 0; for (int i = 1; i < N/2; i++) { if (magnitude[i] > max_magnitude) { max_magnitude = magnitude[i]; peak_index = i; } } double peak_freq = peak_index * Fs / N; printf("Peak frequency: %.2f Hzn", peak_freq); // 输出频谱数据到文件，以便用其他工具绘制 FILE *fp = fopen("spectrum.txt", "w"); if (fp != NULL) { for (int i = 0; i < N/2; i++) { double freq = i * Fs / N; fprintf(fp, "%.2f %.6fn", freq, magnitude[i]); } fclose(fp); } return 0; }

卷积计算

FFT可以用于快速计算线性卷积，这在信号处理中非常有用。下面是一个使用FFT计算卷积的例子：

// 使用FFT计算线性卷积 void fft_convolution(double *x, double *h, double *y, int Nx, int Nh) { int N = 1; while (N < Nx + Nh - 1) { N <<= 1; // 找到大于等于Nx+Nh-1的最小2的幂次 } // 准备FFT输入 Complex *X = (Complex *)malloc(N * sizeof(Complex)); Complex *H = (Complex *)malloc(N * sizeof(Complex)); for (int i = 0; i < N; i++) { X[i] = (i < Nx) ? (x[i] + 0.0 * I) : (0.0 + 0.0 * I); H[i] = (i < Nh) ? (h[i] + 0.0 * I) : (0.0 + 0.0 * I); } // 计算FFT Complex *X_fft = (Complex *)malloc(N * sizeof(Complex)); Complex *H_fft = (Complex *)malloc(N * sizeof(Complex)); fft_iterative(X, X_fft, N); fft_iterative(H, H_fft, N); // 频域相乘 Complex *Y_fft = (Complex *)malloc(N * sizeof(Complex)); for (int i = 0; i < N; i++) { Y_fft[i] = X_fft[i] * H_fft[i]; } // 计算IFFT Complex *Y = (Complex *)malloc(N * sizeof(Complex)); ifft(Y_fft, Y, N); // 提取实部作为卷积结果 for (int i = 0; i < Nx + Nh - 1; i++) { y[i] = creal(Y[i]); } // 释放内存 free(X); free(H); free(X_fft); free(H_fft); free(Y_fft); free(Y); }

滤波器设计与应用

FFT可以用于设计和应用数字滤波器。下面是一个简单的低通滤波器例子：

// 设计和应用低通滤波器 void lowpass_filter(double *x, double *y, int N, double cutoff_freq, double Fs) { // 计算FFT Complex *X = (Complex *)malloc(N * sizeof(Complex)); Complex *x_complex = (Complex *)malloc(N * sizeof(Complex)); for (int i = 0; i < N; i++) { x_complex[i] = x[i] + 0.0 * I; } fft_iterative(x_complex, X, N); // 设计低通滤波器 Complex *H = (Complex *)malloc(N * sizeof(Complex)); int cutoff_index = (int)(cutoff_freq * N / Fs); for (int i = 0; i < N; i++) { if (i <= cutoff_index || i >= N - cutoff_index) { H[i] = 1.0 + 0.0 * I; // 通带 } else { H[i] = 0.0 + 0.0 * I; // 阻带 } } // 应用滤波器 Complex *Y = (Complex *)malloc(N * sizeof(Complex)); for (int i = 0; i < N; i++) { Y[i] = X[i] * H[i]; } // 计算IFFT Complex *y_complex = (Complex *)malloc(N * sizeof(Complex)); ifft(Y, y_complex, N); // 提取实部作为滤波结果 for (int i = 0; i < N; i++) { y[i] = creal(y_complex[i]); } // 释放内存 free(X); free(x_complex); free(H); free(Y); free(y_complex); }

性能优化

查表法优化旋转因子计算

在FFT算法中，旋转因子的计算是一个频繁操作。我们可以预先计算并存储所有旋转因子，避免重复计算：

// 预计算旋转因子表 Complex **precompute_twiddle_factors(int max_N) { int log2_max_N = 0; while ((1 << log2_max_N) < max_N) { log2_max_N++; } Complex **W = (Complex **)malloc((log2_max_N + 1) * sizeof(Complex *)); for (int s = 1; s <= log2_max_N; s++) { int m = 1 << s; W[s] = (Complex *)malloc(m/2 * sizeof(Complex)); for (int k = 0; k < m/2; k++) { W[s][k] = twiddle_factor(m, k); } } return W; } // 使用预计算的旋转因子表的FFT void fft_iterative_optimized(Complex *x, Complex *X, int N, Complex **W) { // 首先进行位反转重排 for (int i = 0; i < N; i++) { X[i] = x[i]; } bit_reverse_reorder(X, N); // 进行FFT计算 for (int s = 1; s <= log2(N); s++) { int m = 1 << s; // 2^s for (int k = 0; k < N; k += m) { for (int j = 0; j < m/2; j++) { Complex t = complex_mul(W[s][j], X[k + j + m/2]); Complex u = X[k + j]; X[k + j] = complex_add(u, t); X[k + j + m/2] = complex_sub(u, t); } } } }

SIMD指令优化

现代CPU支持SIMD（Single Instruction, Multiple Data）指令集，如SSE、AVX等，可以同时处理多个数据。我们可以使用SIMD指令来优化FFT中的复数运算：

#ifdef __SSE2__ #include <emmintrin.h> // 使用SSE指令优化的复数加法 void complex_add_sse(Complex *a, Complex *b, Complex *result, int N) { for (int i = 0; i < N; i++) { __m128d va = _mm_load_pd((double *)&a[i]); __m128d vb = _mm_load_pd((double *)&b[i]); __m128d vresult = _mm_add_pd(va, vb); _mm_store_pd((double *)&result[i], vresult); } } // 使用SSE指令优化的复数乘法 void complex_mul_sse(Complex *a, Complex *b, Complex *result, int N) { for (int i = 0; i < N; i++) { double ar = creal(a[i]); double ai = cimag(a[i]); double br = creal(b[i]); double bi = cimag(b[i]); __m128d va = _mm_set_pd(ai, ar); __m128d vb_brbr = _mm_set_pd(br, br); __m128d vb_bibi = _mm_set_pd(bi, bi); __m128d temp1 = _mm_mul_pd(va, vb_brbr); __m128d temp2 = _mm_mul_pd(va, vb_bibi); __m128d shuf = _mm_shuffle_pd(temp2, temp2, 1); __m128d temp3 = _mm_add_pd(temp1, shuf); __m128d temp4 = _mm_sub_pd(temp1, shuf); __m128d vresult = _mm_shuffle_pd(temp4, temp3, 2); _mm_store_pd((double *)&result[i], vresult); } } #endif

多线程并行计算

FFT算法具有很好的并行性，我们可以使用多线程来加速计算：

#include <pthread.h> typedef struct { Complex *X; int N; int s; int start; int end; Complex **W; } FFTThreadData; // 线程函数 void *fft_thread_func(void *arg) { FFTThreadData *data = (FFTThreadData *)arg; int m = 1 << data->s; for (int k = data->start; k < data->end; k += m) { for (int j = 0; j < m/2; j++) { Complex t = complex_mul(data->W[data->s][j], data->X[k + j + m/2]); Complex u = data->X[k + j]; data->X[k + j] = complex_add(u, t); data->X[k + j + m/2] = complex_sub(u, t); } } return NULL; } // 多线程FFT void fft_iterative_parallel(Complex *x, Complex *X, int N, Complex **W, int num_threads) { // 首先进行位反转重排 for (int i = 0; i < N; i++) { X[i] = x[i]; } bit_reverse_reorder(X, N); // 进行FFT计算 for (int s = 1; s <= log2(N); s++) { int m = 1 << s; // 2^s // 分配线程数据 pthread_t threads[num_threads]; FFTThreadData thread_data[num_threads]; int chunk_size = N / num_threads; // 创建线程 for (int i = 0; i < num_threads; i++) { thread_data[i].X = X; thread_data[i].N = N; thread_data[i].s = s; thread_data[i].start = i * chunk_size; thread_data[i].end = (i == num_threads - 1) ? N : (i + 1) * chunk_size; thread_data[i].W = W; pthread_create(&threads[i], NULL, fft_thread_func, &thread_data[i]); } // 等待线程完成 for (int i = 0; i < num_threads; i++) { pthread_join(threads[i], NULL); } } }

常见问题与解决方案

非2的幂次长度的FFT

标准的基-2 FFT算法要求序列长度是2的幂次。对于非2的幂次长度的序列，有几种处理方法：

零填充：将序列填充到最近的2的幂次长度。这是最简单的方法，但会引入额外的计算和频谱泄漏。

// 零填充到最近的2的幂次 int next_power_of_2(int n) { int power = 1; while (power < n) { power <<= 1; } return power; } // 零填充并计算FFT void fft_with_zero_padding(double *x, Complex *X, int Nx) { int N = next_power_of_2(Nx); Complex *x_padded = (Complex *)malloc(N * sizeof(Complex)); for (int i = 0; i < N; i++) { x_padded[i] = (i < Nx) ? (x[i] + 0.0 * I) : (0.0 + 0.0 * I); } fft_iterative(x_padded, X, N); free(x_padded); }

混合基FFT：使用更通用的FFT算法，如混合基FFT，可以处理任意长度的序列。混合基FFT将序列分解为多个较小长度的DFT，这些长度可以是2、3、4、5等。
素数因子算法：当序列长度可以分解为互素的因数时，可以使用素数因子算法。

频谱泄漏与窗函数

在实际应用中，由于信号的截断，可能会出现频谱泄漏现象。为了减少频谱泄漏，可以使用窗函数对信号进行加权：

// 应用窗函数 typedef enum { RECTANGULAR, HANNING, HAMMING, BLACKMAN } WindowType; // 窗函数 void apply_window(double *x, double *windowed, int N, WindowType type) { for (int i = 0; i < N; i++) { double w; switch (type) { case RECTANGULAR: w = 1.0; break; case HANNING: w = 0.5 - 0.5 * cos(2 * PI * i / (N - 1)); break; case HAMMING: w = 0.54 - 0.46 * cos(2 * PI * i / (N - 1)); break; case BLACKMAN: w = 0.42 - 0.5 * cos(2 * PI * i / (N - 1)) + 0.08 * cos(4 * PI * i / (N - 1)); break; default: w = 1.0; } windowed[i] = x[i] * w; } }

数值精度问题

在FFT计算中，可能会遇到数值精度问题，特别是对于大型FFT。以下是一些解决方案：

使用更高精度的数据类型：例如，使用double而不是float。
归一化：在每一步蝶形运算后进行归一化，防止数值过大或过小。

// 归一化的FFT void fft_iterative_normalized(Complex *x, Complex *X, int N) { // 首先进行位反转重排 for (int i = 0; i < N; i++) { X[i] = x[i]; } bit_reverse_reorder(X, N); // 进行FFT计算 for (int s = 1; s <= log2(N); s++) { int m = 1 << s; // 2^s Complex Wm = twiddle_factor(m, 1); for (int k = 0; k < N; k += m) { Complex W = 1.0 + 0.0 * I; for (int j = 0; j < m/2; j++) { Complex t = complex_mul(W, X[k + j + m/2]); Complex u = X[k + j]; X[k + j] = complex_add(u, t); X[k + j + m/2] = complex_sub(u, t); // 归一化 X[k + j] /= sqrt(2); X[k + j + m/2] /= sqrt(2); W = complex_mul(W, Wm); } } } }

Kahan求和算法：使用Kahan求和算法来减少累加过程中的舍入误差。

// Kahan求和算法 Complex kahan_sum(Complex *arr, int N) { Complex sum = 0.0 + 0.0 * I; Complex c = 0.0 + 0.0 * I; // 补偿 for (int i = 0; i < N; i++) { Complex y = arr[i] - c; Complex t = sum + y; c = (t - sum) - y; sum = t; } return sum; }