经典傅里叶算法小集合 附完整c代码

时间：2024-01-20 11:26:33

前面写过关于傅里叶算法的应用例子。

《基于傅里叶变换的音频重采样算法 (附完整c代码)》

当然也就是举个例子，主要是学习傅里叶变换。

这个重采样思路还有点瑕疵，

稍微改一下，就可以支持多通道，以及提升性能。

当然思路很简单，就是切分，合并。

留个作业哈。

本文不讲过多的算法思路，傅里叶变换的各种变种，

绝大多数是为提升性能，支持任意长度而作。

当然各有所长，

当时提到参阅整理的算法：

https://github.com/cpuimage/StockhamFFT

https://github.com/cpuimage/uFFT

https://github.com/cpuimage/BluesteinCrz

https://github.com/cpuimage/fftw3

例如 :

Stockham 是优化速度,

BluesteinCrz 是支持任意长度，

uFFT是经典实现。

当然，各有利弊，精度也不一。

最近一直对傅里叶算法没放手。

还是想要抽点时间，不依赖第三方库，实现一份不差于fftw的算法，

既要保证精度，又要保证性能，同时还要支持任意长度。

目前还在进行中，目前项目完成了45%左右。

越是学习，看的资料林林总总，越觉得傅里叶变换的应用面很广。

花点时间，采用纯c ，实现了经典的傅里叶算法，

调整代码逻辑，慢慢开始有点清晰了。

前人栽树后人乘凉，为了学习方便，

把本人用纯c实现的经典傅里叶算法开源出来给大家学习。

算法逻辑写得简洁明了，我也是尽力了。

当然，可能还有更好的实现思路，更多的改进算法。

不过，我的目的更多是便于学习和理解算法。

希望能帮助到一些也在学习傅里叶变换算法的同学。

贴上完整算法代码：

#include <stdio.h>

#include <stdbool.h>

#include <math.h>

#include <stddef.h>

#include <string.h>

#include <stdlib.h>

#ifndef M_PI

#define M_PI 3.14159265358979323846f

#endif

typedef struct {

    float real, imag;

} cmplx;

cmplx cmplx_mul_add(const cmplx c, const cmplx a, const cmplx b) {

    const cmplx ret = {

            (a.real * b.real) + c.real - (a.imag * b.imag),

            (a.imag * b.real) + (a.real * b.imag) + c.imag

    };

    return ret;

}

void fft_Stockham(cmplx *input, cmplx *output, size_t n, int flag) {

    size_t half = n >> ;

    cmplx *tmp = (cmplx *) calloc(sizeof(cmplx), n);

    cmplx *y = (cmplx *) calloc(sizeof(cmplx), n);

    memcpy(y, input, sizeof(cmplx) * n);

    for (size_t r = half, l = ; r >= ; r >>= ) {

        cmplx *tp = y;

        y = tmp;

        tmp = tp;

        float factor_w = -flag * M_PI / l;

        cmplx w = {cosf(factor_w), sinf(factor_w)};

        cmplx wj = {, };

        for (size_t j = ; j < l; j++) {

            size_t jrs = j * (r << );

            for (size_t k = jrs, m = jrs >> ; k < jrs + r; k++) {

                const cmplx t = {(wj.real * tmp[k + r].real) - (wj.imag * tmp[k + r].imag),

                                 (wj.imag * tmp[k + r].real) + (wj.real * tmp[k + r].imag)};

                y[m].real = tmp[k].real + t.real;

                y[m].imag = tmp[k].imag + t.imag;

                y[m + half].real = tmp[k].real - t.real;

                y[m + half].imag = tmp[k].imag - t.imag;

                m++;

            }

            const float t = wj.real;

            wj.real = (t * w.real) - (wj.imag * w.imag);

            wj.imag = (wj.imag * w.real) + (t * w.imag);

        }

        l <<= ;

    }

    memcpy(output, y, sizeof(cmplx) * n);

    free(tmp);

    free(y);

}

