[实战] 基8 FFT/IFFT算法原理与实现(完整C代码)

基8 FFT/IFFT算法原理与实现

本文从数学原理开始，不使用任何第三方库，用C语言实现了基于基8的FFT和IFFT算法，效率高，可以移植修改，适合各种场景应用。同时用python验证了C语言实现的正确性，可放心使用。
闲话少说，实战开始：

数学原理

离散傅里叶变换（DFT）定义为：
$$X[k]=\sum _{n=0}^{N-1}x[n]\cdot e^{-j2\pi kn/N},k=0,1,\dots ,N-1$$
逆变换（IDFT）为：
$$x[n]=\frac{1}{N}\sum _{k=0}^{N-1}X[k]\cdot e^{j2\pi kn/N},n=0,1,\dots ,N-1$$

基8算法要求$N=8^m$。将输入序列分解为8个子序列：
$$x_r[m]=x[8m+r],r=0,1,\dots ,7;m=0,1,\dots ,N/8-1$$

DFT可表示为：
X[k]=∑r=07WNrk⋅DFT{xr[m]}[kmod  (N/8)] X[k] = \sum_{r=0}^{7} W_N^{rk} \cdot DFT\{x_r[m]\}[k \mod (N/8)] X[k]=r=0∑7​WNrk​⋅DFT{xr​[m]}[kmod(N/8)]
其中$W_N=e^{-j2\pi /N}$为旋转因子。

算法流程

1. 位反转置换：将输入序列按8进制位反转顺序重排
2. 蝶形运算：分${\log }_{8}N$级迭代：
   • 第$s$级处理$N/8^s$个8点组
   • 每组应用8点DFT矩阵运算
   • 乘以旋转因子$W_N^{rk}$（$r$=子序列索引，$k$=频率索引）
3. IFFT处理：FFT结果取共轭→执行FFT→取共轭→除以$N$

