DIT-FFT算法c语言实现

#include <stdio.h> #include <math.h> #include "fft_lib.h" //复数结构体定义 struct complex{ double real; double imag; }; //复数加法 complex complex_add(complex a1,complex a2) { complex out; out.real=a1.real+a2.real; out.imag=a1.imag+a2.imag; return out; } //复数减法 complex complex_sub(complex a1,complex a2) { complex out; out.real=a1.real-a2.real; out.imag=a1.imag-a2.imag; return out; } //复数乘法 complex complex_multiply(complex a1, complex a2) { complex out; out.real = a1.real * a2.real - a1.imag * a2.imag; out.imag = a1.real * a2.imag + a1.imag * a2.real; return out; } //欧拉公式 :fft里面的旋转因子 complex euler_formula(double P,int N) { double x; complex out; x=2*PI/N*P; out.real=cos(x); out.imag=-sin(x); return out; } /* 说明： m!=0：返回数值为该数值是2的m次方 m=0：输入数值有问题 */ int charge_number(int N) { int m=0; for(;N%2==0;N=N/2){ if(N%2==0) m++; } if(N!=1) m=0; return m; } //倒序 void reverse_order(complex A[],int N) { int i,j; //正序和倒序位置 int lh,n1; complex index; int k; lh=N/2; n1=N-2; j=lh; //重新排序部分 for(i=1;i<n1;i++){ if(i<j){ index=A[j]; A[j]=A[i]; A[i]=index; } k=lh; while(j>=k){ j=j-k; k=k/2; } j=j+k; } } //蝶形运算 void count_fft(complex A[],int M,int N) { int L; //当前在那一级蝶形运算 int B; //跨度 int J,K; //旋转因子以及蝶形运算 int P; //旋转因子系数 complex W; //旋转因子 complex index; //计算过程中的中间值 complex index1; //用来暂存A[K]某一时刻的值 for(L=1;L<=M;L++){ B=(int)pow(2,L-1); for(J=0;J<=B-1;J++){ P=(int)(J*pow(2,M-L)); for(K=J;K<=N-1;K=K+2*B){ W=euler_formula(P,N); index=complex_multiply(A[K+B],W); index1=A[K]; //这个中间变量十分重要，如果不使用中间变量，下一行代码执行后，A[K]值就已经改变 后半段一定就是错的了 A[K]=complex_add(index1,index); A[K+B]=complex_sub(index1,index); } } } } void fft(complex A[],int M,int N) { int i=0; reverse_order(A,N); //倒序 count_fft(A,M,N); //蝶形运算 printf("FFT的结果为：\n"); for(i=0;i<N;i++) printf("第%d组的值为: %-f\t+ %fj \n",i+1,A[i].real,A[i].imag); }