lvbo.cpp

地震勘探滤波程序

#include"stdio.h"
#include"math.h"
#include"stdlib.h"
#define PI 3.1415926
/*双精度型的一维数组，输入（输出）信号的实部和虚部*/
/*m0:2 的次方数，FFT的点数nfft＝2m0*/
/*inv= 1 forward transform, inv = -1 inverse transform*/
void fft(double sr[],double sx[],int m0,int inv)
{
	int i,j,lm,li,k,lmx,np,lix,mm2;
	double t1,t2,c,s,cv,st,ct;
	if(m0<0)
		return;
	lmx = 1;
	for(i=1;i<=m0;++i)
		lmx += lmx;   //form 2**m0
	cv = 2.0*PI/(double)lmx;
	ct = cos(cv);
	st = -inv*sin(cv);
	np = lmx;
	mm2 = m0-2;
	/*fft butterfly numeration*/
	for(i=1;i<=mm2;++i)
	{
		lix = lmx;
		lmx /= 2;
		c =  ct;
		s = st;
		for(li=0;li<np;li+=lix)
		{
			j = li;
			k = j+lmx;
			t1 = sr[j]-sr[k];
			t2 = sx[j]-sx[k];
			sr[j] += sr[k];
			sx[j] += sx[k];
			sr[k] = t1;
			sx[k] = t2;
			++j;
			++k;
			t1 = sr[j]-sr[k];
			t2 = sx[j]-sx[k];
			sr[j] += sr[k];
			sx[j] +=sx[k];
			sr[k] = c*t1-s*t2;
			sx[k] = s*t1+c*t2;
		}
		for(lm=2;lm<lmx;++lm)
        {
			cv = c;
			c = ct*c-st*s;
			s = st*cv+ct*s;
			for(li=0;li<np;li+=lix)
			{
				j = li+lm;
				k = lmx+j;
				t1 = sr[j]-sr[k];
				t2 = sx[j]-sx[k];
				sr[j] +=sr[k];
				sx[j] +=sx[k];
				sr[k] = c*t1-s*t2;
				sx[k] = s*t1+c*t2;
			}
		}
		cv = ct;
		ct = 2.0*ct*ct-1.0;
		st = 2.0*st*cv;
	}
	 /*4 points DFT*/
	if(m0>=2)
		for(li=0;li<np;li+=4)
		{
			j = li;
			k = j+2;
			t1 = sr[j]-sr[k];
			t2 = sx[j]-sx[k];
			sr[j] += sr[k];
			sx[j] += sx[k];
			sr[k] = t1;
			sx[k] = t2;
			++j;
			++k;
			t1 = sr[j]-sr[k];
			t2 = sx[j]-sx[k];
			sr[j] += sr[k];
			sx[j] += sx[k];
			sr[k] = inv*t2;
			sx[k] = -inv*t1;
		}
		/*2 points DFT*/
		for(li=0;li<np;li+=2)
		{
			j = li;
			k = j+1;
			t1 = sr[j]-sr[k];
			t2 = sx[j]-sx[k];
			sr[j] += sr[k];
			sx[j] += sx[k];
            sr[k] = t1;
			sx[k] = t2;
		}
		/*sort according to bit reversal*/
		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[j];
				sr[j] = sr[i];
				sr[i] = t1;
				t1 = sx[j];
				sx[j] = sx[i];
				sx[i] = t1;

			}
		}
		/* if Inverse FFT,multiply 1.0/np */
		if(inv!=-1)
			return;
		t1 = 1.0/np;
		for(i=0;i<np;++i)
		{
			sr[i] *= t1;
			sx[i] *=t1;
		}
}