void fft_radix3(cmplx *in, cmplx *result, size_t n, int flag) {

    if (n < ) {

        memcpy(result, in, sizeof(cmplx) * n);

        return;

    }

    size_t radix = ;

    size_t np = n / radix;

    cmplx *res = (cmplx *) malloc(sizeof(cmplx) * n);

    cmplx *f0 = res;

    cmplx *f1 = f0 + np;

    cmplx *f2 = f1 + np;

    for (size_t i = ; i < np; i++) {

        for (size_t j = ; j < radix; j++) {

            res[i + j * np] = in[radix * i + j];

        }

    }

    fft_radix3(f0, f0, np, flag);

    fft_radix3(f1, f1, np, flag);

    fft_radix3(f2, f2, np, flag);

    float wexp0 = - * (float) M_PI * (flag) / (float) (n);

    cmplx wt = {cosf(wexp0), sinf(wexp0)};

    cmplx w0 = {, };

    for (size_t i = ; i < np; i++) {

        const float w0r = w0.real;

        w0.real = (w0r * wt.real) - (w0.imag * wt.imag);

        w0.imag = (w0.imag * wt.real) + (w0r * wt.imag);

    }

    cmplx w = {, };

    for (size_t j = ; j < radix; j++) {

        cmplx wj = w;

        for (size_t k = ; k < np; k++) {

            result[k + j * np] = cmplx_mul_add(f0[k], cmplx_mul_add(f1[k], f2[k], wj), wj);

            const float wjr = wj.real;

            wj.real = (wjr * wt.real) - (wj.imag * wt.imag);

            wj.imag = (wj.imag * wt.real) + (wjr * wt.imag);

        }

        const float wr = w.real;

        w.real = (wr * w0.real) - (w.imag * w0.imag);

        w.imag = (w.imag * w0.real) + (wr * w0.imag);

    }

    free(res);

}

void fft_radix5(cmplx *x, cmplx *result, size_t n, int flag) {

    if (n < ) {

        memcpy(result, x, sizeof(cmplx) * n);

        return;

    }

    size_t radix = ;

    size_t np = n / radix;

    cmplx *res = (cmplx *) calloc(sizeof(cmplx), n);

    cmplx *f0 = res;

    cmplx *f1 = f0 + np;

    cmplx *f2 = f1 + np;

    cmplx *f3 = f2 + np;

    cmplx *f4 = f3 + np;

    for (size_t i = ; i < np; i++) {

        for (size_t j = ; j < radix; j++) {

            res[i + j * np] = x[radix * i + j];

        }

    }

    fft_radix5(f0, f0, np, flag);

    fft_radix5(f1, f1, np, flag);

    fft_radix5(f2, f2, np, flag);

    fft_radix5(f3, f3, np, flag);

    fft_radix5(f4, f4, np, flag);

    float wexp0 = - * (float) M_PI * (flag) / (float) (n);

    cmplx wt = {cosf(wexp0), sinf(wexp0)};

    cmplx w0 = {, };

    for (size_t i = ; i < np; i++) {

        const float w0r = w0.real;

        w0.real = (w0r * wt.real) - (w0.imag * wt.imag);

        w0.imag = (w0.imag * wt.real) + (w0r * wt.imag);

    }

    cmplx w = {, };

    for (size_t j = ; j < radix; j++) {

        cmplx wj = w;

        for (size_t k = ; k < np; k++) {

            result[k + j * np] = cmplx_mul_add(f0[k], cmplx_mul_add(f1[k], cmplx_mul_add(f2[k],

                                                                                         cmplx_mul_add(f3[k], f4[k],

                                                                                                       wj), wj), wj),

                                               wj);

            const float wjr = wj.real;

            wj.real = (wjr * wt.real) - (wj.imag * wt.imag);

            wj.imag = (wj.imag * wt.real) + (wjr * wt.imag);

        }

        const float wr = w.real;

        w.real = (wr * w0.real) - (w.imag * w0.imag);

        w.imag = (w.imag * w0.real) + (wr * w0.imag);

    }

    free(res);

}

void fft_radix6(cmplx *input, cmplx *output, size_t n, int flag) {

    if (n < ) {

        memcpy(output, input, sizeof(cmplx) * n);

        return;

    }

    size_t radix = ;

    size_t np = n / radix;

    cmplx *res = (cmplx *) calloc(sizeof(cmplx), n);

    cmplx *f0 = res;

    cmplx *f1 = f0 + np;

    cmplx *f2 = f1 + np;

    cmplx *f3 = f2 + np;

    cmplx *f4 = f3 + np;

    cmplx *f5 = f4 + np;

    for (size_t i = ; i < np; i++) {

        for (size_t j = ; j < radix; j++) {

            res[i + j * np] = input[radix * i + j];

        }

    }

    fft_radix6(f0, f0, np, flag);

    fft_radix6(f1, f1, np, flag);

    fft_radix6(f2, f2, np, flag);

    fft_radix6(f3, f3, np, flag);

    fft_radix6(f4, f4, np, flag);

    fft_radix6(f5, f5, np, flag);

    float wexp0 = - * (float) M_PI * (flag) / (float) (n);

    cmplx wt = {cosf(wexp0), sinf(wexp0)};

    cmplx w0 = {, };

    for (size_t i = ; i < np; i++) {

        const float w0r = w0.real;

        w0.real = (w0r * wt.real) - (w0.imag * wt.imag);

        w0.imag = (w0.imag * wt.real) + (w0r * wt.imag);

    }

    cmplx w = {, };

    for (size_t j = ; j < radix; j++) {

        cmplx wj = w;

        for (size_t k = ; k < np; k++) {

            output[k + j * np] = cmplx_mul_add(f0[k], cmplx_mul_add(f1[k], cmplx_mul_add(f2[k],

                                                                                         cmplx_mul_add(f3[k],

                                                                                                       cmplx_mul_add(

                                                                                                               f4[k],

                                                                                                               f5[k],

                                                                                                               wj), wj),

                                                                                         wj), wj), wj);

            const float wjr = wj.real;

            wj.real = (wjr * wt.real) - (wj.imag * wt.imag);

            wj.imag = (wj.imag * wt.real) + (wjr * wt.imag);

        }

        const float wr = w.real;

        w.real = (wr * w0.real) - (w.imag * w0.imag);

        w.imag = (w.imag * w0.real) + (wr * w0.imag);

    }

    free(res);

}

