TTF.rar_地震资料处理

#include "math.h" #include "malloc.h" #include "conio.h" #include "stdio.h" #include "string.h" #define pi (double)3.14159265359 #define m 6 //复数定义 typedef struct { double re; double im; }COMPLEX; //复数加运算 COMPLEX Add(COMPLEX c1,COMPLEX c2) { COMPLEX c; c.re=c1.re+c2.re; c.im=c1.im+c2.im; return c; } //复数乘运算 COMPLEX Mul(COMPLEX c1,COMPLEX c2) { COMPLEX c; c.re=c1.re*c2.re-c1.im*c2.im; c.im=c1.re*c2.im+c2.re*c1.im; return c; } //复数减运算 COMPLEX Sub(COMPLEX c1,COMPLEX c2) { COMPLEX c; c.re=c1.re-c2.re; c.im=c1.im-c2.im; return c; } //快速傅里叶变换 void FFT(COMPLEX *TD,COMPLEX *FD,int power) { int count; int i,j,k,bfsize,p; double angle; COMPLEX *W,*X1,*X2,*X; count=1<<power; W=(COMPLEX *)malloc(sizeof(COMPLEX)*count); X1=(COMPLEX *)malloc(sizeof(COMPLEX)*count); X2=(COMPLEX *)malloc(sizeof(COMPLEX)*count); for (i=0;i<count/2;i++) { angle=-i*pi*2/count; W[i].re=cos(angle); W[i].im=sin(angle); } memcpy(X1,TD,sizeof(COMPLEX)*count); for (k=0;k<power;k++) { for (j=0;j<1<<k;j++) { bfsize=1<<(power-k); for (i=0;i<bfsize/2;i++) { p=j*bfsize; X2[i+p]=Add(X1[i+p],X1[i+p+bfsize/2]); X2[i+p+bfsize/2]=Mul(Sub(X1[i+p],X1[i+p+bfsize/2]),W[i*(1<<k)]); } } X=X1; X1=X2; X2=X; } for (j=0;j<count;j++) { p=0; for (i=0;i<power;i++) { if (j&(1<<i)) p+=1<<(power-i-1); } FD[j]=X1[p]; } free(W); free(X1); free(X2); } //逆变换 void IFFT(COMPLEX *FD,COMPLEX *TD,int power) { int i,count; COMPLEX *x; count=1<<power; x=(COMPLEX *)malloc(sizeof(COMPLEX)*count); memcpy(x,FD,sizeof(COMPLEX)*count); for (i=0;i<count;i++) { x[i].im=-x[i].im; } FFT(x,TD,power); for(i=0;i<count;i++) { TD[i].re/=count; TD[i].im=-TD[i].im/count; } free(x); } void main() { int count; COMPLEX *TV,*FV; int i; count=1<<m; TV=(COMPLEX *)malloc(sizeof(COMPLEX)*count); FV=(COMPLEX *)malloc(sizeof(COMPLEX)*count); for (i=0;i<count;i++) { TV[i].re=i*3/5.0; TV[i].im=0; } FFT(TV,FV,m); IFFT(FV,TV,m); for (i=0;i<count;i++) { printf("x[%d]=%f+i %f\n",i,TV[i].re,TV[i].im); } free(TV); free(FV); getch(); }