?? convolution model.cpp
字號:
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define PI 3.1415926
#define f0 40
void main()
{
int i,j,k,num1,num2;
double dt,*R,wavelet[81],reflect[460],*s,*x;
FILE *fp,*reflection,*record;
/***************************Ricker子波離散化******************************/
dt=0.001; /*離散時間間隔*/
R=wavelet;
k=-40;
for(i=0;k<=40;k++)
{
wavelet[i]=(1-2*(PI*f0*k*dt)*(PI*f0*k*dt))*exp(-(PI*f0*k*dt)*(PI*f0*k*dt));
i++;
}
/****************************輸出數(shù)據(jù)到文件*******************************/
if((fp=fopen("D:\\code\\data\\X_wavelet.txt","w+"))==NULL)
{
printf("cann't open file\n");
return;
}
for(i=0;i<=80;i++)
if (fprintf(fp,"%f\n",wavelet[i])==0)
printf("file write error\n");
fclose(fp);
/****************************輸入反射系數(shù)序列******************************/
num1=81;
num2=231;
for(i=0;i<=230;i++)
reflect[i]=0;
reflect[40]=0.2353;
reflect[75]=0.1871;
reflect[110]=-0.3731;
reflect[140]=0.4247;
reflect[175]=-0.2455;
reflect[230]=0.1871;
s=(double *)malloc((num1+num2-1)*sizeof(double ));
x=(double *)malloc(num2*sizeof(double ));
/**************************輸出反射系數(shù)離散序列到文件**********************/
if((reflection=fopen("D:\\code\\data\\reflection.txt","w+"))==NULL)
{
printf("cann't open file\n");
return;
}
for(i=0;i<=230;i++)
if (fprintf(reflection,"%f\n",reflect[i])==0)
printf("file write error\n");
fclose(reflection);
/*********************wavelet[]與reflect[]褶積的過程***********************/
if(num1>=num2) /*先考慮wavelet數(shù)組長度大于reflect的情況*/
{
for(i=0;i<num1+num2-1;i++)
{
if(i<num2)
{
for(s[i]=0,j=0;j<=i;j++)
s[i]+=wavelet[j]*reflect[i-j];
}
if(num2<=i&&i<num1)
{
for(s[i]=0,j=i-num2+1;j<=i;j++)
s[i]+=wavelet[j]*reflect[i-j];
}
if(i>=num1)
{
for(s[i]=0,j=i-num2+1;j<num1;j++)
s[i]+=wavelet[j]*reflect[i-j];
}
}
}
else /*再考慮reflect數(shù)組長度大于wavelet的情況*/
{
for(i=0;i<num1+num2-1;i++)
{
if(i<num1)
{
for(s[i]=0,j=0;j<=i;j++)
s[i]+=wavelet[j]*reflect[i-j];
}
if(num1<=i&&i<num2-1)
{
for(s[i]=0,j=0;j<num1;j++)
s[i]+=wavelet[j]*reflect[i-j];
}
if(num2-1<=i)
{
for(s[i]=0,j=i-num2+1;j<num1;j++)
s[i]+=wavelet[j]*reflect[i-j];
}
}
}
for(j=0,i=40;i<(num1+num2-1-40);i++,j++) /*輸出褶積結(jié)果*/
x[j]=s[i];
for(j=0;j<num2;j++) /*輸出褶積結(jié)果*/
printf("x[%d]=%f\n ",j,x[j]);
/*******************************輸出結(jié)果離散序列***************************/
if((record=fopen("D:\\code\\data\\record.txt","w+"))==NULL)
{
printf("cann't open file\n");
return;
}
for(i=0;i<num2;i++)
if (fprintf(record,"%f\n",x[i])==0)
printf("file write error\n");
fclose(record);
}
?? 快捷鍵說明
復(fù)制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -