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

https://github.com/cpuimage/StockhamFFT

https://github.com/cpuimage/uFFT

https://github.com/cpuimage/BluesteinCrz

https://github.com/cpuimage/fftw3

Stockham 是优化速度,

BluesteinCrz 是支持任意长度，

uFFT是经典实现。

```#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 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];
}
}
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++) {
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 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];
}
}
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++) {
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 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];
}
}
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++) {
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 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];
}
}
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++) {
f4[k],
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:
break;
case 5:
break;
case 6:
break;
case 7:
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:
break;
case 5:
break;
case 6:
break;
case 7:
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

0 条评论

• ### 不用第三方解码库取得图片宽高 附完整C++算法实现代码

在特定的应用场景下，有时候我们只是想获取图片的宽高， 但不想通过解码图片才取得这个信息。 预先知道图片的宽高信息，进而提速图片加载，预处理等相关操作以提升体验。...

• ### 快速均值模糊算法

前段时间在网上看到一个快速均值模糊算法，性能很不错。 源博客： http://www.lellansin.com/super-fast-blur-%E6%A8%...

• ### 人脸姿态校正算法 附完整C++示例代码

那么假如一张图片只有一个人脸，其实很好判断，通过眼睛的位置的坐标，根据两眼的直线角度，

• ### 乱炖数据之2700余篇“简书交友”专题文章数据的花式玩法

简书上有个“简书交友”专题，经常会有人写些自己的情况、贴贴自己的照片然后投稿到这一专题，有介绍的比较详细的比如下图所示（侵删），较为规整和全面；

• ### Hadoop基础教程-第12章 Hive：进阶（12.4 Hive Metastore）（草稿）

第12章 Hive：进阶 12.4 Hive Metastore 12.4.1 三种配置方式 Hive Metastore有三种配置方式，分别是： Embedd...

• ### 江湖急诏令：腾讯数据库王者挑战赛赏金万两募英豪！

腾讯云数据库王者挑战赛开始了 可以花几分钟参加比赛免费将☟☟抱回家！ MacBook/iPhone 11/AirPods 25台Kindle 8万元腾讯云...

• ### 【报名中】数据库大咖们与你聊聊云上实践的那些事儿

导语：6月29日，深圳腾讯滨海大厦3F多功能厅，云+社区邀您参加《腾讯云数据库行业实战分享会》沙龙活动，腾讯云数据库携手微众银行、销售易、小程序·云开发，为您...

• ### 【手写文字识别】-JavaAPI示例代码

手写文字识别-JavaAPI示例代码 不知不觉手写文字识别百度已经开始邀测了。需要的小伙伴去申请了哦。申请方式加入文字识别群找PM。或者工单提交申请。都要说明自...

• ### IOS开发之自定义Button(集成三种回调模式)

前面在做东西的时候都用到了storyboard，在今天的代码中就纯手写代码自己用封装个Button。这个Button继承于UIView类，在封装的时候用上啦...