Make your own free website on Tripod.com
Periodogram Method for Power Spectrum Estimation

Return to Main Software Page

C
C-----------------------------------------------------------------------
C MAIN PROGRAM: MODIFIED PERIODOGRAM METHOD FOR POWER SPECTRUM
C               ESTIMATION-PMPSE
C AUTHORS:      L R RABINER
C               BELL LABORATORIES, MURRAY HILL, NEW JERSEY 07974
C               R W SCHAFER AND D DLUGOS
C               DEPT OF ELECTRICAL ENGINEERING
C               GEORGIA INSTITUTE OF TECHNOLOGY
C               ATLANTA, GEORGIA 30332
C
C METHOD BASED ON TECHNIQUE DESCRIBED BY P D WELCH, IEEE
C TRANS ON AUDIO AND ELECT, VOL 15, NO 2, PP 70-73, 1967.
C
C INPUT:        M IS THE FFT LENGTH (MUST BE A POWER OF 2)
C                   2 <= M <= MAXM (1024)
C               IWIN IS THE WINDOW TYPE
C                   IWIN = 1  RECTANGULAR WINDOW
C                   IWIN = 2  HAMMING WINDOW
C               L IS THE WINDOW LENGTH
C                   L <= M
C               N IS THE MAXIMUM NUMBER OF SAMPLES AVAILABLE
C                   FOR ANALYSIS
C               FS IS THE SAMPLING FREQUENCY IN HZ
C               IMD REQUESTS ADDITIONAL RUNS
C                   IMD = 1   NEW RUN
C                   IMD = 0   TERMINATES PROGRAM
C-----------------------------------------------------------------------
C
      DIMENSION XA(1024), XFR(513), SXX(513), WD(1024)
      DIMENSION JWIN(2,4)
      DIMENSION ILAG(513)
      COMPLEX X(1024), XMN
      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 IS USER-INTERACTIVE
C        THAT IS - A QUESTION IS WRITTEN ON THE USER
C        TERMINAL (IOUT1) AND THE USER TYPES IN THE ANSWER.
C
C OUTPUT: ALL OUTPUT IS WRITTEN ON THE STANDARD
C         OUTPUT UNIT (IOUT2).
C
      IND = I1MACH(1)
      IOUT1 = I1MACH(4)
      IOUT2 = I1MACH(2)
C
C SET MAXIMUM FFT SIZE MAXM DEPENDING ON DIMENSION WITHIN PROGRAM
C
      MAXM = 1024
      LHM = MAXM/2 + 1
C
C FILL LAG ARRAY FOR PRINTING
C
      DO 10 I=1,LHM
        ILAG(I) = I - 1
  10  CONTINUE
  20  CONTINUE
C
C READ IN ANALYSIS PARAMETERS M,IWIN,L,N,FS
C
      WRITE (IOUT1,9999)
9999  FORMAT (16H FFT LENGTH(I4)=)
      READ (IND,9997) M
      IF (M.GT.MAXM) WRITE (IOUT1,9998)
9998  FORMAT (27H M TOO LARGE--REENTER VALUE)
      IF (M.GT.MAXM) GO TO 20
9997  FORMAT (I4)
      WRITE (IOUT1,9996)
9996  FORMAT (43H WINDOW TYPE(I1)   1=RECTANGULAR, 2=HAMMING)
      READ (IND,9995) IWIN
9995  FORMAT (I1)
      WRITE (IOUT1,9994)
9994  FORMAT (19H WINDOW LENGTH(I4)=)
      READ (IND,9997) L
      WRITE (IOUT1,9993)
9993  FORMAT (40H MAXIMUM NUMBER OF ANALYSIS SAMPLES(I5)=)
      READ (IND,9992) N
9992  FORMAT (I5)
      WRITE (IOUT1,9991)
9991  FORMAT (33H SAMPLING FREQUENCY IN HZ(F10.4)=)
      READ (IND,9990) FS
9990  FORMAT (F10.4)
C
C NSECT = THE TOTAL NUMBER OF ANALYSIS SECTIONS
C    NP = THE TOTAL NUMBER OF SAMPLES ACTUALLY USED
C OVERLAP OF 2 TO 1 IS USED ON ADJACENT ANALYSIS SECTIONS
C    NP = N IF (N-L/2)/(L/2) = AN INTEGER
C
      MHLF1 = M/2 + 1
      NSECT = (N-L/2)/(L/2)
      NP = NSECT*(L/2) + L/2
      WRITE (IOUT2,9989) JWIN(IWIN,1), JWIN(IWIN,2), JWIN(IWIN,3),
     *    JWIN(IWIN,4)
