前面寫過關於傅里葉演算法的應用例子。 《基於傅里葉變換的音頻重採樣演算法 (附完整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這個首碼。
以上,權當拋磚引玉。
若有其他相關問題或者需求也可以郵件聯繫俺探討。
郵箱地址是:
[email protected]