Program for Minimum-p Synthesis of Recursive Digital Filters

Return to Main Software Page

C
C-----------------------------------------------------------------------
C  MAIN PROGRAM: COMBINED AMPLITUDE AND GROUP DELAY SYNTHESIS USING A
C                MINIMUM P ERROR CRITERION AND THE FLETCHER POWELL
C                SUBROUTINE FOR FUNCTION MINIMIZATION
C
C  AUTHOR:       A. G. DECZKY
C                INSTITUT DE MICROTECHNIQUE, UNIVERSITE DE NEUCHATEL
C                RUE DE LA MALADIERE 71,  2000 NEUCHATEL, SWITZERLAND
C
C  REFERENCE:    'SYNTHESIS OF RECURSIVE DIGITAL FILTERS USING THE
C                MINIMUM P ERROR CRITERION', A. G. DECZKY, IEEE TRANS.
C                AUDIO ELECTROACOUST., VOL. AU-20, NO. 5, PP. 257-263,
C                OCTOBER 1972.
C-----------------------------------------------------------------------
C
      LOGICAL DEB, DEBFP, LONG
      INTEGER WEIGHT
      COMMON /CM05/ AS(201), DS(201), WF(201), WG(201)
      COMMON /CM10/ NP, NPD
      COMMON /CM15/ DEB
      COMMON /CM20/ C0(20), C1(20), C2(20), D1(20), D2(20)
      COMMON /CM25/ DEBFP, LONG
      COMMON /CM30/ FV(100), XIT(100), KOUNT
      COMMON /CM35/ FD(201), FDN(201)
      COMMON /CM40/ FIK(10), NPK(10), NK
      COMMON /CM45/ WEIGHT(10)
      COMMON /CM50/ N1, N2, N3, N4, N5, MODE
      COMMON /CM55/ ALFA, BETA, NIT
      COMMON /CM60/ WFK(10), WGK(10)
      COMMON /UNIT/ IND, INP, IOUTD, IOUTP, PI, SPM
      DIMENSION X(40), G(40), IP(10), IQ(10)
C
C INITIALIZE MACHINE DEPENDENT CONSTANTS
C
      IND = I1MACH(1)
      INP = I1MACH(1)
      IOUTD = I1MACH(2)
      IOUTP = I1MACH(3)
      PI = 4.0*ATAN(1.0)
      R1 = R1MACH(1)
      SPM = SQRT(R1)
C
C INPUT DATA
C
C MODE = TYPE OF APPROXIMATION
C = 1 MAGNITUDE ONLY
C = 2 GROUP DELAY ONLY
C = 3 GROUP DELAY EQUALISATION
C = 4 MAGNITUDE AND GROUP DELAY
C
      READ (IND,9999) MODE
C
C NZR = NO OF REAL ZEROS
C NZU = NO OF ZEROS ON THE UNIT CIRCLE
C NZC = NO OF OF COMPLEX ZEROS
C NPR = NO OF REAL POLES
C NPC = NO OF OF COMPLEX POLES
C
      READ (IND,9999) NZR, NZU, NZC, NPR, NPC
      N1 = NZR
      N2 = N1 + NZU
      N3 = N2 + NZC*2
      N4 = N3 + NPR
      N5 = N4 + NPC*2 - 1
      N = N5 + 2
      IF (MODE.EQ.4) N = N + 1
C
C X = INITIAL PARAMETER VECTOR
C ZEROES AND POLES SPECIFIED IN POLAR COORDINATES AS FOLLOWS
C REAL ZEROS COME FIRST, THEN ZEROS ON THE UNIT CIRCLE, THEN COMPLEX
C ZEROS, THEN REAL POLES AND FINALLY COMPLEX POLES
C REAL ZEROES (POLES) ARE SPECIFIED AS RADIUS1,RADIUS2,...
C ZEROES ON THE UNIT CIRCLE ARE SPECIFIED AS ANGLE1,ANGLE2,...
C COMPLEX ZEROES (POLES) ARE SPECIFIED AS ANGLE1,RADIUS1,ANGLE2,
C RADIUS2,...
C WHERE THE ANGLES ARE SPECIFIED IN HZ AS PHI/(2 PI)*FSAMPLE
C AND ONLY ZEROES (POLES) WITH POSITIVE ANGLES ARE SPECIFIED
C THEIR COMPLEX CONJUGATE PAIRS BEING UNDERSTOOD
C THE LAST PARAMETER(S) IS (ARE): K0 IF MODE=1, TAU0 IF MODE=2,3,
C K0,TAU0 IF MODE=4
C
      READ (IND,9998) (X(I),I=1,N)
C
C NK = NO OF INTERVALS
C
      READ (IND,9999) NK
C
C NPK(I) = NO. OF POINTS ON THE INTERVAL FIK(I) - FIK(I+1)
C
      READ (IND,9999) (NPK(I),I=1,NK)
C
C WEIGHT(I) = SPACING OF THE POINTS NPK(I) ON THE I TH INTERVAL
C = 1   POINTS EQUISPACED
C = 2   POINTS SPACED ON A SINE ABSCISSA, HALF CYCLE
C = 3   POINTS SPACED ON A COS  ABSCISSA, HALF CYCLE
C = 4   POINTS SPACED ON A COS  ABSCISSA, FULL CYCLE
C
      READ (IND,9999) (WEIGHT(I),I=1,NK)
C
C WFK(I) = RELATIVE WEIGHTING OF THE I TH INTERVAL FOR THE MAGNITUDE
C
      READ (IND,9998) (WFK(I),I=1,NK)
C WGK(I) = RELATIVE WEIGHTING OF THE I TH INTERVAL FOR THE GROUP DELAY
C
      READ (IND,9998) (WGK(I),I=1,NK)
C
      NK = NK + 1
C
C FIK(I) = STARTING POINT OF THE  I  TH APPROX. INTERVAL
C
      READ (IND,9998) (FIK(I),I=1,NK)
C
C IPK = NO OF SUCCESSIVE INDICES OF APPROXIMATION
C
      READ (IND,9999) IPK
C
C IP(II) = SUCCESSIVE INDICES OF APPROXIMATION FOR THE MAGNITUDE
C
      READ (IND,9999) (IP(II),II=1,IPK)
C
C IQ(II) =  SUCCESSIVE INDICES OF APPROXIMATION FOR THE GROUP DELAY
C
      READ (IND,9999) (IQ(II),II=1,IPK)
C
C INV = NUMBER OF TOTAL CYCLES IF NO CONVERGENCE
C
      READ (IND,9999) INV
C
C F1 = LOWER FREQUENCY LIMIT FOR THE FREQUENCY RESPONSE
C F2 = UPPER FREQUENCY LIMIT FOR THE FREQUENCY RESPONSE
C FSAMPL = SAMPLING FREQUENCY
C IZ  = NUMBER OF POINTS FOR THE FREQUENCY RESPONSE
C (IF IZ LE 0 NO RESPONSE CALCULATED )
C
      READ (IND,9998) F1, F2, FSAMPL
      READ (IND,9999) IZ
      WS = FSAMPL/2.0
C
C ALFA= RELATIVE WEIGHTING OF THE MAGNITUDE
C RELATIVE WEIGHTING OF THE GROUP DELAY = 1.0 - ALFA
C THIS PARAMETER IS ONLY USED IF MODE = 4
C
      READ (IND,9998) ALFA
      BETA = 1.0 - ALFA
C
C LONG = LOGICAL PARAMETER FOR TYPE OF OUTPUT
C = .TRUE. FOR LONG OUTPUT
C = .FALSE. FOR SHORT OUTPUT
C
C
      READ (IND,9997) LONG
C
C LIMIT = NUMBER OF ITERATIONS PER CYCLE BEFORE RESTARTING WITH
C STEEPEST DESCENT
C
      LIMIT = 4*N
C
C EPS = ERROR PARAMETER FOR FMFP SUBROUTINE
C
      EPS = 1.E-3
