polarize.m


"polarize.m" is a complicated times series program that reads in SAC (Seismic Analysis Code) binary files, does some filtering of the data, and then performs an eigenvalue/eigenvector decomposition of the particle motion to investigate estimates of azimuth and incidence angle from arrivals in the three component waveforms.  Plots include the input time series, filtered time series, eigenvalues as a function of time, azimuth and incidence angle calculations from the eigenvectors, and rose plots of the same.  Links are given to functions called in "polarize.m".


 
 
polarize.m
More comments
function polarize(station,delt,ttot,twin,hilb,flp,fhi)
%
%  function polarize(station,delt,ttot,twin,hilb,flp,fhi)
%
%  Program to read in 3 component waveform data
%  Create the covariance matrix for a moving time window
%  Find the principal components and infer polarization
%
% input:
%   station = station name for sacfile prefix
%   delt = sampling interval
%   ttot = total number of seconds to analyze in traces
%   twin = time window length, each time shift will be 1/2 of the
%      window length
%   hilb = 0, no hilbert transform of vertical component
%      = 1, hilbert transform
%   flp = low frequency corner frequency of a 2nd order butterworth
%     filter used to filter the data, if 0, then no filtering
%   fhi = hi frequency corner frequency of the filter

ename=strcat(station,'.EGE.NM.1.sac');  %construct the filename
nname=strcat(station,'.EGN.NM.1.sac');
zname=strcat(station,'.EGZ.NM.1.sac');
npts=ttot/delt + 1;
[t,e] = readsac(ename,npts);     % read the data
[t,n] = readsac(nname,npts);
[t,z] = readsac(zname,npts);
e=dmean(e);           % remove the mean from each
n=dmean(n);
z=dmean(z);
ntest=1;
if ntest == 1;
 else;
 % A synthetic test case
 npts=1000;
 ttot=delt*(npts-1);
 azrayl=30.;
 azlove=45;
 amprayl=2;
 amplove=1.;
 [t,e,n,z]=testchirps(npts,delt,amplove,amprayl,azlove,azrayl);
end;
if hilb ==1;     % hilbert transform the vertical component
 zh=hilbert(z);
 z=-imag(zh);
 else;
end;
f1=figure('name','DATA SEISMOGRAMS');  % plot the data
subplot(3,1,1);
plot(t,e);
xlabel('time sec');
ylabel(strcat('EW Comp at ',station));
subplot(3,1,2);
plot(t,n);
xlabel('time sec');
ylabel(strcat('NS Comp at ',station));
subplot(3,1,3);
plot(t,z);
xlabel('time sec');
ylabel(strcat('Z comp at ',station));
%
% filter the data
%
if flp > 0;
   e1=bandpass(e,flp,fhi,npts,delt);
   n1=bandpass(n,flp,fhi,npts,delt);
   z1=bandpass(z,flp,fhi,npts,delt);
   e=e1;
   n=n1;
   z=z1;
   %
 f2=figure('name','FILTERED SEISMOGRAMS');  % plot the filtered data
 subplot(3,1,1);
 plot(t,e);
 xlabel('time sec');
 ylabel(strcat('EW Comp at ',station));
 subplot(3,1,2);
 plot(t,n);
 xlabel('time sec');
 ylabel(strcat('NS Comp at ',station));
 subplot(3,1,3);
 plot(t,z);
 xlabel('time sec');
 ylabel(strcat('Z comp at ',station));
else;
end;

%
%   Moving window loop
%
npts1=fix(ttot/delt) + 1   %  total number of samples to analyze
nwin=fix(twin/delt) + 1    %  number of samples in a time window
npshift=fix(twin/(2*delt))+1  % number of samples to shift over
kfin=fix((npts1-nwin)/(npshift+1))+1 % number of time windows considered
mxde1=0.;
mxde2=0.;
mxde3=0.;
for k=1:kfin;
   nwinst=(k-1)*(npshift-1)+1;  % start of time window
   nwinfn=nwinst+nwin-1;    % end of time window
   a=csigm(e,n,z,nwinst,nwinfn);  % signal matrix
   %a
   c=a'*a;         % correlation matrix
   %c
   [v1,d1]=eig(c);      % eigenvalue/eigenvectors
   %d1
   %v1
   [v,d]=order(v1,d1);     % order the eigenvalues and eigenvectors
   %d
   %v
   ang1(k)=atan2(v(1,1),v(2,1)) * 180/pi; % azimuth for first  eigenvalue
   ang2(k)=atan2(v(1,2),v(2,2)) * 180/pi; % azimuth for second eigenvalue
   ang3(k)=atan2(v(1,3),v(2,3)) * 180/pi; % azimuth for third  eigenvalue
 de1(k)=d(1);
 de2(k)=d(2);
   de3(k)=d(3);
   mxde1=max(mxde1,de1(k));  % find the maximum values
   mxde2=max(mxde2,de2(k));
   mxde3=max(mxde3,de3(k));
 vang1(k)=acos(abs(v(3,1)))* 180/pi; %angle from the vertical
 vang2(k)=acos(abs(v(3,2)))* 180/pi;
 vang3(k)=acos(abs(v(3,3)))* 180/pi;
 t2(k)=delt*(nwinst-1);    % assign time for this window to the window start
