corrot.m,v

具有特色的地震数据处理源码

head	3.0;
access;
symbols;
locks; strict;
comment	@// @;


3.0
date	2000.06.13.19.19.57;	author gilles;	state Exp;
branches;
next	2.0;

2.0
date	99.05.21.18.45.20;	author mah;	state Exp;
branches;
next	1.5;

1.5
date	99.02.23.15.33.23;	author kay;	state Exp;
branches;
next	1.4;

1.4
date	99.02.22.20.39.44;	author kay;	state Exp;
branches;
next	1.3;

1.3
date	99.02.19.16.17.50;	author kay;	state Exp;
branches;
next	1.2;

1.2
date	99.02.19.15.20.26;	author kay;	state Exp;
branches;
next	1.1;

1.1
date	99.02.19.14.56.17;	author kay;	state Exp;
branches;
next	;


desc
@correlation matrix method
@


3.0
log
@Release 3
@
text
@function [dataout]=corrot(datain,headw1,tint,comp1,comp2)

%corrot -> function to rotate components of 3-C borehole data (DSI etc.)
%into radial and transverse components using matrix eigenvalue algorithm.
%
%function [dataout]=corrot(datain,headw1,tint,comp1,comp2)
%
%
%INPUT VARIABLES
%'datain' must be in official DSI data format
%Each record must represent a component h1, h2, or z.
%This can be achieved using 'sortrec'.
%
%headw1 = 	header word containing first break picks
%tint =     half width of time window to use around first breaks (s)
%comp1 =    record representing one of the components to be rotated
%comp2 =    record number of other component to be rotated
%
%OUTPUT VARIABLES
%Trace header word 4 contains component information.
%This word will be incremented by 3 for components that have
%been rotated.
%
%By convention: h1=>1; h2=>2; z=>3; radial or oriented horizontal=>4;
%transverse horizontal=>5; direct P-arrival (Pd)=>6; orthoganal
%component to Pd=>7 (see Handbook of Geophysical Exploration, sect.1,vol.14B).  
%Components 6 and 7 are the results of rotating 3 and 4; 4 and 5 are
%the results of rotating 1 and 2.
%The largest singular value of the covarience matrix is stored in header
%word 10, and the ratio of sigular value 2 to singular value 1 is stored 
%in word 11.
%
%DSI customized VSP processing software
%by I. Kay and G. Perron (Jan 1998)

%$Id: corrot.m,v 2.0 1999/05/21 18:45:20 mah Exp gilles $
%$Log: corrot.m,v $
%Revision 2.0  1999/05/21 18:45:20  mah
%Release 2
%
%Revision 1.5  1999/02/23 15:33:23  kay
%Using matlab's cov is a bad idea here...
%
%Revision 1.2  1999/01/28 14:07:30  kay
%Added svd(1) to header word 10 and ratio of svd(2)/svd(1) to word 11.
%
%Revision 1.1  1999/01/06 19:09:07  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]=corrot(datain,headw1,tint,comp1,comp2)');

%check to make sure data is separated into components
for i=3:-1:1 %get number of traces in each component
 ntr(i)=datain.th{i}(12,1);
end %for

if (ntr(1)~=ntr(2)) | (ntr(1)~=ntr(3))
 error('check data format - different number of traces in components');
end%if

if length(datain.dat)~=3
 error('data must have only 3 records - one for each of x, y and z');
end %if

%***********************************************************************
%w=180/pi;
tstart=datain.fh{9}; %start time in seconds
int=datain.fh{8}; %sampling interval in seconds
nsamp=datain.fh{7}; %number of points per trace
dataout=datain;
a=comp1;
b=comp2;

ntr=ntr(1);
%increment component trace header word
dataout.th{a}(4,:)=datain.th{a}(4,:)+3;
dataout.th{b}(4,:)=datain.th{b}(4,:)+3;
ca=dataout.th{a}(4,1);
cb=dataout.th{b}(4,1);
rotang=zeros(1,ntr); %initialize variable for angles

%only want to rotate picked traces
traces=find(datain.th{a}(headw1,:)~=0);

for j=traces
	samp1=round((datain.th{a}(headw1,j)-tstart)/int-(tint/int)) +1;
	if samp1<1
	samp1=1;
	end %if
	samp2=samp1+2.*round(tint/int) +1;
% I've decided that cov() is bad.  It removes the mean, so if the
%window of data does not have a zero mean, the relative
%amplitudes/energies between traces will change, biasing the result. 

   %corrmat=cov(datain.dat{a}(samp1:samp2,j),datain.dat{b}(samp1:samp2,j));

%instead, use the raw data...

   v1=datain.dat{a}(samp1:samp2,j);
   v2=datain.dat{b}(samp1:samp2,j);
   corrmat=[ dot(v1,v1) dot(v1,v2); dot(v2,v1) dot(v2,v2)];

%For a linearly polarized waveform with source signature w(t), arriving at
% an angle \phi from H1 (x), the observed data is 
%		u=w(t)[\cos(\phi) \hat{x} + \sin(\phi) \hat{y}]
% The covarience matrix [ <x,x>, <x,y>, <y,x>, <y,y>] is then of the form
% [ (\cos(\phi))^2, \cos(\phi)*\sin(\phi); 
%     \cos(\phi)*\sin(\phi), (\sin(\phi))^2]
% \phi can be estimated from the elements of the covarience matrix as:

   rotang(j)=atan2(corrmat(1,2),corrmat(1,1));

   rotmat=[cos(rotang(j)) -sin(rotang(j)); sin(rotang(j)) cos(rotang(j))];
   rotdata=(rotmat'*[datain.dat{a}(:,j) datain.dat{b}(:,j)]')';

   dataout.dat{a}(:,j)=rotdata(:,1); %rotated h1 component
   dataout.dat{b}(:,j)=rotdata(:,2); %rotated h2 component


