利用c语言实现FFT运算

#include "pch.h" #include "tools.h" complex* x; //同时表示输入xn和输出Xk int N; //表示进行fft的点数 int bits; //log2(N),同时也可以表示fft级数 void fft_init(int fft_point, int xn_length, complex* xn) { int i = 0; N = fft_point; bits = (int)log2(N); x = (complex*)malloc(N * sizeof(complex)); //若序列长度小于fft点数，自动补零 for (i = 0;i < N;i++) { if (i <= xn_length) { x[i].real = xn[i].real; x[i].imaginary = xn[i].imaginary; } else { x[i].real = 0; x[i].imaginary = 0; } } } void ifft_init(int ifft_point, complex* xk) //利用fft算法的对称性实现ifft { int i; N = ifft_point; bits = (int)log2(N); x = (complex*)malloc(N * sizeof(complex)); for (i = 0;i < N;i++) { //利用fft求ifft,获得xk的共轭 x[i].real = xk[i].real; x[i].imaginary = -xk[i].imaginary; } } complex add(complex a, complex b) //复数加法 { complex result; result.real = a.real + b.real; result.imaginary = a.imaginary + b.imaginary; return result; } complex sub(complex a, complex b) //复数减法 { complex result; result.real = a.real - b.real; result.imaginary = a.imaginary - b.imaginary; return result; } complex multiplication(complex a, complex b) //复数乘法 { complex result; result.real = a.real*b.real - a.imaginary*b.imaginary; result.imaginary = a.real*b.imaginary + a.imaginary*b.real; return result; } double model(complex a) //求模值 { double result; double real2 = pow(a.real, 2); double imaginary2 = pow(a.imaginary, 2); result = pow(real2 + imaginary2, 0.5); return result; } complex conjugate(complex a) //求共轭 { complex result; result.real = a.real; result.imaginary = -a.imaginary; return result; } void reverse() //输入倒序 { complex temp; /*int bits = (int)log2(N);*/ char* flag = (char*)malloc(N * sizeof(char)); //用于标记是否被换过，若为1，则已经换过 int i, j, current_n, target_n; for (i = 0;i < N;i++) { flag[i] = 0; } for (i = 0;i < N;i++) { current_n = i; target_n = 0; //获取两个互为倒序的标号，换种思路 for (j = 0;j < bits;j++) { target_n = target_n + (int)((current_n % 2) * pow(2, bits - j - 1)); current_n /= 2; } current_n = i; //对应标号值互换 if (current_n != target_n && flag[current_n] != 1) { temp = x[current_n]; x[current_n] = x[target_n]; x[target_n] = temp; flag[target_n] = 1; } } free(flag); } complex w_builder(int m,int k) //产生旋转因子，m表示级数 { complex W; double base = pow(2, m); W.real = cos(2 * pi*k / base); W.imaginary = -sin(2 * pi*k / base); return W; } void butterfly(int x1_point, int x2_point, complex wn) //蝶形计算单元 { complex result1, result2, T; T = multiplication(x[x2_point], wn); result1 = add(x[x1_point], T); result2 = sub(x[x1_point], T); x[x1_point] = result1; x[x2_point] = result2; } void single_fft(int m) //进行单级的fft运算 { //首先进行分组 int point_distance = (int)pow(2, m - 1); //结点跨距 int group_distance = 2 * point_distance; //分组间隔 int group_count = N / group_distance; //分组数 int group_header = 0; int x1, x2; //参与计算两结点的标号 complex w; int i, j; //循环控制符 //分组 for (i = 0;i < group_count;i++) { //获取一组的标号范围 group_header = i * group_distance; //将每组分成两半，前一半均为x1,后一半均为x2 for (j = 0;j < group_distance / 2;j++) { w = w_builder(m, j); x1 = j + group_header; x2 = x1 + point_distance; butterfly(x1, x2, w); } } } complex* fft(int fft_point, int xn_length, complex* xn) //fft最终封装 { int i; //初始化 fft_init(fft_point, xn_length, xn); //输入倒序 reverse(); //fft运算 for (i = 1;i <= bits;i++) { single_fft(i); } return x; } complex* ifft(int xk_length, complex* xk) //ifft最终封装 { int i; //对输入序列进行处理 ifft_init(xk_length, xk); //套用fft算法 reverse(); for (i = 1;i <= bits;i++) { single_fft(i); } //对输出序列进行处理 for (i = 0;i < N;i++) { x[i] = conjugate(x[i]); //求共轭 x[i].real /= N; x[i].imaginary /= N; } return x; }