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