C
C DEBUG  = .TRUE . FOR DEBUGGING PRINTOUT
C
      DEB = .FALSE.
C
C DEBFP  = .TRUE . FOR DEBUGGING PRINTOUT FROM FMFP  SUBROUTINE
C
      DEBFP = .FALSE.
      READ (IND,9997) DEB, DEBFP
C
9999  FORMAT (I4)
9998  FORMAT (F16.9)
9997  FORMAT (L1)
      CALL COSYFP(X, G, IZ, F1, F2, WS, N, EPS, IP, IQ, LIMIT, INV, IPK)
      STOP
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  COSYFP
C MAIN SUBROUTINE FOR AMPLITUDE AND/OR G.D. SYNTHESIS USING FP
C-----------------------------------------------------------------------
C
      SUBROUTINE COSYFP(X, G, IZ, F1, F2, WS, N, EPS, IP, IQ, LIMIT,
     *    INV, IPK)
      LOGICAL DEB, DEBFP, LONG
      INTEGER WEIGHT
      COMMON /CM05/ AS(201), DS(201), WF(201), WG(201)
      COMMON /CM10/ NP, NPD
      COMMON /CM15/ DEB
      COMMON /CM20/ C0(20), C1(20), C2(20), D1(20), D2(20)
      COMMON /CM25/ DEBFP, LONG
      COMMON /CM30/ FV(100), XIT(100), KOUNT
      COMMON /CM35/ FD(201), FDN(201)
      COMMON /CM40/ FIK(10), NPK(10), NK
      COMMON /CM45/ WEIGHT(10)
      COMMON /CM50/ N1, N2, N3, N4, N5, MODE
      COMMON /CM55/ ALFA, BETA, NIT
      COMMON /CM60/ WFK(10), WGK(10)
      COMMON /CM65/ IPX, IQX
      COMMON /CM70/ A(201), B(201)
      COMMON /UNIT/ IND, INP, IOUTD, IOUTP, PI, SPM
      DIMENSION GDH(201)
      DIMENSION X(1), G(1), H(201), U0(20), V0(20), U(20), V(20),
     *    A0(20), A1(20), A2(20), B0(20), B1(20), B2(20), AMPL(201),
     *    PHI(201), DB(201), GD(201), IP(10), IQ(10)
C
C INITIALIZE CONSTANTS
C
      NPD = 0
      NPA = 0
      NP = 0
      KP = 1
      KOUNT = 0
      NIT = 0
      IER = 0
      NIT = 0
      NKOUNT = 0
      IQX = IQ(1)
      IPX = IP(1)
      NN = N/4
      CONST = 1.0
      DO 10 I=1,20
        A0(I) = 1.0
        B0(I) = 1.0
        A1(I) = 0.0
        A2(I) = 0.0
        B1(I) = 0.0
        B2(I) = 0.0
  10  CONTINUE
C
C INITIALIZE X(N) AND X(N-1) (IF NECESSARY)
C
      GO TO (20, 30, 30, 40), MODE
  20  X(N) = 1.0
      GO TO 50
  30  X(N) = 0.0
      GO TO 50
  40  X(N-1) = 1.0
      X(N) = 0.0
  50  WSPI = PI/WS
C
C COMPUTE FREQUENCY POINTS ON THE VARIOUS INTERVALS
C
      DO 110 K=2,NK
        NPX = NPK(K-1)
        NPX1 = NPX - 1
        FNPX = FLOAT(NPX)
        FI2 = FIK(K)
        FI1 = FIK(K-1)
        IWG = WEIGHT(K-1)
        WGK1 = WGK(K-1)
        WFK1 = WFK(K-1)
C
        IF (NPX.EQ.0) GO TO 110
C
        DO 100 J=1,NPX
          JK = J + NP
          WG(JK) = WGK1
          WF(JK) = WFK1
C
          FLJNP = FLOAT(J-1)/FNPX
          GO TO (60, 70, 80, 90), IWG
  60      FDN(JK) = FLJNP*(FI2-FI1) + FI1
          GO TO 100
  70      FDN(JK) = SIN(FLJNP*PI/2.0)*(FI2-FI1) + FI1
          GO TO 100
  80      FDN(JK) = COS(FLJNP*PI/2.0)*(FI1-FI2) + FI2
          GO TO 100
  90      FDN(JK) = 0.5*(FI1+FI2) + COS(FLJNP*PI)*(FI1-FI2)*0.5
 100    CONTINUE
C
        IF (WFK1.NE.0.) NPA = NPA + NPX
        IF (WGK1.NE.0.) NPD = NPD + NPX
        NP = NP + NPX
 110  CONTINUE
C
      DO 120 J=1,NP
        FD(J) = FDN(J)*WSPI
 120  CONTINUE
C
      N1P = N1 + 1
      N2P = N2 + 1
      N4P = N4 + 1
C
C NORMALISE THE ANGLE OF THE POLES AND ZEROES
C
      IF (N1P.GT.N2) GO TO 140
      DO 130 J=N1P,N2
        X(J) = X(J)/WS
 130  CONTINUE
 140  IF (N2P.GT.N3) GO TO 160
      DO 150 J=N2P,N3,2
        X(J) = X(J)/WS
 150  CONTINUE
 160  IF (N4P.GT.N5) GO TO 180
      DO 170 J=N4P,N5,2
        X(J) = X(J)/WS
 170  CONTINUE
C
C READ IN FILTER COEFFICIENTS
C
 180  IF (MODE.NE.3) GO TO 190
      READ (INP,9999) CONST0, M
9999  FORMAT (E14.7, I4)
      CALL INCOCA(M, CONST0, C0, C1, C2, D1, D2)
      CALL DGSPEC(GDH, M)
C
C SET UP ARRAYS OF DESIRED FUNCTIONS FOR MAGN. AND GD
C
 190  DO 240 J=1,NP
        FJ = FDN(J)
        GO TO (200, 210, 220, 230), MODE
 200    AS(J) = FS(FJ)
        GO TO 240
 210    DS(J) = DGS(FJ)
        GO TO 240
 220    DS(J) = -GDH(J)
        GO TO 240
 230    AS(J) = FS(FJ)
        DS(J) = DGS(FJ)
 240  CONTINUE
C
      CALL FUNCT(N, X, F, G)
C
C COMPUTE INITIAL APPROX. TO X(N) AND X(N-1)
C
      CALL XNXM1(X, N, NP)
C
      GO TO (250, 260, 260, 270), MODE
 250  IF (LONG) WRITE (IOUTD,9983) X(N)
      GO TO 280
 260  IF (LONG) WRITE (IOUTD,9982) X(N)
      GO TO 280
 270  IF (LONG) WRITE (IOUTD,9981) X(N-1), X(N)
 280  CALL FUNCT(N, X, F, G)
      IF (LONG) WRITE (IOUTD,9998) F, IER, KOUNT, NIT
      IF (LONG) WRITE (IOUTD,9994)
      IF (LONG) WRITE (IOUTD,9995) (J,X(J),J,G(J),J=1,N)
C
C MAIN LOOP CALLING THE F.P. SUBROUTINE
C
 290  CALL FMFP(N, X, F, G, EPS, LIMIT, IER, H)
C
      NKOUNT = NKOUNT + KOUNT
      IF (LONG) WRITE (IOUTD,9998) F, IER, KOUNT, NIT
      IF (LONG) WRITE (IOUTD,9995) (J,X(J),J,G(J),J=1,N)
C
      IF (IER.EQ.0) GO TO 310
      IF (INV-1) 310, 310, 300
 300  INV = INV - 1
      GO TO 290
 310  IF (IER.NE.0 .OR. KP.GE.IPK) GO TO 320
      KP = KP + 1
      IQX = IQ(KP)
      IPX = IP(KP)
      GO TO 290
