Make your own free website on Tripod.com
A Coherence and Cross Spectral Estimation Program

Return to Main Software Page

C
C-----------------------------------------------------------------------
C MAIN PROGRAM:  A COHERENCE AND CROSS SPECTRAL ESTIMATION PROGRAM
C AUTHORS:       G. C. CARTER, J. F. FERRIE
C                NAVAL UNDERWATER SYSTEMS CENTER
C                NEW LONDON, CONNECTICUT 06320
C
C INPUT:         NNN IS THE NUMBER OF DATA POINTS PER SEGMENT
C                  4 < NNN < 1025
C                ISR IS THE SAMPLING RATE
C                NDSJP IS THE NUMBER OF DISJOINT SEGMENTS
C                SFX IS THE SCALE FACTOR FOR THE INPUT DATA STORED IN
C                  THE XX ARRAY
C                SFY IS THE SCALE FACTOR FOR THE INPUT DATA STORED IN
C                  THE YY ARRAY
C-----------------------------------------------------------------------
C
C SPECIFICATION AND TYPE STATEMENTS
C
      DIMENSION XX(1024), YY(1024)
      DIMENSION GXX(513), GYY(513), GXYRE(513), GXYIM(513)
      DIMENSION WEGHT(513), PHI(513)
      DIMENSION LINE(50)
      EQUIVALENCE (WEGHT(1),PHI(1))
C
C SET UP MACHINE CONSTANTS
C
      IOIN1 = I1MACH(1)
      IPRTR = I1MACH(2)
      SMALL = R1MACH(1)
C
C READ INPUT CONTROL PARAMETERS FROM COMPUTER DATA CARD
C
      READ (IOIN1,9999) NNN, ISR, NDSJP, SFX, SFY
C NNN IS THE NUMBER OF DATA POINTS PER SEGMENT
C ISR IS THE SAMPLING RATE
C NDSJP IS THE NUMBER OF DISJOINT SEGMENTS
C SFX AND SFY ARE SCALE FACTORS FOR THE INPUT DATA
9999  FORMAT (3I5, 2F10.5)
      NFFTS = NDSJP
C
C PRINT INPUT CONTROL PARAMETERS
C
      WRITE (IPRTR,9998) NNN, ISR, NDSJP, SFX, SFY
