首页 > > 程序设计 > C/C++ >

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

2018-12-04 07:13:17来源：博客园 阅读 ()

新老客户大回馈,云服务器低至5折

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

《基于傅里叶变换的音频重采样算法 (附完整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 >> 1;
    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 = 1; r >= 1; r >>= 1) {
        cmplx *tp = y;
        y = tmp;
        tmp = tp;
        float factor_w = -flag * M_PI / l;
        cmplx w = {cosf(factor_w), sinf(factor_w)};
        cmplx wj = {1, 0};
        for (size_t j = 0; j < l; j++) {
            size_t jrs = j * (r << 1);
            for (size_t k = jrs, m = jrs >> 1; 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 <<= 1;
    }
    memcpy(output, y, sizeof(cmplx) * n);
    free(tmp);
    free(y);
}

void fft_radix3(cmplx *in, cmplx *result, size_t n, int flag) {
    if (n < 2) {
        memcpy(result, in, sizeof(cmplx) * n);
        return;
    }
    size_t radix = 3;
    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 = 0; i < np; i++) {
        for (size_t j = 0; 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 = -2 * (float) M_PI * (flag) / (float) (n);
    cmplx wt = {cosf(wexp0), sinf(wexp0)};
    cmplx w0 = {1, 0};
    for (size_t i = 0; 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 = {1, 0};
    for (size_t j = 0; j < radix; j++) {
        cmplx wj = w;
        for (size_t k = 0; 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 < 2) {
        memcpy(result, x, sizeof(cmplx) * n);
        return;
    }
    size_t radix = 5;
    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 = 0; i < np; i++) {
        for (size_t j = 0; 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 = -2 * (float) M_PI * (flag) / (float) (n);
    cmplx wt = {cosf(wexp0), sinf(wexp0)};
    cmplx w0 = {1, 0};
    for (size_t i = 0; 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 = {1, 0};
    for (size_t j = 0; j < radix; j++) {
        cmplx wj = w;
        for (size_t k = 0; 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 < 2) {
        memcpy(output, input, sizeof(cmplx) * n);
        return;
    }
    size_t radix = 6;
    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 = 0; i < np; i++) {
        for (size_t j = 0; 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 = -2 * (float) M_PI * (flag) / (float) (n);
    cmplx wt = {cosf(wexp0), sinf(wexp0)};
    cmplx w0 = {1, 0};
    for (size_t i = 0; 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 = {1, 0};
    for (size_t j = 0; j < radix; j++) {
        cmplx wj = w;
        for (size_t k = 0; 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 < 2) {
        memcpy(result, x, sizeof(cmplx) * n);
        return;
    }
    size_t radix = 7;
    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 = 0; i < np; i++) {
        for (size_t j = 0; 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 = -2 * (float) M_PI * (flag) / (float) (n);
    cmplx wt = {cosf(wexp0), sinf(wexp0)};
    cmplx w0 = {1, 0};
    for (size_t i = 0; 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 = {1, 0};
    for (size_t j = 0; j < radix; j++) {
        cmplx wj = w;
        for (size_t k = 0; 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 = 1 << ((unsigned int) (ilogbf((float) (2 * n - 1))));
    if (m < 2 * n - 1) {
        m <<= 1;
    }
    cmplx *y = (cmplx *) calloc(sizeof(cmplx), 3 * m);
    cmplx *w = y + m;
    cmplx *ww = w + m;
    float a0 = (float) M_PI / n;
    w[0].real = 1;
    if (flag == -1) {
        y[0].real = input[0].real;
        y[0].imag = -input[0].imag;
        for (size_t i = 1; 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[0].real = input[0].real;
        y[0].imag = input[0].imag;
        for (size_t i = 1; 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, 1);
    fft_Stockham(w, ww, m, 1);
    for (size_t i = 0; 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, -1);
    float scale = 1.0f / m;
    if (flag == -1) {
        for (size_t i = 0; 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 = 0; 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 - 1);
    if (t == 0) {
        return 2;
    }
    for (size_t i = 3; i <= 7; i++) {
        size_t n2 = n;
        while (n2 % i == 0) {
            n2 /= i;
        }
        if (n2 == 1) {
            return i;
        }
    }
    return n;
}

void FFT(cmplx *input, cmplx *output, size_t n) {
    memset(output, 0, sizeof(cmplx) * n);
    if (n < 2) {
        memcpy(output, input, sizeof(cmplx) * n);
        return;
    }
    size_t p = base(n);
    switch (p) {
        case 2:
            fft_Stockham(input, output, n, 1);
            break;
        case 3:
            fft_radix3(input, output, n, 1);
            break;
        case 5:
            fft_radix5(input, output, n, 1);
            break;
        case 6:
            fft_radix6(input, output, n, 1);
            break;
        case 7:
            fft_radix7(input, output, n, 1);
            break;
        default:
            fft_Bluestein(input, output, n, 1);
            break;
    }
}

void IFFT(cmplx *input, cmplx *output, size_t n) {
    memset(output, 0, sizeof(cmplx) * n);
    if (n < 2) {
        memcpy(output, input, sizeof(cmplx) * n);
        return;
    }
    size_t p = base(n);
    switch (p) {
        case 2:
            fft_Stockham(input, output, n, -1);
            break;
        case 3:
            fft_radix3(input, output, n, -1);
            break;
        case 5:
            fft_radix5(input, output, n, -1);
            break;
        case 6:
            fft_radix6(input, output, n, -1);
            break;
        case 7:
            fft_radix7(input, output, n, -1);
            break;
        default: {
            fft_Bluestein(input, output, n, -1);
            break;
        }
    }
    float scale = 1.0f / n;
    for (size_t i = 0; 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 = 513;
    cmplx *input = (cmplx *) calloc(sizeof(cmplx), N);
    cmplx *output = (cmplx *) calloc(sizeof(cmplx), N);
    for (size_t i = 0; i < N; ++i) {
        input[i].real = i;
        input[i].imag = 0;
    }
    for (size_t i = 0; i < N; ++i) {
        printf("(%f %f) \t", input[i].real, input[i].imag);
    }
    for (int i = 0; i < 100; i++) {
        FFT(input, output, N);
    }
    printf("\n");
    IFFT(output, input, N);
    for (size_t i = 0; i < N; ++i) {
        printf("(%f %f) \t", input[i].real, input[i].imag);
    }
    free(input);
    free(output);
    getchar();
    return 0;
}

项目地址:

https://github.com/cpuimage/cpuFFT

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

以上，权当抛砖引玉。

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

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

标签：

版权申明：本站文章部分自网络，如有侵权，请联系：west999com@outlook.com
特别注意：本站所有转载文章言论不代表本站观点，本站所提供的摄影照片，插画，设计作品，如需使用，请与原作者联系，版权归原作者所有

上一篇：C语言二分查找

下一篇：【leetcode 简单】 第三十五题 环形链表

相关文章

IDC资讯： 主机资讯 注册资讯 托管资讯 vps资讯 网站建设

网站运营： 建站经验 策划盈利 搜索优化 网站推广 免费资源

网站联盟： 联盟新闻 联盟介绍 联盟点评 网赚技巧

行业资讯： 搜索引擎 网络游戏 电子商务 广告传媒

网络编程： Asp.Net编程 Asp编程 Php编程 Xml编程 Access Mssql Mysql 其它

服务器技术： Web服务器 Ftp服务器 Mail服务器 Dns服务器 安全防护

软件技巧： 其它软件 Word Excel Powerpoint Ghost Vista QQ空间 QQ FlashGet 迅雷

网页制作： FrontPages Dreamweaver Javascript css photoshop fireworks Flash

程序设计： Java技术 C/C++ VB delphi

网络知识： 网络协议 网络安全 网络管理 组网方案 Cisco技术

操作系统： Win2000 WinXP Win2003 Mac OS Linux FreeBSD

热门词条
最新资讯
热门关注
热门标签