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

您所在的位置:网站首页 傅里叶快速算法 经典傅里叶算法小集合 附完整c代码

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

2024-07-16 02:51| 来源: 网络整理| 查看: 265

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

《基于傅里叶变换的音频重采样算法 (附完整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 #include #include #include #include #include #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; 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


【本文地址】


今日新闻


推荐新闻


CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3