网上找了半天都找不到合适的代码，要么是太重量级了，比如 FFTW，要么就是不可用，写得垃圾。
所以自己写了一个。

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

void *fft_init   (int len);
void  fft_free   (void *c);
void  fft_execute(void *c, float *in, float *out);

typedef struct {
    int    N;
    float *W;
    float *T;
    int   *order;
} FFT_CONTEXT;

// r = c1 + c2
static inline void complex_add(float *r, float *c1, float *c2)
{
    r[0] = c1[0] + c2[0];
    r[1] = c1[1] + c2[1];
}

// r = c1 - c2
static inline void complex_sub(float *r, float *c1, float *c2)
{
    r[0] = c1[0] - c2[0];
    r[1] = c1[1] - c2[1];
}

// r = c1 * c2
static inline void complex_mul(float *r, float *c1, float *c2)
{
    r[0] = c1[0] * c2[0] - c1[1] * c2[1];
    r[1] = c1[1] * c2[0] + c1[0] * c2[1];
}

static int reverse_bits(int n)
{
    n = ((n & 0xAAAAAAAA) >> 1 ) | ((n & 0x55555555) << 1 );
    n = ((n & 0xCCCCCCCC) >> 2 ) | ((n & 0x33333333) << 2 );
    n = ((n & 0xF0F0F0F0) >> 4 ) | ((n & 0x0F0F0F0F) << 4 );
    n = ((n & 0xFF00FF00) >> 8 ) | ((n & 0x00FF00FF) << 8 );
    n = ((n & 0xFFFF0000) >> 16) | ((n & 0x0000FFFF) << 16);
    return n;
}

static void fft_execute_internal(FFT_CONTEXT *ctxt, float *data, int n, int w)
{
    int i;
    for (i=0; i<n/2; i++) {
        // C = (A + B)
        // D = (A - B) * W
        float  A[2] = { data[(0   + i) * 2 + 0], data[(0   + i) * 2 + 1] };
        float  B[2] = { data[(n/2 + i) * 2 + 0], data[(n/2 + i) * 2 + 1] };
        float *C    = &(data[(0   + i) * 2]);
        float *D    = &(data[(n/2 + i) * 2]);
        float *W    = &(ctxt->W[i*w*2]);
        float  T[2];
        complex_add(C, A, B);
        complex_sub(T, A, B);
        complex_mul(D, T, W);
    }

    n /= 2;
    w *= 2;
    if (n > 1) {
        fft_execute_internal(ctxt, data + 0    , n, w);
        fft_execute_internal(ctxt, data + n * 2, n, w);
    }
}

void *fft_init(int n)
{
    int i;

    FFT_CONTEXT *ctxt = (FFT_CONTEXT*)malloc(sizeof(FFT_CONTEXT));
    if (!ctxt) return NULL;

    ctxt->N     = n;
    ctxt->W     = (float*)malloc(n * sizeof(float) * 1);
    ctxt->T     = (float*)malloc(n * sizeof(float) * 2);
    ctxt->order = (int  *)malloc(n * sizeof(int  ) * 1);
    if (!ctxt->W || !ctxt->T || !ctxt->order) {
        fft_free(ctxt);
        return NULL;
    }

    for (i=0; i<ctxt->N/2; i++) {
        ctxt->W[i * 2 + 0] = cos(2 * M_PI * i / ctxt->N);
        ctxt->W[i * 2 + 1] =-sin(2 * M_PI * i / ctxt->N);
    }

    int shift = (int)(32 - log(n)/log(2));
    for (i=0; i<ctxt->N; i++) {
        ctxt->order[i] = (unsigned)reverse_bits(i) >> shift;
    }
    return ctxt;
}

void fft_free(void *c)
{
    FFT_CONTEXT *ctxt = (FFT_CONTEXT*)c;
    if (!ctxt) return;
    if (ctxt->W    ) free(ctxt->W    );
    if (ctxt->T    ) free(ctxt->T    );
    if (ctxt->order) free(ctxt->order);
    free(ctxt);
}

