Ooura FFT算法代码

下面是Ooura FFT算法的C代码，其中包括了初始化、正向变换和反向变换函数：

#include <math.h>

#define PI 3.14159265358979323846264338327L

typedef struct {
    int n;
    int *ip;
    double *w;
} ooura_fft_plan;

void ooura_fft_init(ooura_fft_plan *p, int n)
{
    p->n = n;
    p->ip = (int *)malloc(sizeof(int) * (2 + (int)sqrt(n / 2)));
    p->w = (double *)malloc(sizeof(double) * n / 2);
    int *ip = p->ip;
    double *w = p->w;

    int n1 = n - 1, nd2 = n / 2, i, j, k, n2, n4, n8;
    double c1 = 0.5, c2 = -0.5, t1, t2, t3, t4;

    if (ip[0] == 0) {
        n2 = n / 2;
        n4 = n / 4;
        n8 = n / 8;
        t1 = sin(PI / n);
        t2 = -2.0 * t1 * t1;
        t3 = sqrt(2.0) * t1;
        t4 = 2.0 * t3 * t1;
        w[0] = 1.0;
        w[n4] = c1 - c2;
        w[n2] = 0.0;
        w[n2 + n4] = c1 + c2;
        for (i = 1, j = n1; i < n8; i++, j--) {
            t1 = w[i - 1] * w[n4 + i - 1];
            t2 = w[n4 - i - 1] * w[n2 + i - 1];
            t3 = w[n2 - i - 1] * w[n2 + n4 - i - 1];
            t4 = w[n4 + i - 1] * w[n2 + n4 - i - 1];
            w[i] = t1 + c1 * t2 + c2 * t3;
            w[j] = c1 * t1 + t2 - t3 * c2;
            w[n2 - i] = c2 * t1 + t2 + c1 * t3;
            w[n2 + j] = c1 * t4 - c2 * t2;
            w[n2 + i] = c1 * t3 + t2 + c2 * t1;
            w[n2 + j - 1] = c2 * t4 + c1 * t2;
        }
        ip[0] = 1;
    }
    for (i = 1; i < (int)sqrt(n / 2); i++) {
        if (ip[i] == 0) {
            for (j = 0, k = 2 * i; k < n; j++, k += 2 * i) {
                ip[i + j] = k;
            }
        }
    }
}

void ooura_fft_forward(ooura_fft_plan *p, double *a)
{
    int n = p->n, *ip = p->ip;
    double *w = p->w, t1, t2, r1, r2;
    int j, k, m;

    m = p->n / 2;
    for (j = 1, k = 0; j < n - 1; j++) {
        int i = ip[j];
        if (i > j) {
            t1 = a[j * 2];
            t2 = a[j * 2 + 1];
            a[j * 2] = a[i * 2];
            a[j * 2 + 1] = a[i * 2 + 1];
            a[i * 2] = t1;
            a[i * 2 + 1] = t2;
        }
        k = m;
        while (k <= i) {
            i -= k;
            k /= 2;
        }
        if (k > j) {
            r1 = a[j * 2] + a[j * 2 + 1];
            r2 = a[j * 2] - a[j * 2 + 1];
            a[j * 2] = a[i * 2] + r2 * w[k - m];
            a[j * 2 + 1] = a[i * 2 + 1] - r1 * w[k - m];
            a[i * 2] -= r2 * w[k - m];
            a[i * 2 + 1] += r1 * w[k - m];
        } else {
            r1 = a[j * 2] + a[j * 2 + 1];
            r2 = a[j * 2] - a[j * 2 + 1];
            a[j * 2] = a[i * 2] + r1 * w[k - m];
            a[j * 2 + 1] = a[i * 2 + 1] + r2 * w[k - m];
            a[i * 2] -= r1 * w[k - m];
            a[i * 2 + 1] -= r2 * w[k - m];
        }
    }
    k = 2;
    for (m = n / 2; m > 2; m /= 2) {
        for (j = 0; j < m - 2; j += 2) {
            for (int i = j; i < n; i += k) {
                t1 = a[i * 2 + k] * w[(k / 2) - m];
                t2 = a[i * 2 + k + 1] * w[(k / 2) - m];
                a[i * 2 + k] = a[i * 2] - t1;
                a[i * 2 + k + 1] = a[i * 2 + 1] - t2;
                a[i * 2] += t1;
                a[i * 2 + 1] += t2;
            }
        }
        k *= 2;
    }
}

void ooura_fft_backward(ooura_fft_plan *p, double *a)
{
    int n = p->n, *ip = p->ip;
    double *w = p->w, t1, t2, r1, r2;
    int j, k, m;

    m = p->n / 2;
    for (j = 1, k = 0; j < n - 1; j++) {
        int i = ip[j];
        if (i > j) {
            t1 = a[j * 2];
            t2 = a[j * 2 + 1];
            a[j * 2] = a[i * 2];
            a[j * 2 + 1] = a[i * 2 + 1];
            a[i * 2] = t1;
            a[i * 2 + 1] = t2;
        }
        k = m;
        while (k <= i) {
            i -= k;
            k /= 2;
        }
        if (k > j) {
            r1 = a[j * 2] + a[j * 2 + 1];
            r2 = a[j * 2] - a[j * 2 + 1];
            a[j * 2] = a[i * 2] + r2 * w[k - m];
            a[j * 2 + 1] = a[i * 2 + 1] - r1 * w[k - m];
            a[i * 2] -= r2 * w[k - m];
            a[i * 2 + 1] += r1 * w[k - m];
        } else {
            r1 = a[j * 2] + a[j * 2 + 1];
            r2 = a[j * 2] - a[j * 2 + 1];
            a[j * 2] = a[i * 2] + r1 * w[k - m];
            a[j * 2 + 1] = a[i * 2 + 1] + r2 * w[k - m];
            a[i * 2] -= r1 * w[k - m];
            a[i * 2 + 1] -= r2 * w[k - m];
        }
    }
    k = 2;
    for (m = n / 2; m > 2; m /= 2) {
        for (j = 0; j < m - 2; j += 2) {
            for (int i = j; i < n; i += k) {
                t1 = a[i * 2 + k] * w[(k / 2) - m];
                t2 = a[i * 2 + k + 1] * w[(k / 2) - m];
                a[i * 2 + k] = a[i * 2] - t1;
                a[i * 2 + k + 1] = a[i * 2 + 1] - t2;
                a[i * 2] += t1;
                a[i * 2 + 1] += t2;
            }
        }
        k *= 2;
    }

    for (j = 0; j < n; j++) {
        a[j] /= n;
    }
}

使用方法：

1. 调用`ooura_fft_init`初始化FFT计算所需的内存
2. 将需要进行FFT变换的数据存储到`double`类型的数组中
3. 调用`ooura_fft_forward`进行正向变换，或`ooura_fft_backward`进行反向变换
4. 变换结果将存储在原始数据数组中，可以直接使用

此外，需要注意的是，该实现中的数组下标均从0开始，即第i个元素的下标为i-1，与其他实现可能有所不同。