地球物理频谱分析实验

#include

#include

void fft(float *xr, float *xi, int nr, float T);

int IBIT(int j, int m);

void main ( )

{

float x[500];

int i,j,k,m;

float ra[512],ri[512],a[500];

FILE *fp1,*fp2;

fp1=fopen(\"nose59-500\

fp2=fopen(\"gongkang\

for(i=0; i<25; i++) //选择第24条数据

{

for(m=0;m<512;m++)

{

ra[m]=0; //将ra[512],ri[512]数据设为0

ri[m]=0;

}

fread(&x[0],4,500,fp1);

}

for(j=0;j<500;j++)

{

ra[j]=x[j]; //将数组x的所有值赋于ra中

}

fft(ra,ri,9,1.0);

for(k=0;k<256;k++)

{

a[k]=(float)sqrt(ra[k]*ra[k]+ri[k]*ri[k]);

fprintf(fp2,\"%f\\n\ //将数组a中的值输出并保存在文件fp中

}

fclose(fp1);

fclose(fp2);

}

void fft(float *xr, float *xi, int nr, float T)

{

int i, j, k, l, n, n2, nr1, i1, j1, k2, k1n2;

float c, s, p, tr, ti, trc, tic, ars;

n=(int)(pow(2,nr));

n2=n/2;

nr1=nr-1;

k=0;

for(l=1; l<=nr; l++)

{

loop1: for(j=1; j<=n2; j++)

{

k2=k/(int)(pow(2,nr1));

p=(float)(IBIT(k2,nr));

ars=(float)(6.2831852*p/(float)(n));

c=(float)(cos(ars));

s=(float)(sin(ars));

k1n2=k+n2;

tr=xr[k1n2]*c+xi[k1n2]*s*T;

ti=xi[k1n2]*c-xr[k1n2]*s*T;

xr[k1n2]=xr[k]-tr;

xi[k1n2]=xi[k]-ti;

xr[k]=xr[k]+tr;

xi[k]=xi[k]+ti;

k++;

}

k+=n2;

if(k{ goto loop1; }

else

{

k=0;

nr1=nr1-1;

n2 /=2;

}

}

for(j=1; j<=n; j++)

{

i=IBIT(j-1,nr)+1;

if(i>j)

{

j1=j-1;

i1=i-1;

trc =xr[j1];

tic =xi[j1];

xr[j1]=xr[i1];

xi[j1]=xi[i1];

xr[i1]=trc;

xi[i1]=tic;

}

}

if(T<0.0)

{

for(j=0; j<=n; j++)

{

xr[j]/=n;

xi[j]/=n;

}

}

}

int IBIT(int j, int m)

{

int i, it, j1, j2;

it=0;

j1=j;

for(i=1; i<=m; i++)

{

j2=j1/2;

it=it*2+(j1-2*j2);

j1=j2;

}

return it;

}