Correlation Method for Power Spectrum Estimation

Return to Main Software Page

C
C-----------------------------------------------------------------------
C MAIN PROGRAM: CORRELATION METHOD FOR POWER SPECTRUM ESTIMATION - CMPSE
C AUTHORS:      L R RABINER
C               BELL LABORATORIES, MURRAY HILL, NEW JERSEY 07974
C               R W SCHAFER AND D DLUGOS
C               GEORGIA INSTITUTE OF TECHNOLOGY, ATLANTA, GEORGIA 30332
C
C THIS METHOD IS BASED ON THE TECHNIQUE DESCRIBED BY C M RADER, IN THE
C IEEE TRANS ON AUDIO AND ELECT, VOL 18, NO 4, PP 439-442, 1970.
C
C INPUT:        M IS THE SECTION SIZE(MUST BE A POWER OF 2)
C                   2 <= M <= 512
C               N IS THE NUMBER OF SAMPLES TO BE USED IS THE ANALYSIS
C               MODE IS THE DATA FORMAT TYPE
C                   MODE = 0   AUTO CORRELATION
C                   MODE = 1   CROSS CORRELATION
C                   MODE = 2   AUTO COVARIANCE
C                   MODE = 3   CROSS COVARIANCE
C               FS IS THE SAMPLING FREQUENCY IN HZ
C               IWIN IS THE WINDOW TYPE
C                   IWIN = 1   RECTANGULAR WINDOW
C                   IWIN = 2   HAMMING WINDOW
C               L IS THE NUMBER OF CORRELATION VALUES USED IN
C                   THE SPECTRAL ESTIMATE
C                   2 <= L <= M
C               NFFT IS THE SIZE FFT USED TO GIVE THE SPECTRAL ESTIMATE
C                   L <= NFFT <= MAXM
C               IMD REQUESTS ADDITIONAL RUNS
C                   IMD = 1   NEW RUN
C                   IMD = 0   TERMINATE PROGRAM
C-----------------------------------------------------------------------
C
      DIMENSION XA(512), XFR(257)
      DIMENSION JWIN(2,4)
      DIMENSION ILAG(257)
      COMPLEX X(512), Z(257), XMN, XI, YI
      INTEGER TTI, TTO
      DATA JWIN(1,1), JWIN(1,2), JWIN(1,3), JWIN(1,4) /1HR,1HE,1HC,1HT/
      DATA JWIN(2,1), JWIN(2,2), JWIN(2,3), JWIN(2,4) /1HH,1HA,1HM,1HG/
C
C DEFINE I/O DEVICE CODES
C INPUT: INPUT TO THIS PROGRAM IS USER-INTERACTIVE
C        THAT IS - A QUESTION IS WRITTEN ON THE USER
C        TERMINAL (TTO) AND THE USER TYPES IN THE ANSWER.
C
C OUTPUT: ALL OUTPUT IS WRITTEN ON THE STANDARD
C         OUTPUT UNIT (LPT)
C
      TTI = I1MACH(1)
      TTO = I1MACH(4)
      LPT = I1MACH(2)
      MAXM = 512
      MAXH = MAXM/2 + 1
C
C FILL LAG ARRAY FOR PRINTING
C
      DO 10 I=1,MAXH
        ILAG(I) = I - 1
  10  CONTINUE
C
C READ IN ANALYSIS PARAMETERS M,N
C
  20  WRITE (TTO,9999)
9999  FORMAT (18H SECTION SIZE(I4)=)
      READ (TTI,9998) M
9998  FORMAT (I4)
      IF (M.GT.0 .AND. M.LE.MAXM) GO TO 30
      WRITE (TTO,9997) M
9997  FORMAT (31H ILLEGAL INPUT -- REENTER VALUE)
      GO TO 20
  30  WRITE (TTO,9996)
9996  FORMAT (38H TOTAL NUMBER OF ANALYSIS SAMPLES(I5)=)
      READ (TTI,9995) N
9995  FORMAT (I5)
C
C NSECT IS THE TOTAL NUMBER OF ANALYSIS SECTIONS
C LSHFT IS THE SHIFT BETWEEN ADJACENT ANALYSIS SAMPLES
C
      LSHFT = M/2
      MHLF1 = LSHFT + 1
      NSECT = (FLOAT(N)+FLOAT(LSHFT)-1.)/FLOAT(LSHFT)
C
C READ IN MODE DATA TYPE FORMAT
C
      WRITE (TTO,9994)
9994  FORMAT (10H MODE(I1)=)
      READ (TTI,9993) MODE