9998  FORMAT (/1X, 5HNNN =, I6, 5X, 5HISR =, I7, 5X, 7HNDSJP =, I7,
     *    5X//1X, 5HSFX =, E15.8, 8X, 5HSFY =, E15.8/)
C
C CALCULATE CONSTANTS
C
      TPI = 8.0*ATAN(1.0)
      DEG = 360.0/TPI
      IF (NNN.GT.0 .AND. NNN.LE.1024) GO TO 10
      WRITE (IPRTR,9997)
9997  FORMAT (10X, 9HNNN ERROR)
      STOP
  10  CONTINUE
      VARX = 0.0
      VARY = 0.0
      DT = 1.0/FLOAT(ISR)
      SF = SQRT(ABS(SFX*SFY))
C
C PRINT OUT USER INFORMATION
C
      TIME = FLOAT(NDSJP*NNN)*DT
      WRITE (IPRTR,9996) NDSJP, TIME
9996  FORMAT (10X, 3HTHE, I4, 25H DISJOINT PIECES COMPRISE, F8.2,
     *    16H SECONDS OF DATA)
C
C COMPUTE NEW COMPOSITE NUMBER NNN
C
      CALL HICMP(NNN, NPFFT)
      IF (NPFFT.GT.1024) STOP
      WRITE (IPRTR,9995) NPFFT
9995  FORMAT (10X, 21HNUMBER OF POINT FFT =, I5/)
C
C CALCULATE CONSTANTS
C
      NNNP1 = NNN + 1
      NNND2 = NNN/2
      NND21 = NNND2 + 1
      NP2 = NPFFT + 2
      ND2 = NPFFT/2
      ND2P1 = ND2 + 1
      DF = 1.0/(DT*FLOAT(NPFFT))
      FNYQ = FLOAT(ISR)/2.0
      CONST = 0.25*DT/FLOAT(NNN)
      FLOW = 0.0
      FHIGH = FNYQ
      ISTRT = IFIX(FLOW/DF) + 1
      ISTOP = IFIX(FHIGH/DF) + 1
C
C COMPUTE AND SAVE WEIGHTING FUNCTION
C
      TEMP = TPI/FLOAT(NNN+1)
      SCL = SQRT(2.0/3.0)
      DO 20 I=1,NNND2
        WEGHT(I) = SCL*(1.0-COS(TEMP*FLOAT(I)))
  20  CONTINUE
C
C STORE ZEROS IN THE SUMMING ARRAYS
C
      CALL ZERO(GXX, ND2P1)
      CALL ZERO(GYY, ND2P1)
      CALL ZERO(GXYRE, ND2P1)
      CALL ZERO(GXYIM, ND2P1)
C
C COMPUTE AND SUM NPFFT ESTIMATES
C
      DO 80 KOUNT=1,NFFTS
C
        CALL ZERO(XX, NPFFT)
        CALL ZERO(YY, NPFFT)
C
C LOAD XX AND YY ARRAYS WITH NNN DATA POINTS
C
        CALL LOAD(XX, YY, NNN, KOUNT, ISR)
C
C PRINT OF FIRST 50 INPUT VALUES
C
        IF (KOUNT.NE.1) GO TO 40
        WRITE (IPRTR,9994)
9994    FORMAT (1H1, 9X, 41HPRINTOUT OF FIRST 50 VALUES OF INPUT DATA///
     *      )
        LPMAX = MIN0(NPFFT,50)
        DO 30 I=1,LPMAX
          WRITE (IPRTR,9993) I, XX(I), YY(I)
9993      FORMAT (1X, I5, 1X, 2F15.8, 6X)
  30    CONTINUE
        WRITE (IPRTR,9992)
9992    FORMAT (/1H1)
  40    CONTINUE
C
C REMOVE THE LINEAR TREND AND COMPUTE THE VARIANCE
C IF IS3 = 0 DO NOT REMOVE DC COMPONENT OR SLOPE
C = 1 REMOVE THE DC COMPONENT
C > 1 REMOVE THE DC COMPONENT AND SLOPE
C
        IS3 = 0
        CALL LREMV(XX, NNN, IS3, DX, SX)
        CALL LREMV(YY, NNN, IS3, DY, SY)
        VARXI = 0.0
        VARYI = 0.0
        DO 50 I=1,NNN
          VARXI = VARXI + XX(I)*XX(I)
          VARYI = VARYI + YY(I)*YY(I)
  50    CONTINUE
        VARXI = VARXI/FLOAT(NNN-1)
        VARYI = VARYI/FLOAT(NNN-1)
        WRITE (IPRTR,9991) KOUNT, DX, DY, SX, SY, VARXI, VARYI, IS3
9991    FORMAT (1X, I3, 4H DX=, E12.5, 4H DY=, E12.5, 4H SX=, E12.5,
     *      4H SY=, E12.5/4H VX=, E12.5, 4H VY=, E12.5, I5)
        VARX = VARX + VARXI
        VARY = VARY + VARYI
C
C WEIGHT THE INPUT DATA WITH COSINE WINDOW
C
        DO 60 I=1,NNND2
          ITMP = NNNP1 - I
          XX(I) = XX(I)*WEGHT(I)
          YY(I) = YY(I)*WEGHT(I)
          XX(ITMP) = XX(ITMP)*WEGHT(I)
          YY(ITMP) = YY(ITMP)*WEGHT(I)
  60    CONTINUE
C
C COMPUTE FORWARD FFT
C
        CALL FFT842(0, NPFFT, XX, YY)
C
C COMPUTE SPECTRA
C
        GXX(1) = GXX(1) + 4.0*XX(1)**2
        DO 70 K=2,ND2P1
          J = NP2 - K
          GXX(K) = GXX(K) + (XX(K)+XX(J))**2 + (YY(K)-YY(J))**2
          GYY(K) = GYY(K) + (YY(K)+YY(J))**2 + (XX(J)-XX(K))**2
          GXYRE(K) = GXYRE(K) + XX(K)*YY(J) + XX(J)*YY(K)
          GXYIM(K) = GXYIM(K) + XX(J)**2 + YY(J)**2 - XX(K)**2 -
     *        YY(K)**2
  70    CONTINUE
        GYY(1) = GYY(1) + 4.0*YY(1)**2
        GXYRE(1) = GXYRE(1) + 2.0*(XX(1)*YY(1))
        GXYIM(1) = 0.0
C
C GO BACK FOR NEXT SEGMENT
C
  80  CONTINUE
C
C NORMALIZE ESTIMATES
C
      FNSG = FLOAT(NFFTS)
      OFNSG = 1.0/FNSG
      VARX = VARX*OFNSG
      VARY = VARY*OFNSG
      TEMP1 = CONST*OFNSG*SFX
      TEMP2 = CONST*OFNSG*SFY
      TEMP4 = CONST*OFNSG*SF
      TEMP3 = 2.0*TEMP4
      DO 90 K=1,ND2P1
        GXX(K) = GXX(K)*TEMP1
        GYY(K) = GYY(K)*TEMP2
        GXYRE(K) = GXYRE(K)*TEMP3
        GXYIM(K) = GXYIM(K)*TEMP4
  90  CONTINUE
C
C PRINT OUT VARIANCES
C
      WRITE (IPRTR,9990) VARX, VARY
9990  FORMAT (/10X, 26HAVERAGE VARIANCES ARE, VX=, E12.6, 4H VY=,
     *    E12.6//)
      VARX = 0.0
      VARY = 0.0
      DO 100 K=1,ND2P1
        VARX = VARX + GXX(K)
        VARY = VARY + GYY(K)
 100  CONTINUE
      VARX = VARX*DF*2.0/SFX
      VARY = VARY*DF*2.0/SFY
      WRITE (IPRTR,9989) VARX, VARY
9989  FORMAT (/10X, 29HINTEGRATED VARIANCES ARE, VX=, E12.6, 4H VY=,
     *    E12.6//)
C
C CONVERT GXX TO DB AND PLOT
C
      DO 110 I=1,ND2P1
        XX(I) = GXX(I)
        PHI(I) = 10.0*ALOG10(AMAX1(GXX(I),SMALL))
 110  CONTINUE
      WRITE (IPRTR,9988)
9988  FORMAT (1H1/15X, 25HPLOT OF AUTO SPECTRUM GXX)
      CALL PRPRT(PHI, LINE, DF, FLOW, FHIGH, IPRTR)
C
C COMPUTE AND DISPLAY AUTOCORRELATION FUNCTION OF INPUT SIGNAL XX
C
      CALL ZERO(YY, NPFFT)
      DO 120 K=2,ND2P1
        J = NP2 - K
        XX(J) = XX(K)
 120  CONTINUE
      CALL FFT842(1, NPFFT, XX, YY)
      ONDRO = 1.0/XX(1)
      DO 130 I=1,ND2P1
        XX(I) = XX(I)*ONDRO
 130  CONTINUE
      WRITE (IPRTR,9987)
9987  FORMAT (1H1/20X, 27HPLOT OF AUTOCORRELATION RXX)
      CALL COPLT(XX, ND2P1, DT, 1, IPRTR, LINE)
C
C CONVERT GYY TO DB AND PLOT
C
      DO 140 I=1,ND2P1
        XX(I) = GYY(I)
        PHI(I) = 10.0*ALOG10(AMAX1(GYY(I),SMALL))
 140  CONTINUE
      WRITE (IPRTR,9986)
9986  FORMAT (1H1/15X, 25HPLOT OF AUTO SPECTRUM GYY)
      CALL PRPRT(PHI, LINE, DF, FLOW, FHIGH, IPRTR)
C
C COMPUTE AND DISPLAY AUTOCORRELATION FUNCTION OF INPUT SIGNAL YY
C
      CALL ZERO(YY, NPFFT)
      DO 150 K=2,ND2P1
        J = NP2 - K
        XX(J) = XX(K)
 150  CONTINUE
      CALL FFT842(1, NPFFT, XX, YY)
      ONDRO = 1.0/XX(1)
      DO 160 I=1,ND2P1
        XX(I) = XX(I)*ONDRO
 160  CONTINUE
      WRITE (IPRTR,9985)
9985  FORMAT (1H1/20X, 27HPLOT OF AUTOCORRELATION RYY)
      CALL COPLT(XX, ND2P1, DT, 1, IPRTR, LINE)
C
C COMPUTE AND DISPLAY PHASE FROM AVERAGED GXYRE AND GXYIM SPECTRUM
C
      GXYIM(1) = 0.0
      PHI(1) = 0.0
      DO 200 K=2,ND2P1
        XXK = GXYRE(K)
        IF (XXK) 190, 170, 190
 170    IF (GXYIM(K)) 190, 180, 190
 180    XXK = 1.0
 190    PHI(K) = DEG*ATAN2(GXYIM(K),XXK)
 200  CONTINUE
C
C PLOT PHASE FROM -PHLIM TO PHLIM
C
      PHLIM = 1800.0
      DO 210 K=2,ND2P1
        X = PHI(K) - PHI(K-1)
        PHI(K) = PHI(K) - SIGN(360.,X)*AINT(0.5+ABS(X)/360.0)
        IF (PHI(K).GT.PHLIM) PHI(K) = PHI(K) - PHLIM
        IF (PHI(K).LT.(-PHLIM)) PHI(K) = PHI(K) + PHLIM
 210  CONTINUE
      WRITE (IPRTR,9984)
9984  FORMAT (1H1, 10X, 42HDUMP OF CONTINUOUS PHASE VALUES IN DEGREES/)
      WRITE (IPRTR,9983)
9983  FORMAT (3X, 9HFREQUENCY)
      DO 220 I=ISTRT,ISTOP
        TEMP = DF*FLOAT(I-1)
        WRITE (IPRTR,9982) TEMP, PHI(I)
9982    FORMAT (4X, F8.3, 2X, F10.2)
 220  CONTINUE
C
C COMPUTE CROSS SPECTRUM AND MAGNITUDE SQUARED COHERENCE
C
      DO 230 K=1,ND2P1
        PHI(K) = GXYRE(K)**2 + GXYIM(K)**2
        XX(K) = PHI(K)/(GXX(K)*GYY(K))
 230  CONTINUE
      WRITE (IPRTR,9981)
9981  FORMAT (1H1, 10X, 39HDUMP OF THE MAGNITUDE SQUARED COHERENCE/)
      WRITE (IPRTR,9983)
      DO 240 I=ISTRT,ISTOP
        TEMP = DF*FLOAT(I-1)
        WRITE (IPRTR,9982) TEMP, XX(I)
 240  CONTINUE
C
C CONVERT GXY TO DB AND PLOT
C
      DO 250 I=1,ND2P1
        PHI(I) = 5.0*ALOG10(AMAX1(PHI(I),SMALL))
 250  CONTINUE
      WRITE (IPRTR,9980)
9980  FORMAT (1H1/15X, 26HPLOT OF CROSS SPECTRUM GXY)
      CALL PRPRT(PHI, LINE, DF, FLOW, FHIGH, IPRTR)
C
C COMPUTE SIX GENERALIZED CROSS CORRELATION FUNCTIONS
C
      CALL PRCES(GXX, GYY, GXYRE, GXYIM, NPFFT, XX, YY, IPRTR, DT, LINE)
C
C COMPUTE MODULUS OF TRANSFER FUNCTION IN DB AND PLOT
C
      DO 260 K=1,ND2P1
        TEMP = (GXYRE(K)**2+GXYIM(K)**2)/GXX(K)**2
        PHI(K) = 10.0*ALOG10(AMAX1(TEMP,SMALL))
 260  CONTINUE
      WRITE (IPRTR,9979)
9979  FORMAT (1H1, 10X, 29HPLOT OF THE TRANSFER FUNCTION)
      CALL PRPRT(PHI, LINE, DF, FLOW, FHIGH, IPRTR)
C
C TERMINATE PROGRAM
C
      STOP
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  COPLT
C PLOT A CORRELOGRAM ON A 72 COLUMN PRINTER
C-----------------------------------------------------------------------
C
      SUBROUTINE COPLT(DATA, N, DT, ISWCH, IPRTR, LINE)
C
C INPUT:  DATA = ARRAY OF N VALUES WHICH CONTAIN THE CORRELOGRAM
C         N = NUMBER OF CORRELOGRAM POINTS TO PLOT
C         DT = TIME BETWEEN CORRELOGRAM POINTS
C         ISWCH = 1 FOR AUTO CORRELATION PLOT
C               = 2 FOR CROSS CORRELATION PLOT
C         IPRTR = LOGICAL UNIT NUMBER OF 72 COLUMN PRINTER
C         LINE = INTEGER SCRATCH ARRAY IN CALLING ROUTINE OF AT LEAST
C                45 WORDS
C
      DIMENSION DATA(1), LINE(1)
      DATA ISTAR /1H*/
C
C FIND PEAK AND MINIMUM VALUES OF ARRAY DATA
C
      FMIN = 10000.0
      PEAK = -10000.0
      DO 10 K=1,N
        PEAK = AMAX1(PEAK,DATA(K))
        FMIN = AMIN1(FMIN,DATA(K))
  10  CONTINUE
      WRITE (IPRTR,9999) FMIN, PEAK
9999  FORMAT (///5X, 6HFMIN =, E15.8, 8X, 6HPEAK =, E15.8///1X, 5HINDEX,
     *    4X, 4HTIME, 4X, 5HVALUE/)
C
C PLOT CORRELOGRAM ON A 72 COLUMN PRINTER
C ALL VALUES OF ARRAY DATA ARE SCALED TO FIT ON PRINTER
C
      DO 20 K=1,45
        LINE(K) = ISTAR
  20  CONTINUE
      ND2 = N/2
      ND2P1 = ND2 + 1
      DELTA = 45.0/(PEAK-FMIN)
      DO 50 K=1,N
        IF (ISWCH.EQ.2) GO TO 30
        TIME = DT*FLOAT(K-1)
        J = K
        GO TO 40
  30    TIME = DT*FLOAT(K-ND2)
        J = K - ND2
  40    INDEX = IFIX((DATA(K)-FMIN)*DELTA)
        IF (INDEX.LT.1) INDEX = 1
        IF (INDEX.GT.45) INDEX = 45
        WRITE (IPRTR,9998) J, TIME, DATA(K), (LINE(I),I=1,INDEX)
9998    FORMAT (1X, I5, F9.5, F9.4, 1X, 45A1)
  50  CONTINUE
C
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  HICMP
C THIS SUBROUTINE COMPUTES A NEW COMPOSITE NUMBER
C-----------------------------------------------------------------------
C
      SUBROUTINE HICMP(NNN, NEWNN)
C
C INPUT:  NNN = NUMBER OF DATA POINTS
C OUTPUT: NEWNN = A NEW COMPOSITE NUMBER ( A POWER OF 2 ) > OR = TO NNN
C
      DO 10 I=1,15
        M = I
        NT = 2**I
        IF (NNN.LE.NT) GO TO 20
  10  CONTINUE
C
  20  NEWNN = 2**M
C
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  LOAD
C THIS SUBROUTINE GENERATES TWO CHANNELS OF SYNTHETIC DATA FOR THE
C MAIN PROGRAM. THIS SUBROUTINE CAN BE REPLACED BY A DISC OR A
C MAGNETIC TAPE READ.
C-----------------------------------------------------------------------
C
      SUBROUTINE LOAD(XX, YY, NNN, KOUNT, ISR)
C
C INPUT:  NNN = NUMBER OF DATA POINTS TO BE GENERATED PER CHANNEL
C         KOUNT = NUMBER OF CURRENT SPECTRAL ESTIMATE
C         ISR = INTEGER SAMPLING RATE
C OUTPUT: XX = FIRST CHANNEL OF TIME DATA TO BE PROCESSED. THIS
C              BROADBAND SIGNAL IS GENERATED BY NON LINEARLY DISTORTING
C              A SIGNAL CONSISTING OF THE SUM OF FIVE SINUSOIDS. (THUS
C              GXX CONSISTS NOT ONLY OF THE FIVE SINE WAVES BUT MANY
C              INTERMODULATION PRODUCTS.
C         YY = SECOND CHANNEL OF TIME DATA TO BE PROCESSED. THIS
C              BROADBAND SIGNAL IS DETERMINISTICALLY RELATED TO THE XX
C              ARRAY WITH BOTH A LINEAR AND INCOHERENT COMPONENT
C              ADVANCED (DELAYED) BY ND8 UNITS.
C
      DIMENSION XX(1), YY(1)
      DIMENSION PHASE(5), FREQ(5)
C
      DT = 1.0/FLOAT(ISR)
      TPI = 8.0*ATAN(1.0)
      FREQ(1) = 10.0
      FREQ(2) = 27.0
      FREQ(3) = 43.9
      FREQ(4) = 71.8
      FREQ(5) = 108.31
      TPID = TPI/10.0
      DO 10 K=1,5
        PHASE(K) = FLOAT(K)*TPI*0.2
        FREQ(K) = TPI*FREQ(K)
  10  CONTINUE
C
      ND8 = NNN/8
      NLOOP = NNN + ND8
      DO 30 I=1,NLOOP
        TIME = FLOAT((KOUNT-1)*NNN+I)*DT
        SUM = 0.0
        DO 20 K=1,5
          SUM = SUM + SIN(FREQ(K)*TIME+PHASE(K))
  20    CONTINUE
        IF (SUM.GT.1.0) SUM = 1.0
        IF (SUM.LT.(-1.0)) SUM = -1.0
        IF (I.LE.NNN) XX(I) = SUM
        TEMP = SUM + 2.0*(SUM**2)
        J = I - ND8
        IF (I.GT.ND8) YY(J) = TEMP
  30  CONTINUE
C
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  LREMV
C THIS SUBROUTINE CAN REMOVE THE DC COMPONENT AND SLOPE OF AN ARRAY OF
C DATA IF DESIRED
C-----------------------------------------------------------------------
C
      SUBROUTINE LREMV(XX, NNN, ISWCH, DC, SLOPE)
C
C INPUT:  XX = INPUT DATA ARRAY
C         NNN = NUMBER OF POINTS IN DATA ARRAY
C         ISWCH = 0 DO NOT REMOVE DC COMPONENT OR SLOPE
C               = 1 REMOVE THE DC COMPONENT
C               > 1 REMOVE THE DC COMPONENT AND SLOPE
C OUTPUT: DC = DC COMPONENT OF DATA
C         SLOPE = SLOPE OF DATA
C
      DIMENSION XX(1)
C
C ESTABLISH CONSTANTS
C
      FLN = FLOAT(NNN)
      DC = 0.0
      SLOPE = 0.0
C
      DO 10 I=1,NNN
        DC = DC + XX(I)
        SLOPE = SLOPE + XX(I)*FLOAT(I)
  10  CONTINUE
C
C COMPUTE STATISTICS
C
      DC = DC/FLN
      SLOPE = 12.0*SLOPE/(FLN*(FLN*FLN-1.0)) - 6.0*DC/(FLN-1.0)
C
C DETERMINE KIND OF TREND REMOVAL
C
      IF (ISWCH-1) 60, 40, 20
C
C REMOVE TREND (MEAN AND SLOPE)
C
  20  CONTINUE
      FLN = DC - 0.5*(FLN+1.0)*SLOPE
      DO 30 I=1,NNN
        XX(I) = XX(I) - FLOAT(I)*SLOPE - FLN
  30  CONTINUE
      GO TO 60
C
C REMOVE THE DC COMPONENT
C
  40  CONTINUE
      DO 50 I=1,NNN
        XX(I) = XX(I) - DC
  50  CONTINUE
C
  60  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  PRCES
C THIS SUBROUTINE COMPUTES AND PLOTS SIX GENERALIZED CROSS CORRELATION
C FUNCTIONS
C-----------------------------------------------------------------------
C
      SUBROUTINE PRCES(GXX, GYY, GXYRE, GXYIM, NPFFT, XX, YY, IPRTR,
     *    DT, LINE)
C
C INPUT:  GXX = ARRAY OF AUTO SPECTRAL VALUES OF XX DATA CHANNEL
C         GYY = ARRAY OF AUTO SPECTRAL VALUES OF YY DATA CHANNEL
C         GXYRE = ARRAY REPRESENTING REAL PART OF CROSS SPECTRAL
C                 DENSITY FUNCTION
C         GXYIM = ARRAY REPRESENTING IMAGINARY PART OF CROSS SPECTRAL
C                 DENSITY FUNCTION
C         NPFFT = NUMBER REPRESENTING FOURIER TRANSFORM SIZE
C         XX = SCRATCH ARRAY OF LENGTH NPFFT
C         YY = SCRATCH ARRAY OF LENGTH NPFFT
C         IPRTR = LOGICAL UNIT NUMBER OF 80 COLUMN PRINTER
C         LINE = INTEGER SCRATCH ARRAY IN CALLING ROUTINE OF AT LEAST
C                50 WORDS
C
      DIMENSION GXX(1), GYY(1), GXYRE(1), GXYIM(1), XX(1), YY(1)
      DIMENSION LINE(1)
C
C CALCULATE CONSTANTS
C
      SMALL = R1MACH(1)
      SMALL = AMAX1(0.0001,SMALL)
      NP2 = NPFFT + 2
      ND2P1 = (NPFFT/2) + 1
      ND2 = NPFFT/2
      ND2M1 = ND2 - 1
C
C PROCESS SIX GENERALIZED CROSS CORRELATION FUNCTIONS
C
      DO 90 NTIME=1,6
C
        DO 10 K=1,ND2P1
C
          IF (NTIME.EQ.1) TEMP = 1.0/SQRT(GXX(K)*GYY(K))
          IF (NTIME.EQ.2) TEMP = 1.0/SQRT(GXYRE(K)**2+GXYIM(K)**2)
          IF (NTIME.EQ.3) TEMP = 1.0
          IF (NTIME.EQ.4) TEMP = 1.0/GXX(K)
          GXYMG = SQRT(GXYRE(K)**2+GXYIM(K)**2)
          COHR2 = GXYMG**2/(GXX(K)*GYY(K))
          COHR2 = AMIN1(COHR2,1.0-SMALL)
          IF (NTIME.EQ.5) TEMP = COHR2/((1.-COHR2)*GXYMG)
          TEMP1 = GXX(K) - GXYMG
          H = 1.0
          IF (ABS(TEMP1).LT.SMALL) TEMP1 = SMALL*SIGN(H,TEMP1)
          IF (NTIME.EQ.6) TEMP = GXYMG/(TEMP1**2)
C
          XX(K) = GXYRE(K)*TEMP
          YY(K) = GXYIM(K)*TEMP
C
  10    CONTINUE
C
        DO 20 K=2,ND2P1
          J = NP2 - K
          XX(J) = XX(K)
          YY(J) = -YY(K)
  20    CONTINUE
        YY(ND2P1) = 0.0
C
C COMPUTE INVERSE FFT
C
        CALL FFT842(1, NPFFT, XX, YY)
C
        TEMP = 0.0
        DO 30 K=1,NPFFT
          IF (TEMP.GE.ABS(XX(K))) GO TO 30
          KOFMX = K
          TEMP = ABS(XX(K))
  30    CONTINUE
C
        TEMP = 1.0/TEMP
        DO 40 K=1,NPFFT
          XX(K) = XX(K)*TEMP
  40    CONTINUE
C
        DO 50 I=1,ND2P1
          ITMP1 = I + ND2M1
          YY(ITMP1) = XX(I)
  50    CONTINUE
        DO 60 I=1,ND2M1
          ITMP1 = ND2P1 + I
          YY(I) = XX(ITMP1)
  60    CONTINUE
        NLAG = 100
        ITMP1 = ND2 - NLAG
        ITMP2 = ND2 + NLAG
        IF (ITMP1.GE.1) GO TO 70
        ITMP1 = 1
  70    IF (ITMP2.LE.NPFFT) GO TO 80
        ITMP2 = NPFFT
  80    CONTINUE
        XMIN = -DT*FLOAT(1-NLAG)
        XMAX = DT*FLOAT(1+NLAG)
C
C PLOT GENERALIZED CROSS CORRELATION FUNCTIONS
C
        IF (NTIME.EQ.1) WRITE (IPRTR,9999)
        IF (NTIME.EQ.2) WRITE (IPRTR,9998)
        IF (NTIME.EQ.3) WRITE (IPRTR,9997)
        IF (NTIME.EQ.4) WRITE (IPRTR,9996)
        IF (NTIME.EQ.5) WRITE (IPRTR,9995)
        IF (NTIME.EQ.6) WRITE (IPRTR,9994)
9999    FORMAT (1H1/20X, 21HPLOT OF SCOT FUNCTION)
9998    FORMAT (1H1/20X, 21HPLOT OF PHAT FUNCTION)
9997    FORMAT (1H1/20X, 25HPLOT OF CROSS CORRELATION)
9996    FORMAT (1H1/20X, 24HPLOT OF IMPULSE RESPONSE)
9995    FORMAT (1H1/20X, 23HPLOT OF H-T(I) FUNCTION)
9994    FORMAT (1H1/20X, 23HPLOT OF ECKART FUNCTION)
        CALL COPLT(YY, NPFFT, DT, 2, IPRTR, LINE)
C
  90  CONTINUE
C
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  PRPRT
C THIS SUBROUTINE PLOTS A POWER SPECTRUM ON A 72 COLUMN PRINTER FROM
C AN ARRAY OF DB VALUES
C-----------------------------------------------------------------------
C
      SUBROUTINE PRPRT(POWER, LINE, DF, FLOW, FHIGH, IPRTR)
C
C INPUT:  POWER = AN ARRAY OF POWER SPECTRAL VALUES IN DB
C         LINE = INTEGER SCRATCH ARRAY IN CALLING ROUTINE OF AT LEAST
C                50 WORDS WHICH IS USED TO STORE THE CHARACTER *
C         DF = FREQUENCY RESOLUTION IN HERTZ
C         FLOW = STARTING FREQUENCY OF SIGNAL TO BE PLOTTED
C                MINIMUM VALUE = 0.0 HERTZ
C         FHIGH = ENDING FREQUENCY OF SIGNAL TO BE PLOTTED
C                 MAXIMUM VALUE = FLOAT(ISR/2) HERTZ
C         IPRTR = LOGICAL UNIT NUMBER OF 72 COLUMN PRINTER
C
      DIMENSION POWER(1), LINE(1)
      DATA ISTAR /1H*/
C
C FIND PEAK AND MINIMUM DB VALUES OF ARRAY POWER BETWEEN FLOW AND FHIGH
C
      ISTRT = IFIX(FLOW/DF) + 1
      ISTOP = IFIX(FHIGH/DF) + 1
      FMIN = 10000.0
      PEAK = -10000.0
      DO 10 K=ISTRT,ISTOP
        PEAK = AMAX1(PEAK,POWER(K))
        FMIN = AMIN1(FMIN,POWER(K))
  10  CONTINUE
      WRITE (IPRTR,9999) FMIN, PEAK
9999  FORMAT (///5X, 6HFMIN =, F7.2, 3H DB, 4X, 6HPEAK =, F7.2, 3H DB//
     *    1X, 5HINDEX, 4X, 4HFREQ, 5X, 2HDB/)
C
C PLOT SPECTRUM ON PRINTER
C
      DO 20 K=1,50
        LINE(K) = ISTAR
  20  CONTINUE
C
      FBEG = FLOAT(IFIX(FLOW/DF))*DF
      DO 30 K=ISTRT,ISTOP
        FREQ = FBEG + DF*FLOAT(K-ISTRT)
        INDEX = IFIX(POWER(K)-FMIN)/2
        IF (INDEX.LT.1) INDEX = 1
        IF (INDEX.GT.50) INDEX = 50
        WRITE (IPRTR,9998) K, FREQ, POWER(K), (LINE(I),I=1,INDEX)
9998    FORMAT (I6, F8.3, F7.2, 1X, 50A1)
  30  CONTINUE
C
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  ZERO
C THIS SUBROUTINE STORES ZEROES IN A FLOATING POINT ARRAY
C-----------------------------------------------------------------------
C
      SUBROUTINE ZERO(ARRAY, NUMBR)
C
C INPUT:  ARRAY = AN ARRAY OF FLOATING POINT VALUES TO BE
C                 ZERO FILLED
C         NUMBR = NUMBER OF ARRAY VALUES
C
      DIMENSION ARRAY(1)
C
      DO 10 K=1,NUMBR
        ARRAY(K) = 0.0
  10  CONTINUE
C
      RETURN
      END
Return to Main Software Page