1、把3個聲明函數的地方的註釋去掉了
2、把函數定義時的參數名字由m改成m1(因為m已經被定義成常量)
下面是修改後的源代碼:
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define N 1199
#define m 81
#define NN 1500
void main()
{
FILE *fp1,*fp2,*fp3;
double s[NN],b[NN],rxx[2397],y[NN],y1[NN],a;
int i,k;
void conv(double w[],int m1,double r[],int n,double s[],int l);
void cor(double *x,int n,double *h,int m1,double *y);
int lvs(double t[], int n,double b[],double x[]);
if((fp1=fopen("E:\\Seis.txt","r"))==NULL) //讀地震合成記錄
printf("Can't open file.\n");
else
{
for(i=0;i<N;i++)
{ fscanf(fp1,"%f",&a);
s[i]=a;}
}
for(i=0;i<m;i++) //右端項
{
if(i==0) b[i]=1.0;
else b[i]=0.0;
}
//調用相關函數,計算地震記錄的自相關rxx[l]
cor(s,N,s,N,rxx);
//調用萊文森算法,計算反濾波因子,並輸出
k=lvs(rxx,m,b,y);
printf("%d\n",k);
fp2=fopen("E:\\fanzibo.txt","w");
if(fp2==NULL) {printf("can't open file\n"); exit(0);}
for(i=0;i<m;i++)
{ fprintf(fp2,"%f\n",y[i]);}
fclose(fp2);
//調用褶積函數,計算反射系數
conv(y,m,s,N,y1,m+N-1);
fp3=fopen("E:\\fanzheji.txt","w");
if(fp3==NULL) {printf("can't open file\n"); exit(0);}
for(i=0;i<m+N-1;i++)
{ fprintf(fp3,"%f\n",y1[i]); }
fclose(fp3);
}
void conv(double w[],int m1,double r[],int n,double s[],int l) //褶積
{ int k,i;
for(k=0;k<l;k++)
{ s[k]=0.0; //初始化
for(i=0;i<m;i++)
{if(k-i>=0&&k-i<=n-1)
s[k]=s[k]+w[i]*r[k-i];}
}
}
//相關計算(計算x[n]與h[m]的互相關,結果為y[n])
void cor(double *x,int n,double *h,int m1,double *y)
{
int i,j;
for(i=0;i<n;i++)
{
y[i]=0.0;
for(j=0;j<m;j++)
{
if(i<=j)
y[i]+=x[j]*h[j-i];
}
}
}
int lvs(double t[], int n,double b[],double x[])
{
int i,j,k;
double a,beta,q,c,h,*y,*s;
s=(double *)malloc(n*sizeof(double));
y=(double *)malloc(n*sizeof(double));
a=t[0];
if(fabs(a)+1.0==1.0)
{
free(s);free(y);
printf("fail\n");
return(-1);
}
y[0]=1.0;x[0]=b[0]/a;
for(k=1;k<=n-1;k++)
{
beta=0.0;q=0.0;
for(j=0;j<=k-1;j++)
{
beta=beta+y[j]*t[j+1];
q=q+x[j]*t[k-j];
}
if(fabs(a)+1.0==1.0)
{
free(s);free(y);
printf("fail\n");
return(-1);
}
c=-beta/a;s[0]=c*y[k-1];y[k]=y[k-1];
if(k!=1)
for(i=1;i<=k-1;i++) s[i]=y[i-1]+c*y[k-i-1];
a=a+c*beta;
if(fabs(a)+1.0==1.0)
{
free(s);free(y);
printf("fail\n");
return(-1);
}
h=(b[k]-q)/a;
for(i=0;i<=k-1;i++)
{
x[i]=x[i]+h*s[i];
y[i]=s[i];
}
x[k]=h*y[k];
}
free(s);free(y);
return(1);
}