[实战] 基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点DFT
                dft8(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=64
    complex 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);
    
    // 执行FFT
    for (int i = 0; i < N; i++) fft_out[i] = input[i];
    fft_base8(fft_out, N);
    save_bin("fft_c.bin", fft_out, N);
    
    // 执行IFFT
    for (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 plt

def 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-06
    IFFT 最大误差: 1.459477e-06 平均误差:6.917030e-07
    IFFT_原始信号 最大误差: 1.513257e-06 平均误差:6.714861e-07
   

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

关键特性

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

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

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

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