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