#include<stdio.h>
#include <math.h>
#include<malloc.h>
#include<string.h>
#include<stdlib.h>
#define R 61
#define S 61
#define N 195
#define dt 0.04
#define PI 3.141593
int main()
{
int r,s,n,k,nfft;
float w,t;
float ***p,***p_re,***p_im,***p0,***p0_re,***p0_im,***m,***m_re,***m_im;
FILE *fp1,*fp2;
void fft(double sr[61][195][61],double sx[61][195][61],int m0,int inv);
float***create_3_array(int m,int n,int t);
p=create_3_array(r,n,s);
p_re=create_3_array(r,n,s);
p_im=create_3_array(r,n,s);
p0=create_3_array(r,n,s);
p0_re=create_3_array(r,n,s);
p0_im=create_3_array(r,n,s);
m=create_3_array(r,n,s);
m_re=create_3_array(r,n,s);
m_im=create_3_array(r,n,s);
/*讀資料到p[][][]*/
if (fp1=fopen("shot_fullwave_new44_use-6_new_sun.dat", "rb")==NULL)
{
printf("Open failed\n");
exit(0);
}
else
for(s=0;s<=S; s++)
{
for(r=0;r<=R;r++)
{
for(n=0;n<=N;n++)
{
fread(&p[r][n][s],sizeof(float),1,fp1);
}
}
putchar('\n');
}
fclose(fp1);
/*計算fft次數*/
k=log(N)/log(2);
if(N>pow(2,k))
k=k+1;
nfft=pow(2,k);
printf("nfft=%d k=%d\n",nfft,k);
/*為p_re、p_im賦初值*/
for(s=0;s<=S; s++)
{
for(r=0;r<=R;r++)
{
for(n=0;n<=N;n++)
{
p_re[r][n][s]=p[r][n][s];
p_im[r][n][s]=0.0;
}
}
}
/*對p_re、p_im對時間維度做fft、ifft*/
fft(p_re,p_im,k,1);
fft(p_re,p_im,k,-1);
/*讀入資料*/
fp2=fopen("output_m.dat","wb");
if(fp2==NULL)
printf("NULL2");
else
{
for(s=0;s<S;s++)
{
for(r=0;r<R;r++)
{
for(n=0;n<N;n++)
fwrite(&p[r][n][s],sizeof(float),1,fp2);
}
}
}
return 0;
}
/*fft ifft*/
void fft(double sr[61][195][61],double sx[61][195][61],int m0,int inv)
{
int r,s,i,j,lm,li,k,lmx,np,lix,mm2;
double t1,t2,c,f,cv,st,ct;
if(m0<0)
return;
lmx=1;
for(i=1;i<=m0;++i)
lmx+=lmx;
cv=2.0*PI/(double)lmx;
ct=cos(cv);
st=-inv*sin(cv);
np=lmx;
mm2=m0-2;
for(s=0;s<S;s++)
{
for(r=0;r<R;r++)
{
for(i=1;i<=mm2;++i)
{
lix=lmx;
lmx/=2;
c=ct;
f=st;
for(li=0;li<np;li+=lix)
{
j=li;
k=j+lmx;
t1=sr[r][j][s]-sr[r][k][s];
t2=sx[r][j][s]-sx[r][k][s];
sr[r][j][s]+=sr[r][k][s];
sx[r][j][s]+=sx[r][k][s];
sr[r][k][s]=t1;
sx[r][k][s]=t2;
++j;
++k;
t1=sr[r][j][s]-sr[r][k][s];
t2=sx[r][j][s]-sx[r][k][s];
sr[r][j][s]+=sr[r][k][s];
sx[r][j][s]+=sx[r][k][s];
sr[r][k][s]=c*t1-f*t2;
sx[r][k][s]=f*t1+c*t2;
}
for(lm=2;lm<lmx;++lm)
{
cv=c;
c=ct*c-st*f;
f=st*cv+ct*f;
for(li=0;li<np;li+=lix)
{
j=li+lm;
k=lmx+j;
t1=sr[r][j][s]-sr[r][k][s];
t2=sx[r][j][s]-sx[r][k][s];
sr[r][j][s]+=sr[r][k][s];
sx[r][j][s]+=sx[r][k][s];
sr[r][k][s]=c*t1-f*t2;
sx[r][k][s]=f*t1+c*t2;
}
}
cv=ct;
ct=2.0*ct*ct-1.0;
st=2.0*st*cv;
}
if(m0>=2)
for(li=0;li<np;li+=4)
{
j=li;
k=j+2;
t1=sr[r][j][s]-sr[r][k][s];
t2=sx[r][j][s]-sx[r][k][s];
sr[r][j][s]+=sr[r][k][s];
sx[r][j][s]+=sx[r][k][s];
sr[r][k][s]=t1;
sx[r][k][s]=t2;
++j;
++k;
t1=sr[r][j][s]-sr[r][k][s];
t2=sx[r][j][s]-sx[r][k][s];
sr[r][j][s]+=sr[r][k][s];
sx[r][j][s]+=sx[r][k][s];
sr[r][k][s]=inv*t2;
sx[r][k][s]=-inv*t1;
}
for(li=0;li<np;li+=2)
{
j=li;
k=j+1;
t1=sr[r][j][s]-sr[r][k][s];
t2=sx[r][j][s]-sx[r][k][s];
sr[r][j][s]+=sr[r][k][s];
sx[r][j][s]+=sx[r][k][s];
sr[r][k][s]=t1;
sx[r][k][s]=t2;
}
lmx=np/2;
j=0;
for(i=1;i<np-1;++i)
{
k=lmx;
while(k<=j)
{
j-=k;
k/=2;
}
j+=k;
if(i<j)
{
t1=sr[r][j][s];
sr[r][j][s]=sr[r][i][s];
sr[r][i][s]=t1;
t1=sx[r][i][s];
sx[r][j][s]=sx[r][i][s];
sx[r][i][s]=t1;
}
}
if(inv!=-1)
return;
t1=1.0/np;
for(i=0;i<np;++i)
{
sr[r][i][s]*=t1;sx[r][i][s]*=t1;
}
}
}
}
float***create_3_array(int m,int n,int t)
{
float***array=NULL;
int i,j;
array=(float***)malloc(sizeof(float**)*m);
for(i=0;i<m;i++)
{
array[i]=(float**)malloc(sizeof(float*)*n);
for(j=0;j<n;j++)
array[i][j]=(float*)malloc(sizeof(float)*t);
}
return array;
}
uj5u.com熱心網友回復:
int main()
{
int r,s,n,k,nfft;
float w,t;
float ***p,***p_re,***p_im,***p0,***p0_re,***p0_im,***m,***m_re,***m_im;
FILE *fp1,*fp2;
void fft(double sr[61][195][61],double sx[61][195][61],int m0,int inv);
float***create_3_array(int m,int n,int t);
p=create_3_array(r,n,s);
p_re=create_3_array(r,n,s);
p_im=create_3_array(r,n,s);
p0=create_3_array(r,n,s);
p0_re=create_3_array(r,n,s);
p0_im=create_3_array(r,n,s);
m=create_3_array(r,n,s);
m_re=create_3_array(r,n,s);
m_im=create_3_array(r,n,s);
/*讀資料到p[][][]*/
//if (fp1=fopen("shot_fullwave_new44_use-6_new_sun.dat", "rb")==NULL)
if ((fp1=fopen("shot_fullwave_new44_use-6_new_sun.dat", "rb"))==NULL)
{
printf("Open failed\n");
exit(0);
}
//else
int i, j;
p = (float ***)malloc(sizeof(float **) * (S+1));
for (i = 0; i <= S; i++) {
p[i] = (float **)malloc(sizeof(float *) * (R+1));
for (j = 0; j <= R; j++) {
p[i][j] = (float *)malloc(sizeof(float) * N);
}
}
for(s=0;s<=S; s++)
{
for(r=0;r<=R;r++)
{
for(n=0;n<=N;n++)
{
fread(&p[r][n][s],sizeof(float),1,fp1);
}
}
putchar('\n');
}
fclose(fp1);
/*計算fft次數*/
k=log(N)/log(2);
if(N>pow(2,k))
k=k+1;
nfft=pow(2,k);
printf("nfft=%d k=%d\n",nfft,k);
/*為p_re、p_im賦初值*/
for(s=0;s<=S; s++)
{
for(r=0;r<=R;r++)
{
for(n=0;n<=N;n++)
{
p_re[r][n][s]=p[r][n][s];
p_im[r][n][s]=0.0;
}
}
}
/*對p_re、p_im對時間維度做fft、ifft*/
fft(p_re,p_im,k,1);
fft(p_re,p_im,k,-1);
/*讀入資料*/
fp2=fopen("output_m.dat","wb");
if(fp2==NULL)
printf("NULL2");
else
{
for(s=0;s<S;s++)
{
for(r=0;r<R;r++)
{
for(n=0;n<N;n++)
fwrite(&p[r][n][s],sizeof(float),1,fp2);
}
}
}
return 0;
}
供參考~
打開檔案有問題,會導致例外,使用野指標也會導致例外,運行之后如果還出現問題,建議繼續查查;
uj5u.com熱心網友回復:
檔案打開失敗,會報錯,所以//exit(0), 其他問題見注釋,供參考:#include<stdlib.h>
#include<stdio.h>
#include <math.h>
#include<malloc.h>
#include<string.h>
#define R 61
#define S 61
#define N 195
#define dt 0.04
#define PI 3.141593
//void fft(double sr[61][195][61],double sx[61][195][61],int m0,int inv);
void fft(float ***sr,float ***sx,int m0,int inv);
float ***create_3_array(int m,int n,int t);
int main()
{
int r,s,n,k,nfft;
float w,t;
float ***p,***p_re,***p_im;//***p0,***p0_re,***p0_im,***m,***m_re,***m_im;
FILE *fp1,*fp2;
p = create_3_array(R,N,S);//(r,n,s);
p_re=create_3_array(R,N,S);//(r,n,s);
p_im=create_3_array(R,N,S);//(r,n,s);
//這些沒用到
//p0=create_3_array(R,N,S);//(r,n,s);
//p0_re=create_3_array(R,N,S);//(r,n,s);
//p0_im=create_3_array(R,N,S);//(r,n,s);
//m=create_3_array(R,N,S);//(r,n,s);
//m_re=create_3_array(R,N,S);//(r,n,s);
//m_im=create_3_array(R,N,S);//(r,n,s);
/*讀資料到p[][][]*/
fp1=fopen("shot_fullwave_new44_use-6_new_sun.dat", "rb");
if (fp1 == NULL)
{
printf("Open failed\n");
//exit(0);
}
else
for(r=0;r<=R; r++)//for(s=0;s<=S; s++)
{
for(n=0;n<=N;n++)//for(r=0;r<=R;r++)
{
for(s=0;s<=S; s++)//for(n=0;n<=N;n++)
{
fread(&p[r][n][s],sizeof(float),1,fp1);
}
}
putchar('\n');
}
fclose(fp1);
/*計算fft次數*/
k=log(N)/log(2);
if(N>pow(2,k))
k=k+1;
nfft=pow(2,k);
printf("nfft=%d k=%d\n",nfft,k);
/*為p_re、p_im賦初值*/
for(r=0;r<=R;r++)//for(s=0;s<=S; s++)
{
for(n=0;n<=N;n++)//for(r=0;r<=R;r++)
{
for(s=0;s<=S; s++)//for(n=0;n<=N;n++)
{
p_re[r][n][s]=p[r][n][s];
p_im[r][n][s]=0.0;
}
}
}
/*對p_re、p_im對時間維度做fft、ifft*/
fft(p_re,p_im,k,1);
fft(p_re,p_im,k,-1);
/*讀入資料*/
fp2=fopen("output_m.dat","wb");
if(fp2==NULL)
printf("NULL2");
else
{
for(r=0;r<R;r++)//for(s=0;s<S;s++)
{
for(n=0;n<N;n++)//for(r=0;r<R;r++)
{
for(s=0;s<S;s++)//for(n=0;n<N;n++)
fwrite(&p[r][n][s],sizeof(float),1,fp2);
} //這里保存的是 p[r][n][s] ?????
}
}
system("pause");
return 0;
}
/*fft ifft*/
void fft(float ***sr,float ***sx,int m0,int inv)
{
int r,s,i,j,lm,li,k,lmx,np,lix,mm2;
double t1,t2,c,f,cv,st,ct;
if(m0<0)
return;
lmx=1;
for(i=1;i<=m0-1;++i)//for(i=1;i<=m0;++i)陣列越界
lmx+=lmx;
cv=2.0*PI/(double)lmx;
ct=cos(cv);
st=-inv*sin(cv);
np=lmx;
mm2=m0-2;
for(s=0;s<S;s++)
{
for(r=0;r<R;r++)
{
for(i=1;i<=mm2;++i)
{
lix=lmx;
lmx/=2;
c=ct;
f=st;
for(li=0;li<np;li+=lix)
{
j=li;
k=j+lmx;
t1=sr[r][j][s]-sr[r][k][s];
t2=sx[r][j][s]-sx[r][k][s];
sr[r][j][s]+=sr[r][k][s];
sx[r][j][s]+=sx[r][k][s];
sr[r][k][s]=t1;
sx[r][k][s]=t2;
++j;
++k;
t1=sr[r][j][s]-sr[r][k][s];
t2=sx[r][j][s]-sx[r][k][s];
sr[r][j][s]+=sr[r][k][s];
sx[r][j][s]+=sx[r][k][s];
sr[r][k][s]=c*t1-f*t2;
sx[r][k][s]=f*t1+c*t2;
}
for(lm=2;lm<lmx;++lm)
{
cv=c;
c=ct*c-st*f;
f=st*cv+ct*f;
for(li=0;li<np;li+=lix)
{
j=li+lm;
k=lmx+j;
t1=sr[r][j][s]-sr[r][k][s];
t2=sx[r][j][s]-sx[r][k][s];
sr[r][j][s]+=sr[r][k][s];
sx[r][j][s]+=sx[r][k][s];
sr[r][k][s]=c*t1-f*t2;
sx[r][k][s]=f*t1+c*t2;
}
}
cv=ct;
ct=2.0*ct*ct-1.0;
st=2.0*st*cv;
}
if(m0>=2)
for(li=0;li<np;li+=4)
{
j=li;
k=j+2;
t1=sr[r][j][s]-sr[r][k][s];
t2=sx[r][j][s]-sx[r][k][s];
sr[r][j][s]+=sr[r][k][s];
sx[r][j][s]+=sx[r][k][s];
sr[r][k][s]=t1;
sx[r][k][s]=t2;
++j;
++k;
t1=sr[r][j][s]-sr[r][k][s];
t2=sx[r][j][s]-sx[r][k][s];
sr[r][j][s]+=sr[r][k][s];
sx[r][j][s]+=sx[r][k][s];
sr[r][k][s]=inv*t2;
sx[r][k][s]=-inv*t1;
}
for(li=0;li<np;li+=2)
{
j=li;
k=j+1;
t1=sr[r][j][s]-sr[r][k][s];
t2=sx[r][j][s]-sx[r][k][s];
sr[r][j][s]+=sr[r][k][s];
sx[r][j][s]+=sx[r][k][s];
sr[r][k][s]=t1;
sx[r][k][s]=t2;
}
lmx=np/2;
j=0;
for(i=1;i<np;++i)
{
k=lmx;
while(k<=j)
{
j-=k;
k/=2;
}
j+=k;
if(i<j)
{
t1=sr[r][j][s];
sr[r][j][s]=sr[r][i][s];
sr[r][i][s]=t1;
t1=sx[r][i][s];
sx[r][j][s]=sx[r][i][s];
sx[r][i][s]=t1;
}
}
if(inv!=-1)
return;
t1=1.0/np;
for(i=0;i<np;++i)
{
sr[r][i][s]*=t1;
sx[r][i][s]*=t1;
}
}
}
}
float ***create_3_array(int m,int n,int t)
{
float***array=NULL;
int i,j;
array=(float***)malloc(sizeof(float**)*(m+1));
for(i=0;i<=m;i++)
{
array[i]=(float**)malloc(sizeof(float*)*(n+1));
for(j=0;j<=n;j++)
array[i][j]=(float*)malloc(sizeof(float)*(t+1));
}
return array;
}
轉載請註明出處,本文鏈接:https://www.uj5u.com/houduan/274203.html
標籤:C語言
上一篇:求助大佬