end %for j=trace
firsttr=traces(1);
pilotsamp1=round((datain.th{a}(headw1,firsttr)-tstart)/int-(tint/int)) +1;
if pilotsamp1<1
	pilotsamp1=1;
end %if
pilotsamp2=pilotsamp1+2.*round(tint/int)+1;
pilot=dataout.dat{a}(pilotsamp1:pilotsamp2,1);

for j=traces
  samp1=round((datain.th{a}(headw1,j)-tstart)/int-(tint/int)) +1;
  if samp1<1
     samp1=1;
  end %if
  samp2=samp1+(pilotsamp2-pilotsamp1);
  corrmat=cov(pilot, dataout.dat{a}(samp1:samp2,j));
  %v2=dataout.dat{a}(samp1:samp2,j);
  %corrmat=[ dot(pilot,pilot) dot(pilot, v2); dot(v2, pilot) dot(v2,v2)];
  if  corrmat(1,2) < 0.
     rotang(j)=rotang(j)+pi;
     dataout.dat{a}(:,j)=dataout.dat{a}(:,j).*-1. ;
     dataout.dat{b}(:,j)=dataout.dat{b}(:,j).*-1. ;
  end %if
end %for j=traces
dataout.th{a}(5,:)=mod(rotang*180/pi+360,360);
dataout.th{b}(5,:)=mod(rotang*180/pi+360,360);
@


2.0
log
@Release 2
@
text
@d36 1
a36 1
%$Id: corrot.m,v 1.5 1999/02/23 15:33:23 kay Exp mah $
d38 3
@


1.5
log
@Using matlab's cov is a bad idea here...
@
text
@d36 1
a36 1
%$Id: corrot.m,v 1.2 1999/01/28 14:07:30 kay Exp $
d38 3
@


1.4
log
@First version in which I have some real confidence...
@
text
@d15 3
a17 4
%tint =         time window to use around first breaks (s)
%(window starts 'tint' sec. before pick and ends 'tint' sec. after pick
%comp1 =        record representing one of the components to be rotated
%comp2 =        record number of other component to be rotated
d108 25
a132 10
   samp1=round((datain.th{a}(headw1,j)-tstart)/int-(tint/int)) +1;
   if samp1<1
      samp1=1;
   end %if
   samp2=samp1+2.*round(tint/int) +1;
   %v1=datain.dat{a}(samp1:samp2,j);
   %v2=datain.dat{b}(samp1:samp2,j);
   %corrmat=[ dot(v1,v1) dot(v1,v2); dot(v2,v1) dot(v2,v2)];
   corrmat=cov(datain.dat{a}(samp1:samp2,j),datain.dat{b}(samp1:samp2,j));
 
d134 1
@


1.3
log
@added a cross correlation of the first trace with all other traces.
If the cross correlation is negative, then multiply the trace by -1 and
add 180 to the rotation angle.
@
text
@d89 1
a89 1
w=180/pi;
d114 8
d123 2
a124 1
  corrmat=cov(datain.dat{a}(samp1:samp2,j),datain.dat{b}(samp1:samp2,j));
a125 3
  rotang(j)=-atan2(corrmat(1,2),corrmat(1,1));
  rotmat=[cos(rotang(j)) sin(rotang(j)); -sin(rotang(j)) cos(rotang(j))];
  rotdata=(rotmat'*[datain.dat{a}(:,j) datain.dat{b}(:,j)]')';
d127 3
a129 7
  dataout.dat{a}(:,j)=rotdata(:,1); %rotated h1 component
  dataout.dat{b}(:,j)=rotdata(:,2); %rotated h2 component


end %for j=traces

pilotsamp1=round((datain.th{a}(headw1,j)-tstart)/int-(tint/int)) +1;
d137 8
a144 8
   samp1=round((datain.th{a}(headw1,j)-tstart)/int-(tint/int)) +1;
   if samp1<1
      samp1=1;
   end %if
   samp2=samp1+(pilotsamp2-pilotsamp1);
  % size(pilot)
   %size(datain.dat{a}(samp1:samp2,j))
  corrmat=cov(pilot, datain.dat{a}(samp1:samp2,j));
d146 3
a148 3
  	rotang(j)=-rotang(j);
	dataout.dat{a}(:,j)=dataout.dat{a}(:,j).*-1. ;
	dataout.dat{b}(:,j)=dataout.dat{b}(:,j).*-1. ;
d151 2
a152 3

dataout.th{a}(5,:)=rotang;
dataout.th{b}(5,:)=rotang;
@


1.2
log
@*** empty log message ***
@
text
@d4 1
a4 1
%into radial and transverse components using crosscorrelation matrix.
d37 1
a37 1
%$Id: corrot.m,v 1.1 1999/02/19 14:56:17 kay Exp $
d39 4
a42 1
%Revision 1.1  1999/02/19 14:56:17  kay
a44 1

d117 2
a118 2
  rotang(j)=atan2(corrmat(1,2),corrmat(1,1));
  rotmat=[cos(rotang(j)) -sin(rotang(j)); sin(rotang(j)) cos(rotang(j))];
d124 24
a147 1
 
@


1.1
log
@Initial revision
@
text
@d37 5
a41 2
%$Id: $
%$Log:$
d115 2
a116 2
  rotang(j)=-atan2(corrmat(1,2),corrmat(1,1));
  rotmat=[cos(rotang(j)) sin(rotang(j)); -sin(rotang(j)) cos(rotang(j))];
@