CTITLESAMXPH -- MIXED PHASE WAVELETS 00000010 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CA AUTHOR D. NYMAN 00000020 CA DESIGNER D. NYMAN 00000030 CA LANGUAGE S/370 FORTRAN H 00000040 CA SYSTEM IBM AND CRAY CA WRITTEN 11/13/78 00000050 C REVISED MO/DA/YR BY PROGRAMMER FOR REASON 00000060 C REVISED 10/15/90 MAA. MODIFIED FOR CRAY USE. 00000070 CA 00000080 CA 00000090 CA CALL SAMXPH 00000100 CA (WAV, LENG, SMPINT, PHASE, U, NP, NX, V, NUSM, KPPRNT, IER) 00000110 CA 00000120 CA INPUT WAV = MINIMUM PHASE WAVELET ARRAY R4 00000130 CA INPUT LENG = NUMBER OF POINTS IN ARRAY WAV I4 00000140 CA INPUT SMPINT = SAMPLE INTERVAL I4 00000150 CA INPUT PHASE = ARRAY OF PHASE VALUES FOR MIXED PHASE I4 00000160 CA WAVELETS 00000170 CA OUTPUT U = WORK ARRAY OF LENGTH NP+2 FOR FOURIER R4 00000180 CA TRANSFORM 00000190 CA INPUT NP = FOURIER TRANSFORM LENGTH, A POWER OF 2 I4 00000200 CA INPUT NX = FOURIER TRANSFORM MAGNITUDE: NP=2**NX I4 00000210 CA OUTPUT V = WORK ARRAY R4 00000220 CA INPUT NUSM = LENGTH OF V I4 00000230 CA INPUT KPPRNT = FORTRAN LOGICAL UNIT NUMBER FOR PRINT I4 00000240 CA OUTPUT IER = 0: NORMAL TERMINATION I4 00000250 CA 1: FOURIER TRANSFORM SIN/COS TABLE INADEQUATE 00000260 CA 2: WORK ARRAY V TOO SHORT 00000270 CA 00000280 CA 00000290 CA SAMXPH CALCULATES MIXED PHASE WAVELETS FROM THE MINIMUM PHASE 00000300 CA WAVELET IN ARRAY WAV. THE MIXED PHASE WAVELETS ARE A MIXTURE OF 00000310 CA THE GIVEN MINIMUM PHASE WAVELET AND THE ZERO PHASE WAVELET WITH 00000320 CA THE SAME POWER SPECTRUM, WHERE THE PROPORTION OF EACH IS GIVEN BY 00000330 CA THE VALUES IN PHASE. THE MIXED PHASE WAVELET VALUES ARE PRINTED 00000340 CA AND THE WAVELETS ARE DISPLAYED VIA A PRINTER PLOT. PHASE VALUES 00000350 CA ARE IN THE RANGE -1.0 TO 1.0, WHERE -1.0 DENOTES MINIMUM PHASE, 00000360 CA 0.0 DENOTES ZERO PHASE, AND 1.0 DENOTES MAXIMUM PHASE. ANY 00000370 CA NUMBER OF PHASE VALUES MAY BY SPECIFIED. THE PROCESS TERMINATES 00000380 CA WHEN A PHASE VALUE OUT OF RANGE IS DETECTED. 00000390 CA 00000400 CA 00000410 C=======================================================================00000420 C EJECT 00000430 C=======================================================================00000440 C 00000450 C LOCAL VARIABLES AND CONSTANTS 00000460 C 00000470 C IC = INDEX TO START OF NEXT MIXED PHASE WAVELET IN V I4 00000480 C IC0 = INDEX TO START OF FIRST MIXED PHASE WAVELET IN V I4 00000490 C KPLT = NUMBER OF MIXED PHASE WAVELETS IN V I4 00000500 C LENV = LENGTH OF EACH MIXED PHASE WAVELET I4 00000510 C P = CONTINUOUS PHASE VS. FREQUENCY OF INPUT WAVELET R4 00000520 C PHS = PHASES OF WAVELETS IN V R4 00000530 C SN = SIGN OF ZERO FREQUENCY FOURIER TRANSFORM VALUE R4 00000540 C T1 = TIME OF FIRST VALUE OF MIXED PHASE WAVELETS I4 00000550 C 00000560 C=======================================================================00000570 C EJECT 00000580 C 00000590 SUBROUTINE SAMXPH (WAV, LENG, SMPINT, PHASE, U, 00000600 * NP, NX, V, NUSM, KPPRNT, IER) 00000610 C 00000620 EXTERNAL S1ATP 00000630 INTEGER SMPINT,T1 00000640 DIMENSION WAV(1),PHASE(1),U(1),V(1),PHS(3) 00000650 C 00000660 C SET MIXED PHASE WAVELET LENGTH AND START TIME 00000670 C 00000680 LENV=LENG*2-1 00000690 T1=(1-LENG)*SMPINT 00000700 C 00000710 C FOURIER TRANSFORM INPUT WAVELET 00000720 C 00000730 CALL ARMVE (WAV,U,LENG) 00000740 CALL ARSET (U(LENG+1),NP-LENG,0) 00000750 CALL S2DFT2 (NX,U,*60) 00000760 C 00000770 C SET SN, CHANGE SIGN OF FOURIER TRANSFORM IF SN < 0 00000780 C 00000790 SN=SIGN(1.,U(1)) 00000800 IF (SN.LT.0.) CALL ARREVF (U,U,NP) 00000810 C 00000820 C CALCULATE FOURIER AMPLITUDE AND CONTINUOUS PHASE 00000830 C 00000840 P0=0. 00000850 C 00000860 DO 10 J=3,NP,2 00000870 P=ATAN2(U(J+1),U(J)) 00000880 IF (ABS(P-P0).GT.3.14159265) P0=P0+SIGN(6.28318531,P-P0) 00000890 U(J)=SQRT(U(J)*U(J)+U(J+1)*U(J+1))*SN 00000900 U(J+1)=U(J-1)+P-P0 00000910 10 P0=P 00000920 C 00000930 C CALCULATE AND DISPLAY MIXED PHASE WAVELETS 00000940 C 00000950 IC0=NP+3 00000960 IC=IC0 00000970 KPLT=0 00000980 C 00000990 DO 40 I=1,1000 00001000 P0=-PHASE(I) 00001010 IF (ABS(P0).GT.1.) GO TO 30 00001020 IF (IC+LENV.GT.NUSM) GO TO 70 00001030 C 00001040 C RECONSTRUCT FOURIER TRANSFORM, BUT WITH 00001050 C A PHASE OF -PHASE(I) * ORIGINAL PHASE 00001060 C 00001070 DO 20 K=1,NP,2 00001080 P=U(K+1)*P0 00001090 V(K)=U(K)*COS(P) 00001100 20 V(K+1)=U(K)*SIN(P) 00001110 C 00001120 V(NP+1)=0. 00001130 V(NP+2)=0. 00001140 C 00001150 C INVERSE FOURIER TRANSFORM AND MOVE RESULT TO V(IC) 00001160 C 00001170 CALL S2DFI2 (NX,V,*60) 00001180 CALL ARMVE (V(NP+2-LENG),V(IC),LENG-1) 00001190 CALL ARMVE (V,V(IC+LENG-1),LENG) 00001200 C 00001210 KPLT=KPLT+1 00001220 JC=IC+LENV-1 00001230 IF (KPLT.EQ.1) WRITE (KPPRNT,2010) 00001240 WRITE (KPPRNT,2020) PHASE(I),(V(K),K=IC,JC) 00001250 PHS(KPLT)=PHASE(I) 00001260 IF (KPLT.LT.3) GO TO 40 00001270 C 00001280 C WHEN THREE WAVELETS HAVE BEEN 00001290 C GENERATED CALL PRINTER PLOT 00001300 C 00001310 30 IF (KPLT.EQ.0) GO TO 50 00001320 WRITE (KPPRNT,2030) (PHS(K),K=1,KPLT) 00001330 CALL USPPLT (V(IC0),LENV,KPLT,LENV,T1,SMPINT,1,3,KPPRNT) 00001340 IC=IC0-LENV 00001350 KPLT=0 00001360 IF (ABS(P0).GT.1.) GO TO 50 00001370 C 00001380 40 IC=IC+LENV 00001390 C 00001400 C NORMAL TERMINATION 00001410 C 00001420 50 IER=0 00001430 GO TO 100 00001440 C 00001450 C FOURIER TRANSFORM SIN/COS TABLE INADEQUATE 00001460 C 00001470 60 IER=1 00001480 GO TO 100 00001490 C 00001500 C WORK ARRAY V TOO SHORT 00001510 C 00001520 70 IER=2 00001530 C 00001540 100 RETURN 00001550 C 00001560 2010 FORMAT ('1MIXED PHASE WAVELETS -- TIME ZERO IS CENTER VALUE') 00001570 C 00001580 2020 FORMAT (//' WAVELET FOR RELATIVE PHASE OF'F7.3/(1X,10F10.4)) 00001590 C 00001600 2030 FORMAT (//3(F25.3,' PHASE'11X)) 00001610 C 00001620 END 00001630