前面写过关于傅里叶算法的应用例子.
基于傅里叶变换的音频重采样算法 (附完整 c 代码)
当然也就是举个例子, 主要是学习傅里叶变换.
这个重采样思路还有点瑕疵,
稍微改一下, 就可以支持多通道, 以及提升性能.
当然思路很简单, 就是切分, 合并.
留个作业哈.
本文不讲过多的算法思路, 傅里叶变换的各种变种,
绝大多数是为提升性能, 支持任意长度而作.
当然各有所长,
当时提到参阅整理的算法:
- https://github.com/cpuimage/StockhamFFT
- https://github.com/cpuimage/uFFT
- https://github.com/cpuimage/BluesteinCrz
- https://github.com/cpuimage/fftw3
例如 :
https://github.com/cpuimage/StockhamFFT 是优化速度,
https://github.com/cpuimage/BluesteinCrz 是支持任意长度,
https://github.com/cpuimage/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 这个前缀.
以上, 权当抛砖引玉.
来源: https://www.cnblogs.com/cpuimage/p/9452536.html