end;
%
f3=figure('name','Eigenvalues and Inferred Azimuth');
subplot(3,1,1);
plot(t2,de1,'-or',t2,de2,'-dg',t2,de3,'-+b');
xlabel('time sec');
ylabel('eigenvalues');
subplot(3,1,2);
plot(t2,ang1,'-or',t2,ang2,'-dg',t2,ang3,'-+b');
xlabel('time sec');
ylabel('Azimuth ');
subplot(3,1,3);
plot(t2,vang1,'-or',t2,vang2,'-dg',t2,vang3,'-+b');
xlabel('time sec');
ylabel('incidence angle ');
%
%  Rose plot
f4=figure('name','Azimuth Distribution');
subplot(2,3,1);
title('Azimuth - Largest Eigenvalue');
rose(ang1*pi/180,100);
subplot(2,3,2);
title('Azimuth - Intermediate Eigenvalue');
rose(ang2*pi/180,100);
subplot(2,3,3);
title('Azimuth - Smallest Eigenvalue');
rose(ang3*pi/180,100);
%
%  Compass plots
%f5=figure('name','Compass Plots');
subplot(2,3,4);
title('Azimuth - Largest Eigenvalue');
compass(de1.*cos(ang1*pi/180),de1.*sin(ang1*pi/180));
subplot(2,3,5);
title('Azimuth - Intermediate Eigenvalue');
compass(de2.*cos(ang2*pi/180),de2.*sin(ang2*pi/180));
subplot(2,3,6);
title('Azimuth - Smallest Eigenvalue');
compass(de3.*cos(ang3*pi/180),de3.*sin(ang3*pi/180));
% Do the same thing for the highest amplitudes for each
% eigenvalue
%
nskip=1;
if nskip == 1;
   else;
neig1=0;
neig2=0;
neig3=0;
for k=1:kfin;
   if de1(k) >= 0.5*mxde1;
      neig1=neig1+1;
      angm1(neig1)=ang1(k);
   else;
   end;
   if de2(k) >= 0.5*mxde2;
      neig2=neig2+1;
      angm2(neig2)=ang2(k);
   else;
   end;
 if de3(k) >= 0.5*mxde3;
      neig3=neig3+1;
      angm3(neig3)=ang3(k);
   else;
   end;
end;
%
subplot(2,3,4);
title('Azimuth - Largest Eigenvalue,50% Threshold');
rose(angm1*pi/180,100);
subplot(2,3,5);
title('Azimuth - Intermediate Eigenvalue,50% Threshold');
rose(angm2*pi/180,100);
subplot(2,3,6);
title('Azimuth - Smallest Eigenvalue,50% Threshold');
rose(angm3*pi/180,100);
end;

 


 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

Only the first part of the  filename is given.  This part reconstructs the rest.
 
 

readsac is a routine written by Roy Greenfield and Tom Battenhouse.
 

Remove the mean.
 
 
 
 

This was used when checking the program results.
 
 
 
 
 
 
 
 
 
 
 
 
 

Plot the raw data.  We are going to plot three seismograms in the window.
 
 
 
 
 
 
 
 
 
 
 

Filter the data between flp and fhi.
 
 
 
 
 
 
 

Plot the filtered seismograms.
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

Apply a moving window analysis to the 3 component data.
 
 
 
 
 
 
 
 
 
 
 
 

Load a time portion of the signals into the signal matrix.

Construct the correlation matrix.

Find the eigenvalues and eigenvectors of the signal matrix.
 

Order the resulting eigenvectors/eigenvalues in decending order.
 

Calculate azimuths from the eigenvector components.
 
 
 
 
 
 
 
 

Calculate the apparent incidence angles from the eigenvectors.
 
 
 
 
 

Plot this stuff.
 
 
 
 
 
 
 
 
 
 
 
 

Now make a Rose diagram of azimuths.
The top set is directly from the eigenvectors.
 
 
 
 
 
 
 
 
 

This set of "compass" plots is weighted by the eigenvalues.
 
 
 
 
 
 
 
 
 
 
 
 
 

An experimental portion not invoked by the example shown here.
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

 



Results of this module using the following input parameters:

>>polarize('LOT1',0.005,30.0,0.5,1,1.0,5.0)

(Note:  The data come from the Gujarat aftershock deployment!)





Matlab WebSite: http://www.mathworks.com/



 

Back