filcalc.m

德国人开发的地震处理分析软件

function filcalc

global A;
global numero;
global deltat;
global archivo;
global letras;
global linutil;
global vacio;
global MaxA;
global fildeltaf;
global Afil;
global nfilgraf;
global nnfilgraf;
global numerofil;
global Ampfil;
global fcifil;
global fcffil;
global Nestfil;
global estNfil;
global espefil;
global densfil;
global velofil;
global amorfil;
global filAmin;
global filAmax;
global filFmin;
global filFmax;
global SenXfil;
global SenYfil;
global tipofilF;
global filnombre;
global filruta;
global tipofilS;
global tipofilA;
global directorio;
global directorio1;

izq=numero(SenXfil);
der=numero(SenXfil+1)-1;
Nn=der-izq+1;

Nfft=2^(fix(log(der-izq+1)/log(2))+1);
df=1/deltat(SenXfil)/Nfft;
FFi=df;
FFf=Nfft/2*df;
Nf=Nfft/2;
i=sqrt(-1);

x = A(izq:der);
xfft = fft(x,Nfft);

if Nestfil==1;
   haskel=1;
else;
   cd (directorio1);

   fid=fopen('haskel.dat','wt');
   fprintf(fid,'%10.3f%10i%10i%10.3f\n',FFf,Nf,Nestfil,MaxA(SenXfil));
   for J=1:Nestfil;
       fprintf(fid,'%10.3f%10.3f%10.3f%10.3f\n',velofil(J),densfil(J),espefil(J),amorfil(J)); 
   end;
   if tipofilA==1;
      directorio2=[filruta,filnombre];
      fprintf(fid,'%80s\n',directorio2);
      fprintf(fid,'%10.3f\n',abs(xfft(1:Nf)));
      fclose(fid);
      !haskel2
   else;
      fclose(fid);
      !haskel
   end

   load haskel.trf;

   delete haskel.trf
   delete haskel.dat

   haskel=haskel(:,1).*exp(i*haskel(:,2));
   haskel2=conj(flipud(haskel));
   haskel=[1 ; haskel; haskel2(2:Nf)];
   if tipofilS==1;haskel=1./haskel;end;
end;

if fcifil==0;
   Ncerosi=0;
else;
   Ncerosi=fix(fcifil/df);
end;

if fcffil>=FFf;
   Ncerosf=0;
else;
   Ncerosf=fix((FFf-fcffil)/df);
end;

filtro=[zeros(Ncerosi,1);Ampfil*ones(Nf-Ncerosi-Ncerosf,1);zeros(Ncerosf,1)];
if tipofilF==1;filtro=abs(filtro-1);end;
filtro2=flipud(filtro);

filtro=haskel.*[0;filtro;filtro2(2:Nf)];

xfft=xfft.*filtro;

z=ifft(xfft);
z=real(z(1:Nn));
Amax = max(max(z),abs(min(z)));

lineas=0;
letras = [letras;length(SenYfil)];
nombre=[SenYfil vacio(length(SenYfil)+1:80)];
delta=deltat(SenXfil);

numero   =  [numero;numero(length(numero))+length(z)];
linutil  =  [linutil;lineas];
deltat   =  [deltat;delta];
archivo  =  [archivo;nombre];
MaxA     =  [MaxA;Amax];
A        =  [A;z];

fildeltaf=[fildeltaf;df];
Afil=[Afil;abs(filtro(2:Nf+1))];
nfilgraf=[nfilgraf length(deltat)];
numerofil=[numerofil;numerofil(length(numerofil))+Nf];
nnfilgraf=length(deltat);

cd (directorio);

filgra;