gusucode.com > ​DSISoft是由加拿大地质调查局发布的用于垂直地震剖面(VSP)数据处理的免费软件包 > dsisoftv3/main/autocor.m

    function [dataout]=autocor(datain)

%[dataout]=autocor(datain)
%
%Autocorrelation
%Performs autocorrelation on each trace.
%Uses signal processing toolbox but contains commented out code that works
%without the toolbox.
%
%Similar to INSIGHT module 'auto'.
%
%DSI customized VSP processing software
%written by Kristen Beaty December, 1997

%$Id: autocor.m,v 3.0 2000/06/13 19:19:45 gilles Exp $
%$Log: autocor.m,v $
%Revision 3.0  2000/06/13 19:19:45  gilles
%Release 3
%
%Revision 2.0  1999/05/21 18:45:09  mah
%Release 2
%
%Revision 1.1  1999/01/06 19:09:01  kay
%Initial revision
%
%
%Copyright (C) 1998 Seismology and Electromagnetic Section/
%Continental Geosciences Division/Geological Survey of Canada
%
%This library is free software; you can redistribute it and/or
%modify it under the terms of the GNU Library General Public
%License as published by the Free Software Foundation; either
%version 2 of the License, or (at your option) any later version.
%
%This library is distributed in the hope that it will be useful,
%but WITHOUT ANY WARRANTY; without even the implied warranty of
%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
%Library General Public License for more details.
%
%You should have received a copy of the GNU Library General Public
%License along with this library; if not, write to the
%Free Software Foundation, Inc., 59 Temple Place - Suite 330,
%Boston, MA  02111-1307, USA.
%
%DSI Consortium
%Continental Geosciences Division
%Geological Survey of Canada
%615 Booth St.
%Ottawa, Ontario
%K1A 0E9
%
%email: dsi@cg.nrcan.gc.ca

disp('[dataout]=autocor(datain)');

dataout=datain;

nrec=datain.fh{12}; %get number of records
samples=datain.fh{7}; %number of points per trace
fpoints=2^nextpow2(3*samples); %number of points to be used in fft

for COUNT=1:nrec %loop over records
 ntpr=datain.th{COUNT}(12,1); %number of traces this record

 baseline=mean(datain.dat{COUNT},1);
 basegr=meshgrid(baseline,1:samples);
 dataout.dat{COUNT}=datain.dat{COUNT}-basegr;
 %each trace must have mean value of 0

  %----find the auto correlation functions in the Fourier domain---------

  %take fft of flipped data before zero padding to use for autocorrelation
  %au_fft=fft(flipud(dataout.dat{COUNT}),fpoints);
  %pilot=zeros(samples*3,ntpr); %initialize pilot
  %pilot(samples+1:2*samples,:)=dataout.dat{COUNT};
  %pilot_fft=fft(pilot,fpoints);

  %multiply magnitudes and add phases to get autocorrelation function fft
  %pfabs=abs(pilot_fft);
  %pfang=angle(pilot_fft);

  %perform autocorrelation
  %mag=abs(au_fft).*pfabs;
  %phase=angle(au_fft)+pfang;
  %aucorr_fft=mag.*exp(phase.*i); %convert to complex
  %aucorr=ifft(aucorr_fft,fpoints); %inverse fft
  %normgr=meshgrid(max(real(aucorr)),1:fpoints);
  %aucorr=real(aucorr)./normgr;
  %dataout.dat{COUNT}=flipud(aucorr(samples+1:2*samples,:));
%-----------------use signal processing toolbox xcor below------------

 aucorr=zeros(samples*2+1:datain.th{COUNT}(12,1));
 for tr=1:datain.th{1}(12,1)
  aucorr(:,tr)=xcorr(dataout.dat{COUNT}(:,tr));
 end %loop over traces
 normgr=meshgrid(max(aucorr),1:size(aucorr,1));
 dataout.dat{COUNT}=aucorr./normgr;
 dataout.dat{COUNT}=dataout.dat{COUNT}(samples:size(dataout.dat{COUNT},1),:);

end %loop over records

%end of function autocor