C C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CTITLESAMNPHZ -- MINIMUM PHASE SPECTRUM OF GIVEN AMPLITUDE SPECTRUM 00020001 CA CA AUTHOR B. S. BOK 00030001 CA DESIGNER B. S. BOK 00040001 CA LANGUAGE FORTRAN 00050001 CA SYSTEM IBM / CRAY 00060001 CA WRITTEN AUGUST, 1990 00070001 C REVISED 7/03-91 BSB CHANGED WHITE NOISE FACTOR C REVISED 12-21-91 JJC MODIFIED TO MEET EDP STANDARDS. CA 00090001 CA CALL SAMNPHZ (A, MA, P, MP, NN, NL, IOPT, DLAY, WK, NWK, IER) 0000000 CA 00110001 CA IN/OUT A = AMPLITUDE SPECTRUM (MA*NN) R4 00120001 CA INPUT MA = COLUMN LENGTH OF A I4 00130001 CA OUTPUT P = COMPUTED MINIMUM PHASE SPECTRUM (MP*NN) R4 CA INPUT MP = COLUMN LENGTH OF P TO STORE PHASE SPECTRUMI4 00140001 CA INPUT NN = FFT LENGTH OF 2**NL I4 CA INPUT NL = FFT LENGTH OF MAGNITUDE I4 00150001 CA INPUT IOPT = >0, IF TO TAKE OUT TIME DELAY (SLOPE) I4 00160001 CA = 0 , IF NOT TO TAKE OUT SLOPE FROM PHASE CA OUTPUT DLAY = TIME DELAY R4 CA (SLOPE OF PHASE SPECTRUM TIMES (NYQ-1)/PI) CA IN/OUT WK = WORK ARRAY OF LENGTH AT LEAST 2*NN + 2 R4 CA IN/OUT NWK = SIZE OF WORK ARRAY AVAILABLE I4 CA INPUT IER = 0; IF WORK ARRAY SIZE ENOUGH I4 00170001 CA CA THIS SUBROUTINE COMPUTES THE MINIMUM PHASE SPECTRUM OF A GIVEN CA AMPLITUDE SPECTRUM. CA C REMARKS: INPUT AMPLITUDE SPECTRUM WILL BE SLIGHTLY ALTERED AFTER C MINPHZ CALL TO MINIMIZE ARTIFACTS IN WAVELET C************************************************************** C * C SUBROUTINES AND FUNCTIONS CALLED FROM THIS ROUTINE * C * C SAHILBT MFORSP SAPOLAR SAUNWRP SAPHZFT SABWRAP * C * C************************************************************** SUBROUTINE SAMNPHZ (A,MA, P,MP, NN,NL, IOPT,DLAY, WK, NWK, IER) 0000000 C IMPLICIT INTEGER (A-Z) C DIMENSION A(MA,NN), P(MP,NN), WK(NWK) 0000000 C REAL A REAL ALOG REAL AMAX REAL AMEAN REAL DLAY REAL HOLD REAL P REAL PHI REAL SUM REAL WK C C CHECK FOR WORK ARRAY SIZE 0000000 C 0000000 IF (NWK .GE. 2*NN+2) GO TO 100 IER = 99 RETURN C C ***************** READY FOR MIN PHASE CALCULATION *************** C 100 CONTINUE IER = 0 NYQ = NN / 2 + 1 C C GET AVERAGE AMPLITUDE IN EACH FREQUENCY C AMEAN = 0.0 DO 120 N = 1, NYQ AMEAN = AMEAN + A(1, N) ** 2 120 CONTINUE AMEAN = AMEAN / NYQ + 1.0E-30 C C TAKE LOGARITHM OF A, GET THE PHASE AND WRAP INTO (-PI,PI) 0000000 C 0000000 DO 140 N = 1, NN WK(N) = ALOG ( A(1, N) + 0.001 * AMEAN ) 140 CONTINUE CALL SAHILBT ( WK, WK(NN + 1), NN, NL ) C C ZERO OUT PHASE FOR AMP BELOW 0.1 % OF MAX AMP C AT BOTH ENDS OF SPECTRUM C AMAX = 0.0 DO 160 N = 2, NYQ AMAX = AMAX1 ( A(1, N), AMAX ) 160 CONTINUE HOLD = 0.001 * AMAX DO 180 N = 2, NYQ IF (A(1, N) .GT. HOLD) GO TO 200 WK(NN + N) = 0.0 180 CONTINUE 200 CONTINUE DO 220 N = NYQ - 1, 2, -1 IF (A(1, N) .GT. HOLD) GO TO 240 WK(NN + N) = 0.0 220 CONTINUE 240 CONTINUE C 0000000 C MOVE PHASE INTO P AND MAKE ANTISYMMETRIC AROUND NYQUIST FREQUENCY 0000000 C 0000000 P(1, 1) = WK(NN + 1) DO 260 N = 2, NYQ P(1, NN - N + 2) = -WK(NN + N) P(1, N) = WK(NN + N) 260 CONTINUE C C REFINE AMP/PHASE TO GET WELL-BEHAVING WAVELET C FIRST GET CRUDE MIN-PHASE WAVELET C DO 280 N = 1, NN WK(N) = A(1, N) WK(NN + N) = P(1, N) 280 CONTINUE DO 300 N = 1, NN HOLD = WK(N) * COS ( WK(NN + N) ) WK(NN + N) = -WK(N) * SIN ( WK(NN + N) ) WK(N) = HOLD 300 CONTINUE CJJ CALL BBFFT ( WK(1), WK(NN + 1), NN, NL, 1 ) CALL MFORSP (NL, WK(1), WK(1+NN), 1 ) C C NEXT TRUNCATE TRACE AT 99 PERCENT OF TOTAL ENERGY C SUM = 0.0 DO 320 N = 1, NYQ SUM = SUM + WK(N) * WK(N) WK(NN + N) = SUM 320 CONTINUE SUM = 0.99 * WK(NN + NYQ) DO 340 N = 1, NN IF (WK(NN + N) .GE. SUM) GO TO 360 340 CONTINUE N = NYQ 360 CONTINUE DO 380 M = N, NYQ WK(M) = 0.0 380 CONTINUE C C AND ZERO OUT PHASE PART BEFORE INVERSE FFT C DO 400 N = 1, NN WK(NN + N) = 0.0 400 CONTINUE C C FINALLY GET REFINED AMP-PHASE SPECTRUM C CJJ CALL BBFFT ( WK(1), WK(NN + 1), NN, NL, -1 ) CALL MFORSP(NL, WK(1), WK(1+NN), -1 ) C CALL POLAR ( WK(1), WK(NN + 1), NN, WK(1), WK(NN + 1) ) CALL SAPOLAR ( WK(1), 1, WK(NN + 1), 1, NN, WK(1), 1, 1 WK(NN + 1), 1 ) C 0000000 C TAKE OUT INTEGER TIME SAMPLE DELAY (SLOPE IN PHASE SPECTRUM) 0000000 C FROM PHASE SPECTRUM IN WK(NN+I) C 0000000 IF (IOPT .GT. 0) THEN CALL SAUNWRP ( WK(NN + 1), NYQ ) CALL SAPHZFT ( WK, WK(NN + 1), NYQ, 0.5, 0, PHI, DLAY, WK(NN + 1 NYQ + 1) ) ILAY = IFIX ( DLAY + 0.5 ) DO 420 N = 1, NYQ WK(NN + N) = WK(NN + N) + ILAY * (N - 1 ) * 3.141592654 / ( 1 NYQ - 1 ) 420 CONTINUE CALL SABWRAP ( WK(NN + 1), NYQ ) ENDIF C C MOVE AMP/P AND MAKE ANTISYMMETRIC AROUND NYQUIST FREQUENCY 0000000 C 0000000 A(1, 1) = WK(1) P(1, 1) = -WK(NN + 1) DO 440 N = 2, NYQ A(1, NN - N + 2) = WK(N) A(1, N) = WK(N) P(1, NN - N + 2) = WK(NN + N) P(1, N) = -WK(NN + N) 440 CONTINUE C C 0000000 RETURN END