C
C COMPUTE QUADRATIC FACTORS, POLES AND ZEROS FROM ARG. VECTOR
C
 320  CALL FACTOR(X, A1, A2, B1, B2, U0, V0, U, V, CONST, JA, JB, JZ,
     *    JP, WS)
C
      NN = MAX0(JA-1,JB-1)
      JZ1 = JZ - 1
      JP1 = JP - 1
C
      IF (MODE.NE.3) GO TO 350
      JA = JA + M
      JB = JB + M
      NNP = NN + 1
      NM = NN + M
C
      IF (NNP.GT.NM) GO TO 340
      DO 330 J=NNP,NM
        JM = J - NN
        A1(J) = C1(JM)
        A2(J) = C2(JM)
        B1(J) = D1(JM)
        B2(J) = D2(JM)
 330  CONTINUE
 340  NN = NM
C
 350  IF (DEB) WRITE (IOUTD,9980) JA, JB
      IF (DEB) WRITE (IOUTD,9979) (A1(J),A2(J),B1(J),B2(J),J=1,NN)
C
      GO TO (360, 390, 370, 380), MODE
 360  CONST = X(N)
      GO TO 390
 370  CONST = CONST*CONST0
      GO TO 390
 380  CONST = X(N-1)
C
C OUTPUT SECTION
C
 390  WRITE (IOUTD,9996)
      WRITE (IOUTD,9986) IPX, IQX
      WRITE (IOUTD,9998) F, IER, NKOUNT, NIT
      WRITE (IOUTD,9993)
      WRITE (IOUTD,9995) (J,X(J),J,G(J),J=1,N)
      IF (JZ1.GT.0) WRITE (IOUTD,9989) (U0(J),V0(J),J=1,JZ1)
      IF (JP1.GT.0) WRITE (IOUTD,9988) (U(J),V(J),J=1,JP1)
      WRITE (IOUTD,9987) CONST
      WRITE (IOUTD,9985) (A0(J),A1(J),A2(J),J=1,NN)
      WRITE (IOUTD,9984) (B0(J),B1(J),B2(J),J=1,NN)
      IF (IZ.LE.0) GO TO 400
      CALL FREDIC(NN, CONST, A0, A1, A2, B1, B2, F1, F2, WS, IZ, AMPL,
     *    PHI, DB, GD, FD)
      WRITE (IOUTD,9997) (FD(J),AMPL(J),DB(J),PHI(J),GD(J),J=1,IZ)
C
 400  WRITE (IOUTP,9992) NN
      WRITE (IOUTP,9991) CONST, NN
      WRITE (IOUTP,9990) (A0(J),A1(J),A2(J),B1(J),B2(J),J=1,NN)
      RETURN