9993  FORMAT (I1)
      WRITE (TTO,9992)
9992  FORMAT (33H SAMPLING FREQUENCY IN HZ(F10.4)=)
      READ (TTI,9991) FS
9991  FORMAT (F10.4)
      WRITE (LPT,9990) M, N, MODE, FS
9990  FORMAT (3H M=, I4, 4H  N=, I5, 7H  MODE=, I1, 16H  SAMPLING FREQU,
     *    5HENCY=, F10.4)
      IF (MODE.LT.2) GO TO 80
C
C SS IS GENERATOR SAMPLE NUMBER
C NRD IS NUMBER OF SAMPLES OF GENERATOR OUTPUT TO BE COMPUTED
C
      SS = 1.
      NRD = LSHFT
      XSUM = 0.
      YSUM = 0.
C
C LOOP TO CALCULATE MEANS OF X AND Y DATA
C USE GETX TO READ NRD SAMPLES FROM X GENERATOR STARTING AT SAMPLE SS
C USE GETY TO READ NRD SAMPLES FROM Y GENERATOR IF CROSS VARIANCE
C
      DO 70 K=1,NSECT
        IF (K.EQ.NSECT) NRD = N - (K-1)*NRD
        CALL GETX(XA, NRD, SS)
        DO 40 I=1,NRD
          XSUM = XSUM + XA(I)
  40    CONTINUE
        IF (MODE.EQ.2) GO TO 60
        CALL GETY(XA, NRD, SS)
        DO 50 I=1,NRD
          YSUM = YSUM + XA(I)
  50    CONTINUE
  60    SS = SS + FLOAT(NRD)
  70  CONTINUE
      XMEAN = XSUM/FLOAT(N)
      YMEAN = YSUM/FLOAT(N)
      IF (MODE.EQ.2) YMEAN = XMEAN
      XMN = CMPLX(XMEAN,YMEAN)
      WRITE (LPT,9989)
