An Optimization Program for the Design of Digital Filter Transfer Functions
Return to Main Software Page
C
C-----------------------------------------------------------------------
C MAIN PROGRAM: OPTIMIZATION PROGRAM FOR THE DESIGN OF
C DIGITAL FILTER TRANSFER FUNCTIONS
C AUTHORS: M. T. DOLAN AND J. F. KAISER
C BELL LABORATORIES, MURRAY HILL, NEW JERSEY 07974
C
C INPUT: NORD IS THE ORDER, MUST BE EVEN, .LE. 20
C T IS THE SAMPLING TIME
C XK IS THE GAIN FACTOR
C X IS THE COEFFICIENT VECTOR
C NC IS THE NUMBER OF COEFS TO REMAIN CONSTANT
C NL IS THE ARRAY OF INDICES OF FIXED COEFS
C IF NC=0, DO NOT INPUT NL ARRAY
C NG IS THE GAIN OPTION
C NF IS THE NUMBER OF FREQUENCIES FOR OPTIMIZATION
C FREQ IS THE ARRAY OF FREQUENCIES
C N1 IS THE INDEX FOR TOP RANGE OF LOWER SPECS
C SPEC1 IS THE ARRAY OF LOWER SPECS
C TH1 IS THE ARRAY OF WEIGHTS FOR LOWER SPECS
C N2 IS THE LOWER INDEX FOR RANGE OF UPPER SPECS
C N3 IS THE UPPER INDEX FOR RANGE OF UPPER SPECS
C SPEC2 IS THE ARRAY OF UPPER SPECS
C TH2 IS THE ARRAY OF WEIGHTS FOR UPPER SPECS
C NGID IS THE OPTIMIZATION GUIDE OPTION
C NOUT IS THE OUTPUT OPTION
C NN IS THE INITIAL STEP OPTION
C VBND IS THE PARAMETER USED BY INITIAL STEP OPTIONS
C NPEN IS THE PENALTY FUNCTION CHOICE
C KP DETERMINES STOP CRITERION FOR MINOR ITERATIONS
C EPS IS A PARAMETER USED BY 2 MINOR STOP OPTIONS
C KTIME IS A PARAMETER USED BY 2 MINOR STOP OPTIONS
C NFINR DETERMINES STOP CRITERION FOR MAJOR ITERATIONS
C PERR IS THE PERCENTAGE VALUE USED FOR MAJOR STOP
C OPTION 1 WITH INTERIOR PENALTY FUNCTIONS
C TAZ IS THE SUM OF CONSTRAINTS-SQUARED USED FOR MAJOR
C STOP OPTION 1 WITH EXTERIOR PENALTY FUNCTION
C NR IS THE NUMBER OF R VALUES
C INIR IS THE R GUIDE OPTION - 1 TO GUIDE
C - 0 FOR AUTOMATIC
C RINI IS THE INITIAL R VALUE SET BY ROUTINE ZR IF INIR=0
C RCHNG IS THE R CHANGE FACTOR SET BY MAIN IF INIR=0
C-----------------------------------------------------------------------
C
DIMENSION X(42), S(42), XB(42), H(42,42), G(42), Y(42), GB(42)
DIMENSION SPEC1(42), SPEC2(42), TH1(42), TH2(42), FREQ(42), W(42)
COMMON X, F, DEL, VBND, EPS, KTIME, NIT, NFE, NGE, NQ, K
COMMON /SUBC/ S, XB, GB, FB, H, G, IPANIC
COMMON /FGVARY/ T, XK, FREQ, W, SPEC1, SPEC2, TH1, TH2, R, NX, N,
* NF, N1, N2, N3
COMMON /GRAD/ TAZ, NC, NL(42), NOUT, NVI, NPEN, NFINR
COMMON /RVALUE/ RINI, INIR
COMMON /WRITE/ IW
C
C IR AND IW ARE MACHINE DEPENDENT READ AND WRITE DEVICE NUMBERS.
C THE WRITE DEVICE NUMBER, IW, IS PASSED TO THE SUBROUTINES
C WHICH NEED IT THROUGH LABELLED COMMON.
C
IR = I1MACH(1)
IW = I1MACH(2)
C
C READ ORDER, SAMPLING TIME, AND GAIN FACTOR
C
READ (IR,9999) NORD, T, XK
9999 FORMAT (I3, 2E16.8)
C
N = NORD
NW = 2*N
C
C READ THE COEFFICIENT VECTOR
C
READ (IR,9998) (X(I),I=1,NW)
9998 FORMAT (4E16.8)
C
C READ NUMBER OF COEFS TO REMAIN CONSTANT
C
READ (IR,9997) NC
9997 FORMAT (10I3)
IF (NC.EQ.0) GO TO 10
C
C READ ARRAY OF INDICES OF FIXED COEFS IF NC .NE. 0
C
READ (IR,9997) (NL(I),I=1,NC)
C
C READ GAIN OPTION AND NUMBER OF FREQS FOR OPTIMIZATION
C
10 READ (IR,9997) NG, NF
C
C READ ARRAY OF FREQUENCIES
C
READ (IR,9998) (FREQ(I),I=1,NF)
C
C READ INDEX FOR TOP RANGE OF LOWER SPECS
C
READ (IR,9997) N1
C
C READ ARRAY OF LOWER SPECS
C
READ (IR,9998) (SPEC1(I),I=1,N1)
C
C READ ARRAY OF WEIGHTS FOR LOWER SPECS
C
READ (IR,9998) (TH1(I),I=1,N1)
C
C READ LOWER AND UPPER INDICES FOR RANGE OF UPPER SPECS
C
READ (IR,9997) N2, N3
C
N23 = N3 - N2 + 1
C
C READ ARRAY OF UPPER SPECS
C
READ (IR,9998) (SPEC2(I),I=1,N23)
C
C READ ARRAY OF WEIGHTS FOR UPPER SPECS
C
READ (IR,9998) (TH2(I),I=1,N23)
C
C NGID IS THE OPTIMIZATION GUIDE OPTION. INPUT 0
C FOR AUTOMATIC SELECTION OF PARAMETERS AND NO FURTHER
C INPUT CARDS ARE NECESSARY
C READ OPTIMIZATION GUIDE OPTION
C
READ (IR,9997) NGID
IF (NGID.EQ.0) GO TO 20
C
C READ ADDITIONAL CARDS ONLY FOR NGID=1
C TO SIMPLIFY MATTERS, ALL DATA IS ENTERED.
C USERS SHOULD INPUT 0 FOR PARAMETERS THAT WILL NOT BE USED.
C
C IF NGID=1, READ OUTPUT OPTION, INITIAL STEP OPTION,
C PARAMETER FOR STEP OPTION, PENALTY FUNCTION CHOICE,
C STOP CRITERION FOR MINOR ITERATIONS, TWO PARAMETERS
C FOR STOP OPTION, STOP CRITERION FOR MAJOR ITERATIONS,
C THREE PARAMETERS FOR STOP OPTION, R-GUIDE OPTION,
C INITIAL R VALUE, AND THE R CHANGE FACTOR.
C
READ (IR,9996) NOUT, NN, VBND, NPEN, KP, EPS, KTIME
READ (IR,9995) NFINR, PERR, TAZ, NR, INIR, RINI, RCHNG
GO TO 30
9996 FORMAT (2I3, E12.4, 2I3, E12.4, I3)
9995 FORMAT (I3, 2E12.4, 2I3, 2E12.4)
C
C AUTOMATIC SELECTION FOR NGID=0
C FOR CONSISTENCY, PARAMETERS NOT USED WILL BE SET TO 0
C
20 NOUT = 0
NN = 2
VBND = 0.
NPEN = 1
KP = 2
C
C EPS AND KTIME WILL NOT BE USED
C
EPS = 0.
KTIME = 0
NFINR = 2
C
C PERR AND TAZ WILL NOT BE USED
C
PERR = 0.
TAZ = 0.
NR = 4
INIR = 0
C
C RINI AND RCHNG WILL BE SET BY THE PROGRAM
C
RINI = 0.
RCHNG = 0.
C
30 NX = NW + 2
NVI = 0
IF (INIR.EQ.0) RCHNG = 10.
NXM1 = NX - 1
X(NXM1) = SQRT(ABS(XK))
IF (NG.NE.0) GO TO 40
NC = NC + 1
NL(NC) = NXM1
40 TPI = 8.*ATAN(1.)
NQ = NX
DO 50 I=1,NF
W(I) = TPI*FREQ(I)
50 CONTINUE
XK = X(NXM1)**2
CALL ZR(X)
60 NIT = 0
NFE = 0
NGE = 0
K = KP
OLDXNX = X(NX)
70 CALL FLPWL
IF (K) 100, 100, 80
80 CALL STEP(NN)
CALL QUAD
IF (K) 100, 100, 90
90 CALL FINISH
100 CONTINUE
IF (K) 110, 110, 70
110 WRITE (IW,9994) R
WRITE (IW,9993)
9994 FORMAT (26H0FINAL COEFFICIENTS FOR R=, E16.8)
9993 FORMAT (54H A2 A1 B2 ,
* 4H B1)
WRITE (IW,9998) (X(I),I=1,NW)
XK = X(NXM1)**2
WRITE (IW,9992) XK, X(NX), F
9992 FORMAT (6H GAIN=, E15.8, 7H ALPHA=, E15.8, 10H FUNCTION=, E15.8)
CALL RINV(X, Y, NX)
WRITE (IW,9991)
9991 FORMAT (22H0MODIFIED COEFFICIENTS)
WRITE (IW,9998) (Y(I),I=1,NW)
DXK = Y(NXM1)**2
WRITE (IW,9990) DXK
9990 FORMAT (15H MODIFIED GAIN=, E15.8)
WRITE (IW,9989)
9989 FORMAT (10H0GRADIENTS)
WRITE (IW,9998) (G(I),I=1,NW)
WRITE (IW,9986) NFE, NGE
WRITE (IW,9988)
WRITE (IW,9987)
9988 FORMAT (54H0XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX,
* 18HXXXXXXXXXXXXXXXXXX)
9987 FORMAT (1H0)
9986 FORMAT (1H , I5, 21H FUNCTION EVALUATIONS, I10, 14H GRADIENT EVAL,
* 7HUATIONS)
IF (NFINR.EQ.2) GO TO 120
IF (ABS(X(NX)-OLDXNX)-ABS(PERR*.01*X(NX))) 160, 160, 130
120 NR = NR - 1
IF (NR.EQ.0) STOP
130 GO TO (140, 140, 150), NPEN
140 R = R/RCHNG
GO TO 60
150 R = R*RCHNG
GO TO 60
160 STOP
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: ZR
C CALCULATES X(NX)=ALPHA AND INITIAL R
C-----------------------------------------------------------------------
C
SUBROUTINE ZR(X)
C
DIMENSION X(42), G(42), FF(42), BOUND1(42), BOUND2(42),
* DF(42,42), HM(42), HG(42,42)
DIMENSION SPEC1(42), SPEC2(42), TH1(42), TH2(42), FREQ(42), W(42)
COMMON /FGVARY/ T, XK, FREQ, W, SPEC1, SPEC2, TH1, TH2, R, NX, N,
* NF, N1, N2, N3
COMMON /GRAD/ TAZ, NC, NL(42), NOUT, NVI, NPEN, NFINR
COMMON /RVALUE/ RINI, INIR
COMMON /WRITE/ IW
C
NXM1 = NX - 1
XK = X(NXM1)**2
CALL CASC(N, X, W, NF, T, XK, HM, FF, 0, HG, DF)
DO 10 J=1,N1
IF (TH1(J).EQ.0.) GO TO 10
BOUND1(J) = (FF(J)-SPEC1(J))/TH1(J)
10 CONTINUE
JJ = 0
DO 20 J=N2,N3
JJ = JJ + 1
IF (TH2(JJ).EQ.0.) GO TO 20
BOUND2(JJ) = (-FF(J)+SPEC2(JJ))/TH2(JJ)
20 CONTINUE
DO 30 J=1,N1
IF (BOUND1(J).LT.0.) GO TO 100
30 CONTINUE
JJ = 0
DO 40 J=N2,N3
JJ = JJ + 1
IF (BOUND2(JJ).LT.0.) GO TO 100
40 CONTINUE
C
C ALL POSITIVE. FIND MINIMUM HIT
C
DO 50 J=1,N1
IF (BOUND1(J).EQ.0.) GO TO 50
X(NX) = BOUND1(J)
GO TO 70
50 CONTINUE
JJ = 0
DO 60 J=N2,N3
JJ = JJ + 1
IF (BOUND2(JJ).EQ.0.) GO TO 60
X(NX) = BOUND2(JJ)
GO TO 70
60 CONTINUE
70 DO 80 J=1,N1
IF (BOUND1(J).EQ.0.) GO TO 80
IF (X(NX).LE.BOUND1(J)) GO TO 80
X(NX) = BOUND1(J)
80 CONTINUE
JJ = 0
DO 90 J=N2,N3
JJ = JJ + 1
IF (BOUND2(JJ).EQ.0.) GO TO 90
IF (X(NX).LE.BOUND2(JJ)) GO TO 90
X(NX) = BOUND2(JJ)
90 CONTINUE
X(NX) = -X(NX)*.9999
ZETA = X(NX)
GO TO 150
C
C AT LEAST ONE BOUND IS NEGATIVE. FIND THE MAXIMUM MISS.
C
100 X(NX) = 0.
DO 120 J=1,N1
IF (BOUND1(J)) 110, 110, 120
110 TXNX = ABS(BOUND1(J))
IF (TXNX.LT.X(NX)) GO TO 120
X(NX) = TXNX
120 CONTINUE
JJ = 0
DO 140 J=N2,N3
JJ = JJ + 1
IF (BOUND2(JJ)) 130, 130, 140
130 TXNX = ABS(BOUND2(JJ))
IF (TXNX.LT.X(NX)) GO TO 140
X(NX) = TXNX
140 CONTINUE
X(NX) = X(NX)*1.00001
ZETA = X(NX)
150 AZ = 0.
DO 220 J=1,N1
BOUND1(J) = FF(J) + ZETA*TH1(J) - SPEC1(J)
IF (NPEN.EQ.3) GO TO 190
IF (BOUND1(J)) 320, 320, 160
160 GO TO (170, 180), NPEN
170 BOUND1(J) = 1./BOUND1(J)
AZ = AZ + BOUND1(J)
GO TO 220
180 AZ = AZ - ALOG(BOUND1(J))
GO TO 220
190 IF (BOUND1(J)) 200, 220, 210
200 AZ = AZ + BOUND1(J)**2
GO TO 220
210 BOUND1(J) = 0.
220 CONTINUE
JJ = 0
DO 290 J=N2,N3
JJ = JJ + 1
BOUND2(JJ) = -FF(J) + ZETA*TH2(JJ) + SPEC2(JJ)
IF (NPEN.EQ.3) GO TO 260
IF (BOUND2(JJ)) 320, 320, 230
230 GO TO (240, 250), NPEN
240 BOUND2(JJ) = 1./BOUND2(JJ)
AZ = AZ + BOUND2(JJ)
GO TO 290
250 AZ = AZ - ALOG(BOUND2(JJ))
GO TO 290
260 IF (BOUND2(JJ)) 270, 290, 280
270 AZ = AZ + BOUND2(JJ)**2
GO TO 290
280 BOUND2(JJ) = 0.
290 CONTINUE
IF (INIR.EQ.1) GO TO 310
IF (AZ.EQ.0.) GO TO 300
R = ABS(ZETA/AZ)
IF (R.EQ.0.) R = 1.
RETURN
300 R = 1.
RETURN
310 R = RINI
RETURN
320 WRITE (IW,9999)
9999 FORMAT (54H0TROUBLE - CHECK THAT SPECS ARE MET AT FREQUENCIES WIT,
* 11HH 0 WEIGHTS)
STOP
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: FLPWL
C COMPUTES DIRECTION VECTOR
C-----------------------------------------------------------------------
C
SUBROUTINE FLPWL
C
DIMENSION X(42), S(42), XB(42), H(42,42), G(42)
DIMENSION SIGMA(42), GB(42), Y(42), YTH(42), HY(42)
COMMON X, F, DEL, VBND, EPS, KTIME, NIT, NFE, NGE, NX, K
COMMON /SUBC/ S, XB, GB, FB, H, G, IPANIC
COMMON /WRITE/ IW
C
C FIND GRADIENT USING NEWEST X-ARRAY
C
CALL FG
NGE = NGE + 1
IF (IPANIC) 20, 20, 10
10 WRITE (IW,9999)
9999 FORMAT (24H0 FATAL IPANIC IN FLPWL)
STOP
C
C IF THIS IS THE FIRST ITERATION, SET H = UNIT MATRIX
C
20 IF (NIT) 30, 30, 90
30 DO 70 I=1,NX
DO 60 J=1,NX
IF (I-J) 40, 50, 40
40 H(I,J) = 0.
GO TO 60
50 H(I,J) = 1.
60 CONTINUE
70 CONTINUE
C
C IN THIS CASE, SLOPE = -GRADIENT
C G NOW BECOMES THE OLD G STORE IT IN GB
C
DO 80 I=1,NX
GB(I) = G(I)
S(I) = -G(I)
80 CONTINUE
GO TO 210
C
C IF THIS IS NOT THE FIRST ITERATION, CONTROL REACHES HERE.
C PREPARE TO UPDATE H
C
90 STY = 0.
DO 100 I=1,NX
C
C SIGMA VECTOR REPRESENTS CHANGE IN THE X-ARRAY
C
SIGMA(I) = X(I) - XB(I)
C
C Y-VECTOR REPRESENTS CHANGE IN THE GRADIENT
C
Y(I) = G(I) - GB(I)
C
C ONCE THE NEW G-VECTOR IS USED TO FIND Y, IT BECOMES THE OLD G
C
GB(I) = G(I)
STY = SIGMA(I)*Y(I) + STY
100 CONTINUE
IF (STY) 110, 180, 110
C
C A-MATRIX = (SIGMA-VECTOR TIMES SIGMA-TRANSPOSE) DIVIDED BY
C (SIGMA-TRANSPOSE TIMES Y-VECTOR)
C B-MATRIX = -(H-MATRIX TIMES Y-VECTOR) TIMES
C (Y-TRANSPOSE TIMES H-MATRIX) ALL DIVIDED BY
C (Y-TRANSPOSE TIMES H-MATRIX TIMES Y-VECTOR)
C
110 DO 130 I=1,NX
HY(I) = 0.
YTH(I) = 0.
DO 120 J=1,NX
HY(I) = H(I,J)*Y(J) + HY(I)
YTH(I) = H(J,I)*Y(J) + YTH(I)
120 CONTINUE
130 CONTINUE
C
YTHY = 0.
DO 140 I=1,NX
YTHY = YTH(I)*Y(I) + YTHY
140 CONTINUE
IF (YTHY) 150, 180, 150
C
C UPDATE H
C H-MATRIX = H-MATRIX + A-MATRIX + B-MATRIX
C
150 DO 170 I=1,NX
DO 160 J=1,NX
H(I,J) = -HY(I)*HY(J)/YTHY + H(I,J) + (SIGMA(I)*SIGMA(J))/STY
160 CONTINUE
170 CONTINUE
C
C CALCULATE THE SLOPE VECTOR S-VECTOR=-(H-MATRIX TIMES G-VECTOR)
C
180 DO 200 I=1,NX
S(I) = 0.
DO 190 J=1,NX
S(I) = -H(J,I)*G(J) + S(I)
190 CONTINUE
200 CONTINUE
C
210 RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: STEP
C SETS INITIAL STEP
C-----------------------------------------------------------------------
C
SUBROUTINE STEP(N)
C
DIMENSION X(42), S(42), XB(42), H(42,42), G(42), GB(42)
COMMON X, F, DEL, VBND, EPS, KTIME, NIT, NFE, NGE, NX, K
COMMON /SUBC/ S, XB, GB, FB, H, G, IPANIC
COMMON /GRAD/ TAZ, NC, NL(42), NOUT, NVI, NPEN, NFINR
COMMON /WRITE/ IW
C
IF (N-2) 10, 120, 160
C
C OPTION 1
C USER PROVIDES VMAX OR VMAX=.05*(MAX ABS X)
C STEP=MIN(2.,(VMAX/(MAX ABS S))
C
10 VMAX = VBND
IF (VMAX) 60, 20, 60
20 XMAX = ABS(X(1))
DO 40 I=2,NX
IF (ABS(X(I))-XMAX) 40, 40, 30
30 XMAX = ABS(X(I))
40 CONTINUE
VMAX = .05*XMAX
IF (VMAX) 60, 50, 60
50 VMAX = .01
60 SMAX = ABS(S(1))
DO 80 I=2,NX
IF (ABS(S(I))-SMAX) 80, 80, 70
70 SMAX = ABS(S(I))
80 CONTINUE
IF (SMAX) 90, 90, 110
90 DEL = 2.
100 IF (NOUT.EQ.0) RETURN
WRITE (IW,9999) DEL
9999 FORMAT (14H0INITIAL STEP=, E16.8)
RETURN
110 FAC2 = VMAX/SMAX
DEL = AMIN1(2.,FAC2)
GO TO 100
C
C OPTION 2
C USER PROVIDES ESTIMATE OF FUNCTION MINIMUM CALLED VMIN.
C STEP=MIN(2.,-2.*(F-VMIN))/(G-TRANSPOSE*S-VECTOR)
C
120 VMIN = VBND
DEN = 0.
DO 130 I=1,NX
DEN = G(I)*S(I) + DEN
130 CONTINUE
IF (DEN) 140, 90, 140
140 FAC2 = -2.*(F-VMIN)/DEN
IF (FAC2) 90, 90, 150
150 DEL = AMIN1(2.,FAC2)
GO TO 100
C
C OPTION 3 - IF THIS IS THE FIRST TIME THROUGH, USE OPTION 1.
C IF THIS IS NOT THE FIRST TIME, DIVIDE PRESENT STEP BY 2.
C
C OPTION 4 - DO NOTHING USER PROVIDED STEP
C
160 IF (N-4) 170, 190, 190
170 IF (NIT) 180, 10, 180
180 DEL = DEL/2.
GO TO 100
190 RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: QUAD
C PERFORMS QUADRATIC INTERPOLATION
C-----------------------------------------------------------------------
C
SUBROUTINE QUAD
C
DIMENSION X(42), S(42), XB(42), H(42,42), G(42), GB(42)
DIMENSION X0(42), X1(42), X2(42)
COMMON X, F, DEL, VBND, EPS, KTIME, NIT, NFE, NGE, NX, K
COMMON /SUBC/ S, XB, GB, FB, H, G, IPANIC
COMMON /WRITE/ IW
C
C EVALUATE THE FUNCTION WITH THE GIVEN X-ARRAY AND CALL IT F0.
C AT THE SAME TIME, SAVE THIS ARRAY IN XB AS THE OLD X-ARRAY
C AND SAVE ITS FTN VALUE IN FB AS THE OLD FTN VALUE. LATER
C WE WILL STORE THE NEW ARRAY IN X AND ITS FTN VALUE IN F.
C
CALL FUNC(X, F)
NFE = NFE + 1
DO 10 I=1,NX
XB(I) = X(I)
X0(I) = X(I)
10 CONTINUE
FB = F
F0 = F
C
C GENERATE X1-ARRAY AND F1
C
20 DO 30 I=1,NX
X1(I) = X0(I) + DEL*S(I)
30 CONTINUE
CALL FUNC(X1, F1)
NFE = NFE + 1
IF (IPANIC) 50, 50, 40
40 DEL = DEL/5.
GO TO 20
C
C CHECK IF F1 IS SMALLER THAN F0
C
50 IF (F1-F0) 60, 60, 120
60 DO 70 I=1,NX
X2(I) = X1(I) + DEL*S(I)
70 CONTINUE
CALL FUNC(X2, F2)
NFE = NFE + 1
IF (IPANIC) 90, 90, 80
80 DEL = DEL/5.
GO TO 60
90 IF (F2-F1) 100, 170, 170
C
C IF F2 IS SMALLER THAN F1, DISCARD F0, SHIFT ALL RESULTS LEFT,
C DOUBLE THE INCREMENT, AND GENERATE A NEW X2-ARRAY AND F2
C
100 DO 110 I=1,NX
X0(I) = X1(I)
X1(I) = X2(I)
110 CONTINUE
F0 = F1
F1 = F2
DEL = 2.*DEL
GO TO 60
C
C IF CONTROL REACHES HERE, F1 WAS LARGER AND F0.
C CHECK IF IT WAS MORE THAN 2.*F0
C
120 TF0 = 2.*F0
IF (F1-TF0) 130, 160, 160
130 DO 140 I=1,NX
X2(I) = X1(I)
140 CONTINUE
F2 = F1
DEL = DEL/2.
DO 150 I=1,NX
X1(I) = X0(I) + DEL*S(I)
150 CONTINUE
CALL FUNC(X1, F1)
NFE = NFE + 1
IF (F1-F0) 170, 170, 130
C
C F1 WAS MORE THAN 2.*F0, DIVIDE INCREMENT BY 7 AND
C GENERATE NEW X1-ARRAY AND F2
C
160 DEL = DEL/7.
GO TO 20
C
C INTERPOLATE
170 DO 190 I=1,NX
IF (S(I)) 180, 190, 180
180 B = (X1(I)-X0(I))/S(I)
GO TO 200
190 CONTINUE
GO TO 270
200 C = B + DEL
B2 = B**2
C2 = C**2
FAC1 = .5*(F0*(C2-B2)-C2*F1+B2*F2)
IF (FAC1) 210, 260, 210
210 FAC2 = F0*(C-B) - C*F1 + B*F2
IF (FAC2) 220, 260, 220
220 AL = .5*(F0*(C2-B2)-C2*F1+B2*F2)/(F0*(C-B)-C*F1+B*F2)
DO 230 I=1,NX
X0(I) = X0(I) + AL*S(I)
230 CONTINUE
CALL FUNC(X0, F0)
NFE = NFE + 1
C
C COMPARE RESULT WITH F1 AND RETURN THE BETTER ONE
C
IF (F1-F0) 270, 240, 240
240 DO 250 I=1,NX
X(I) = X0(I)
250 CONTINUE
F = F0
GO TO 290
260 K = 0
WRITE (IW,9999)
9999 FORMAT (31H0UNDERFLOW DURING INTERPOLATION)
270 DO 280 I=1,NX
X(I) = X1(I)
280 CONTINUE
F = F1
290 RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: FUNC
C EVALUATES THE PENALTY FUNCTION
C-----------------------------------------------------------------------
C
SUBROUTINE FUNC(X, F)
C
DIMENSION SPEC1(42), SPEC2(42), TH1(42), TH2(42), FREQ(42), W(42)
DIMENSION HG(42,42), DF(42,42), S(42), XB(42), H(42,42), GB(42)
DIMENSION X(42), G(42), FF(42), BOUND1(42), BOUND2(42), HM(42)
COMMON /FGVARY/ T, XK, FREQ, W, SPEC1, SPEC2, TH1, TH2, R, NX, N,
* NF, N1, N2, N3
COMMON /SUBC/ S, XB, GB, FB, H, G, IPANIC
COMMON /GRAD/ TAZ, NC, NL(42), NOUT, NVI, NPEN, NFINR
C
NXM1 = NX - 1
XK = X(NXM1)**2
CALL CASC(N, X, W, NF, T, XK, HM, FF, 0, HG, DF)
AZ = 0.
IPANIC = 0
ZETA = X(NX)
DO 70 J=1,N1
BOUND1(J) = FF(J) + ZETA*TH1(J) - SPEC1(J)
IF (NPEN.EQ.3) GO TO 40
IF (BOUND1(J)) 150, 150, 10
10 GO TO (20, 30), NPEN
20 BOUND1(J) = 1./BOUND1(J)
AZ = AZ + BOUND1(J)
GO TO 70
30 AZ = AZ - ALOG(BOUND1(J))
GO TO 70
40 IF (BOUND1(J)) 50, 70, 60
50 AZ = AZ + BOUND1(J)**2
GO TO 70
60 BOUND1(J) = 0.
70 CONTINUE
JJ = 0
DO 140 J=N2,N3
JJ = JJ + 1
BOUND2(JJ) = -FF(J) + ZETA*TH2(JJ) + SPEC2(JJ)
IF (NPEN.EQ.3) GO TO 110
IF (BOUND2(JJ)) 150, 150, 80
80 GO TO (90, 100), NPEN
90 BOUND2(JJ) = 1./BOUND2(JJ)
AZ = AZ + BOUND2(JJ)
GO TO 140
100 AZ = AZ - ALOG(BOUND2(JJ))
GO TO 140
110 IF (BOUND2(JJ)) 120, 140, 130
120 AZ = AZ + BOUND2(JJ)**2
GO TO 140
130 BOUND2(JJ) = 0.
140 CONTINUE
F = ZETA + R*AZ
GO TO 160
150 IPANIC = 1
NVI = NVI + 1
160 RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: FG
C EVALUATES BOTH THE PENALTY FUNCTION AND ITS GRADIENT
C-----------------------------------------------------------------------
C
SUBROUTINE FG
C
C EVALUATES THE PENALTY FUNCTION AND ITS GRADIENT FOR A FUNCTION
C WITH VALUES FF(J) AND DERIVATIVE VALUES DF(K,J) AT POINTS J FOR
C ELEMENT NUMBER K. THE ELEMENT VALUES ARE STORED AS X(K), THE
C PENALTY FUNCTION AS F, AND THE GRADIENT AS G(K).
C
DIMENSION X(42), G(42), FF(42), BOUND1(42), BOUND2(42),
* DF(42,42), HM(42), HG(42,42)
DIMENSION SPEC1(42), SPEC2(42), TH1(42), TH2(42), FREQ(42), W(42)
DIMENSION S(42), XB(42), H(42,42), GG(42), GB(42)
COMMON /FGVARY/ T, XK, FREQ, W, SPEC1, SPEC2, TH1, TH2, R, NX, N,
* NF, N1, N2, N3
COMMON /GRAD/ TAZ, NC, NL(42), NOUT, NVI, NPEN, NFINR
COMMON /SUBC/ S, XB, GB, FB, H, G, IPANIC
COMMON X, F, DEL, VBND, EPS, KTIME, NIT, NFE, NGE, NQ, KDUM
COMMON /WRITE/ IW
C
CONS = -20./ALOG(10.)
NXM1 = NX - 1
NXM2 = NX - 2
XK = X(NXM1)**2
CALL CASC(N, X, W, NF, T, XK, HM, FF, 1, HG, DF)
AZ = 0.
ZETA = X(NX)
IPANIC = 0
DO 70 J=1,N1
BOUND1(J) = FF(J) + ZETA*TH1(J) - SPEC1(J)
IF (NPEN.EQ.3) GO TO 40
IF (BOUND1(J)) 340, 340, 10
10 GO TO (20, 30), NPEN
20 BOUND1(J) = 1./BOUND1(J)
AZ = AZ + BOUND1(J)
GO TO 70
30 AZ = AZ - ALOG(BOUND1(J))
GO TO 70
40 IF (BOUND1(J)) 50, 70, 60
50 AZ = AZ + BOUND1(J)**2
GO TO 70
60 BOUND1(J) = 0.
70 CONTINUE
JJ = 0
DO 140 J=N2,N3
JJ = JJ + 1
BOUND2(JJ) = -FF(J) + ZETA*TH2(JJ) + SPEC2(JJ)
IF (NPEN.EQ.3) GO TO 110
IF (BOUND2(JJ)) 340, 340, 80
80 GO TO (90, 100), NPEN
90 BOUND2(JJ) = 1./BOUND2(JJ)
AZ = AZ + BOUND2(JJ)
GO TO 140
100 AZ = AZ - ALOG(BOUND2(JJ))
GO TO 140
110 IF (BOUND2(JJ)) 120, 140, 130
120 AZ = AZ + BOUND2(JJ)**2
GO TO 140
130 BOUND2(JJ) = 0.
140 CONTINUE
F = ZETA + R*AZ
IF ((NPEN.EQ.3) .AND. (NFINR.EQ.1)) GO TO 150
GO TO 160
150 IF (AZ.GT.TAZ) GO TO 160
KDUM = 0
C
C THIS NEXT SECTION DETERMINES THE GRADIENT OF THE PENALTY
C FUNCTION WITH RESPECT TO THE VARIABLES X(K), BUT NOT INCLUDING
C IS X(NX). NOR DOES IT INCLUDE XK WHICH IS X(NX-1)
C
160 NXM2 = NX - 2
GO TO (170, 200, 230), NPEN
170 DO 180 J=1,N1
BOUND1(J) = BOUND1(J)*BOUND1(J)
180 CONTINUE
N4 = N3 - N2 + 1
DO 190 J=1,N4
BOUND2(J) = BOUND2(J)*BOUND2(J)
190 CONTINUE
GO TO 260
200 DO 210 J=1,N1
BOUND1(J) = 1./BOUND1(J)
210 CONTINUE
N4 = N3 - N2 + 1
DO 220 J=1,N4
BOUND2(J) = 1./BOUND2(J)
220 CONTINUE
GO TO 260
C
C SIGN CHANGES BELOW SO CODING AT 2 CAN BE USED
C
230 DO 240 J=1,N1
BOUND1(J) = -2.*BOUND1(J)
240 CONTINUE
N4 = N3 - N2 + 1
DO 250 J=1,N4
BOUND2(J) = -2.*BOUND2(J)
250 CONTINUE
260 DO 290 K=1,NXM2
AZ = 0.0
DO 270 J=1,N1
AZ = AZ - DF(K,J)*BOUND1(J)
270 CONTINUE
JJ = 0
DO 280 J=N2,N3
JJ = JJ + 1
AZ = AZ + DF(K,J)*BOUND2(JJ)
280 CONTINUE
G(K) = R*AZ
290 CONTINUE
C
C THE PARTIAL OF THE FTN WRT XK - SAME FOR RECIP LOG & EXT
C BECAUSE OF WHAT IS CURRENTLY STORED IN BOUNDS FOR EACH
C
AZ = 0.
DO 300 J=1,N1
AZ = AZ - BOUND1(J)
300 CONTINUE
JJ = 0
DO 310 J=N2,N3
JJ = JJ + 1
AZ = AZ + BOUND2(JJ)
310 CONTINUE
G(NXM1) = R*AZ*CONS*2./X(NXM1)
C
C THE GRADIENT OF THE PENALTY FUNCTION WITH RESPECT TO ZETA
C SAME CODING FOR ALL PENALTY FTNS
C
AZ = 0.0
DO 320 J=1,N1
AZ = AZ + TH1(J)*BOUND1(J)
320 CONTINUE
JJ = 0
DO 330 J=N2,N3
JJ = JJ + 1
AZ = AZ + TH2(JJ)*BOUND2(JJ)
330 CONTINUE
G(NX) = 1. - R*AZ
GO TO 350
340 IPANIC = 1
WRITE (IW,9999)
9999 FORMAT (29H0 CONSTRAINTS VIOLATED IN FG)
350 IF (NIT.EQ.0) GO TO 370
IF (NOUT.EQ.0) GO TO 380
WRITE (IW,9998) NVI
9998 FORMAT (23H0 CONSTRAINTS VIOLATED, I5, 22H TIMES DURING QUADRATI,
* 5HC FIT)
NVI = 0
WRITE (IW,9997)
9997 FORMAT (13H0COEFFICIENTS)
360 WRITE (IW,9996)
9996 FORMAT (54H A2 A1 B2 ,
* 4H B1)
WRITE (IW,9995) (X(I),I=1,NXM2)
9995 FORMAT (1X, 4E16.8)
XK = X(NXM1)**2
WRITE (IW,9994) XK, X(NX), F
9994 FORMAT (6H GAIN=, E15.8, 7H ALPHA=, E15.8, 10H FUNCTION=, E15.8)
GO TO 380
370 WRITE (IW,9993) R
9993 FORMAT (28H0INITIAL COEFFICIENTS FOR R=, E16.8)
GO TO 360
380 IF (NC.EQ.0) GO TO 400
DO 390 I=1,NC
IFIX = NL(I)
G(IFIX) = 0.
390 CONTINUE
400 RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: CASC
C CALCULATES THE FREQUENCY RESPONSE AND ITS GRADIENT
C-----------------------------------------------------------------------
C
SUBROUTINE CASC(N, AL, W, NF, T, XK, HM, HL, JJ, HG, HGL)
C
C NOTES FOR USE OUTSIDE THE PACKAGE:
C CASC CAN BE USED TO CALCULATE THE FREQUENCY RESPONSE AND ITS
C GRADIENT FOR ANY NUMBER OF FREQUENCY VALUES, NF.
C THE 42 OF THE CURRENT DIMENSION STATEMENT LIMITS THE TRANSFER
C FUNCTION TO 10 SECTIONS EACH HAVING 4 COEFFICIENTS.
C THE EXTRA 2 LOCATIONS ARE USED BY THE PACKAGE AND ARE NOT
C NECESSARY IF CASC IS USED INDEPENDENTLY. TO USE CASC FOR
C A TRANSFER FUNCTION WITH M SECTIONS WHERE M IS GREATER
C THAN 10, CHANGE THE 42 TO 4*M AND CHANGE THE SIZE OF AH TO M.
C THERE ARE 4 COEFFICIENTS PER SECTION BECAUSE THE CONSTANTS MUST
C BE 1.
C N - INTEGER INPUT - ORDER OF THE FILTER .LE. 20
C AL - REAL INPUT - ARRAY OF COEFFICIENTS OF THE
C TRANSFER FUNCTION. SEE "DESCRIPTION OF
C REQUIRED DATA".
C W - REAL INPUT - ARRAY OF FREQUENCY VALUES MULTIPLIED
C BY 2 PI
C NF - INTEGER INPUT - NUMBER OF FREQUENCY VALUES
C T - REAL INPUT - SAMPLING TIME
C XK - REAL INPUT - GAIN FACTOR
C HM - REAL OUTPUT - ARRAY OF MAGNITUDE VALUES
C HL - REAL OUTPUT - ARRAY OF INSERTION LOSS VALUES
C JJ - INTEGER INPUT OPTION FOR CALCULATION OF THE
C GRADIENT JJ=0 FOR NO GRADIENT
C JJ=1 FOR GRADIENT
C HG - REAL OUTPUT MATRIX - GRADIENT OF THE MAGNITUDE
C WHERE HG(I,J) IS THE DERIVATIVE AT FREQUENCY J
C FOR COEFFICIENT I
C HGL - REAL OUTPUT MATRIX - GRADIENT OF THE INSERTION LOSS
C WHERE HGL(I,J) IS THE DERIVATIVE AT FREQUENCY J
C FOR COEFFICIENT I
C
DIMENSION AL(42), W(NF), HM(NF), HL(NF), HG(42,NF), HGL(42,NF),
* AH(10)
CONS = -20./ALOG(10.)
NU = N/2
C
DO 100 I=1,NF
C
C CALCULATE QUANTITIES CONSTANT WITH FREQUENCY
C
WT = W(I)*T
CWT = COS(WT)
TWT = 2.*WT
CTWT = COS(TWT)
SWT = SIN(WT)
STWT = SIN(TWT)
C
C INITIALIZE MAGNITUDE OF H
C
HM(I) = 1.
C
DO 30 K=1,NU
C
C CALCULATE QUANTITIES CONSTANT WITH ALPHA(K)
C
K4 = 4*K
K4M3 = K4 - 3
K4M2 = K4 - 2
K4M1 = K4 - 1
A = AL(K4M3)
B = AL(K4M2)
C = AL(K4M1)
D = AL(K4)
XM1 = A*CTWT + B*CWT + 1.
XM2 = C*CTWT + D*CWT + 1.
XMU1 = A*STWT + B*SWT
XMU2 = C*STWT + D*SWT
XM12 = XM1*XM2
XMU12 = XMU1*XMU2
XM1MU2 = XM1*XMU2
XM2MU1 = XM2*XMU1
DW = XM2**2 + XMU2**2
C
C SAVE MAGNITUDE OF H(K)
RE = (XM12+XMU12)/DW
XIM = (XM1MU2-XM2MU1)/DW
AH(K) = SQRT(RE**2+XIM**2)
IF (JJ) 20, 20, 10
C
C IF NECESSARY, MAKE PARTIAL CALCULATION OF GRADIENT
C
10 DW2 = DW**2
CST = XM2*CTWT + XMU2*STWT
CS = XM2*CWT + XMU2*SWT
XXR = XM12 + XMU12
XXI = XM1MU2 - XM2MU1
PRA = CST/DW
PRB = CS/DW
PRC = (XM1*CTWT+XMU1*STWT)/DW - 2.*XXR*CST/DW2
PRD = (XM1*CWT+XMU1*SWT)/DW - 2.*XXR*CS/DW2
PIA = (XMU2*CTWT-XM2*STWT)/DW
PIB = (XMU2*CWT-XM2*SWT)/DW
PIC = (XM1*STWT-XMU1*CTWT)/DW - 2.*XXI*CST/DW2
PID = (XM1*SWT-XMU1*CWT)/DW - 2.*XXI*CS/DW2
HG(K4M3,I) = (RE*PRA+XIM*PIA)/AH(K)
HG(K4M2,I) = (RE*PRB+XIM*PIB)/AH(K)
HG(K4M1,I) = (RE*PRC+XIM*PIC)/AH(K)
HG(K4,I) = (RE*PRD+XIM*PID)/AH(K)
C
C MULTIPLY MAGNITUDE OF H BY LATEST MAGNITUDE OF H(K)
C
20 HM(I) = HM(I)*AH(K)
30 CONTINUE
C
HM(I) = HM(I)*XK
C
C FIND INSERTION LOSS
C
HL(I) = -20.*ALOG10(HM(I))
C
IF (JJ) 100, 100, 40
C
C IF NECESSARY, COMPLETE CALCULATION OF GRADIENT
C
40 DO 70 K=1,NU
K4 = 4*K
DO 60 L=1,NU
IF (L-K) 50, 60, 50
50 HG(K4-3,I) = HG(K4-3,I)*AH(L)
HG(K4-2,I) = HG(K4-2,I)*AH(L)
HG(K4-1,I) = HG(K4-1,I)*AH(L)
HG(K4,I) = HG(K4,I)*AH(L)
60 CONTINUE
70 CONTINUE
DO 80 K=1,NU
K4 = 4*K
HG(K4-3,I) = HG(K4-3,I)*XK
HG(K4-2,I) = HG(K4-2,I)*XK
HG(K4-1,I) = HG(K4-1,I)*XK
HG(K4,I) = HG(K4,I)*XK
80 CONTINUE
C
C FIND GRADIENT OF INSERTION LOSS
C
CH = CONS/HM(I)
DO 90 K=1,NU
K4 = 4*K
HGL(K4-3,I) = CH*HG(K4-3,I)
HGL(K4-2,I) = CH*HG(K4-2,I)
HGL(K4-1,I) = CH*HG(K4-1,I)
HGL(K4,I) = CH*HG(K4,I)
90 CONTINUE
C
100 CONTINUE
C
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: RINV
C INVERTS POLES OUTSIDE UNIT CIRCLE
C-----------------------------------------------------------------------
C
SUBROUTINE RINV(X, Y, NX)
C
DIMENSION X(42), Y(42)
C
C SQRT OF CORRECTION ON GAIN IS USED BECAUSE GAIN
C IS SQUARED BY CALLING PROGRAM. OPTIMIZER WORKS
C WITH SQRT OF GAIN TO KEEP IT POSITIVE
NW = NX - 2
NXM1 = NX - 1
DO 10 I=1,NX
Y(I) = X(I)
10 CONTINUE
DO 100 KK=3,NW,4
IF (Y(KK).NE.0.) GO TO 20
GO TO 40
20 DISC = Y(KK+1)**2 - 4.*Y(KK)
IF (DISC.GE.0.) GO TO 50
IF (Y(KK).GT.1.) GO TO 30
GO TO 100
30 SAVE = Y(KK)
Y(KK) = 1./Y(KK)
Y(KK+1) = Y(KK+1)/SAVE
Y(NXM1) = Y(NXM1)/SQRT(ABS(SAVE))
GO TO 100
40 IF (ABS(Y(KK+1)).LE.1.) GO TO 100
C
C FOR REAL ROOTS MULT GAIN BY ROOT BEFORE IT'S FLIPPED
C
SAVE = Y(KK+1)
Y(KK+1) = 1./Y(KK+1)
Y(NXM1) = Y(NXM1)/SQRT(ABS(SAVE))
GO TO 100
50 R1 = (-Y(KK+1)+SQRT(DISC))/(2.*Y(KK))
R2 = (-Y(KK+1)-SQRT(DISC))/(2.*Y(KK))
IF (ABS(R1).LT.1.) GO TO 60
GO TO 70
60 Y(NXM1) = Y(NXM1)*SQRT(ABS(R1))
R1 = 1./R1
70 IF (ABS(R2).LT.1.) GO TO 80
GO TO 90
80 Y(NXM1) = Y(NXM1)*SQRT(ABS(R2))
R2 = 1./R2
C
C RECONSTRUCT THE QUADRATIC
C
90 Y(KK) = 1./(R1*R2)
Y(KK+1) = -(R1+R2)/(R1*R2)
100 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: FINISH
C CHECKS THE STOP CRITERION
C-----------------------------------------------------------------------
C
SUBROUTINE FINISH
C
COMMON X, F, DEL, VBND, EPS, KTIME, NIT, NFE, NGE, NX, N
COMMON /SUBC/ S, XB, GB, FB, H, G, IPANIC
DIMENSION X(42), S(42), XB(42), H(42,42), G(42), XX(42), GB(42)
C
C OPTIONS 1,2,4,5 INITIALIZE THE COUNTER KT IN OPTION 3 IN
C CASE THERE IS SWITCHING TO AND FROM OPTION 3
C
C IF N=0 ON RETURN FROM THIS SUBROUTINE, THEN THE STOP
C CRITERION HAS BEEN SATISFIED. NIT, THE NUMBER OF
C ITERATIONS IS RESET TO ZERO.
C
C N NEGATIVE ON RETURN IS A WARNING THAT OVER 1000 ITERATIONS
C HAVE OCCURRED.
C
C N POSITIVE ON RETURN IS NORMAL AND INDICATES THAT THE
C ITERATIONS ARE TO CONTINUE.
C
GO TO (10, 50, 60, 110, 200), N
C
C OPTION 1
C STOP WHEN GRADIENT-TRANSPOSE TIMES GRADIENT (GTG) IS LESS
C THAN SOME NUMBER GIVEN BY USER OR SET TO GTG/1.E+06
C
10 KT = 0
GTG = 0.
DO 20 I=1,NX
GTG = G(I)*G(I) + GTG
20 CONTINUE
IF (EPS) 30, 30, 40
30 EPS = GTG/1.E+06
40 IF (GTG-EPS) 210, 220, 220
C
C OPTION 2
C STOP WHEN PRESENT FUNCTION VALUE MINUS PREVIOUS FUNCTION
C VALUE IS LESS THAN .0001 TIMES PRESENT VALUE
C
50 KT = 0
IF (ABS(F-FB)-ABS(.0001*F)) 210, 210, 220
C
C OPTION 3
C STOP WHEN OPTION 2 IS SATISFIED A GIVEN NUMBER OF TIMES.
C IF NUMBER IS NOT SPECIFIED IT IS SET TO 4.
C
60 IF (KTIME) 70, 70, 80
70 KTIME = 4
80 IF (ABS(F-FB)-ABS(.0001*F)) 100, 100, 90
90 KT = 0
GO TO 220
100 KT = KT + 1
IF (KT-KTIME) 220, 210, 210
C
C OPTION 4
C STOP WHEN MAX ABS (X-NEW MINUS X-OLD) IS LESS THAN A GIVEN
C NUMBER. IF NOT SPECIFIED, THE NUMBER IS SET TO .0005 MAX X
C
110 KT = 0
DELTA = EPS
IF (DELTA) 120, 120, 160
120 XMAX = ABS(X(1))
DO 140 I=2,NX
IF (ABS(X(I))-XMAX) 140, 140, 130
130 XMAX = ABS(X(I))
140 CONTINUE
DELTA = .0005*XMAX
IF (DELTA) 150, 150, 160
150 DELTA = .0005
160 DO 170 I=1,NX
XX(I) = X(I) - XB(I)
170 CONTINUE
DMAX = ABS(XX(1))
DO 190 I=2,NX
IF (ABS(XX(I))-DMAX) 190, 190, 180
180 DMAX = ABS(XX(I))
190 CONTINUE
IF (DMAX-DELTA) 210, 210, 220
C
C OPTION 5
C STOP AFTER A GIVEN NUMBER OF ITERATIONS
C
200 KT = 0
IF (NIT-KTIME) 220, 210, 210
C
210 N = 0
NIT = 0
RETURN
220 NIT = NIT + 1
IF (NIT-1000) 240, 230, 230
230 N = -1
240 RETURN
END
Return to Main Software Page