C
C FORMATS
C
9998  FORMAT (/3H F=, E12.5, 5X, 4HIER=, I2, 5X, 17HNO OF ITERATIONS=,
     *    I3//1X, 24HNO OF FUNCT EVALUATIONS=, I4)
9997  FORMAT (1H1//19H FREQUENCY RESPONSE//1X, 9HFREQUENCY, 6X,
     *    9HMAGNITUDE, 5X, 10HLOSS IN DB, 8X, 5HPHASE, 7X, 9HGROUP DEL,
     *    2HAY//(E11.4, 4E15.5))
9996  FORMAT (1H1//49H SYNTHESIS OF DIGITAL FILTER WITH SPEC. MAGNITUDE,
     *    20H AND/OR GROUP DELAY /)
9995  FORMAT (/(10X, 2HX(, I2, 2H)=, E15.8, 10X, 2HG(, I2, 2H)=, E15.8))
9994  FORMAT (//33H INITIAL ARG. VECTOR AND GRADIENT)
9993  FORMAT (//34H COMPUTED ARG. VECTOR AND GRADIENT)
9992  FORMAT (52H COEFFICIENTS OF DIGITAL FILTER WITH SPEC. AMPLITUDE,
     *    19H AND/OR GROUP DELAY, / 3H N=, I2)
9991  FORMAT (1X, E14.7, I4)
9990  FORMAT (1X, 5E14.7)
9989  FORMAT (/24X, 13HZ PLANE ZEROS//12X, 12HANGLE/2PI*FS, 16X,
     *    6HRADIUS//(2E25.8))
9988  FORMAT (/24X, 13HZ PLANE POLES//12X, 12HANGLE/2PI*FS, 16X,
     *    6HRADIUS//(2E25.8))
9987  FORMAT (/1X, 6HCONST=, E15.8)
9986  FORMAT (/36H THE ERROR CRITERION USED IS MINIMUM, I2, 9H FOR THE ,
     *    13HMAGNITUDE AND//29X, 7HMINIMUM, I2, 14H FOR THE GROUP,
     *    6H DELAY)
9985  FORMAT (//24X, 22HNUMERATOR COEFFICIENTS///14X, 2HA0, 20X, 2HA1,
     *    20X, 2HA2//(3E22.8))
9984  FORMAT (//23X, 24HDENOMINATOR COEFFICIENTS///14X, 2HB0, 20X, 2HB1,
     *    20X, 2HB2//(3E22.8))
9983  FORMAT (1H1//5H K0 =, E15.8)
9982  FORMAT (1H1//5H T0 =, E15.8)
9981  FORMAT (1H1//5H K0 =, E15.8, 10X, 5H T0 =, E15.8)
9980  FORMAT (1H1//2I10//)
9979  FORMAT (//(4E18.5))
C
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  XNXM1
C SUBROUTINE TO CALCULATE INITIAL APPROX. TO X(N),X(N-1)
C-----------------------------------------------------------------------
C
      SUBROUTINE XNXM1(X, N, NP)
      DIMENSION X(N)
      COMMON /CM05/ AS(201), DS(201), WF(201), WG(201)
      COMMON /CM50/ N1, N2, N3, N4, N5, MODE
      COMMON /CM70/ A(201), B(201)
      COMMON /UNIT/ IND, INP, IOUTD, IOUTP, PI, SPM
C
      AGD = 0.0
      AGN = 0.0
      AD = 0.0
      AN = 0.0
      DO 40 J=1,NP
        GO TO (10, 20, 20, 30), MODE
  10    WFJ = WF(J)
        IF (WFJ.LT.SPM) WFJ = SPM
        HJ = A(J)/WFJ + AS(J)
        WHJ = HJ*WFJ
        AN = AN + WHJ*AS(J)
        AD = AD + HJ*WHJ
        GO TO 40
  20    AGN = AGN + B(J)
        AGD = AGD + WG(J)
        GO TO 40
  30    WFJ = WF(J)
        IF (WFJ.LT.SPM) WFJ = SPM
        HJ = A(J)/WFJ + AS(J)
        WHJ = HJ*WFJ
        AN = AN + WHJ*AS(J)
        AD = AD + HJ*WHJ
        AGN = AGN + B(J)
        AGD = AGD + WG(J)
  40  CONTINUE
C
      GO TO (50, 60, 60, 70), MODE
  50  X(N) = AN/AD
      GO TO 80
  60  X(N) = AGN/AGD
      GO TO 80
  70  X(N) = AGN/AGD
      X(N-1) = AN/AD
  80  RETURN
      END
C
C-----------------------------------------------------------------------
C FUNCTION:  RANGE
C FUNCTION TO LIMIT THE RANGE OF ANGLES (X) TO 0 TO PI
C-----------------------------------------------------------------------
C
      FUNCTION RANGE(X)
C
      X = ABS(X)
      IX = IFIX(X)
      MX = 2*(IX/2)
      RANGE = X - FLOAT(MX)
      IF (RANGE.GT.1.0) RANGE = 2.0 - RANGE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  DGSPEC
C TO CALCULATE THE GROUP DELAY OF THE NON-EQUALISED FILTER
C-----------------------------------------------------------------------
C
      SUBROUTINE DGSPEC(GDH, M)
      DIMENSION GDH(201)
      COMMON /CM10/ NP, NPD
      COMMON /CM20/ C0(20), C1(20), C2(20), D1(20), D2(20)
      COMMON /CM35/ FD(201), FDN(201)
      COMMON /UNIT/ IND, INP, IOUTD, IOUTP, PI, SPM
      COMPLEX NUMC, NUMD, DENC, DEND
C
      DO 20 J=1,NP
        FJ = FD(J)
        GDH(J) = 0.
        CS = COS(FJ)
        SN = SIN(FJ)
        CS2 = CS + CS
        SN2 = SN + SN
C
        DO 10 I=1,M
          C11 = C0(I) + C2(I)
          C12 = C1(I)
          C22 = C0(I) - C2(I)
          D11 = 1. + D2(I)
          D12 = D1(I)
          D22 = 1. - D2(I)
          NUMC = CMPLX(CS2+C12,SN2)
          NUMD = CMPLX(CS2+D12,SN2)
          DENC = CMPLX(C11*CS+C12,C22*SN)
          DEND = CMPLX(D11*CS+D12,D22*SN)
          GDH(J) = GDH(J) + REAL(NUMD/DEND) - REAL(NUMC/DENC)
  10    CONTINUE
C
  20  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  INCOCA
C READ COEFFICIENTS FOR ONE DIGITAL CASCADE FILTER
C E14.7 FORMAT, FILE ON UNIT INP
C-----------------------------------------------------------------------
C
      SUBROUTINE INCOCA(NTERM, CONST, A0, A1, A2, B1, B2)
      COMMON /UNIT/ IND, INP, IOUTD, IOUTP, PI, SPM
      REAL A0(NTERM), A1(NTERM), A2(NTERM), B1(NTERM), B2(NTERM)
      DO 10 N=1,NTERM
        READ (INP,9999) A0(N), A1(N), A2(N), B1(N), B2(N)
  10  CONTINUE
9999  FORMAT (5E14.7)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  FREDIC
C SUBROUTINE TO COMPUTE THE FREQUENCY RESPONSE OF A DIGITAL FILTER IN
C CASCADE FORM
C-----------------------------------------------------------------------
C
      SUBROUTINE FREDIC(NTERM, CONST, A0, A1, A2, B1, B2, FU, FO, WS,
     *    IZ, AMPL, PHI, DB, GD, FD)
      COMMON /UNIT/ IND, INP, IOUTD, IOUTP, PI, SPM
      COMPLEX N, D, GWS, NN, DD, N1, N2
      DIMENSION A0(NTERM), A1(NTERM), A2(NTERM), B1(NTERM), B2(NTERM)
      DIMENSION AMPL(IZ), PHI(IZ), DB(IZ), FD(IZ), GD(IZ)
      DIMENSION AA(40), AB(40), BA(40), BB(40)
C
      WSXPI = PI/WS
      WU = FU
      Z = (FO-FU)/FLOAT(IZ-1)
C
      DO 10 I=1,NTERM
        AA(I) = A0(I) + A2(I)
        AB(I) = A0(I) - A2(I)
        BA(I) = 1.0 + B2(I)
        BB(I) = 1.0 - B2(I)
  10  CONTINUE
C
      DO 40 J=1,IZ
        WT = WU + Z*FLOAT(J-1)
        FD(J) = WT
        WT = WT*WSXPI
        CWT = COS(WT)
        SWT = SIN(WT)
        D = CMPLX(1.0,0.)
        N = D
        QL = 0.0
        QQ = QL
C
        DO 20 I=1,NTERM
          NN = CMPLX(AA(I)*CWT+A1(I),AB(I)*SWT)
          DD = CMPLX(BA(I)*CWT+B1(I),BB(I)*SWT)
          N2 = CMPLX(2.0*CWT+B1(I),2.0*SWT)
          N1 = CMPLX(2.0*A0(I)*CWT+A1(I),2.0*A0(I)*SWT)
          N = N*NN
          D = D*DD
C
          IF (CABS(DD).LT.SPM) DD = CMPLX(SPM,SPM)
          Q2 = REAL(N2/DD)
C
          IF (CABS(NN).LT.SPM) NN = CMPLX(SPM,SPM)
          Q1 = REAL(N1/NN)
          IF (ABS(Q1-QL).GT.50.) Q1 = 1.0
          QL = Q1
          QQ = QQ + Q2 - Q1
  20    CONTINUE
C
        GD(J) = QQ
C
        IF (CABS(D).LT.SPM) D = CMPLX(SPM,SPM)
        GWS = N/D
C
        AMPL(J) = CONST*CABS(GWS)
        GABS = AMPL(J)
        IF (GABS.NE.0.) GO TO 30
C
        PHI(J) = 0.0
        DB(J) = 1000.
        GO TO 40
C
  30    GWSI = AIMAG(GWS)
        GWSR = REAL(GWS)
        PHI(J) = ATAN(GWSI/GWSR)
        DB(J) = -ALOG10(GABS)*20.0
  40  CONTINUE
C
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  FACTOR
C SUBROUTINE TO RECALCULATE QUADRATIC FACTORS
C AND POLES AND ZEROS FROM PARAMETER VECTOR X
C-----------------------------------------------------------------------
C
      SUBROUTINE FACTOR(X, A1, A2, B1, B2, U0, V0, U, V, CONST, JA, JB,
     *    JZ, JP, WS)
      COMMON /CM50/ N1, N2, N3, N4, N5, MODE
      COMMON /UNIT/ IND, INP, IOUTD, IOUTP, PI, SPM
      DIMENSION X(1), A1(1), A2(1), B1(1), B2(1), U0(1), V0(1)
      DIMENSION U(1), V(1)
C
      J = 0
      JA = 1
      JB = 1
      JZ = 1
      CONST = 1.0
  10  IF (J.GE.N1) GO TO 40
      J = J + 1
      R = X(J)
      GO TO (20, 30, 30, 30), MODE
  20  A2(JA) = 0.0
      A1(JA) = -R
      U0(JZ) = 0.
      V0(JZ) = R
      JA = JA + 1
      JZ = JZ + 1
      GO TO 80
  30  A1(JA) = -R - 1.0/R
      A2(JA) = 1.0
      U0(JZ) = 0.0
      V0(JZ) = R
      U0(JZ+1) = 0.0
      V0(JZ+1) = 1.0/R
      JA = JA + 1
      JZ = JZ + 2
      GO TO 80
C
  40  IF (J.GE.N2) GO TO 50
      J = J + 1
      XJ = X(J)
      XJ = RANGE(XJ)
      A1(JA) = -2.0*COS(XJ*PI)
      A2(JA) = 1.0
      V0(JZ) = 1.0
      U0(JZ) = XJ*WS
      JA = JA + 1
      JZ = JZ + 1
      GO TO 80
C
  50  IF (J.GE.N3) GO TO 90
      J = J + 2
      XJM = X(J-1)
      XJM = RANGE(XJM)
      CSP = COS(XJM*PI)
      R = X(J)
      GO TO (60, 70, 70, 70), MODE
  60  A2(JA) = R*R
      A1(JA) = -(R+R)*CSP
      U0(JZ) = XJM*WS
      V0(JZ) = R
      JA = JA + 1
      JZ = JZ + 1
      GO TO 80
  70  A2(JA) = R*R
      A1(JA) = -(R+R)*CSP
      A2(JA+1) = 1.0/A2(JA)
      A1(JA+1) = -2.0/R*CSP
      V0(JZ) = R
      V0(JZ+1) = 1.0/R
      U0(JZ) = XJM*WS
      U0(JZ+1) = XJM*WS
      JA = JA + 2
      JZ = JZ + 2
  80  IF (J.LT.N3) GO TO 10
C
  90  NN = MAX0(N4,N5)
      JP = 1
 100  IF (J.GE.N4) GO TO 140
      J = J + 1
      R = X(J)
      GO TO (110, 120, 120, 110), MODE
 110  B2(JB) = 0.0
      B1(JB) = -R
      U(JP) = 0.0
      V(JP) = R
      GO TO 130
 120  B2(JB) = 0.0
      B1(JB) = -R
      A1(JA) = -1.0/R
      A2(JA) = 0.0
      CONST = CONST*R
      U(JP) = 0.0
      V(JP) = R
      U0(JZ) = 0.0
      V0(JZ) = R
      JA = JA + 1
      JZ = JZ + 1
 130  JB = JB + 1
      JP = JP + 1
      GO TO 180
C
 140  IF (J.GT.N5) GO TO 190
      J = J + 2
      R = X(J)
      XJM = X(J-1)
      XJM = RANGE(XJM)
      GO TO (150, 160, 160, 150), MODE
 150  B2(JB) = R*R
      B1(JB) = -(R+R)*COS(XJM*PI)
      U(JP) = XJM*WS
      V(JP) = R
      GO TO 170
 160  B2(JB) = R*R
      B1(JB) = -(R+R)*COS(XJM*PI)
      A1(JA) = B1(JB)/B2(JB)
      A2(JA) = 1.0/B2(JB)
      CONST = CONST*B2(JB)
      U(JP) = XJM*WS
      V(JP) = R
      U0(JZ) = XJM*WS
      V0(JZ) = 1.0/R
      JA = JA + 1
      JZ = JZ + 1
 170  JB = JB + 1
      JP = JP + 1
 180  IF (J.LT.NN) GO TO 100
C
 190  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  FUNCT
C SUBROUTINE TO CALCULATE THE FUNCTION TO BE MINIMIZED
C FROM THE PARAMETER VECTOR X AS A FUNCTION OF FREQUENCY
C-----------------------------------------------------------------------
C
      SUBROUTINE FUNCT(N, X, F, G)
      REAL NUM, NUMM, NUMP
      LOGICAL DEB
      COMMON /CM05/ AS(201), DS(201), WF(201), WG(201)
      COMMON /CM10/ NP, NPD
      COMMON /CM15/ DEB
      COMMON /CM35/ FD(201), FDN(201)
      COMMON /CM50/ N1, N2, N3, N4, N5, MODE
      COMMON /CM55/ ALFA, BETA, NIT
      COMMON /CM65/ IP, IQ
      COMMON /CM70/ A(201), B(201)
      COMMON /UNIT/ IND, INP, IOUTD, IOUTP, PI, SPM
      DIMENSION X(N), G(N), GZR(20), GZU(20), GPR(20), GZCP(20),
     *    GZCR(20), GPCP(20), GPCR(20), GPRD(20), GPCPD(20), GPCRD(20)
C
      RSPM = SQRT(SPM)
      F = 0.0
C
      DO 10 J=1,N
        G(J) = 0.0
  10  CONTINUE
      XN1 = X(N-1)
      XN = X(N)
      NIT = NIT + 1
      IF (DEB) WRITE (IOUTD,9999) NIT
9999  FORMAT (//27H NO. OF FUNCT. EVALUATIONS=, I4//)
C
      DO 400 I=1,NP
C
        FI = FD(I)
        CS = COS(FI)
        CS2 = 2.0*CS
        H = XN
        IF (MODE.EQ.4) H = XN1
        HD = 1.0
        HN = 1.0
        DG = 0.0
C
        NN = MAX0(N4,N5)
C
        J = 0
  20    J = J + 1
        IF (J.GT.N1) GO TO 50
        R = X(J)
        GO TO (30, 190, 190, 40), MODE
  30    NUM = 1. + R*R - R*CS2
        IF (NUM.LT.SPM) NUM = SPM
        HN = HN*NUM
        GZR(J) = (R-CS)/NUM
        GO TO 190
  40    NUM = 1.0/R + R - CS2
        IF (NUM.LT.SPM) NUM = SPM
        HN = HN*NUM**2
        GZR(J) = (1.0-1.0/R**2)/NUM
        GO TO 190
C
  50    IF (J.GT.N2) GO TO 70
        GO TO (60, 190, 190, 60), MODE
  60    P = X(J)*PI
        NUMP = 2.0*(1.0-COS(FI+P))
        NUMM = 2.0*(1.0-COS(FI-P))
        IF (NUMM.LT.SPM) NUMM = SPM
        IF (NUMP.LT.SPM) NUMP = SPM
        HN = HN*NUMM*NUMP
        GZU(J) = -SIN(FI-P)/NUMM + SIN(FI+P)/NUMP
        GO TO 190
C
  70    IF (J.GE.N3) GO TO 100
        R = X(J+1)
        P = X(J)*PI
        GO TO (80, 190, 190, 90), MODE
  80    RP = 1.0 + R*R
        R2 = R + R
        CSNM = COS(FI-P)
        CSNP = COS(FI+P)
        NUMM = RP - R2*CSNM
        NUMP = RP - R2*CSNP
        IF (NUMP.LT.SPM) NUMP = SPM
        IF (NUMM.LT.SPM) NUMM = SPM
        HN = HN*NUMM*NUMP
        GZCP(J) = R*(-SIN(FI-P)/NUMM+SIN(FI+P)/NUMP)
        GZCR(J) = (R-CSNM)/NUMM + (R-CSNP)/NUMP
        J = J + 1
        GO TO 190
  90    CSNP = 2.0*COS(FI+P)
        CSNM = 2.0*COS(FI-P)
        R2 = R + 1.0/R
        NUMM = R2 - CSNM
        NUMP = R2 - CSNP
        IF (NUMP.LT.RSPM) NUMP = RSPM
        IF (NUMM.LT.RSPM) NUMM = RSPM
        HN = HN*(NUMM*NUMP)**2
        GZCP(J) = 2.0*(-SIN(FI-P)/NUMM+SIN(FI+P)/NUMP)
        GZCR(J) = (1.0-1.0/R**2)*(1.0/NUMP+1.0/NUMM)
        J = J + 1
        GO TO 190
C
 100    IF (J.GT.N4) GO TO 140
        R = X(J)
        RSQ = R**2
        R2 = 2.0*R
        DEN = 1.0 + RSQ - R2*CS
        IF (DEN.LT.RSPM) DEN = RSPM
        GO TO (110, 120, 120, 130), MODE
 110    HD = HD*DEN
        GPR(J) = (R-CS)/DEN
        GO TO 190
 120    DG = DG + 2.0*(1.0-R*CS)/DEN
        GPRD(J) = 2.0*((1.0+RSQ)*CS-R2)/DEN**2
        GO TO 190
 130    HD = HD*DEN
        GPR(J) = (R-CS)/DEN
        DG = DG + (1.0-R*CS)/DEN
        GPRD(J) = ((1.0+RSQ)*CS-R2)/DEN**2
        GO TO 190
C
 140    P = X(J)*PI
        R = X(J+1)
        R2 = R + R
        RP = 1.0 + R*R
        CSNP = COS(FI+P)
        CSNM = COS(FI-P)
        SNDP = SIN(FI+P)
        SNDM = SIN(FI-P)
        DENP = RP - R2*CSNP
        DENM = RP - R2*CSNM
        IF (ABS(DENP).LT.RSPM) DENP = RSPM
        IF (ABS(DENM).LT.RSPM) DENM = RSPM
        GO TO (150, 160, 160, 170), MODE
 150    HD = HD*DENM*DENP
        GPCR(J) = (R-CSNM)/DENM + (R-CSNP)/DENP
        GPCP(J) = R*(-SNDM/DENM+SNDP/DENP)
        GO TO 180
 160    DENP2 = DENP**2
        DENM2 = DENM**2
        DG = DG + 2.0*((1.0-R*CSNM)/DENM+(1.0-R*CSNP)/DENP)
        GPCPD(J) = R2*(2.0-RP)*(SNDM/DENM2-SNDP/DENP2)
        GPCRD(J) = 2.0*((RP*CSNM-R2)/DENM2+(RP*CSNP-R2)/DENP2)
        GO TO 180
 170    HD = HD*DENM*DENP
        GPCR(J) = (R-CSNM)/DENM + (R-CSNP)/DENP
        GPCP(J) = R*(-SNDM/DENM+SNDP/DENP)
        DENP2 = DENP**2
        DENM2 = DENM**2
        DG = DG + (1.0-R*CSNM)/DENM + (1.0-R*CSNP)/DENP
        GPCPD(J) = R*(2.0-RP)*(SNDM/DENM2-SNDP/DENP2)
        GPCRD(J) = (RP*CSNM-R2)/DENM2 + (RP*CSNP-R2)/DENP2
 180    J = J + 1
C
 190    IF (J.LT.NN) GO TO 20
        XNP = FLOAT(NP)
        XIP = FLOAT(IP)
        XNPD = FLOAT(NPD)
        XIQ = FLOAT(IQ)
        GO TO (200, 210, 210, 220), MODE
 200    H = H*SQRT(HN/HD)
        E = H - AS(I)
        IF (ABS(E).LT.SPM) E = SPM
        A(I) = E*WF(I)
        EP = E**IP/XNP
        E2H = XIP*(EP/E)*H*WF(I)
        E2P = E2H*PI
        GO TO 230
C
 210    ED = DG - DS(I) - XN
        IF (ABS(ED).LT.SPM) ED = SPM
        B(I) = ED*WG(I)
        EDP = ED**IQ/XNPD
        ED2H = XIQ*(EDP/ED)*WG(I)
        ED2P = ED2H*PI
        GO TO 230
C
 220    H = H*SQRT(HN/HD)
        E = H - AS(I)
        IF (ABS(E).LT.SPM) E = SPM
        A(I) = E*WF(I)
        EP = E**IP/XNP
        E2H = XIP*(EP/E)*H*WF(I)*ALFA
        E2P = E2H*PI
        ED = DG - DS(I) - XN
        IF (ABS(ED).LT.SPM) ED = SPM
        B(I) = ED*WG(I)
        EDP = ED**IQ/XNPD
        ED2H = XIQ*(EDP/ED)*WG(I)*BETA
        ED2P = ED2H*PI
C
 230    J = 0
 240    J = J + 1
        IF (J.GT.N1) GO TO 250
        G(J) = G(J) + E2H*GZR(J)
        GO TO 360
C
 250    IF (J.GT.N2) GO TO 260
        G(J) = G(J) + E2P*GZU(J)
        GO TO 360
C
 260    IF (J.GE.N3) GO TO 270
        JJ = J + 1
        G(J) = G(J) + E2P*GZCP(J)
        G(JJ) = G(JJ) + E2H*GZCR(J)
        J = JJ
        GO TO 360
C
 270    IF (J.GT.N4) GO TO 310
        GO TO (280, 290, 290, 300), MODE
 280    G(J) = G(J) - E2H*GPR(J)
        GO TO 360
 290    G(J) = G(J) + ED2H*GPRD(J)
        GO TO 360
 300    G(J) = G(J) - E2H*GPR(J) + ED2H*GPRD(J)
        GO TO 360
C
 310    JJ = J + 1
        GO TO (320, 330, 330, 340), MODE
 320    G(J) = G(J) - E2P*GPCP(J)
        G(JJ) = G(JJ) - E2H*GPCR(J)
        GO TO 350
 330    G(J) = G(J) + ED2P*GPCPD(J)
        G(JJ) = G(JJ) + ED2H*GPCRD(J)
        GO TO 350
 340    G(J) = G(J) - E2P*GPCP(J) + ED2P*GPCPD(J)
        G(JJ) = G(JJ) - E2H*GPCR(J) + ED2H*GPCRD(J)
 350    J = JJ
 360    IF (J.LT.NN) GO TO 240
C
        GO TO (370, 380, 380, 390), MODE
 370    F = F + EP*WF(I)
        G(N) = G(N) + E2H/XN
        GO TO 400
 380    F = F + EDP*WG(I)
        G(N) = G(N) - ED2H
        GO TO 400
 390    F = F + EP*WF(I)*ALFA + EDP*WG(I)*BETA
        G(N-1) = G(N-1) + E2H/XN1
        G(N) = G(N) - ED2H
C
 400  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  FMFP
C TO FIND THE LOCAL MINIMUM OF A FUNCTION OF SEVERAL VARIABLES
C BY THE METHOD OF FLETCHER AND POWELL
C-----------------------------------------------------------------------
C
      SUBROUTINE FMFP(N, X, F, G, EPS, LIMIT, IER, H)
      COMMON /CM25/ DEB, LONG
      COMMON /CM30/ FV(100), XIT(100), KOUNT
      COMMON /CM50/ M1, M2, M3, M4, M5, MODE
      COMMON /UNIT/ IND, INP, IOUTD, IOUTP, PI, SPM
      LOGICAL DEB, LONG
      DIMENSION H(201), X(N), G(N), Y(40)
C
C COMPUTE FUNCTION VALUE AND GRADIENT VECTOR FOR INITIAL ARGUME
C
      CALL FUNCT(N, X, F, G)
C
C RESET ITERATION COUNTER AND GENERATE IDENTITY MATRIX
C
      KOUNT = 0
      IER = 0
      N2 = N + N
      N3 = N2 + N
      N31 = N3 + 1
  10  K = N31
      IF (DEB) WRITE (IOUTD,9982)
      DO 30 J=1,N
        NJ = N - J
        H(K) = 1.0
        IF (NJ.LE.0) GO TO 40
        DO 20 L=1,NJ
          KL = K + L
          H(KL) = 0.0
  20    CONTINUE
        K = KL + 1
  30  CONTINUE
C
C START ITERATION LOOP, UPDATE COUNTER AND SAVE FUNCTION VALUE
C
  40  KOUNT = KOUNT + 1
      FV(KOUNT) = F
      XIT(KOUNT) = KOUNT
C
C SAVE FUNCTION VALUE, ARGUMENT VECTOR AND GRADIENT VECTOR
C
      IF (DEB .OR. LONG) WRITE (IOUTD,9999) F, KOUNT
      IF (DEB) WRITE (IOUTD,9997)
      IF (DEB) WRITE (IOUTD,9998) (X(J),G(J),J=1,N)
      OLDF = F
C
C DETERMINE DIRECTION VECTOR H
C
      DO 70 J=1,N
        K = N + J
        H(K) = G(J)
        K = K + N
        H(K) = X(J)
        T = 0.0
        K = J + N3
        DO 60 L=1,N
          T = T - G(L)*H(K)
          IF (L.GE.J) GO TO 50
          K = K + N - L
          GO TO 60
  50      K = K + 1
  60    CONTINUE
        H(J) = T
  70  CONTINUE
      IF (DEB) WRITE (IOUTD,9981)
      IF (DEB) WRITE (IOUTD,9996) (J,H(J),J=1,N)
C
C CHECK WHETHER FUNCTION WILL DECREASE STEPPING ALONG H
C
      GNRM = 0.0
      HNRM = 0.0
      DY = 0.0
C
C CALCULATE DIRECTIONAL DERIVATIVE AND TEST VALUES FOR DIRECTION
C VECTOR H AND GRADIENT VECTOR G.
C
      DO 80 J=1,N
        HNRM = HNRM + ABS(H(J))
        GNRM = GNRM + ABS(G(J))
        DY = DY + H(J)*G(J)
  80  CONTINUE
      IF (DEB) WRITE (IOUTD,9995) HNRM, GNRM
      IF (DEB) WRITE (IOUTD,9992) DY, F
C
C REPEAT SEARCH IN DIRECTION OF STEEPEST DESCENT IF DIRECTIONAL
C DERIVATIVE APPEARS TO BE POSITIVE OR ZERO OR IF H IS
C
C
      IF (DY.GE.0.0 .OR. (HNRM/GNRM).LE.EPS) GO TO 500
C
C SEARCH MINIMUM ALONG DIRECTION H
C SEARCH ALONG H FOR POSITIVE DIRECTIONAL DERIVATIVE
C
      FY = F
      HN = 1.0
      IF (HNRM.GT.1.0) HN = 1.0/HNRM
      EPSC = 2.*R1MACH(3)
      UC = 1.0 - EPSC
C
C FIND MAX. STEP SIZE TO ENSURE STABILITY
C
      IP1 = M3 + 1
      IP2 = M4
      IP3 = M4 + 2
      IP4 = M5 + 1
C
C CONSIDER REAL POLES
C
      IF (M3.EQ.M4) GO TO 100
      DO 90 J=IP1,IP2
        IF (ABS(X(J)+H(J)*HN).LT.UC) GO TO 90
        IF (H(J).LT.0.0) HN = (-UC-X(J))/H(J)
        IF (H(J).GT.0.0) HN = (UC-X(J))/H(J)
  90  CONTINUE
C
C CONSIDER COMPLEX POLES
C
 100  IF (M4.EQ.M5+1) GO TO 120
      DO 110 J=IP3,IP4,2
        R = X(J) + H(J)*HN
        IF (R.GT.0. .AND. R.LT.UC) GO TO 110
        IF (H(J).GT.0.0) HN = (UC-X(J))/H(J)
        IF (H(J).LT.0.0) HN = -X(J)/H(J)
 110  CONTINUE
 120  IF (DEB) WRITE (IOUTD,9994) HN
      DO 130 J=1,N
        Y(J) = X(J) + HN*H(J)
 130  CONTINUE
      AMBDA = HN
      CALL FUNCT(N, Y, F, G)
      DYN = DY*HN
      DYY = 0.0
      DO 140 J=1,N
        DYY = DYY + H(J)*G(J)
 140  CONTINUE
      IF (DEB) WRITE (IOUTD,9991) DYY, F
      IF (DYY.NE.0.0) GO TO 160
      DO 150 J=1,N
        X(J) = Y(J)
 150  CONTINUE
      GO TO 390
 160  CONTINUE
      DEN = F - FY - DYN
      IF (DEN.EQ.0.0) GO TO 170
      ALFA = -0.5*DYN/DEN*HN
      IF (DEB) WRITE (IOUTD,9990) ALFA, HN
C
C USE ESTIMATED STEP SIZE IF POSITIVE OTHERWISE USE HN
C
      IF (ALFA.LE.0.0) GO TO 170
      IF (F.GT.FY .AND. DYY.LT.0.) GO TO 170
      IF (DYY*(ALFA-HN).GE.0.0) GO TO 170
      AMBDA = ALFA
 170  ALFA = 0.0
C
C SAVE FUNCTION AND DERIVATIVE VALUES FOR OLD ARGUMENT
C
 180  DX = DY
      FX = FY
C
C STEP ARGUMENT ALONG H
C FIND MAX. STEP SIZE TO ENSURE STABILITY
C CONSIDER REAL POLES FIRST
C
      IF (M3.EQ.M4) GO TO 200
      DO 190 J=IP1,IP2
        IF (ABS(X(J)+H(J)*AMBDA).LT.UC) GO TO 190
        IF (H(J).LT.0.0) AMBDA = (-UC-X(J))/H(J)
        IF (H(J).GT.0.0) AMBDA = (UC-X(J))/H(J)
 190  CONTINUE
C
C CONSIDER COMPLEX POLES
C
 200  IF (M4.EQ.M5+1) GO TO 220
      DO 210 J=IP3,IP4,2
        R = X(J) + H(J)*AMBDA
        IF (R.GT.0. .AND. R.LT.UC) GO TO 210
        IF (H(J).GT.0.0) AMBDA = (UC-X(J))/H(J)
        IF (H(J).LT.0.0) AMBDA = -X(J)/H(J)
 210  CONTINUE
 220  DO 230 I=1,N
        X(I) = X(I) + AMBDA*H(I)
 230  CONTINUE
      IF (DEB) WRITE (IOUTD,9994) AMBDA
      IF (DEB) WRITE (IOUTD,9993) (X(I),I=1,N)
C
C COMPUTE FUNCTION VALUE AND GRADIENT FOR NEW ARGUMENT
C
      CALL FUNCT(N, X, F, G)
      FY = F
C
C COMPUTE DIRECTIONAL DERIVATIVE DY FOR NEW ARGUMENT. TERMINATE
C SEARCH IF DY IS POSITIVE. IF DY IS ZERO THE MINIMUM IS FOUND
C
      DY = 0.0
      DO 240 I=1,N
        DY = DY + H(I)*G(I)
 240  CONTINUE
      IF (DEB) WRITE (IOUTD,9989) DY, F
      DYMN = AMIN1(-DX,DY)
      OF = F
      IF (DY) 250, 400, 260
C
C TERMINATE SEARCH ALSO IF THE FUNCTION VALUE INDICATES THAT
C A MINIMUM HAS BEEN PASSED
C
 250  IF (FY.GE.FX) GO TO 260
      IF (ABS(DX-DY).LE.EPSC) GO TO 520
C
C REPEAT SEARCH AND DOUBLE STEP SIZE FOR FURTHER SEARCHES
C
      AMBDA = AMBDA + ALFA
      ALFA = AMBDA
C
C END OF SEARCH LOOP
C TERMINATE IF THE CHANGE IN ARGUMENT GETS VERY LARGE
C
      IF (HNRM*AMBDA.LE.1.0) GO TO 180
C
C LINEAR SEARCH TECHNIQUE INDICATES THAT NO MINIMUM EXISTS
C
      IER = 2
      RETURN
C
C INTERPOLATE CUBICALLY IN THE INTERVAL DEFINED BY THE SEARCH
C ABOVE AND COMPUTE THE ARGUMENT X FOR WHICH THE INTERPOLATION
C POLYNOMIAL IS MINIMIZED
C
 260  T = 0.0
 270  CONTINUE
      IF (DEB) WRITE (IOUTD,9987) FX, FY
      IF (DEB) WRITE (IOUTD,9986) DX, DY
      IF (AMBDA.EQ.0.0) GO TO 390
      Z = 3.0*(FX-FY)/AMBDA + DX + DY
      DALFA = Z*Z - DX*DY
      IF (DALFA.LT.0.0) GO TO 500
      W = SQRT(DALFA)
      ALFA = DY - DX + W + W
      IF (ALFA.EQ.0.0) GO TO 280
      ALFA = (DY-Z+W)/ALFA
      GO TO 290
 280  ALFA = (Z+DY-W)/(Z+DX+Z+DY)
 290  ALFA = ALFA*AMBDA
      DO 300 I=1,N
        X(I) = X(I) + (T-ALFA)*H(I)
 300  CONTINUE
C
C TERMINATE IF THE VALUE OF THE ACTUAL FUNCTION AT X IS LESS
C THAN THE FUNCTION VALUES AT THE INTERVAL ENDS. OTHERWISE REDUCE
C THE INTERVAL BY CHOOSING ONE END POINT EQUAL TO X AND REPEAT
C THE INTERPOLATION. WHICH END POINT IS CHOSEN DEPENDS ON THE
C VALUE OF THE FUNCTION AND ITS GRADIENT AT X.
C
      CALL FUNCT(N, X, F, G)
      DALFA = 0.0
      DO 310 I=1,N
        DALFA = DALFA + H(I)*G(I)
 310  CONTINUE
      IF (DEB) WRITE (IOUTD,9988) DALFA, F
      FYE = 1.0 - FY/F - EPS
      FXE = 1.0 - FX/F - EPS
      IF (DEB) WRITE (IOUTD,9985) FXE, FYE
      IF (FXE.GT.0.0) GO TO 320
      IF (FYE.LE.0.0) GO TO 390
 320  DALFA = 0.0
      DO 330 I=1,N
        DALFA = DALFA + H(I)*G(I)
 330  CONTINUE
      IF (DALFA.GE.0.0) GO TO 360
      IF (FXE) 350, 340, 360
 340  IF (DX.EQ.DALFA) GO TO 390
 350  FX = F
      AMBDA = ALFA
      T = AMBDA
      DX = DALFA
      GO TO 270
 360  IF (FYE) 380, 370, 380
 370  IF (DY.EQ.DALFA) GO TO 390
 380  AMBDA = AMBDA - ALFA
      DY = DALFA
      FY = F
      GO TO 260
C
C IF DIR. DER. AT X NOT << DIR. DER. AT INTERVAL ENDS, REPEAT INTERP
C
 390  DER = ABS(DALFA/DYMN)
      IF (DER.LE.1.0E-2 .OR. ABS(1.0-OF/F).LE.EPS) GO TO 400
      OF = F
      IF (DEB) WRITE (IOUTD,9983) DER, OLDF
      IF (DALFA) 350, 350, 380
 400  CONTINUE
C
C TERMINATE IF FUNCTION HAS NOT BEEN REDUCED DURING LAST ITERATION
C
      IF (OLDF+EPS.LT.F) GO TO 500
C
C COMPUTE DIFFERENCE VECTORS OF ARGUMENT AND GRADIENT FROM
C TWO CONSECUTIVE ITERATIONS
C
      DO 410 J=1,N
        K = N + J
        H(K) = G(J) - H(K)
        K = N + K
        H(K) = X(J) - H(K)
 410  CONTINUE
      IF (DEB) WRITE (IOUTD,9988) DALFA, F
C
C TEST LENGTH OF ARGUMENT DIFFERENCE VECTOR AND DIRECTION VECTOR
C IF AT LEAST N ITERATIONS HAVE BEEN EXECUTED. TERMINATE IF
C BOTH ARE LESS THAN EPS
C
      IER = 0
      IF (KOUNT.LT.N) GO TO 430
      T = 0.0
      DO 420 J=1,N
        K = N2 + J
        T = T + ABS(H(K))
 420  CONTINUE
      IF (HNRM.GT.EPS) GO TO 430
      IF (T.LE.EPS) GO TO 540
C
C TERMINATE IF NUMBER OF ITERATIONS WOULD EXCEED LIMIT
C
 430  IF (KOUNT.GE.LIMIT) GO TO 490
C
C PREPARE UPDATING OF MATRIX H
C
      Z = 0.0
      ALFA = Z
      DO 460 J=1,N
        W = 0.0
        K = J + N3
        DO 450 L=1,N
          KL = N + L
          W = W + H(KL)*H(K)
          IF (L.GE.J) GO TO 440
          K = K + N - L
          GO TO 450
 440      K = K + 1
 450    CONTINUE
        K = N + J
        KN = K + N
        ALFA = ALFA + W*H(K)
        Z = Z + H(K)*H(KN)
        H(J) = W
 460  CONTINUE
      IF (DEB) WRITE (IOUTD,9984) Z, ALFA
C
C REPEAT SEARCH IN DIRECTION OF STEEPEST DESCENT IF RESULTS
C ARE NOT SATISFACTORY
C
      IF (Z*ALFA.EQ.0.0) GO TO 10
C
C UPDATE MATRIX H
C
      K = N31
      DO 480 L=1,N
        KL = N2 + L
        DO 470 J=L,N
          NJ = N2 + J
          H(K) = H(K) + H(KL)*H(NJ)/Z - H(L)*H(J)/ALFA
          K = K + 1
 470    CONTINUE
 480  CONTINUE
      GO TO 40
C
C END OF ITERATION LOOP
C NO CONVERGENCE AFTER LIMIT ITERATIONS
C
 490  IER = 1
      RETURN
C
C RESTORE OLD VALUES OF FUNCTION AND ARGUMENTS
C
 500  DO 510 J=1,N
        K = N2 + J
        X(J) = H(K)
 510  CONTINUE
      CALL FUNCT(N, X, F, G)
C
C REPEAT SEARCH IN DIRECTION OF STEEPEST DESCENT IN DERIVATIVE
C FAILS TO BE SUFFICIENTLY SMALL
C
      IF (GNRM.LE.EPS) GO TO 530
C
C TEST FOR REPEATED FAILURE OF ITERATION
C
 520  IF (IER.LT.0) GO TO 540
      IER = -1
      GO TO 10
 530  IER = 0
 540  RETURN
C
C FORMATS
C
9999  FORMAT (8H OLDF = , E11.4, 10X, 13HITERATION NO., I3)
9998  FORMAT (//(2E30.5))
9997  FORMAT (//31H FUNCTION AND GRADIENT VECTORS-)
9996  FORMAT (3(5H   H(, I3, 4H) = , E13.6))
9995  FORMAT (//9H HNRM   =, E12.5, 10X, 7HGNRM  =, E12.5)
9994  FORMAT (//11H STEP SIZE=, E12.5)
9993  FORMAT (16H NEW ARG. VECTOR/(4E15.7))
9992  FORMAT (//9H DY(XO) =, E12.5, 10X, 7HF(XO) =, E12.5)
9991  FORMAT (//10H DY(XO+H)=, E12.5, 10X, 8HF(XO+H)=, E12.5)
9990  FORMAT (//6H ALFA=, E12.5, 10X, 3HHN=, E12.5)
9989  FORMAT (//12H DY(XO+A H)=, E12.5, 10X, 10HF(XO+A H)=, E12.5)
9988  FORMAT (//7H DYMIN=, E14.7, 10X, 4HFMN=, E14.7)
9987  FORMAT (//4H FX=, E14.7, 10X, 3HFY=, E14.7)
9986  FORMAT (//4H DX=, E14.7, 10X, 3HDY=, E14.7)
9985  FORMAT (//5H FXE=, E14.7, 10X, 4HFYE=, E14.7)
9984  FORMAT (//3H Z=, E14.7, 10X, 5HALFA=, E14.7)
9983  FORMAT (/5H DER=, E14.7, 10X, 6H OLDF=, E14.7)
9982  FORMAT (//38H ************** START STEEPEST DESCENT)
9981  FORMAT (//30X, 20H DIRECTION VECTOR H //)
      END
C
C-----------------------------------------------------------------------
C FUNCTION:  FS
C FUNCTION FOR SPECIFYING THE REQUIRED MAGNITUDE
C-----------------------------------------------------------------------
C
      FUNCTION FS(FI)
      COMMON /CM40/ FIK(10), NPK(10), NK
      COMMON /UNIT/ IND, INP, IOUTD, IOUTP, PI, SPM
C
      DO 60 K=2,NK
        K1 = K - 1
C
        IF (FI.LT.FIK(K1) .OR. FI.GE.FIK(K)) GO TO 60
C
        GO TO (10, 20, 30, 40, 50), K1
  10    FS = 1.0
        GO TO 60
  20    FS = 0.0
        GO TO 60
  30    FS = 0.0
        GO TO 60
  40    FS = 0.0
        GO TO 60
  50    FS = 0.0
  60  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C FUNCTION:  DGS
C FUNCTION FOR SPECIFYING THE REQUIRED GROUP DELAY
C-----------------------------------------------------------------------
C
      FUNCTION DGS(FI)
      COMMON /CM40/ FIK(10), NPK(10), NK
      COMMON /UNIT/ IND, INP, IOUTD, IOUTP, PI, SPM
C
      DO 60 K=2,NK
        K1 = K - 1
C
        IF (FI.LT.FIK(K1) .OR. FI.GE.FIK(K)) GO TO 60
C
        GO TO (10, 20, 30, 40, 50), K1
  10    DGS = 4.0
        GO TO 60
  20    DGS = 4.0
        GO TO 60
  30    DGS = 0.0
        GO TO 60
  40    DGS = 0.0
        GO TO 60
  50    DGS = 0.0
  60  CONTINUE
      RETURN
      END
Return to Main Software Page