1周学FFT——第7天 Cooley-Tukey算法的C语言实现

Cooley-Tukey算法以发明者J. W. Cooley和John Tukey命名。Cooley-Tukey算法是最著名的FFT算法。它可以与其他DFT算法合并混用，比如将Cooley-Tukey算法与Rader算法或Bluestein算法合并使用，可以处理含有大质因数的情况（而不是填零凑基-2）。

Cooley-Tukey算法是一种递归式算法，最早由著名的数学小王子高斯发明（很难想象高斯在什么样的背景下展开对这一问题的讨论，还是仅仅出于数学考量），Cooley和Tukey在160年之后再次独立的发现了它。

以下是Cooley-Tukey算法的伪代码以及C语言实现，实现了第5天所描述的时间抽选奇偶分解基-2FFT算法。

伪代码

X0,...,N−1 ← ditfft2(x, N, s):                         基-2 FFT，x是采样序列，N是点数，s是步长
    if N = 1 then
        X0 ← x0                                        1点DFT，递归已经到底了
    else
        X0,...,N/2−1 ← ditfft2(x, N/2, 2s)             一半点数的偶数DFT（x0, x2s, x4s, ...）
        XN/2,...,N−1 ← ditfft2(x+s, N/2, 2s)           一半点数的奇数DFT（xs, xs+2s, xs+4s, ...）
        for k = 0 to N/2−1 do                          用蝶形运算，将两个子DFT合并为完整的DFT
            t ← Xk
            Xk ← t + exp(−2πi k/N) Xk+N/2
            Xk+N/2 ← t − exp(−2πi k/N) Xk+N/2
        end for
    end if

C语言实现

根据伪代码，可以写出N点基-2FFT算法的C语言实现。

首先定义和复数的数据结构MyComplex。

#define PI acos(-1)

typedef struct _MyComplex
{
    float re;
    float im;
}MyComplex;

然后定义几个复数运算函数，包括：

• complex_add()，complex_sub()和complex_mul()，这三个函数参与运算，以后需要对其做进一步的优化；
• complex_show()，complex_abs()不参与主要运算，可不进行优化。

void complex_add(MyComplex* cr, MyComplex* c1, MyComplex* c2)
{
    cr->re = c1->re + c2->re;
    cr->im = c1->im + c2->im;
}

void complex_sub(MyComplex* cr, MyComplex* c1, MyComplex* c2)
{
    cr->re = c1->re - c2->re;
    cr->im = c1->im - c2->im;
}

void complex_mul(MyComplex* cr, MyComplex* c1, MyComplex* c2)
{
    cr->re = c1->re * c2->re - c1->im * c2->im;
    cr->im = c1->re * c2->im + c1->im * c2->re;
}

void complex_show(char* name, MyComplex* c)
{
    printf("%s = %f + j%f\n", name, c->re, c->im);
}
void complex_abs(float* ca, MyComplex* c)
{
    *ca = c->re * c->re + c->im * c->im;
    *ca = sqrt(*ca);
}

然后是时间抽选奇偶分解基-2 FFT算法的主体：

void* dit_r2_fft(float* xn, int N, int stride, MyComplex* Xk)
{
    MyComplex X1k[N/2];
    MyComplex X2k[N/2];

    if (N==1)
    {
        // 递归终止
        Xk[0].re = xn[0];
        Xk[0].im = 0;
    }
    else
    {
        // 偶数N/2点DFT
        dit_r2_fft(xn, N/2, 2*stride, X1k);
        // 奇数N/2点DFT
        dit_r2_fft((float*)(xn+stride), N/2, 2*stride, X2k);

        // 蝶形运算
        for (int k=0;k<=N/2-1;k++)
        {
            MyComplex t;
            MyComplex WNk;

            // 这种求WNk的方法会消耗大量时间，必须优化
            WNk.re = cos(2*PI/N*k);
            WNk.im = -sin(2*PI/N*k);

            complex_mul(&t, &WNk, &X2k[k]);
            complex_add(&Xk[k], &X1k[k], &t);
            complex_sub(&Xk[k+N/2], &X1k[k], &t);
        }
    }
}

调用并且测试结果：

void main(void)
{
    int N = 16;       // FFT点数
    float fs = 16;    // 采样频率
    float dt = 1/fs;  // 采样间隔（周期）

    float xn[N];      // 采样信号序列
    MyComplex Xk[N];  // FFT结果，频谱序列

    float Xabs[N];    // 频率序列绝对值（模）

    // 制作采样序列，实际当中由采样电路和采样程序获得
    for (int i=0;i<N;i++)
    {
        // 1Hz与2Hz信号混合
        xn[i] = 5*sin(2*PI*1*dt*i) + 6*sin(2*PI*2*dt*i);
    }
    for (int i = 0; i < N; i++)
    {
        printf("xn[%d] = %f ", i, xn[i]);
    }
    printf("\n");

    // FFT
    dit_r2_fft(xn, N, 1, Xk);

    // 对频率序列Xk求绝对值
    for (int i = 0; i < N; i++)
    {
        complex_abs(&Xabs[i], &Xk[i]);
        Xabs[i] = Xabs[i]/N*2;
        printf("Xabs[%d] = %f\n", i, Xabs[i]);
    }

    // 简单绘制频谱
    for (int i = 0; i < N; i++)
    {
        printf("%f Hz\t", i*fs/N);
        for (int j=0; j <= Xabs[i]; j++)
        {
            printf("*");
        }
        printf("\n");
    }
}

运行结果如下图所示：

FFT运行结果 在龙芯1C+RT-Thread上运行的效果如图

习题：

1. 用C语言实现Cooley-Tukey算法；
2. 选择合适的点数与采样频率，利用所写的Cooley-Tukey算法，获取信号的频谱信息。