9989  FORMAT (13H WINDOW TYPE=, 4A1)
      WRITE (IOUT2,9987)
      WRITE (IOUT2,9988) M, NP, L, FS
9988  FORMAT (3H M=, I4, 5H  NP=, I5, 4H  L=, I4, 18H  SAMPLING FREQUEN,
     *    3HCY=, F10.4)
C
C CALCULATE MEAN OF DATA.
C
      SS = 1.
      XSUM = 0.
      NS1 = NSECT + 1
      L1 = L/2
      DO 40 K=1,NS1
        CALL GETX(XA, L, SS)
        DO 30 I=1,L1
          XSUM = XSUM + XA(I)
  30    CONTINUE
        SS = SS + FLOAT(L1)
  40  CONTINUE
      XMEAN = XSUM/FLOAT(NP)
      XMN = CMPLX(XMEAN,XMEAN)
      WRITE (IOUT2,9987)
9987  FORMAT (//)
      WRITE (IOUT2,9986) XMEAN
9986  FORMAT (7H XMEAN=, E14.5)
C
C GENERATE WINDOW
C
      U = FLOAT(L)
      IF (IWIN.NE.2) GO TO 60
      U = 0.
      FL = FLOAT(L-1)
      TWOPI = 8.*ATAN(1.0)
      DO 50 I=1,L
        FI = FLOAT(I-1)
        WD(I) = .54 - .46*COS(TWOPI*FI/FL)
        U = U + WD(I)*WD(I)
  50  CONTINUE
  60  CONTINUE
C
C LOOP TO ACCUMULATE SPECTRA 2 AT A TIME
C
      SS = 1.
      DO 70 I=1,MHLF1
        SXX(I) = 0.
  70  CONTINUE
C
C READ L/2 SAMPLES TO INITIALIZE BUFFER
C
      NRD = L/2
      L2 = L/2 + 1
      CALL GETX(XA(L2), NRD, SS)
      SS = SS + FLOAT(NRD)
      IMN = L/2 + 1
      KMX = (NSECT+1)/2
      NSECTP = (NSECT+1)/2
      NRD = L
      DO 190 K=1,KMX
C
C MOVE DOWN UPPER HALF OF XA BUFFER
C
        DO 80 I=1,L1
          J = L1 + I
          X(I) = CMPLX(XA(J),0.)
  80    CONTINUE
        IF (K.NE.KMX .OR. NSECTP.EQ.NSECT) GO TO 100
        DO 90 I=IMN,NRD
          XA(I) = 0.0
  90    CONTINUE
        NRD = L/2
 100    CALL GETX(XA, NRD, SS)
        DO 110 I=1,L1
          J = I + L1
          X(J) = CMPLX(XA(I),XA(J)) - XMN
          X(I) = CMPLX(REAL(X(I)),XA(I)) - XMN
 110    CONTINUE
        IF (K.NE.KMX .OR. NSECTP.EQ.NSECT) GO TO 130
C
C AN ODD NUMBER OF SECTIONS--ZERO OUT THE SECOND PART
C
        DO 120 I=1,L
          X(I) = CMPLX(REAL(X(I)),0.)
 120    CONTINUE
 130    CONTINUE
        SS = SS + FLOAT(NRD)
        IF (IWIN.NE.2) GO TO 150
        DO 140 I=1,L
          X(I) = X(I)*WD(I)
 140    CONTINUE
 150    CONTINUE
        IF (L.EQ.M) GO TO 170
        LP1 = L + 1
        DO 160 I=LP1,M
          X(I) = (0.,0.)
 160    CONTINUE
 170    CONTINUE
        CALL FFT(X, M, 0)
        DO 180 I=2,MHLF1
          J = M + 2 - I
          SXX(I) = SXX(I) + REAL(X(I)*CONJG(X(I))+X(J)*CONJG(X(J)))
 180    CONTINUE
        SXX(1) = SXX(1) + REAL(X(1)*CONJG(X(1)))*2.
 190  CONTINUE
C
C NORMALIZE SPECTRAL ESTIMATE AND OBTAIN CORRELATION FUNCTION
C USING INVERSE FFT OF POWER SPECTRUM
C
      FNORM = 2.*U*FLOAT(NSECT)
      SXX(1) = SXX(1)/FNORM
      X(1) = CMPLX(SXX(1),0.)
      DO 200 I=2,MHLF1
        SXX(I) = SXX(I)/FNORM
        X(I) = CMPLX(SXX(I),0.)
        J = M + 2 - I
        X(J) = X(I)
 200  CONTINUE
      CALL FFT(X, M, 1)
      DO 210 I=1,MHLF1
        XA(I) = REAL(X(I))
 210  CONTINUE
C
C CORRELATION ESTIMATE IS IN XA FROM 1 TO MHLF1
C COMPUTE LOG OF POWER SPECTRUM ESTIMATE
C
      XFS = FS/FLOAT(M)
      DO 220 I=1,MHLF1
        XFR(I) = FLOAT(I-1)*XFS
        TMP = ALOG10(SXX(I))
        SXX(I) = 20.*TMP
 220  CONTINUE
C
C LOG POWER SPECTRUM (DB) IS IN ARRAY SXX
C IF DESIRED, THE USER MAY INSERT CODE AT THIS POINT TO PLOT THE LOG
C POWER SPECTRUM
C
      WRITE (IOUT2,9987)
      WRITE (IOUT2,9985)
9985  FORMAT (19H LOG POWER SPECTRUM)
      WRITE (IOUT2,9987)
      WRITE (IOUT2,9984)
9984  FORMAT (5X, 4HFREQ, 7X, 2HDB, 5X, 4HFREQ, 7X, 2HDB, 5X, 4HFREQ,
     *    7X, 2HDB, 5X, 4HFREQ, 7X, 2HDB)
      WRITE (IOUT2,9983) (XFR(I),SXX(I),I=1,MHLF1)
9983  FORMAT (4(1X, F8.2, 1X, F8.3))
C
C CORRELATION FUNCTION IS IN ARRAY XA
C IF DESIRED, THE USER MAY INSERT CODE AT THIS POINT TO PLOT THE
C CORRELATION FUNCTION
C
      WRITE (IOUT2,9987)
      WRITE (IOUT2,9982)
9982  FORMAT (21H CORRELATION FUNCTION)
      WRITE (IOUT2,9987)
      WRITE (IOUT2,9981)
9981  FORMAT (1X, 3HLAG, 2X, 4HCORR, 5X, 3HLAG, 2X, 4HCORR, 5X, 3HLAG,
     *    2X, 4HCORR, 5X, 3HLAG, 2X, 4HCORR, 5X, 3HLAG, 2X, 4HCORR)
      WRITE (IOUT2,9980) (ILAG(I),XA(I),I=1,MHLF1)
9980  FORMAT (5(1X, I3, E10.3))
      WRITE (IOUT2,9987)
      WRITE (IOUT2,9979)
9979  FORMAT (////)
      WRITE (IOUT1,9978)
9978  FORMAT (23H MORE DATA(1=YES,0=NO)=)
      READ (IND,9995) IMD
      IF (IMD.EQ.1) GO TO 20
      STOP
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  GETX
C READ DATA ROUTINE GENERATE X(N) FOR A SINE INPUT OF FREQUENCY FS/10.
C WHERE FS IS THE SAMPLING FREQUENCY IN 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 SINE WAVE FREQUENCY IS 1000 HZ WITH AN ASSUMED SAMPLING FREQUENCY
C OF 10000 HZ
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:  FFT
C JIM COOLEY'S SIMPLE FFT PROGRAM--USES DECIMATION IN TIME ALGORITHM
C X IS AN N=2**M POINT COMPLEX ARRAY THAT INITIALLY CONTAINS THE INPUT
C AND ON OUTPUT CONTAINS THE TRANSFORM
C THE PARAMETER INV SPECIFIED DIRECT TRANSFORM IF 0 AND INVERSE IF 1
C-----------------------------------------------------------------------
C
      SUBROUTINE FFT(X, N, INV)
      COMPLEX X(1), U, W, T, CMPLX
C
C   X = COMPLEX ARRAY OF SIZE N--ON INPUT X CONTAINS
C       THE SEQUENCE TO BE TRANSFORMED
C       ON OUTPUT X CONTAINS THE DFT OF THE INPUT
C   N = SIZE OF FFT TO BE COMPUTED--N=2**M FOR 1.LE.M.LE.15
C INV = PARAMETER TO DETERMINE WHETHER TO DO A DIRECT TRANSFORM (INV=0)
C       OR AN INVERSE TRANSFORM (INV=1)
C
      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.*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