用前篇日志所写的互相关cross-correlation程序进行编写.

//---------------------------------------------------

InOutData: 一维数组,可读写, 保存有二维剖面数据.按道存取,NX道,NT采样点.

NX: 数据总道数.

NT: 每道采样点数.

StartP: 计算互相关的起始点,[0,NT-1]

EndP: 计算互相关的终止点,(StartP,NT-1]

maxDelayP: 最大延时或起前点数.

reftrace: 参考道号,[0,NT-1]

fillzero: 校正时,不存在的时间是否充零, true:充零, false:充该道的平均值.

//------------------------------------------------------------------------------
void NmoCorrectionfun(double **InOutData,int NX,int NT,int StartP,int EndP,int maxDelayP,int reftrace,bool fillzero)
{

/*
coded by yiling
email: yiling@email.jlu.edu.cn
*/

int n;
int maxposition;
double *x;
double *y;
double themean;
double themax;
int pro;
char err[256];
double *rxy;
int * Nmo;
double *tempdata;
double *data;
bool rr;

n=EndP-StartP+1;
data=*InOutData;
rxy=new double[maxDelayP*2+1];
x=new double[n];
y=new double[n];
Nmo=new int[NX];
tempdata=new double[NX*NT];

for(int jj=0;jj<n;jj++)
        x[jj]=data[reftrace*NT+StartP+jj];

for(int ii=0;ii<NX;ii++)
        {
        for(int jj=0;jj<n;jj++)
                y[jj]=data[ii*NT+StartP+jj];

        rr=CrossCorrelationMatlab(x,y,rxy,n,maxDelayP, &pro,err);

        themax=0;
        for(int jj=0;jj<(maxDelayP*2+1);jj++)
                if(themax<rxy[jj])
                        {
                        themax=rxy[jj];
                        Nmo[ii]=jj;
                        }
        }

for(int ii=0;ii<NX;ii++)
        {
        themean=0;
        for(int jj=0;jj<NT;jj++)
                themean +=data[ii*NT+jj];
        themean /=NT;

        if(fillzero)
                {
                 for(int jj=0;jj<NT;jj++)
                        tempdata[ii*NT+jj]=0;
                }
        else
                {
                 for(int jj=0;jj<NT;jj++)
                        tempdata[ii*NT+jj]=themean;
                }

        Nmo[ii]=Nmo[ii]-maxDelayP;       //=Nmo[ii]+StartP-n+1;

        if(Nmo[ii]<0)       //y延时
                {
                 for(int jj=-Nmo[ii];jj<NT;jj++)
                        tempdata[ii*NT+jj+Nmo[ii]]=data[ii*NT+jj];
                }
        else
                {
                 for(int jj=Nmo[ii];jj<NT;jj++)
                        tempdata[ii*NT+jj]=data[ii*NT+jj-Nmo[ii]];
                }

        }

(*InOutData)=tempdata;
delete[] rxy;
delete[] x;
delete[] y;
delete[] Nmo;
delete[] data;

}