1楼
#include<stdio.h>
#include <math.h>
void dft(double x[],double y[],double a[],double b[],int n,int sign);
void main()
{FILE *fp;
int i,j;
float dt=0.002;
float r[200],t;
float ri[200];
float x[200];
double f0,zb[31],Br[31],Bi[31],Ar[31],Ai[31],za[31],temp[31],pi;
double alfa;
int n,m;
f0=35.0;
pi=3.14159;
if((fp=fopen("c_deltat.txt","r"))==NULL)
printf("Can't open file.\n");
else
{
for(i=0;i<200;i++)
{ fscanf(fp,"%f",&r[i]);
}
}
//--------------------------------------------//
alfa=(f0*f0*log(1.8))*2.0;
printf("%lf",alfa);
for(i=0;i<31;i++)
{zb[i]=sin(2*pi*f0*i*dt)*exp(-alfa*(i*dt)*(i*dt));
//printf("%f\n",b[i]);
}
for(i=0;i<200;i++)
x[i]=0;
//计算卷积
for(n=0;n<198;n++)
{for(m=0;m<31;m++)
if((n-m>=0)&&(n-m<198)) x[n]=x[n]+zb[m]*r[n-m];
}
for(i=0;i<31;i++) temp[i]=0.0;
for(i=0;i<31;i++) Ar[i]=0.0;
for(i=0;i<31;i++) Ai[i]=0.0;
for(i=0;i<31;i++) Br[i]=0.0;
for(i=0;i<31;i++) Bi[i]=0.0;
for(i=0;i<31;i++) za[i]=0.0;
dft(zb,temp,Br,Bi,31,1);
for(i=0;i<31;i++) Ar[i]=1.0/Br[i];
for(i=0;i<31;i++) Ai[i]=1.0/Bi[i];
//计算傅立叶反变换
dft(Ar,Ai,za,temp,31,-1);
//计算反褶积
for(i=0;i<200;i++)
ri[i]=0;
for(n=0;n<200;n++)
{for(m=0;m<31;m++)
if((n-m>=0)&&(n-m<200)) ri[n]=ri[n]+za[m]*x[n-m];
}
if((fp=fopen("hecheng.txt","w"))==NULL)
printf("Can't open file.\n");
else
{for(i=0;i<200;i++)
fprintf(fp,"%f\n",x[i]);
}
fclose(fp);
if((fp=fopen("zibo.txt","w"))==NULL)
printf("Can't open file.\n");
else
{for(i=0;i<31;i++)
fprintf(fp,"%f\n",za[i]);
}
fclose(fp);
//输出反卷积结果
if((fp=fopen("fanzheji.txt","w"))==NULL)
printf("Can't open file.\n");
else
{for(i=0;i<200;i++)
fprintf(fp,"%f\n",ri[i]);
}
fclose(fp);
}
/************************************************************************
* 离散傅立叶变换与反变换
* 输入: x--要变换的数据的实部
* y--要变换的数据的虚部
* a--变换结果的实部
* b--变换结果的虚部
* n--数据长度
* sign--sign=1时，计算离散傅立叶正变换;sign=-1时;计算离散傅立叶反变换
************************************************************************/
void dft(double x[],double y[],double a[],double b[],int n,int sign)
{
int i,k;
double c,d,q,w,s;
q=6.28318530718/n;
for(k=0;k<n;k++)
{
w=k*q;
a[k]=b[k]=0.0;
for(i=0;i<n;i++)
{
d=i*w;
c=cos(d);
s=sin(d)*sign;
a[k]+=c*x[i]+s*y[i];
b[k]+=c*y[i]-s*x[i];
}
}
if(sign==-1)
{
c=1.0/n;
for(k=0;k<n;k++)
{
a[k]=c*a[k];
b[k]=c*b[k];
}
}
}
