Linear Prediction Analysis Programs

C
C-----------------------------------------------------------------------
C  MAIN PROGRAM: TEST PROGRAM FOR COVAR AND AUTO
C  AUTHORS:      A. H. GRAY, JR.  AND J. D. MARKEL
C      GRAY -    UNIV. OF CALIF., SANTA BARBARA, CA 93109
C      BOTH -    SIGNAL TECHNOLOGY, INC., 15 W. DE LA GUERRA STREET
C                SANTA BARBARA, CA 93101
C-----------------------------------------------------------------------
C
DIMENSION X(100), A(21), RC(21)
IOUTD = I1MACH(2)
X(1) = 1.
X(2) = 1.44
X(3) = 1.44*X(2) - 1.26
DO 10 J=4,40
X(J) = 1.44*X(J-1) - 1.26*X(J-2) + .81*X(J-3)
10  CONTINUE
9999  FORMAT (1X, I3, 2E15.6)
9998  FORMAT (1X, E12.6/)
N = 40
DO 20 M=2,3
CALL COVAR(N, X, M, A, ALPHA, RC)
MP = M + 1
WRITE (IOUTD,9999) (I,A(I),RC(I),I=1,MP)
WRITE (IOUTD,9998) ALPHA
CALL AUTO(N, X, M, A, ALPHA, RC)
WRITE (IOUTD,9999) (I,A(I),RC(I),I=1,MP)
WRITE (IOUTD,9998) ALPHA
20  CONTINUE
DO 30 JB=1,40
J = 44 - JB
X(J) = X(J-3)
30  CONTINUE
DO 40 J=1,3
JB = 43 + J
X(J) = 0.
X(JB) = 0.
40  CONTINUE
M = 3
N = 46
CALL COVAR(N, X, M, A, ALPHA, RC)
MP = M + 1
WRITE (IOUTD,9999) (I,A(I),RC(I),I=1,MP)
WRITE (IOUTD,9998) ALPHA
CALL AUTO(N, X, M, A, ALPHA, RC)
WRITE (IOUTD,9999) (I,A(I),RC(I),I=1,MP)
WRITE (IOUTD,9998) ALPHA
STOP
END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  COVAR
C A SUBROUTINE FOR IMPLEMENTING THE COVARIANCE
C METHOD OF LINEAR PREDICTION ANALYSIS
C-----------------------------------------------------------------------
C
SUBROUTINE COVAR(N, X, M, A, ALPHA, GRC)
C
C   INPUTS:   N - NO. OF DATA POINTS
C             X(N) - INPUT DATA SEQUENCE
C             M - ORDER OF FILTER (M<21, SEE NOTE*)
C  OUTPUTS:   A - FILTER COEFFICIENTS
C             ALPHA - RESIDUAL "ENERGY"
C             A - FILTER COEFFICIENTS
C             GRC - "GENERALIZED REFLECTION COEFFICIENTS",
C
C *PROGRAM LIMITED TO M=20, BECAUSE OF THE DIMENSIONS
C B(M*(M-1)/2), BETA(M), AND CC(M+1)
C
DIMENSION X(1), A(1), GRC(1)
DIMENSION B(190), BETA(20), CC(21)
MP = M + 1
MT = (MP*M)/2
DO 10 J=1,MT
B(J) = 0.
10  CONTINUE
ALPHA = 0.
CC(1) = 0.
CC(2) = 0.
DO 20 NP=MP,N
NP1 = NP - 1
ALPHA = ALPHA + X(NP)*X(NP)
CC(1) = CC(1) + X(NP)*X(NP1)
CC(2) = CC(2) + X(NP1)*X(NP1)
20  CONTINUE
B(1) = 1.
BETA(1) = CC(2)
GRC(1) = -CC(1)/CC(2)
A(1) = 1.
A(2) = GRC(1)
ALPHA = ALPHA + GRC(1)*CC(1)
MF = M
DO 130 MINC=2,MF
DO 30 J=1,MINC
JP = MINC + 2 - J
N1 = MP + 1 - JP
N2 = N + 1 - MINC
N3 = N + 2 - JP
N4 = MP - MINC
CC(JP) = CC(JP-1) + X(N4)*X(N1) - X(N2)*X(N3)
30    CONTINUE
CC(1) = 0.
DO 40 NP=MP,N
N1 = NP - MINC
CC(1) = CC(1) + X(N1)*X(NP)
40    CONTINUE
MSUB = (MINC*MINC-MINC)/2
MM1 = MINC - 1
N1 = MSUB + MINC
B(N1) = 1.
DO 80 IP=1,MM1
ISUB = (IP*IP-IP)/2
IF (BETA(IP)) 150, 150, 50
50      GAM = 0.
DO 60 J=1,IP
N1 = ISUB + J
GAM = GAM + CC(J+1)*B(N1)
60      CONTINUE
GAM = GAM/BETA(IP)
DO 70 JP=1,IP
N1 = MSUB + JP
N2 = ISUB + JP
B(N1) = B(N1) - GAM*B(N2)
70      CONTINUE
80    CONTINUE
BETA(MINC) = 0.
DO 90 J=1,MINC
N1 = MSUB + J
BETA(MINC) = BETA(MINC) + CC(J+1)*B(N1)
90    CONTINUE
IF (BETA(MINC)) 150, 150, 100
100    S = 0.
DO 110 IP=1,MINC
S = S + CC(IP)*A(IP)
110    CONTINUE
GRC(MINC) = -S/BETA(MINC)
DO 120 IP=2,MINC
M2 = MSUB + IP - 1
A(IP) = A(IP) + GRC(MINC)*B(M2)
120    CONTINUE
A(MINC+1) = GRC(MINC)
S = GRC(MINC)*GRC(MINC)*BETA(MINC)
ALPHA = ALPHA - S
IF (ALPHA) 150, 150, 130
130  CONTINUE
140  RETURN
150  CONTINUE
C
C WARNING - SINGULAR MATRIX
C
IOUTD = I1MACH(2)
WRITE (IOUTD,9999)
9999  FORMAT (34H WARNING - SINGULAR MATRIX - COVAR)
GO TO 140
END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  AUTO
C A SUBROUTINE FOR IMPLEMENTING THE AUTOCORRELATION
C METHOD OF LINEAR PREDICTION ANALYSIS
C-----------------------------------------------------------------------
C
SUBROUTINE AUTO(N, X, M, A, ALPHA, RC)
C
C   INPUTS:   N - NO. OF DATA POINTS
C             X(N) - INPUT DATA SEQUENCE
C             M - ORDER OF FILTER (M<21, SEE NOTE*)
C  OUTPUTS:   A - FILTER COEFFICIENTS
C             ALPHA - RESIDUAL "ENERGY"
C             RC - REFLECTION COEFFICIENTS
C
C *PROGRAM LIMITED TO M<21 BECAUSE OF DIMENSIONS OF R(.)
C
DIMENSION X(1), A(1), RC(1)
DIMENSION R(21)
MP = M + 1
DO 20 K=1,MP
R(K) = 0.
NK = N - K + 1
DO 10 NP=1,NK
N1 = NP + K - 1
R(K) = R(K) + X(NP)*X(N1)
10    CONTINUE
20  CONTINUE
RC(1) = -R(2)/R(1)
A(1) = 1.
A(2) = RC(1)
ALPHA = R(1) + R(2)*RC(1)
DO 50 MINC=2,M
S = 0.
DO 30 IP=1,MINC
N1 = MINC - IP + 2
S = S + R(N1)*A(IP)
30    CONTINUE
RC(MINC) = -S/ALPHA
MH = MINC/2 + 1
DO 40 IP=2,MH
IB = MINC - IP + 2
AT = A(IP) + RC(MINC)*A(IB)
A(IB) = A(IB) + RC(MINC)*A(IP)
A(IP) = AT
40    CONTINUE
A(MINC+1) = RC(MINC)
ALPHA = ALPHA + RC(MINC)*S
IF (ALPHA) 70, 70, 50
50  CONTINUE
60  RETURN
70  CONTINUE
C
C WARNING - SINGULAR MATRIX
C
IOUTD = I1MACH(2)
WRITE (IOUTD,9999)
9999  FORMAT (33H WARNING - SINGULAR MATRIX - AUTO)
GO TO 60
END
C