DIT基-2FFT的C语言实现

Zuko 2021-09-01 原文

在学习完FFT原理后，尝试用C语言写的程序，主要用于加深理解，没有考虑实用性

使用倒序输入，同址运算，顺序输出的方法

为了省事，直接用了C99的complex.h库

经测试与Matlab运算结果一致

#include <complex.h>

#define _USE_MATH_DEFINES
#include <math.h>

void fft(double complex x[], int N);
void ButterflyDiag(double complex * base, double complex * x1, double complex * x2, int layers, int presentlayer, int N);    // 蝶形结运算 
void Inversion(double complex x[], int N);   // 原始数据倒序 


void fft(double complex x[], int N)
{
    int layers = (int) log2(N);   // 迭代层数，默认点数为2的整数幂 
    int points, groups, distant, presentlayer;
    int i, j;
    double complex * base;
    
    Inversion(x, N);
    
    points = 2; // 每一层每个小组内的点数 
    for (presentlayer = 1; presentlayer <= layers; presentlayer++) // 当前层 
    {
        groups = N / points; // 每一层的分组数 
        distant = points / 2; // 蝶形运算两点的距离，也是每一组包含的蝶形运算数量 
        for (i = 0; i < groups; i++)
        {
            base = x + i * points; // 每一组的起始位置 
            for (j = 0; j < distant; j++)
                ButterflyDiag(base, base + j, base + j + distant, layers, presentlayer, N);
        }
        points *= 2;
    }
}

void ButterflyDiag(double complex * base, double complex * x1, double complex * x2, int layers, int presentlayer, int N)
{
    double complex tmp1 = *x1;
    double complex tmp2 = *x2;
    int r = (x1 - base) << (layers - presentlayer);
    double complex wnr= cexp(-I * 2.0 * M_PI / N * r);
    
    *x1 = tmp1 + tmp2 * wnr; //同址蝶形运算 
    *x2 = tmp1 - tmp2 * wnr;
    
}

void Inversion(double complex x[], int N)
{
    int pos, i, dst;
    int bit = (int) log2(N);
    double complex tmp;
    
    for(pos = 0; pos < N; pos++)
    {
        dst = 0;
        for(i = 0; i < bit; i++)
            dst += ((pos >> i) & 1) << (bit - 1 - i); // 地址下标的二进制形式取倒序 
        
        if (dst > pos)
        {
            tmp = x[pos];
            x[pos] = x[dst];
            x[dst] = tmp;
        }
    }
    
}