9989  FORMAT (//)
      WRITE (LPT,9988) XMEAN, YMEAN
9988  FORMAT (7H XMEAN=, E14.5, 8H  YMEAN=, E14.5)
C
C LOOP TO ACCUMULATE CORRELATIONS
C
  80  SS = 1.
      NRDY = M
      NRDX = LSHFT
      DO 90 I=1,MHLF1
        Z(I) = (0.,0.)
  90  CONTINUE
      DO 190 K=1,NSECT
        NSECT1 = NSECT - 1
        IF (K.LT.NSECT1) GO TO 110
        NRDY = N - (K-1)*LSHFT
        IF (K.EQ.NSECT) NRDX = NRDY
        IF (NRDY.EQ.M) GO TO 110
        NRDY1 = NRDY + 1
        DO 100 I=NRDY1,M
          X(I) = (0.,0.)
 100    CONTINUE
C
C READ NRDY SAMPLES FROM X GENERATOR STARTING AT SAMPLE SS
C
 110    CALL GETX(XA, NRDY, SS)
        DO 120 I=1,NRDY
          X(I) = CMPLX(XA(I),XA(I))
 120    CONTINUE
        IF (MODE.EQ.0 .OR. MODE.EQ.2) GO TO 140
C
C READ NRDY SAMPLES FROM Y GENERATOR IF CROSS COR. OR CROSS COV.
C
        CALL GETY(XA, NRDY, SS)
        DO 130 I=1,NRDY
          X(I) = CMPLX(REAL(X(I)),XA(I))
 130    CONTINUE
 140    IF (MODE.LT.2) GO TO 160
        DO 150 I=1,NRDY
          X(I) = X(I) - XMN
 150    CONTINUE
 160    NRDX1 = NRDX + 1
        DO 170 I=NRDX1,M
          X(I) = CMPLX(0.,AIMAG(X(I)))
 170    CONTINUE
C
C CORRELATE X AND Y SECTIONS
C DO EVEN-ODD SEPARATION AND ACCUMULATE  CONJG(X)*Y
C
        CALL FFT(X, M, 0)
        DO 180 I=2,LSHFT
          J = M + 2 - I
          XI = (X(I)+CONJG(X(J)))*.5
          YI = (X(J)-CONJG(X(I)))*.5
          YI = CMPLX(AIMAG(YI),REAL(YI))
          Z(I) = Z(I) + CONJG(XI)*YI
 180    CONTINUE
        XI = X(1)
        Z(1) = Z(1) + CMPLX(REAL(XI)*AIMAG(XI),0.)
        XI = X(MHLF1)
        Z(MHLF1) = Z(MHLF1) + CMPLX(REAL(XI)*AIMAG(XI),0.)
        SS = SS + FLOAT(LSHFT)
 190  CONTINUE
C
C INVERSE DFT TO GIVE CORRELATION
C
      DO 200 I=2,LSHFT
        J = M + 2 - I
        X(I) = Z(I)
        X(J) = CONJG(Z(I))
 200  CONTINUE
      X(1) = Z(1)
      X(MHLF1) = Z(MHLF1)
      CALL FFT(X, M, 1)
      FN = FLOAT(N)
      DO 210 I=1,MHLF1
        XA(I) = REAL(X(I))/FN
 210  CONTINUE
C
C IF DESIRED, THE USER MAY INSERT CODE AT THIS POINT TO PLOT
C THE CORRELATION FUNCTION WHICH IS IN THE ARRAY XA.
C
C PRINT THE CORRELATION FUNCTION
C
      WRITE (LPT,9989)
      WRITE (LPT,9987)
9987  FORMAT (21H CORRELATION FUNCTION)
      WRITE (LPT,9989)
      WRITE (LPT,9986)
9986  FORMAT (1X, 3HLAG, 2X, 4HCORR, 5X, 3HLAG, 2X, 4HCORR, 5X, 3HLAG,
     *    2X, 4HCORR, 5X, 3HLAG, 2X, 4HCORR, 5X, 3HLAG, 2X, 4HCORR)
      WRITE (LPT,9985) (ILAG(I),XA(I),I=1,MHLF1)
9985  FORMAT (5(1X, I3, E10.3))
      WRITE (LPT,9989)
C
C WINDOW CORRELATION USING L VALUES TO GIVE SPECTRAL ESTIMATE
C CREATE SYMMETRICAL ARRAY IF MODE=0
C READ IN WINDOW TYPE AND WINDOW LENGTH
C NOTE SPECTRAL ESTIMATE MAY NOT BE MEANINGFUL IF X NOT EQUAL TO Y
C
      WRITE (TTO,9984)
9984  FORMAT (43H WINDOW TYPE(I1)-  1=RECTANGULAR, 2=HAMMING)
      READ (TTI,9983) IWIN
9983  FORMAT (I1)
      WRITE (TTO,9982)
9982  FORMAT (35H NO OF CORRELATION VALUES USED(I4)=)
      READ (TTI,9981) L
9981  FORMAT (I4)
      WRITE (TTO,9980)
9980  FORMAT (14H FFT SIZE(I4)=)
      READ (TTI,9981) NFFT
      NHLF1 = NFFT/2 + 1
      WRITE (LPT,9979) JWIN(IWIN,1), JWIN(IWIN,2), JWIN(IWIN,3),
     *    JWIN(IWIN,4), L, NFFT
9979  FORMAT (13H WINDOW TYPE=, 4A1, 22H  NO OF WINDOW VALUES=, I4,
     *    11H  FFT SIZE=, I4)
C
C WINDOW CORRELATION FUNCTION--BEWARE IF X NOT EQUAL TO Y
C
      PI = 4.0*ATAN(1.0)
      DO 230 I=2,L
        IF (IWIN.EQ.1) GO TO 220
        XA(I) = XA(I)*(0.54+0.46*COS(PI*FLOAT(I-1)/FLOAT(L-1)))
 220    IF (MODE.EQ.1 .OR. MODE.EQ.3) GO TO 230
        J = NFFT + 2 - I
        XA(J) = XA(I)
 230  CONTINUE
      NLAST = NFFT + 1 - L
      IF (MODE.EQ.1 .OR. MODE.EQ.3) NLAST = NFFT
      L1 = L + 1
      DO 240 I=L1,NLAST
        XA(I) = 0.
 240  CONTINUE
      DO 250 I=1,NFFT
        X(I) = CMPLX(XA(I),0.)
 250  CONTINUE
      CALL FFT(X, NFFT, 0)
C
C OBTAIN LOG POWER SPECTRUM  IN DB
C
      XFS = FS/FLOAT(NFFT)
      NHF = NFFT/2
      NHF1 = NHF + 1
      DO 260 I=1,NHF1
        XFR(I) = FLOAT(I-1)*XFS
        T = ALOG10(CABS(X(I)))
        XA(I) = 20.*T
 260  CONTINUE
C
C LOG POWER SPECTRUM (DB) IS IN XA
C IF DESIRED, THE USER MAY INSERT CODE AT THIS POINT TO
C PLOT LOG POWER SPECTRUM AT THIS POINT
C
      WRITE (LPT,9989)
      WRITE (LPT,9978)
9978  FORMAT (19H LOG POWER SPECTRUM)
      WRITE (LPT,9989)
      WRITE (LPT,9977)
9977  FORMAT (5X, 4HFREQ, 7X, 2HDB, 5X, 4HFREQ, 7X, 2HDB, 5X, 4HFREQ,
     *    7X, 2HDB, 5X, 4HFREQ, 7X, 2HDB)
      WRITE (LPT,9976) (XFR(I),XA(I),I=1,NHLF1)
9976  FORMAT (4(1X, F8.2, 1X, F8.3))
      WRITE (LPT,9975)
9975  FORMAT (////)
      WRITE (TTO,9974)
9974  FORMAT (23H MORE DATA(1=YES,0=NO)=)
      READ (TTI,9993) IMD
      IF (IMD.EQ.1) GO TO 20
      STOP
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: GETX
C GENERATE X(N) FOR A SINE INPUT OF FREQUENCY 1000 HZ WITH AN ASSUME
C SAMPLING FREQUENCY OF 10000 HZ
C-----------------------------------------------------------------------
C
      SUBROUTINE GETX(X, NRD, SS)
      DIMENSION X(1)
C
C   X = ARRAY OF SIZE NRD TO HOLD GENERATOR OUTPUT DATA
C NRD = NUMBER OF SAMPLES TO BE CREATED
C  SS = STARTING SAMPLE OF GENERATOR OUTPUT
C
      TPI = 8.*ATAN(1.0)
      CF = 1000./10000.
      DO 10 I=1,NRD
        XSMP = (SS-1.) + FLOAT(I-1)
        X(I) = COS(TPI*CF*XSMP)
  10  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: GETY
C GENERATE Y(N) FOR A SINE INPUT OF FREQUENCY 1000 HZ WITH AN
C ASSUMED SAMPLING FREQUENCY OF 10000 HZ
C-----------------------------------------------------------------------
C
      SUBROUTINE GETY(Y, NRD, SS)
      DIMENSION Y(1)
C
C   Y = ARRAY OF SIZE NRD TO HOLD GENERATOR OUTPUT DATA
C NRD = NUMBER OF SAMPLES TO BE CREATED
C  SS = STARTING SAMPLE OF GENERATOR OUTPUT
C
      TPI = 8.*ATAN(1.0)
      CF = 1000./10000.
      DO 10 I=1,NRD
        XSMP = (SS-1.) + FLOAT(I-1)
        Y(I) = SIN(TPI*CF*XSMP)
  10  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: FFT
C JIM COOLEY'S SIMPLE FFT PROGRAM USING DECIMATION IN TIME ALGORITHM
C-----------------------------------------------------------------------
C
      SUBROUTINE FFT(X, N, INV)
C
C   X = 2**M COMPLEX ARRAY THAT INITIALLY CONTAINS INPUT
C       AND ON RETURN CONTAINS TRANSFORM
C   N = 2**M POINTS
C INV = 0, DIRECT TRANSFORM
C INV = 1, INVERSE TRANSFORM
C
      COMPLEX X(1), U, W, T, CMPLX
      M = ALOG(FLOAT(N))/ALOG(2.) + .1
      NV2 = N/2
      NM1 = N - 1
      J = 1
      DO 40 I=1,NM1
        IF (I.GE.J) GO TO 10
        T = X(J)
        X(J) = X(I)
        X(I) = T
  10    K = NV2
  20    IF (K.GE.J) GO TO 30
        J = J - K
        K = K/2
        GO TO 20
  30    J = J + K
  40  CONTINUE
      PI = 4.0*ATAN(1.0)
      DO 70 L=1,M
        LE = 2**L
        LE1 = LE/2
        U = (1.0,0.0)
        W = CMPLX(COS(PI/FLOAT(LE1)),-SIN(PI/FLOAT(LE1)))
        IF (INV.NE.0) W = CONJG(W)
        DO 60 J=1,LE1
          DO 50 I=J,N,LE
            IP = I + LE1
            T = X(IP)*U
            X(IP) = X(I) - T
            X(I) = X(I) + T
  50      CONTINUE
          U = U*W
  60    CONTINUE
  70  CONTINUE
      IF (INV.EQ.0) RETURN
      DO 80 I=1,N
        X(I) = X(I)/CMPLX(FLOAT(N),0.)
  80  CONTINUE
      RETURN
      END
Return to Main Software Page