void fft_radix7(cmplx *x, cmplx *result, size_t n, int flag) {

    if (n < ) {

        memcpy(result, x, sizeof(cmplx) * n);

        return;

    }

    size_t radix = ;

    size_t np = n / radix;

    cmplx *res = (cmplx *) calloc(sizeof(cmplx), n);

    cmplx *f0 = res;

    cmplx *f1 = f0 + np;

    cmplx *f2 = f1 + np;

    cmplx *f3 = f2 + np;

    cmplx *f4 = f3 + np;

    cmplx *f5 = f4 + np;

    cmplx *f6 = f5 + np;

    for (size_t i = ; i < np; i++) {

        for (size_t j = ; j < radix; j++) {

            res[i + j * np] = x[radix * i + j];

        }

    }

    fft_radix7(f0, f0, np, flag);

    fft_radix7(f1, f1, np, flag);

    fft_radix7(f2, f2, np, flag);

    fft_radix7(f3, f3, np, flag);

    fft_radix7(f4, f4, np, flag);

    fft_radix7(f5, f5, np, flag);

    fft_radix7(f6, f6, np, flag);

    float wexp0 = - * (float) M_PI * (flag) / (float) (n);

    cmplx wt = {cosf(wexp0), sinf(wexp0)};

    cmplx w0 = {, };

    for (size_t i = ; i < np; i++) {

        const float w0r = w0.real;

        w0.real = (w0r * wt.real) - (w0.imag * wt.imag);

        w0.imag = (w0.imag * wt.real) + (w0r * wt.imag);

    }

    cmplx w = {, };

    for (size_t j = ; j < radix; j++) {

        cmplx wj = w;

        for (size_t k = ; k < np; k++) {

            result[k + j * np] = cmplx_mul_add(f0[k], cmplx_mul_add(f1[k], cmplx_mul_add(f2[k],

                                                                                         cmplx_mul_add(f3[k],

                                                                                                       cmplx_mul_add(

                                                                                                               f4[k],

                                                                                                               cmplx_mul_add(

                                                                                                                       f5[k],

                                                                                                                       f6[k],

                                                                                                                       wj),

                                                                                                               wj), wj),

                                                                                         wj), wj), wj);

            const float wjr = wj.real;

            wj.real = (wjr * wt.real) - (wj.imag * wt.imag);

            wj.imag = (wj.imag * wt.real) + (wjr * wt.imag);

        }

        const float wr = w.real;

        w.real = (wr * w0.real) - (w.imag * w0.imag);

        w.imag = (w.imag * w0.real) + (wr * w0.imag);

    }

    free(res);

}

