CTITLESABWAVE - COMPUTE THE AUTOCORRELATION OF A BANDPASS WAVELET.
C***********************************************************************
C                 COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991            *
C***********************************************************************
C
C
CA    AUTHOR            H. W. SWAN
CA    DESIGNER          H. W. SWAN
CA    SYSTEM            IBM/CRAY
CA    LANGUAGE          VS FORTRAN VERSION 2.3
CA    WRITTEN           12-09-88
CA    LAST REVISED      01-03-90
CA
C     REVISED    HWS    01-03-90   ADDED DOCUMENTATION
C     REVISED    JJC    05-23-91   CHANGED IMPLICIT NONE TO (A-Z).
C     REVISED    JJC    12-17-91   MODIFIED TO MEET SPARC STANDARDS.
CA
CA    CALLING SEQUENCE:
CA
CA      CALL SABWAVE(FFT, NS, WAVELT, TWVLTP, MAXH, NH, FS, WPARS)
CA
CA IN/OUT ARGUMENT TYPE DESCRIPTION
CA ------ -------- ---- -----------
CA
CA  OUT     FFT    (R4) A SCRATCH ARRAY OF LENGTH 'NS+2', WHERE
CA                      'NS+1' IS A POWER OF 2.  USED FOR FFT.
CA
CA   IN     NS      I4  DETERMINES THE LENGTH OF THE FFT ARRAY.
CA
CA  OUT    WAVELT  (R4) THE AUTOCORRELATION OF THE WAVELET.
CA                      (DIMENSIONED 0:2*MAXH)
CA
CA  OUT    TWVLTP  (R4) THE CROSS CORRELATION OF W(T) AND TW'(T).
CA                      (DIMENSIONED -MAXH:MAXH)
CA
CA   IN    MAXH     I4  DETERMINES THE DIMENSIONED LENGTH OF THE
CA                      CORRELATION ARRAYS WAVELT AND TWVLTP.
CA
CA   IN     NH      R4  THE NUMBER OF CORRELATION LAG VALUES TO
CA                      ACTUALLY COMPUTE.
CA
CA   IN     FS      R4  THE SAMPLING FREQUENCY (HZ)
CA
CA   IN     WPARS  (R4) THE FOUR BANDPASS WAVELET CUTOFF FREQUENCIES:
CA                      THE WAVELET SPECTRA WILL HAVE GAUSSIAN TAPERS
CA                      WHOSE CUTOFF FREQUENCIES BEING 30 DB DOWN
CA                      FROM THE PASS BAND.   (UNITS ARE HZ.)
CA                            WPARS(1) = LOW CUTOFF FREQUENCY
CA                            WPARS(2) = LOW PASS FREQUENCY
CA                            WPARS(3) = HIGH PASS FREQUENCY
CA                            WPARS(4) = HIGH CUTOFF FREQUENCY
CA
CA  SUBROUTINES CALLED:
CA
CA      ARMVE   -  ARRAY MOVER
CA      S1FMAG  -  DETERMINE LOG2 OF THE FFT LENGTH.
CA      S2DFI2  -  COMPLEX TO REAL INVERSE FFT.
CA
CA
C
C  NOTE:                    WAVLET <-- W(T) * W(-T);
C                           TWVLTP <-- W(T) * {-TW'(-T)}.
C
        SUBROUTINE SABWAVE(FFT, NS, WAVELT, TWVLTP, MAXH, NH, FS,
     *      WPARS)
C
CJJ     IMPLICIT  NONE
        IMPLICIT INTEGER (A-Z)
C
        REAL      FFT, WAVELT, TWVLTP, WPARS, FS
        INTEGER   NS, MAXH, NH, LFOUR
        INTEGER   I, I2, MAG, IEXP, NPNTS
        DIMENSION FFT(0:NS+1), WAVELT(0:2*MAXH)
        DIMENSION TWVLTP(-MAXH:MAXH), WPARS(4)
        REAL      LC, LP, HP, HC, X0, F, RAT, ANS, Q, THRESH
        REAL      VALUE
        PARAMETER (IEXP=2)
        PARAMETER (THRESH=7.0)
        EXTERNAL  ARMVE, S1FMAG, S2DFI2
C
C  PICK UP THE WAVELET PARAMETERS.
C
        LC = WPARS(1)
        LP = WPARS(2)
        HP = WPARS(3)
        HC = WPARS(4)
C
C  FIRST, COMPUTE THE WAVELET.
C
        X0 = SQRT(1.5 * ALOG(10.0))
        DO 100 I=0, NS/2
            I2 = I * 2
            F  = I * FS / NS
            RAT = 0.0
            IF(F .LE. LP) THEN
                RAT = IEXP * X0 * (F - LP) / (LC - LP)
                ANS = 0.0
            ELSE IF(F .GT. HP) THEN
                RAT = IEXP * X0 * (F - HP) / (HC - HP)
                ENDIF
            ANS = 0.0
            IF(RAT .LE. THRESH) ANS = EXP(-RAT**2)
            FFT(I2) = ANS
            FFT(I2+1) = 0.0
 100        CONTINUE

        NPNTS = 2*MIN0(NH, MAXH) + 1
        CALL S1FMAG(NS, MAG, LFOUR)
        CALL S2DFI2(MAG, FFT, *1000)
        CALL ARMVE(FFT, WAVELT, NPNTS)
C
C  NEXT, COMPUTE THE WAVELET DERIVATIVE.
C
        DO 300 I=0, NS/2
            I2 = I * 2
            F  = I * FS / NS
            RAT = 0.0
            Q   = 0.0
            IF(F .LE. LP) THEN
                RAT = IEXP * X0 * (F - LP) / (LC - LP)
                Q   = 2.0*F*RAT*X0 / (LC - LP) / IEXP
                ANS = 0.0
            ELSE IF(F .GE. HP) THEN
                RAT = IEXP * X0 * (F - HP) / (HC - HP)
                Q   = 2.0*F*RAT*X0 / (HC - HP) / IEXP
                ENDIF
            ANS = 0.0
            IF(RAT .LE. THRESH)
     *          ANS = -(1 - Q) * EXP(-RAT**2)
            FFT(I2) = ANS
            FFT(I2+1) = 0.0
 300        CONTINUE

        CALL S2DFI2(MAG, FFT, *1000)

        NPNTS = MIN0(NH, MAXH)
        DO 400 I=0, NPNTS
            VALUE = FFT(I)
            TWVLTP(I) = VALUE
 400        TWVLTP(-I) = VALUE

        RETURN

 1000   WRITE (*,*) 'FFT LENGTH PROBLEM.'
        STOP
        END