void lvbo(double *xt,int np,float dt,float fc1,float fc2)
{
	double *xr,*xi;
	float *H;
	int i,nfft,k;
	float t,df,f,z;
	FILE *fpar,*fp1,*fp2,*fp3,*fp4,*fp5,*fp6;
	char fil1[80],fil2[80],fil3[80],fil4[80],fil5[80],fil6[80];

//	fpar = fopen("filter.par","r");
//	fscanf(fpar,"%s",fil1);
//	fscanf(fpar,"%s",fil2);
//    fscanf(fpar,"%s",fil3);
///	fscanf(fpar,"%s",fil4);
//	fscanf(fpar,"%s",fil5);
//	fscanf(fpar,"%s",fil6);
//	printf("%s\n",fil1);
//	printf("%s\n",fil2);
//	printf("%s\n",fil3);
//	printf("%s\n",fil4);
//	printf("%s\n",fil5);
//	printf("%s\n",fil6);
/*	fscanf(fpar,"%d",&np);
	printf("np=%d\n",np);
	fscanf(fpar,"%f",&dt);
	printf("dt=%8.3fms\n",dt);
	dt = dt/1000.0;         //**参数文件中dt以ms为单位*
	fscanf(fpar,"%f%f",&fc1,&fc2);
	printf("fc1=%5.1f fc2=%5.1f\n",fc1,fc2);*/
    //计算FFT点数
	k = log(np)/log(2);
    if(np>pow(2,k))
		k = k+1;
	nfft = pow(2,k);
    //计算频率采样间隔
	df = 1.0/(nfft*dt);
//	printf("nfft=%d k=%d\n",nfft,k);
//	printf("dt=%8.4f df=%8.4f\n",dt,df);
    //动态分配数组内存
	xr = (double*)calloc(nfft,sizeof(double));
	xi = (double*)calloc(nfft,sizeof(double));
	H = (float*)calloc(nfft,sizeof(float));
	//生成频率域带通滤波器（零相位）
//	fp1 = fopen(fil1,"w");
	for(i=0;i<=nfft/2;i++)
	{
		f = i*df;
        if(f<=fc2 && f>=fc1)
			H[i] = 1.0;
		else
			H[i] = 0.0;
//		fprintf(fp1,"%4d%10.4f%12.4\n",i,f,H[i]);
	}
    //滤波器是对称的
	for(i=nfft/2+1;i<nfft;i++)
	{
		f = i*df;
		H[i] = H[nfft-1];
//		fprintf(fp1,"%4d %10.4f %12.4\n",i,f,H[i]);
	}
//	fclose(fp1);
//	fp2 = fopen(fil2,"w");
    //傅氏反变换，生成时间域滤波器
	for(i=0;i<nfft;i++)
	{
		xr[i] = H[i];
		xi[i] = 0.0;
	}
	fft(xr,xi,k,-1);
	for(i=-np;i<=-1;i++)
	{
		t = i*dt;
//		fprintf(fp2,"%10.4f%12.4f\n",t,xr[-i]/dt);
	}
	for(i=0;i<=np;i++)
	{
		t = i*dt;
//		fprintf(fp2,"%10.4f%12.4f\n",t,xr[i]/dt);
	}
//	fclose(fp2);

	for(i=0;i<np;i++)
		xr[i] = xt[i];
/*	fp3 = fopen(fil3,"r");
    //输入原始信号
	for(i=0;i<np;i++)
	{
		fscanf(fp3,"%f%f",&t,&z);
		xr[i] = z;
	}
	fclose(fp3);*/
    //点数不够，补零
	if(np<nfft)
	{
        for(i=np;i<nfft;i++)
			xr[i] = 0;
	}
	for(i=0;i<nfft;i++)
	    xi[i] = 0.0;
    //对信号做傅氏正变换，得到振幅谱
	fft(xr,xi,k,1);
//	fp4 = fopen(fil4,"w");
	for(i=0;i<nfft;i++)
	{
        f = i*df;
//		fprintf(fp4,"%4d  %10.4f  %12.4\n",i,f,sqrt(xr[i]*xr[i]+xi[i]*xi[i])*dt);
	}
//	fclose(fp4);
    //频率域滤波Y(f) = X(f)H(f)
//	fp5 = fopen(fil5,"w");
	for(i=0;i<nfft;i++)
	{
		xr[i] = xr[i]*H[i];
        xi[i] = xi[i]*H[i];
		f = i*df;
//		fprintf(fp5,"%4d  %10.4f  %12.4f\n",i,f,sqrt(xr[i]*xr[i]+xi[i]*dt));
	}
//	fclose(fp5);
    //做傅氏反变换，得到时间域滤波结果
	fft(xr,xi,k,-1);
//	fp6 = fopen(fil6,"w");
	for(i=0;i<=np;i++)
	{
		t = i*dt;
		xt[i] = xr[i];
//		fprintf(fp6,"%10.4f%12.4f\n",t,xr[i]);
	}
//	fclose(fp6);
	//释放内存
	free(xr);
	free(xi);
	free(H);
}