C PROGRAM SYNT - SYNTHETIC SEISMOGRAM COMPUTATION C C ************************************************************** C character*80 TITLE,mprint,iprint DIMENSION NDIS(100) COMPLEX S(2048),REF(2048),SS(2048) COMMON/SPECTR/S,REF,DF,NI,NT,NT2,NTM,ILS,IRS,NFS,NFILT,NWIN COMMON/SIG/F(2048),TO,DT,NPTS,NO,NSIG,NPL1,NPL2,NPL3, 1NPL4,IPR1,nder,nint COMMON/WIND/FLO,FLEFT,FRO,FRIGHT,FEXP,ILF,IRF,NFW COMMON/AUXIL/FF(2048),GG(2048),IS(2048),AUX(4),IAUX(4) COMMON/PLOT1/ARESP(6),ASPECT(6),ASYNT(6),JRESP(4),JSPECT(4),JSYNT( 14),NOO,RMIN,RMAX,GLR,TRMIN,TRMAX,GLTR,IR,ITR,NRED,NSH,AREDUC,TITLE 2,NPLOT C c *** mode=0 call serv(mode,lin,lou,llu2,llu3,llu4) if(mode.eq.0)lin=5 if(mode.eq.0)lou=6 if(mode.eq.0)call plots(ldum1,ldum2,7) c *** c READ(lin,110)TITLE READ(lin,101) IPR1,IPR2,LU8,LU9 WRITE(lou,111)TITLE WRITE(lou,101)IPR1,IPR2,LU8,LU9 if(mode.eq.1)lu8=llu2 if(mode.eq.1)lu9=llu3 IF(LU8.EQ.0)LU8=10 IF(LU9.EQ.0)LU9=11 REWIND LU8 REWIND LU9 ANULL=0. NULL=0 tmin=0. tmax=10000. C C READ INPUT PARAMETERS C READ(lin,101)NSIG,NPTS,NT,NWIN,NFILT,NPLOT,NSTOP,NSH, 1nder,nint WRITE(lou,101)NSIG,NPTS,NT,NWIN,NFILT,NPLOT,NSTOP,NSH, 1nder,nint IF(NSH.EQ.0)NSH=7 SHIFT=FLOAT(NSH) IF(NPLOT.EQ.0)GO TO 2 READ(lin,101)NNPLOT,(NDIS(I),I=1,NNPLOT) WRITE(lou,101)NNPLOT,(NDIS(I),I=1,NNPLOT) 2 CONTINUE IF(NT.EQ.0)NT=1024 NT2=NT/2 NTM=NT2+1 READ(LU8,110)MPRINT READ(LU8,110)IPRINT READ(LU8,152)XSOUR,ysour,ZSOUR,TSOUR,RSTEP,FL,DF READ(LU8,102)NDST,NFS,ILS,MCOMP,iloc DT=1./(DF*FLOAT(NT)) WRITE(lou,111)MPRINT WRITE(lou,111)IPRINT WRITE(lou,152)XSOUR,ysour,ZSOUR,TSOUR,RSTEP,FL,DF WRITE(lou,101)NDST,NFS,ILS,MCOMP WRITE(LU9,110)MPRINT WRITE(LU9,110)IPRINT WRITE(LU9,110)TITLE WRITE(LU9,152)XSOUR,ysour,ZSOUR,TSOUR,RSTEP,DT,DF WRITE(LU9,101)NDST,NT,MCOMP,iloc IRS=ILS+NFS IRSS=IRS+1 FR=FLOAT(IRS-1)*DF FMIN=0. NI=10 IF(NT.EQ.2048)NI=11 IF(NT.EQ.512)NI=9 IF(NT.EQ.256)NI=8 IF(NT.EQ.128)NI=7 READ(lin,106)ARESP,JRESP READ(lin,106)ASPECT,JSPECT READ(lin,106)ASYNT,JSYNT WRITE(lou,106)ARESP,JRESP WRITE(lou,106)ASPECT,JSPECT WRITE(lou,106)ASYNT,JSYNT IF(NPLOT.EQ.0)GO TO 7 CALL PLOT(5.,5.,-3) C C INPUT SIGNAL C 7 CALL SIGNAL(lin,lou) C C SPECTRUM OF THE INPUT SIGNAL C CALL FT(lin,lou) IF(NSTOP.EQ.1.and.nplot.ne.0)call plot(0.,0.,999) IF(NSTOP.EQ.1)STOP C C LOOP FOR RECEIVERS C DO 40 IDIST=1,NDST NPR1=0 NPR2=0 NPR3=0 READ(LU8,160)DST,AA WRITE(lou,160)DST,AA NFF=2*NFS DO 49 I=1,NTM FF(I)=0. SS(I)=(0.,0.) 49 REF(I)=(0.,0.) READ(LU8,161)(IS(I),I=1,NFF) DO 41 I=1,NFS A=AA*FLOAT(IS(2*(I-1)+1))/99999.1 B=AA*FLOAT(IS(2*I))/99999.1 FF(I+ILS+1)=SQRT(A*A+B*B) 41 REF(I+ILS+1)=CMPLX(A,B) IF(NPLOT.EQ.0)GO TO 42 DO 43 I=1,NNPLOT IF(IDIST.EQ.NDIS(I))GO TO 45 43 CONTINUE GO TO 42 45 CONTINUE READ(lin,101)NPR1,NPR2,NPR3 WRITE(lou,170)NPR1,NPR2,NPR3 CALL PLOT(5.,0.,-3) CALL SYMBOL(-3.,0.,0.3,'X(KM)=',90.,6) CALL NUMBER(-3.,2.,0.3,DST,90.,2) 42 CONTINUE C C PRINTING AND PLOTTING OF FREQUENCY RESPONSE C IF(IPR2.LE.1.AND.NPR1.EQ.0)GO TO 50 CALL REDUC(AREDUC,1,IRSS,FF,GG,IS) IF(IPR2.LE.1)GO TO 48 WRITE(lou,114)DST,DF,AREDUC,IRSS WRITE(lou,107)(IS(I),I=1,IRSS) 48 IF(NPR1.EQ.0)GO TO 50 NOO=1 nred=npr1 IF(NPR1.EQ.1)CALL ZPLOT(0.,DF,ILS,IRSS,GG) IF(NPR1.EQ.2)CALL ZPLOT(0.,DF,ILS,IRSS,FF) 50 CONTINUE C C AMPLITUDE SPECTRUM OF SYNTHETIC SEISMOGRAM C DO 18 I=1,NTM 18 FF(I)=0. DO 11 I=ILF,IRF SS(I)=REF(I)*S(I) A=REAL(SS(I)) B=AIMAG(SS(I)) 11 FF(I)=SQRT(A*A+B*B) DO 12 I=2,NT2 IJ=NT-I+2 12 SS(IJ)=CONJG(SS(I)) IF(IPR2.LE.2.AND.NPR2.EQ.0)GO TO 60 CALL REDUC(AREDUC,1,IRSS,FF,GG,IS) IF(IPR2.LE.2)GO TO 52 WRITE(lou,115)DST,DF,AREDUC,IRSS WRITE(lou,107)(IS(I),I=1,IRSS) 52 IF(NPR2.EQ.0)GO TO 60 NOO=5 NRED=NPR2 IF(NPR2.EQ.1)CALL ZPLOT(0.,DF,ILS,IRSS,GG) IF(NPR2.EQ.2)CALL ZPLOT(0.,DF,ILS,IRSS,FF) 60 CONTINUE C C GOING BACK TO TIME DOMAIN. SYNTHETIC SEISMOGRAMS C CALL FCOOLR(NI,SS,-1.0) DEL=1.0/FLOAT(NT) DO 20 I=1,NT 20 FF(I)=DEL*REAL(SS(I)) CALL REDUC(AREDUC,1,NT,FF,GG,IS) WRITE(LU9,162)DST,TO,AREDUC,NT WRITE(LU9,107)(IS(I),I=1,NT) C IF(IPR2.EQ.0.AND.NPR3.EQ.0)GO TO 70 IF(IPR2.EQ.0)GO TO 65 WRITE(lou,116)DST,TO,DT,AREDUC,NT WRITE(lou,108)(IS(I),I=1,NT) 65 IF(NPR3.EQ.0)GO TO 70 NOO=6 NRED=NPR3 IF(AREDUC.EQ.0.)GO TO 66 TSTART=TMINO-to IPOM=IFIX(TSTART/DT)+1 67 IF(IPOM.GT.0)GO TO 61 IPOM=IPOM+NT GO TO 67 61 IF(IPOM.LE.NT)GO TO 62 IPOM=IPOM-NT GO TO 61 62 CONTINUE IIAUX=NT-IPOM DO 63 J=1,NT IF(J.GE.IPOM)JAUX=J-IPOM+1 IF(J.LT.IPOM)JAUX=J+IIAUX+1 63 FF(JAUX)=GG(J)*AREDUC TLEN=asynt(2)-asynt(1) NFIN=IFIX(TLEN/DT)+1 CALL REDUC(AREDUC,1,NFIN,FF,GG,IS) 66 CONTINUE IF(NPR3.EQ.1)CALL ZPLOT(TMINO,DT,1,NFIN,GG) IF(NPR3.EQ.2)CALL ZPLOT(TMINO,DT,1,NFIN,FF) 70 CONTINUE C 40 CONTINUE IF(NPLOT.NE.0)CALL PLOT(0.,0.,999) C C END OF THE LOOP FOR RECEIVERS C 101 FORMAT(16I5) 102 FORMAT(26i3) 106 FORMAT(6F10.5,4I5) 107 FORMAT(20I4) 108 FORMAT(1X,20I4) 110 FORMAT(A) 111 FORMAT(1X,A) 114 FORMAT(/1X,'FREQUENCY RESPONSE-AMPLITUDES X(KM)=',F10.5/,2X, 1'DF=',F10.5,2X,'MAX=',E16.7,2X,'N=',I5) 115 FORMAT(/1X,'AMPLITUDE SPECTRUM OF SYNTHETIC SEISMOGRAM X(KM)=', 1F10.5/2X,'DF=',F10.5,2X,'MAX=',E16.8,2X,'N=',I5) 116 FORMAT(/1X,'SYNTHETIC SEISMOGRAM X(KM)=',F10.5/2X,'T0=', 1F10.5,2X,'DT=',F10.5,2X,'MAX=',E16.8,2X,'N=',I5) 152 FORMAT(5F10.5,2E15.7) 154 FORMAT(3E16.8) 160 FORMAT(F10.3,E12.5) 161 FORMAT(12I6) 162 FORMAT(2F10.3,1E12.5,I5) 170 FORMAT(1X,3I5) REWIND LU8 REWIND LU9 STOP END C C SUBROUTINE BORD character*80 TITLE COMMON/PLOT1/ARESP(6),ASPECT(6),ASYNT(6),JRESP(4),JSPECT(4),JSYNT( 14),NOO,RMIN,RMAX,GLR,TRMIN,TRMAX,GLTR,IR,ITR,NRED,NSH,AREDUC,TITLE 2,NPLOT IF(NOO.GT.1)GO TO 1 RMIN=ARESP(1) RMAX=ARESP(2) GLR=ARESP(3) TRMIN=ARESP(4) TRMAX=ARESP(5) GLTR=ARESP(6) IR=JRESP(1) ITR=JRESP(2) ndecx=jresp(3) ndecy=jresp(4) GO TO 3 1 IF(NOO.EQ.2.OR.NOO.GE.6)GO TO 2 RMIN=ASPECT(1) RMAX=ASPECT(2) GLR=ASPECT(3) TRMIN=ASPECT(4) TRMAX=ASPECT(5) GLTR=ASPECT(6) IR=JSPECT(1) ITR=JSPECT(2) ndecx=jspect(3) ndecy=jspect(4) GO TO 3 2 RMIN=ASYNT(1) RMAX=ASYNT(2) GLR=ASYNT(3) TRMIN=ASYNT(4) TRMAX=ASYNT(5) GLTR=ASYNT(6) IR=JSYNT(1) ITR=JSYNT(2) ndecx=jsynt(3) ndecy=jsynt(4) 3 CONTINUE CALL SYMBOL(-1.,0.,.5,TITLE,90.,80) CALL PLOT(0.,0.,3) if(trmax.ge.0.)CALL NUMBER(0.2,-1.5,0.4,TRMAX,90.,ndecy) if(trmax.lt.0.)CALL NUMBER(0.2,-2.,0.4,TRMAX,90.,ndecy) CALL PLOT(0.,0.,3) ZR=(RMAX-RMIN)/FLOAT(IR) ZTR=(TRMAX-TRMIN)/FLOAT(ITR) A=GLR/FLOAT(IR) B=GLTR/FLOAT(ITR) DO 100 J=1,ITR F=B*FLOAT(J) TG=TRMAX-ZTR*FLOAT(J) CALL PLOT(F,0.,2) CALL PLOT(F,0.2,2) if(tg.ge.0.)CALL NUMBER(F+0.2,-1.5,0.4,TG,90.,ndecy) if(tg.lt.0.)CALL NUMBER(F+0.2,-2.,0.4,TG,90.,ndecy) CALL PLOT(F,0.,3) 100 CONTINUE CALL NUMBER(F+0.8,-.2,0.4,RMIN,90.,ndecx) IF(NOO.EQ.1)CALL SYMBOL(F+2.1,0.,0.4,'FREQUENCY RESPONSE',90.,18) IF(NOO.EQ.2)CALL SYMBOL(F+2.1,0.,0.4,'INPUT SIGNAL',90.,12) IF(NOO.EQ.3)CALL SYMBOL(F+2.1,0.,0.4,'AMPLITUDE SPECTRUM OF INPUT 1SIGNAL',90.,34) IF(NOO.EQ.4)CALL SYMBOL(F+2.1,0.,0.4, 1'AMPLITUDE SPECTRUM, WINDOW APPLIED',90.,34) IF(NOO.EQ.5)CALL SYMBOL(F+2.1,0.,0.4,'AMPLITUDE SPECTRUM',90.,18) IF(NOO.EQ.5)CALL SYMBOL(F+2.7,0.,0.4, 1'OF SYNTHETIC SEISMOGRAM',90.,23) IF(NOO.EQ.6)CALL SYMBOL(F+2.1,0.,0.4,'SYNTHETIC SEISMOGRAM',90., 120) IF(NOO.EQ.7)CALL SYMBOL(F+2.1,0.,0.4, 1'INPUT SIGNAL AFTER WINDOWING',90.,28) IF(NOO.EQ.7)NOO=2 CALL PLOT(F,0.,3) DO 200 J=1,IR C=A*FLOAT(J) CALL PLOT(F,C,2) CALL PLOT(F-0.2,C,2) ZG=RMIN+ZR*FLOAT(J) CALL NUMBER(F+0.8,C-.6,0.4,ZG,90.,ndecx) CALL PLOT(F,C,3) 200 CONTINUE IF(NOO.NE.2.AND.NOO.NE.6)CALL SYMBOL(F+1.4,C-2.,0.4, 1'F(HZ)',90.,5) IF(NOO.EQ.2.OR.NOO.EQ.6)CALL SYMBOL(F+1.4,C-2.4,0.4, 1'TIME(S)',90.,7) CALL PLOT(F,C,3) DO 300 J=1,ITR D=B*FLOAT(ITR-J) CALL PLOT(D,C,2) CALL PLOT(D,C-0.2,2) 300 CALL PLOT(D,C,2) IF(NRED.EQ.1) 1CALL SYMBOL(D-0.4,0.,0.3,'REDUCTION FACTOR= ',90.,18) IF(NRED.EQ.1) 1CALL NUMBER(D-0.4,5.4,0.3,AREDUC,90.,5) CALL PLOT(D,C,3) DO 400 J=1,IR E=A*FLOAT(IR-J) CALL PLOT(D,E,2) CALL PLOT(D+0.2,E,2) 400 CALL PLOT(D,E,2) RETURN END C C SUBROUTINE ZPLOT(FMIN,DF,NMIN,NMAX,F) character*80 TITLE DIMENSION F(1) COMMON/PLOT1/ARESP(6),ASPECT(6),ASYNT(6),JRESP(4),JSPECT(4),JSYNT( 14),NOO,RMIN,RMAX,GLR,TRMIN,TRMAX,GLTR,IR,ITR,NRED,NSH,AREDUC,TITLE 2,NPLOT SHIFT=FLOAT(NSH) CALL BORD GMR=GLR/(RMAX-RMIN) GMTR=GLTR/(TRMAX-TRMIN) C DO 1 II=NMIN,NMAX VO=GMR*(FMIN-RMIN+DF*FLOAT(II-1)) IF((VO+0.0001).GE.0.)GO TO 2 1 CONTINUE GO TO 99 C 2 IMAX1=0 WO=GMTR*(TRMAX-F(II)) VOS=VO WOS=WO IF(WO.GT.0.)GO TO 3 IMAX1=1 CALL PLOT(0.,VO,3) GO TO 4 3 CALL PLOT(WO,VO,3) 4 CONTINUE II=II+1 IMAX2=0 C DO 10 I=II,NMAX VO=GMR*(FMIN-RMIN+DF*FLOAT(I-1)) WO=GMTR*(TRMAX-F(I)) IF(VO.GT.(GLR+0.0001))GO TO 11 IF(WO.LT.0.)GO TO 12 IF(IMAX1.EQ.0)GO TO 13 VV=VOS-WOS*(VO-VOS)/(WO-WOS) CALL PLOT(0.,VV,2) CALL PLOT(WO,VO,2) IMAX1=0 GO TO 9 13 CALL PLOT(WO,VO,2) IMAX1=0 GO TO 9 12 IF(IMAX1.EQ.0)GO TO 15 CALL PLOT(0.,VO,2) IMAX1=1 GO TO 9 15 VV=VOS-WOS*(VO-VOS)/(WO-WOS) CALL PLOT(0.,VV,2) CALL PLOT(0.,VO,2) IMAX1=1 GO TO 9 11 IF(IMAX1.EQ.0)GO TO 16 CALL PLOT(0.,RMAX,2) GO TO 99 16 IF(WO.LT.0.)GO TO 18 17 WW=WOS+(WO-WOS)*(GLR-VOS)/(VO-VOS) CALL PLOT(WW,GLR,2) GO TO 99 18 VV=VOS-WOS*(VO-VOS)/(WO-WOS) IF(VV.GT.GLR)GO TO 17 CALL PLOT(0.,VV,2) CALL PLOT(0.,GLR,2) GO TO 99 9 VOS=VO WOS=WO 10 CONTINUE CALL PLOT(GMTR*TRMAX,VOS,2) CALL PLOT(GMTR*TRMAX,GLR,2) C 99 continue CALL PLOT(GLTR+FLOAT(NSH),0.,-3) RETURN END C C SUBROUTINE SIGNAL(lin,lou) character*80 TITLE COMPLEX S(2048),REF(2048) COMMON/SIG/F(2048),TO,DT,NPTS,NO,NSIG,NPL1,NPL2,NPL3, 1NPL4,IPR,nder,nint COMMON/SPECTR/S,REF,DF,NI,NT,NT2,NTM,ILS,IRS,NFS,NFILT,NWIN COMMON/AUXIL/FF(2048),GG(2048),IS(2048),AUX(4),IAUX(4) COMMON/PLOT1/ARESP(6),ASPECT(6),ASYNT(6),JRESP(4),JSPECT(4),JSYNT( 14),NOO,RMIN,RMAX,GLR,TRMIN,TRMAX,GLTR,IR,ITR,NRED,NSH,AREDUC,TITLE 2,NPLOT C NPL1=0 NPL2=0 NPL3=0 NPL4=0 IF(NPLOT.EQ.0)GO TO 3 READ(lin,100) NPL1,NPL2,NPL3,NPL4 WRITE(lou,100)NPL1,NPL2,NPL3,NPL4 3 READ(lin,104)IAUX,AUX,TO WRITE(lou,104)IAUX,AUX,TO DO 1 I=1,NPTS F(I)=0. 1 GG(I)=0. c IF(NSIG)30,20,30 30 READ(lin,100)(IS(I),I=1,NPTS) WRITE(lou,120)NSIG WRITE(lou,102)(IS(I),I=1,NPTS) A=10.**(-NSIG) DO 32 I=1,NPTS F(I)=A*FLOAT(IS(I)) IS(I)=0 32 CONTINUE GO TO 60 20 CALL SFUN 60 CONTINUE C CALL REDUC(AREDUC,1,NPTS,F,GG,IS) IF(IPR.EQ.0)GO TO 61 WRITE(lou,121)TO,DT,AREDUC,NPTS WRITE(lou,107)(IS(I),I=1,NPTS) 61 IF(NPL1.EQ.0)GO TO 10 NOO=2 NRED=NPL1 IF(NPL1.EQ.1)CALL ZPLOT(TO,DT,1,NPTS,GG) IF(NPL1.EQ.2)CALL ZPLOT(TO,DT,1,NPTS,F) C 10 CONTINUE RETURN 100 FORMAT(16I5) 102 FORMAT(2X,16I5) 104 FORMAT(4I5,5F10.5) 107 FORMAT(20I4) 120 FORMAT(2X,'SIGNAL. NSIG=',I3) 121 FORMAT(/,1X,'INPUT SIGNAL T0=',F10.5,2X,'DT=',F10.5,2X, 1'MAX=',E16.8,2X,'N=',I5) END C C SUBROUTINE FT(lin,lou) character*80 TITLE COMPLEX S(2048),REF(2048),saux COMMON/SIG/F(2048),TO,DT,NPTS,NO,NSIG,NPL1,NPL2,NPL3, 1NPL4,IPR,nder,nint COMMON/SPECTR/S,REF,DF,NI,NT,NT2,NTM,ILS,IRS,NFS,NFILT,NWIN COMMON/AUXIL/FF(2048),GG(2048),IS(2048),AUX(4),IAUX(4) COMMON/WIND/FLO,FLEFT,FRO,FRIGHT,FEXP,ILF,IRF,NFW COMMON/PLOT1/ARESP(6),ASPECT(6),ASYNT(6),JRESP(4),JSPECT(4),JSYNT( 14),NOO,RMIN,RMAX,GLR,TRMIN,TRMAX,GLTR,IR,ITR,NRED,NSH,AREDUC,TITLE 2,NPLOT C IF(NWIN.EQ.0)GO TO 1 READ(lin,101)FLO,FLEFT,FRIGHT,FRO,FEXP WRITE(lou,101)FLO,FLEFT,FRIGHT,FRO,FEXP 1 CONTINUE C C PREPARE FOR FCOOLR DO 10 I=1,NPTS 10 S(I)=CMPLX(F(I),0.0) IF(NT.EQ.NPTS)GO TO 11 NPTS1=NPTS+1 DO 12 I=NPTS1,NT 12 S(I)=(0.0,0.0) 11 CONTINUE CALL FCOOLR(NI,S,1.0) DO 21 I=1,NT A=REAL(S(I)) B=AIMAG(S(I)) 21 FF(I)=SQRT(A*A+B*B) c c derivative (for nder.ne.0) or integral (for nint.ne.0) c of the input signal c if(nder.ne.0.or.nint.ne.0)then do 22 i=2,nt aux1=float(i-1)*df if(nder.ne.0)saux=cmplx(0.,-aux1) if(nint.ne.0)saux=cmplx(0.,1./aux1) s(i)=s(i)*saux 22 continue end if C C PLOTTING AMPLITUDE SPECTRUM OF THE INPUT SIGNAL C CALL REDUC(AREDUC,1,IRS,FF,GG,IS) IF(IPR.EQ.0)GO TO 25 WRITE(lou,120)DF,AREDUC,IRS WRITE(lou,107)(IS(I),I=1,IRS) 25 CONTINUE IF(NPL2.EQ.0)GO TO 28 NOO=3 NRED=NPL2 IF(NPL2.EQ.1)CALL ZPLOT(0.,DF,1,IRS,GG) IF(NPL2.EQ.2)CALL ZPLOT(0.,DF,1,IRS,FF) 28 CONTINUE C C WINDOWING OF AMPLITUDE SPECTRUM C ILF=ILS+1 IRF=IRS-1 NFW=IRF-ILF+1 IF(NWIN.EQ.0.AND.NFILT.EQ.0)GO TO 30 IF(NWIN.NE.0)CALL FENSTR IF(NFILT.NE.0)CALL FILTER DO 31 I=1,IRS A=REAL(S(I)) B=AIMAG(S(I)) 31 FF(I)=SQRT(A*A+B*B) CALL REDUC(AREDUC,1,IRS,FF,GG,IS) IF(IPR.EQ.0)GO TO 32 WRITE(lou,121)DF,AREDUC,IRS WRITE(lou,107)(IS(I),I=1,IRS) 32 CONTINUE IF(NPL3.EQ.0)GO TO 30 NOO=4 NRED=NPL3 IF(NPL3.EQ.1)CALL ZPLOT(0.,DF,1,IRS,GG) IF(NPL3.EQ.2)CALL ZPLOT(0.,DF,1,IRS,FF) 30 CONTINUE C C INPUT SIGNAL, FREQUENCY WINDOW AND FILTER APPLIED C IF(IPR.EQ.0.AND.NPL4.EQ.0)GO TO 50 DO 40 I=1,NT 40 REF(I)=S(I) DEL=1./FLOAT(NT) CALL FCOOLR(NI,REF,-1.0) DO 41 I=1,NT 41 FF(I)=DEL*REAL(REF(I)) CALL REDUC(AREDUC,1,NPTS,FF,GG,IS) IF(IPR.EQ.0)GO TO 42 WRITE(lou,122)TO,DT,AREDUC,NPTS WRITE(lou,107)(IS(I),I=1,NPTS) 42 CONTINUE IF(NPL4.EQ.0)GO TO 50 NRED=NPL4 NOO=7 IF(NPL4.EQ.1)CALL ZPLOT(TO,DT,1,NPTS,GG) IF(NPL4.EQ.2)CALL ZPLOT(TO,DT,1,NPTS,F) 50 CONTINUE RETURN C 101 FORMAT(8F10.5) 107 FORMAT(20I4) 120 FORMAT(/1X,'AMPLITUDE SPECTRUM OF THE INPUT SIGNAL'/, 12X,'DF=',F10.5,2X,'MAX=',E16.8,2X,'N=',I5) 121 FORMAT(/1X,'AMPLITUDE SPECTRUM OF INPUT SIGNAL,WINDOW APPLIED'/ 1,2X,'DF=',F10.5,2X,'MAX=',E16.8,2X,'N=',I5) 122 FORMAT(/1X,'INPUT SIGNAL, FREQUENCY WINDOW APPLIED'/2X,'T0=', 1F10.5,2X,'DT=',F10.5,2X,'MAX=',E16.8,2X,'N=',I5) END C C C C SUBROUTINE FENSTR COMPLEX S(2048),REF(2048) COMMON/SPECTR/S,REF,DF,NI,NT,NT2,NTM,ILS,IRS,NFS,NFILT,NWIN COMMON/WIND/FLO,FLEFT,FRO,FRIGHT,FEXP,ILF,IRF,NFW C ILO=FLO/DF+1.1 IRO=FRO/DF+1.1 IF(IRO.GT.IRS)IRO=IRS IF(ILO.LT.ILS)ILO=ILS IF(IRO.GT.NTM)IRO=NTM ILEFT=FLEFT/DF+1.1 IRIGHT=FRIGHT/DF+1.1 IF(IRIGHT.GT.IRO)IRIGHT=IRO IF(ILEFT.LT.ILO)ILEFT=ILO IF(IRIGHT.GT.NTM)IRIGHT=NTM NLEFT=ILEFT-ILO NRIGHT=IRO-IRIGHT IF(NLEFT.GT.0)DLEFT=3.14159/FLOAT(NLEFT) IF(NRIGHT.GT.0)DRIGHT=3.14159/FLOAT(NRIGHT) C DO 1 I=1,NTM J=NT+2-I FIF=I-IRIGHT FAF=I-ILO IF(I.GE.ILEFT.AND.I.LE.IRIGHT)GO TO 1 IF(I.LE.ILO.OR.I.GE.IRO)S(I)=(0.0,0.0) IF(I.GT.ILO.AND.I.LT.ILEFT) 1S(I)=S(I)*(0.5*(1.-COS(DLEFT*FAF)))**FEXP IF(I.GT.IRIGHT.AND.I.LT.IRO) 1S(I)=S(I)*(0.5*(COS(DRIGHT*FIF)+1.))**FEXP IF(I.EQ.1.OR.I.EQ.NTM)GO TO 1 S(J)=CONJG(S(I)) 1 CONTINUE ILF=ILO+1 IRF=IRO-1 NFW=IRF-ILF+1 RETURN END C C SUBROUTINE SFUN COMMON/AUXIL/FF(2048),GG(2048),IS(2048),AUX(4),IAUX(4) COMMON/SIG/F(2048),TO,DT,NPTS,NO,NSIG,NPL1,NPL2,NPL3, 1NPL4,IPR,nder,nint C IF(IAUX(1).NE.1)GO TO 10 C C GABOR SIGNAL C DO 1 I=1,NPTS T=TO-AUX(4)+DT*FLOAT(I-1) G=6.2831853*T*AUX(1) G1=G/AUX(2) G2=G1*G1 F(I)=0. IF(G2.LT.20.)F(I)=EXP(-G2)*COS(G+AUX(3)) 1 CONTINUE RETURN 10 CONTINUE C IF(IAUX(1).NE.2)GO TO 20 C C BERLAGE SIGNAL C DO 2 I=1,NPTS T=TO-AUX(4)+DT*FLOAT(I-1) G=6.2831853*T*AUX(1) G1=AUX(2)*T F(I)=0. IF(G1.LT.20..AND.T.GT.0.)F(I)=EXP(-G1)*SIN(G)*T**AUX(3) 2 CONTINUE RETURN 20 CONTINUE IF(IAUX(1).NE.3)GO TO 30 C C MULLER SIGNAL C G1=FLOAT(IAUX(2)) G2=G1/(G1+2) G3=G1+2 G1=3.14159*G1 G3=3.14159*G3 DO 3 I=1,NPTS T=TO-AUX(2)+DT*FLOAT(I-1) F(I)=0. IF(T.GT.0.AND.T.LE.AUX(1)) 1F(I)=SIN(G1*T/AUX(1))-G2*SIN(G3*T/AUX(1)) 3 CONTINUE RETURN 30 CONTINUE IF(IAUX(1).NE.4)GO TO 40 C C RICKER SIGNAL C DO 4 I=1,NPTS T=TO-AUX(2)+DT*FLOAT(I-1) G1=AUX(1)*T G2=2.*G1 G3=G1*G1 F(I)=0. IF(G3.LT.20.)F(I)=(1.-0.5*G2*G2)*EXP(-G3) 4 CONTINUE RETURN 40 CONTINUE IF(IAUX(1).NE.5)GO TO 50 C C BOX-CAR FUNCTION C DO 5 I=1,NPTS T=TO+DT*FLOAT(I-1) F(I)=0. IF(T.GE.AUX(1).AND.T.LE.AUX(2))F(I)=AUX(3) 5 CONTINUE RETURN 50 CONTINUE IF(IAUX(1).NE.6) GO TO 60 C C RAMP FUNCTION C DO 6 I=1,NPTS F(I)=0. T=TO+DT*FLOAT(I-1) IF(T.LT.AUX(3).AND.T.GE.AUX(2))F(I)=AUX(4) IF(T.LE.AUX(2).AND.T.GE.AUX(1)) 1F(I)=AUX(4)*(T-AUX(1))/(AUX(2)-AUX(1)) 6 CONTINUE RETURN 60 CONTINUE IF(IAUX(1).NE.7)GO TO 70 C C TRIANGLE FUNCTION C DO 7 I=1,NPTS F(I)=0. T=TO+DT*FLOAT(I-1) IF(T.LE.AUX(2).AND.T.GE.AUX(1)) 1F(I)=AUX(4)*(T-AUX(1))/(AUX(2)-AUX(1)) IF(T.LE.AUX(3).AND.T.GE.AUX(2)) 1F(I)=AUX(4)*(AUX(3)-T)/(AUX(3)-AUX(2)) 7 CONTINUE RETURN 70 CONTINUE IF(IAUX(1).NE.8)GO TO 80 C C ONE SAMPLE FUNCTION C DO 8 I=1,NPTS F(I)=0. IF(I.EQ.IAUX(2))F(I)=AUX(1) 8 CONTINUE RETURN 80 CONTINUE C RETURN END C C C SUBROUTINE FCOOLR(K,D,SN) C FAST FOURIER TRANSFORM OF N = 2**K COMPLEX DATA POINTS C REPARTS HELD IN D(1,3,...2N-1), IMPARTS IN D(2,4,...2N). C SUBROUTINE FCOOLR(K,D,SN) DIMENSION INU(15),D(4096) LX=2**K Q1=LX IL=LX SH=SN*6.28318530718/Q1 DO 10 I=1,K IL=IL/2 10 INU(I)=IL NKK=1 DO 40 LA=1,K NCK=NKK NKK=NKK+NKK LCK=LX/NCK L2K=LCK+LCK NW=0 DO 40 ICK=1,NCK FNW=NW AA=SH*FNW W1=COS(AA) W2=SIN(AA) LS=L2K*(ICK-1) DO 20 I=2,LCK,2 J1=I+LS J=J1-1 JH=J+LCK JH1=JH+1 Q1=D(JH)*W1-D(JH1)*W2 Q2=D(JH)*W2+D(JH1)*W1 D(JH)=D(J)-Q1 D(JH1)=D(J1)-Q2 D(J)=D(J)+Q1 20 D(J1)=D(J1)+Q2 DO 29 I=2,K ID=INU(I) IL=ID+ID IF(NW-ID-IL*(NW/IL)) 40,30,30 30 NW=NW-ID 29 CONTINUE 40 NW=NW+ID NW=0 DO 6 J=1,LX IF(NW-J) 8,7,7 7 JJ=NW+NW+1 J1=JJ+1 JH1=J+J JH=JH1-1 Q1=D(JJ) D(JJ)=D(JH) D(JH)=Q1 Q1=D(J1) D(J1)=D(JH1) D(JH1)=Q1 8 DO 9 I=1,K ID=INU(I) IL=ID+ID IF(NW-ID-IL*(NW/IL)) 6,5,5 5 NW=NW-ID 9 CONTINUE 6 NW=NW+ID RETURN END C C SUBROUTINE REDUC(AREDUC,N1,N2,F,G,IS) DIMENSION F(1),G(1),IS(1) DO 1 I=N1,N2 IS(I)=0 1 G(I)=0. AREDUC=0. DO 2 I=N1,N2 2 IF(ABS(F(I)).GT.AREDUC)AREDUC=ABS(F(I)) IF(AREDUC.EQ.0.)GO TO 5 B=1./AREDUC C=999.1*B GO TO 6 5 B=0. C=0. 6 CONTINUE DO 3 I=N1,N2 G(I)=B*F(I) 3 IS(I)=IFIX(C*F(I)) RETURN END C C SUBROUTINE FILTER COMPLEX S(2048),REF(2048) COMMON/SPECTR/S,REF,DF,NI,NT,NT2,NTM,ILS,IRS,NFS,nfilt,nwin RETURN END