psdcal.m

地震资料处理的matlab源码软件包

function psdcal

global A;
global Apsd;
global Atf;
global Acor;
global Acorxy;
global psdgrafX;
global psdgrafY;
global psdN;
global npsdgraf;
global nnpsdgraf;
global ntfgraf;
global nntfgraf;
global ncorgraf;
global nncorgraf;
global ncorxygraf;
global nncorxygraf;
global psddeltaf;
global tfdeltaf;
global deltat;
global numero;
global numeropsd;
global numerotf;
global numerocor;
global numerocorxy;
global Limiteizq;
global Limiteder;

psdNN=psdN/2;
df=1/deltat(psdgrafX)/psdN;

Nox=numero(psdgrafX);
Nnx=numero(psdgrafX+1)-Nox;

if Limiteizq==Limiteder;
   Nizq=0;
   Nder=Nnx;
else;
   Nizq=fix(Limiteizq/deltat(psdgrafX));
   Nder=fix(Limiteder/deltat(psdgrafX));
end;

if Nizq>Nder;Nizq=Mud;Nizq=Nder;Nder=Mud;end;

if Nder>=Nnx;
   Nfx=Nox+Nnx-1;
else;
   Nfx=Nox+Nder-1;
   Nnx=Nder-Nizq;
end;

if Nnx > psdN;
  
   if psdgrafY==[];

      Bpsd=spectrum(A(Nox+Nizq:Nfx),psdN);
      ind=0;
      j=0;
      for i=npsdgraf;
          j=j+1;
          if i==psdgrafX;
             psddeltaf(j)=df;
             Apsd=[Apsd(1:numeropsd(j)-1,:);Bpsd;...
                   Apsd(numeropsd(j+1):length(Apsd),:)];
             dif=psdNN-numeropsd(j+1)+numeropsd(j);
             if j==length(npsdgraf);
                numeropsd=[numeropsd(1:j);numeropsd(j)+psdNN];
             else; 
                numeropsd=[numeropsd(1:j);numeropsd(j)+psdNN;...
                          (numeropsd(j+2:length(numeropsd))+dif)];
             end;
             ind=1;
             break; 
          end;
      end;

      if ind==0;
         Bcor=xcorr(A(Nox+Nizq:Nfx),'coeff');
         Bcor=Bcor(Nnx:2*Nnx-1);
         psddeltaf=[psddeltaf;df];
         Apsd=[Apsd;Bpsd];
         Acor=[Acor;Bcor];
         npsdgraf=[npsdgraf psdgrafX];
         numeropsd=[numeropsd;numeropsd(length(numeropsd))+psdNN];
         ncorgraf=[ncorgraf psdgrafX];
         numerocor=[numerocor;numerocor(length(numerocor))+Nnx];
      end;

      nnpsdgraf=psdgrafX;
      nncorgraf=nnpsdgraf;

   else;

      Noy=numero(psdgrafY);
      Nfy=Noy+Nder-1;  

      Bpsd=spectrum(A(Nox+Nizq:Nfx),...
                    A(Noy+Nizq:Nfy),psdN);
      ind=0;
      j=0;
      for i=npsdgraf;
          j=j+1;
          if i==psdgrafX;
             psddeltaf(j)=df;
             Apsd=[Apsd(1:numeropsd(j)-1,:);Bpsd(:,[1 6]);...
                   Apsd(numeropsd(j+1):length(Apsd),:)];
             dif=psdNN-numeropsd(j+1)+numeropsd(j);
             if j==length(npsdgraf);
                numeropsd=[numeropsd(1:j);numeropsd(j)+psdNN];
             else; 
                numeropsd=[numeropsd(1:j);numeropsd(j)+psdNN;...
                          (numeropsd(j+2:length(numeropsd))+dif)];
             end;
             ind=1;
             break; 
          end;
      end;
      if ind==0;
         Bcor=xcorr(A(Nox+Nizq:Nfx),'coeff');
         Bcor=Bcor(Nnx:2*Nnx-1);
         psddeltaf=[psddeltaf;df];
         Apsd=[Apsd;Bpsd(:,[1 6])];
         Acor=[Acor;Bcor];
         npsdgraf=[npsdgraf psdgrafX];
         numeropsd=[numeropsd;numeropsd(length(numeropsd))+psdNN];
         ncorgraf=[ncorgraf psdgrafX];
         numerocor=[numerocor;numerocor(length(numerocor))+Nnx];
      end;

      ind=0;
      j=0;
      for i=npsdgraf;
          j=j+1;
          if i==psdgrafY;
             psddeltaf(j)=df;
             Apsd=[Apsd(1:numeropsd(j)-1,:);Bpsd(:,[2 7]);...
                   Apsd(numeropsd(j+1):length(Apsd),:)];
             dif=psdNN-numeropsd(j+1)+numeropsd(j);
             if j==length(npsdgraf);
                numeropsd=[numeropsd(1:j);numeropsd(j)+psdNN];
             else; 
                numeropsd=[numeropsd(1:j);numeropsd(j)+psdNN;...
                          (numeropsd(j+2:length(numeropsd))+dif)];
             end;
             ind=1;
             break; 
          end;
      end;
      if ind==0;
         Bcor=xcorr(A(Noy+Nizq:Nfy),'coeff');
         Bcor=Bcor(Nnx:2*Nnx-1);
         psddeltaf=[psddeltaf;df];
         Apsd=[Apsd;Bpsd(:,[2 7])];
         Acor=[Acor;Bcor];
         npsdgraf=[npsdgraf psdgrafY];
         numeropsd=[numeropsd;numeropsd(length(numeropsd))+psdNN];
         ncorgraf=[ncorgraf psdgrafY];
         numerocor=[numerocor;numerocor(length(numerocor))+Nnx];
      end;

      ind=0;
      j=0;
      if ntfgraf~=[0;0];
         for i=ntfgraf(1,:);
             j=j+1;
             if i==psdgrafX & ntfgraf(2,j)==psdgrafY;
                tfdeltaf(j)=df;
                Atf=[Atf(1:numerotf(j)-1,:);Bpsd(:,[3 4 5 8]);...
                     Atf(numerotf(j+1):length(Atf),:)];
                dif=psdNN-numerotf(j+1)+numerotf(j);
                if j==length(ntfgraf(1,:));
                   numerotf=[numerotf(1:j);numerotf(j)+psdNN];
                else; 
                   numerotf=[numerotf(1:j);numerotf(j)+psdNN;...
                             (numerotf(j+2:length(numerotf))+dif)];
                end;
                ind=1;
                break; 
             end;
         end;
      else;ntfgraf=[];
      end;

      if ind==0;
         Bcor=xcorr(A(Nox+Nizq:Nfx),...
                    A(Noy+Nizq:Nfy),'coeff');
         Bcor=Bcor(Nnx:2*Nnx-1);
         tfdeltaf=[tfdeltaf;df];
         Atf=[Atf;Bpsd(:,[3 4 5 8])];
         Acorxy=[Acorxy;Bcor];
         ntfgraf=[ntfgraf [psdgrafX;psdgrafY]];
         numerotf=[numerotf;numerotf(length(numerotf))+psdNN];
         ncorxygraf=[ncorxygraf [psdgrafX;psdgrafY]];
         numerocorxy=[numerocorxy;numerocorxy(length(numerocorxy))+Nnx];
      end;

      nnpsdgraf=[psdgrafX psdgrafY];   
      nncorgraf=nnpsdgraf;
      nntfgraf=[psdgrafX;psdgrafY];
      nncorxygraf=nntfgraf;

   end;
end;

pxxgra;