FFT相位差算法的C语言实现.docVIP

/*此算法为FFT相位差算法的C语言实现，经验证可用（网上有很多此类程序无法运行）；前半部分为FFT算法的子程序，也可用作其它用途*/ #include stdio.h #include stdlib.h #include math.h typedef struct {double real,imag;} complex; complex cexp(complex a) { complex z; z.real=exp(a.real)*cos(a.imag); z.imag=exp(a.real)*sin(a.imag); return(z); } void mcmpfft(complex x[],int n,int isign) { complex t,z,ce; double pisign; int mr,m,l,j,i,nn; for(i=1;i=16;i++) /*n must be power of 2 */ { nn=(int)pow(2,i); if(n==nn) break; } if(i16) { printf( N is not a power of 2 ! \n); return; } z.real=0.0; pisign=4*isign*atan(1.); /*pisign的值为+180度或-180度*/ mr=0; for(m=1;mn;m++) {l=n; while(mr+l=n) l=l/2; mr=mr%l+l; if(mr=m) continue; t.real=x[m].real;t.imag=x[m].imag; x[m].real=x[mr].real;x[m].imag=x[mr].imag; x[mr].real=t.real;x[mr].imag=t.imag; } l=1; while(1) { if(l=n) { if(isign==-1) /*isign=-1 For Forward Transform*/ return; for(j=0;jn;j++) /* Inverse Transform */ { x[j].real=x[j].real/n; x[j].imag=x[j].imag/n; } return; } for(m=0;ml;m++) /*完成当前级所有蝶形运算 */ { for(i=m;in;i=i+2*l)/*完成当前级相同W因子的所有蝶形运算*/ { z.imag=m*pisign/l; ce=cexp(z); t.real=x[i+l].real*ce.real-x[i+l].imag*ce.imag; t.imag=x[i+l].real*ce.imag+x[i+l].imag*ce.real; x[i+l].real=x[i].real-t.real; /*原位运算*/ x[i+l].imag=x[i].imag-t.imag; x[i].real=x[i].real+t.real; x[i].imag=x[i].imag+t.imag; } } l=2*l; /*确定下一级蝶形运算中W因子个数*/ } } void main() { int i,N=1024,k; float pi=3x[10000],A[10000],fi=3,f0=89,fs=450,max=0,o,oo,fi1,f00; complex s[2000]; for (i=0;iN;i++) { x[i]=1*sin(2*pi*f0*i