void fft_execute(void *c, float *in, float *out)
{
    int i;
    FFT_CONTEXT *ctxt = (FFT_CONTEXT*)c;
    memcpy(ctxt->T, in, sizeof(float) * 2 * ctxt->N);
    fft_execute_internal(ctxt, ctxt->T, ctxt->N, 1);
    for (i=0; i<ctxt->N; i++) {
        out[ctxt->order[i] * 2 + 0] = ctxt->T[i * 2 + 0];
        out[ctxt->order[i] * 2 + 1] = ctxt->T[i * 2 + 1];
    }
}

#if 1
static void dft(float *in, float *out, int n)
{
    int   i, j;
    float tmp[2];
    float *w = (float*)malloc(n*2*sizeof(float));
    for (i=0; i<n; i++) {
        w[i * 2 + 0] = cos(2 * M_PI * i / n);
        w[i * 2 + 1] =-sin(2 * M_PI * i / n);
    }
    for (j=0; j<n; j++) {
        out[0] = out[1] = 0;
        for (i=0; i<n; i++) {
            // out += in[j] * w[(j*i)%n]
            complex_mul(tmp, &in[i * 2], &w[((i * j) % n) * 2]);
            complex_add(out, out, tmp);
        }
        out += 2;
    }
    free(w);
}

int main(void)
{
    int i;
    float in [16 * 2] = { 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
    float out[16 * 2];
    void *fft = NULL;

    printf("\nin:\n");
    for (i=0; i<16; i++) {
        printf("( %9.6f, %9.6f )\n", in [i*2+0], in [i*2+1]);
    }
    printf("\n");

    fft = fft_init(16);
    fft_execute(fft, in, out);
    fft_free(fft);

    printf("\nout:\n");
    for (i=0; i<16; i++) {
        printf("( %9.6f, %9.6f )\n", out[i*2+0], out[i*2+1]);
    }
    printf("\n");
    
    dft(in, out, 16);
    printf("\nout:\n");
    for (i=0; i<16; i++) {
        printf("( %9.6f, %9.6f )\n", out[i*2+0], out[i*2+1]);
    }
    printf("\n");
}
#endif

[此贴子已经被作者于2016-5-26 16:09编辑过]

基二频率抽取法，采用了递归，核心代码很简单，核心就是频率抽取的蝶形运算：
C = (A + B)
D = (A - B) * W

fft_execute_internal 采用了递归，将当前输入序列用蝶形运算计算完后，然后再分为两段，递归的调用 fft_execute_internal 进行处理，由此可见递归的简洁高效。

fft_init 中主要对 W 和逆序表 order 进行了计算和保存，以空间换时间。

测试程序中 dft 函数是离散傅里叶变换的数学公式算法，相对 fft 当然是速度慢，但是简单容易理解，也是用来对 fft 进行对比测试验证。

所有的代码，除了内部采用了 FFT_CONTEXT 这样自定义的数据类型，外部看来，全部都是 C 语言的基本数据类型，接口清晰简单，容易理解，也方便移植。

对外部暴露的接口足够简单：
void *fft_init   (int len);
void  fft_free   (void *c);
void  fft_execute(void *c, float *in, float *out);

学到了什么：
算法和理论是工程实践的指导，但是算法和理论必须转换为代码才有实际意义
尽量使用 C 语言原生的数据类型，能不自定义数据类型就不要定义，一旦定义了新了数据类型，就增加了理解、阅读和编写代码的负担（我绝对记得住 float，但是如果是自定义的“复数类型”，到底是 complex 还是 Complex 还是 COMPLEX 还是 complet_t 呢，用别人的库，库里面自定义的数据类型，真是考验你的记性）
尽量对外暴露简单直接的接口，方便移植和使用
有时候自己造轮子还是值得的，好处多多，相对于使用 FFTW 库，我这个代码足够简洁，而且性能也能满足需求

[此贴子已经被作者于2016-5-27 09:45编辑过]