void fft_Bluestein(cmplx *input, cmplx *output, size_t n, int flag) {

    size_t m =  << ((unsigned int) (ilogbf((float) ( * n - ))));

    if (m <  * n - ) {

        m <<= ;

    }

    cmplx *y = (cmplx *) calloc(sizeof(cmplx),  * m);

    cmplx *w = y + m;

    cmplx *ww = w + m;

    float a0 = (float) M_PI / n;

    w[].real = ;

    if (flag == -) {

        y[].real = input[].real;

        y[].imag = -input[].imag;

        for (size_t i = ; i < n; i++) {

            const float wexp = a0 * i * i;

            w[i].real = cosf(wexp);

            w[i].imag = sinf(wexp);

            w[m - i] = w[i];

            y[i].real = (input[i].real * w[i].real) - (input[i].imag * w[i].imag);

            y[i].imag = (-input[i].imag * w[i].real) - (input[i].real * w[i].imag);

        }

    } else {

        y[].real = input[].real;

        y[].imag = input[].imag;

        for (size_t i = ; i < n; i++) {

            const float wexp = a0 * i * i;

            w[i].real = cosf(wexp);

            w[i].imag = sinf(wexp);

            w[m - i] = w[i];

            y[i].real = (input[i].real * w[i].real) + (input[i].imag * w[i].imag);

            y[i].imag = (input[i].imag * w[i].real) - (input[i].real * w[i].imag);

        }

    }

    fft_Stockham(y, y, m, );

    fft_Stockham(w, ww, m, );

    for (size_t i = ; i < m; i++) {

        const float r = y[i].real;

        y[i].real = (r * ww[i].real) - (y[i].imag * ww[i].imag);

        y[i].imag = (y[i].imag * ww[i].real) + (r * ww[i].imag);

    }

    fft_Stockham(y, y, m, -);

    float scale = 1.0f / m;

    if (flag == -) {

        for (size_t i = ; i < n; i++) {

            output[i].real = ((y[i].real * w[i].real) + (y[i].imag * w[i].imag)) * scale;

            output[i].imag = -((y[i].imag * w[i].real) - (y[i].real * w[i].imag)) * scale;

        }

    } else {

        for (size_t i = ; i < n; i++) {

            output[i].real = ((y[i].real * w[i].real) + (y[i].imag * w[i].imag)) * scale;

            output[i].imag = ((y[i].imag * w[i].real) - (y[i].real * w[i].imag)) * scale;

        }

    }

    free(y);

}

size_t base(size_t n) {

    size_t t = n & (n - );

    if (t == ) {

        return ;

    }

    for (size_t i = ; i <= ; i++) {

        size_t n2 = n;

        while (n2 % i == ) {

            n2 /= i;

        }

        if (n2 == ) {

            return i;

        }

    }

    return n;

}

void FFT(cmplx *input, cmplx *output, size_t n) {

    memset(output, , sizeof(cmplx) * n);

    if (n < ) {

        memcpy(output, input, sizeof(cmplx) * n);

        return;

    }

    size_t p = base(n);

    switch (p) {

        case :

            fft_Stockham(input, output, n, );

            break;

        case :

            fft_radix3(input, output, n, );

            break;

        case :

            fft_radix5(input, output, n, );

            break;

        case :

            fft_radix6(input, output, n, );

            break;

        case :

            fft_radix7(input, output, n, );

            break;

        default:

            fft_Bluestein(input, output, n, );

            break;

    }

}

void IFFT(cmplx *input, cmplx *output, size_t n) {

    memset(output, , sizeof(cmplx) * n);

    if (n < ) {

        memcpy(output, input, sizeof(cmplx) * n);

        return;

    }

    size_t p = base(n);

    switch (p) {

        case :

            fft_Stockham(input, output, n, -);

            break;

        case :

            fft_radix3(input, output, n, -);

            break;

        case :

            fft_radix5(input, output, n, -);

            break;

        case :

            fft_radix6(input, output, n, -);

            break;

        case :

            fft_radix7(input, output, n, -);

            break;

        default: {

            fft_Bluestein(input, output, n, -);

            break;

        }

    }

    float scale = 1.0f / n;

    for (size_t i = ; i < n; i++) {

        output[i].real = output[i].real * scale;

        output[i].imag = output[i].imag * scale;

    }

}

int main() {

    printf("Fast Fourier Transform\n");

    printf("blog: http://cpuimage.cnblogs.com/\n");

    printf("A Simple and Efficient FFT Implementation in C");

    size_t N = ;

    cmplx *input = (cmplx *) calloc(sizeof(cmplx), N);

    cmplx *output = (cmplx *) calloc(sizeof(cmplx), N);

    for (size_t i = ; i < N; ++i) {

        input[i].real = i;

        input[i].imag = ;

    }

    for (size_t i = ; i < N; ++i) {

        printf("(%f %f) \t", input[i].real, input[i].imag);

    }

    for (int i = ; i < ; i++) {

        FFT(input, output, N);

    }

    printf("\n");

    IFFT(output, input, N);

    for (size_t i = ; i < N; ++i) {

        printf("(%f %f) \t", input[i].real, input[i].imag);

    }

    free(input);

    free(output);

    getchar();

    return ;

}

项目地址:

https://github.com/cpuimage/cpuFFT

想了好久都没想到取啥名字好，最后还是选择了cpu这个前缀。

以上，权当抛砖引玉。

若有其他相关问题或者需求也可以邮件联系俺探讨。

邮箱地址是: 
gaozhihan@vip.qq.com