C语言实现

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdint.h>typedef float real;
typedef struct { real real, imag; } complex;#define PI 3.14159265358979323846f// 复数运算
complex cmplx(real r, real i) { return (complex){r, i}; }
complex cadd(complex a, complex b) { return cmplx(a.real+b.real, a.imag+b.imag); }
complex csub(complex a, complex b) { return cmplx(a.real-b.real, a.imag-b.imag); }
complex cmul(complex a, complex b) {return cmplx(a.real*b.real - a.imag*b.imag, a.real*b.imag + a.imag*b.real);
}// 位反转置换 (基8)
void bit_reverse(complex *x, int N) {int m = (int)(logf(N)/logf(8) + 0.5f);complex *temp = malloc(N * sizeof(complex));for (int i = 0; i < N; i++) {int rev = 0, n = i;for (int j = 0; j < m; j++) {rev = (rev << 3) | (n & 7);n >>= 3;}temp[rev] = x[i];}for (int i = 0; i < N; i++) x[i] = temp[i];free(temp);
}// 8点DFT核函数
void dft8(complex *x) {complex temp[8];for (int k = 0; k < 8; k++) {temp[k] = cmplx(0, 0);for (int n = 0; n < 8; n++) {real angle = -2*PI*k*n/8;complex w = cmplx(cosf(angle), sinf(angle));temp[k] = cadd(temp[k], cmul(x[n], w));}}for (int i = 0; i < 8; i++) x[i] = temp[i];
}// 基8 FFT
void fft_base8(complex *x, int N) {bit_reverse(x, N);int stages = (int)(logf(N)/logf(8));for (int stage = 0; stage < stages; stage++) {int group_size = 1;for (int i = 0; i <= stage; i++) group_size *= 8;int half_group = group_size / 8;for (int group = 0; group < N; group += group_size) {for (int pos = 0; pos < half_group; pos++) {complex data[8];// 加载数据并应用旋转因子for (int r = 0; r < 8; r++) {int idx = group + r*half_group + pos;real angle = -2*PI*r*pos/(group_size);complex w = cmplx(cosf(angle), sinf(angle));data[r] = cmul(x[idx], w);}// 执行8点DFTdft8(data);// 存储结果for (int r = 0; r < 8; r++) {x[group + r*half_group + pos] = data[r];}}}}
}// 基8 IFFT
void ifft_base8(complex *x, int N) {// 取共轭for (int i = 0; i < N; i++) x[i].imag = -x[i].imag;fft_base8(x, N);// 取共轭并缩放real scale = 1.0f/N;for (int i = 0; i < N; i++) {x[i].real *= scale;x[i].imag = -x[i].imag * scale;}
}// 文件操作
void save_bin(const char *fname, complex *data, int N) {FILE *f = fopen(fname, "wb");fwrite(data, sizeof(complex), N, f);fclose(f);
}void load_bin(const char *fname, complex *data, int N) {FILE *f = fopen(fname, "rb");fread(data, sizeof(complex), N, f);fclose(f);
}int main() {const int N = 64;  // 8^2=64complex input[N], fft_out[N], ifft_out[N];// 生成测试信号: 2个正弦波叠加for (int i = 0; i < N; i++) {real t = 2*PI*i/N;input[i] = cmplx(sinf(3*t) + 0.5f*cosf(7*t), 0);}// 保存输入数据save_bin("input_c.bin", input, N);// 执行FFTfor (int i = 0; i < N; i++) fft_out[i] = input[i];fft_base8(fft_out, N);save_bin("fft_c.bin", fft_out, N);// 执行IFFTfor (int i = 0; i < N; i++) ifft_out[i] = fft_out[i];ifft_base8(ifft_out, N);save_bin("ifft_c.bin", ifft_out, N);return 0;
}

Python验证脚本

import numpy as np
import matplotlib.pyplot as pltdef read_bin(fname, n):return np.fromfile(fname, dtype=np.complex64, count=n)def compare(data1, data2, name):diff = np.abs(data1 - data2)print(f"{name} 最大误差: {np.max(diff):.6e}")plt.figure()plt.plot(diff, 'r-')plt.title(f"{name}误差分析")plt.ylabel("误差幅度")plt.xlabel("采样点")plt.grid()plt.savefig(f"{name}_diff.png")N = 64# 读取C语言生成的数据
input_c = read_bin("input_c.bin", N)
fft_c = read_bin("fft_c.bin", N)
ifft_c = read_bin("ifft_c.bin", N)# 使用NumPy计算参考值
fft_ref = np.fft.fft(input_c)
ifft_ref = np.fft.ifft(fft_ref)# 结果对比
compare(fft_c, fft_ref, "FFT")
compare(ifft_c, ifft_ref, "IFFT")
compare(ifft_c, input_c, "IFFT_原始信号")# 保存参考结果
fft_ref.astype(np.complex64).tofile("fft_py.bin")
ifft_ref.astype(np.complex64).tofile("ifft_py.bin")

验证方法

1. 数据生成：

   • C程序生成测试信号并保存为input_c.bin
   • 执行FFT保存为fft_c.bin
   • 执行IFFT保存为ifft_c.bin

2. Python验证：

   • 读取C生成的二进制文件
   • 使用NumPy计算FFT/IFFT作为参考
   • 计算并报告最大误差
   • 生成误差分析图

FFT结果：

IFFT结果：

IFFT原始数据：

3. 误差分析：

    FFT 最大误差: 6.492365e-05 平均误差:4.865580e-06IFFT 最大误差: 1.459477e-06 平均误差:6.917030e-07IFFT_原始信号 最大误差: 1.513257e-06 平均误差:6.714861e-07
   

   误差主要来源于浮点精度差异，在可接受范围内

关键特性

1. 原位计算：FFT/IFFT均原位操作减少内存占用
2. 精确位反转：8进制位反转确保正确数据访问模式
3. 模块化设计：8点DFT作为独立核函数
4. 完整流程：包含数据生成→变换→存储→验证全流程

此实现完整展示了基8 FFT算法的原理和实现，并通过严格的数值验证确保了正确性。

---

注意：因为基8算法的特性，本文代码只能计算8^(2n)点数的FFT和IFFT

---

研究学习不易，点赞易。
工作生活不易，收藏易，点收藏不迷茫 ：）

---