Program for the Design of Recursive Digital Filters
Return to Main Software Page
C
C-----------------------------------------------------------------------
C MAIN PROGRAM: D O R E D I
C
C DESIGN AND OPTIMIZATION OF RECURSIVE DIGITAL FILTERS
C
C BUTTERWORTH, CHEBYSHEV, AND ELLIPTIC APPROXIMATION;
C OPTIMIZATION OF THE COEFFICIENT WORDLENGTH AND THE OUTPUT NOISE;
C ANALYSIS OF THE DESIGNED FILTER
C
C VERSION: V005-3 4-15-81
C
C AUTHOR: DR.-ING. GUENTER F. DEHNER
C
C INSTITUT FUER NACHRICHTENTECHNIK
C UNIVERSITAET ERLANGEN-NUERNBERG
C
C CAUERSTRASSE 7
C D-8520 ERLANGEN, GERMANY
C
C INPUT: MANUAL, SECTION 2
C-----------------------------------------------------------------------
C
DOUBLE PRECISION DPI, DOMI
C
COMPLEX ZP, ZPS
C
COMMON /CONTR/ IPRUN, IPCON, NINP, NOUT, NDOUT, LSPOUT, NSPOUT
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
COMMON /CONST2/ DPI, DOMI
COMMON /CPLOT/ IPLOT, LNOR, IPAG, OMLO, OMUP, RMAX, RMIN
C
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /TOLSN/ VSN, VD, A
COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
COMMON /RESZ/ ZFA, ZM(18,4), ZZR(16), ZZI(16), ZPR(16), ZPI(16)
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
COMMON /SFILT/ NBS, SFACT, SB2(16), SB1(16), SB0(16), SC1(16),
* SC0(16)
COMMON /FILTRE/ IRCO(5), IECO(16,5), IDCO(16,5), IECOM
COMMON /CGRID/ GR(64), NGR(12)
COMMON /CRECO/ JRCO, JECO, JDCO, JJDCO(5), JMAXV, JTRB2, LREF
COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
* ISTAB, IDEPSL, IDEPSD
COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
COMMON /COPTST/ LOPTS, ISTOR
COMMON /CRENO/ LCNO, ICNO(16,5)
COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
* LPOT2, JSTRU
COMMON /CREST1/ JSTRUS, JSTRUD
COMMON /CREST2/ KSEQ(16,2)
COMMON /CNOISE/ RI, RIN, RE, REN, FAC
C
COMMON /SCRAT/ ADUM(32)
COMMON /CPOL/ ZP(16,2), ZPS(16,2)
COMMON /COPST2/ PN(100,2), TF(100,2), TFA(100,2)
COMMON /CFFUNC/ PHI(5), BF2(5), BF1(5), BF0(5), IBB(5), ICOR
COMMON /CNFUNC/ AQC(5), BN2(5), BN1(5), BN0(5)
COMMON /CPOW/ PNU, PNC, AND, ITCORP
COMMON /OUTDAT/ IP, PRE, PIM, IZ, ZRE, ZIM
COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
C
CALL DEFIN0
C
NDDEL = -1
CALL HEAD(KA2)
C
IPRUN = 0
IPCON = 0
IPLOT = 0
10 CALL DORINP
C
IF (NDOUT.NE.(-1) .OR. NDDEL.NE.(-1)) GO TO 20
CALL HEAD(KA4)
NDDEL = 0
20 IWLRR = IWLR
GO TO (30, 40, 30, 80), IPRUN
C
30 NGR(8) = 0
CALL SYNTHE
NINP = 1
IF (IWLRR.EQ.100) IWLR = IWL
IF (IPRUN.LT.0) GO TO 50
IF (IPRUN.NE.3) GO TO 90
C
40 IF (NINP.EQ.0) GO TO 100
IF (NINP.GE.3) NGR(8) = 0
C
CALL SEQSTA
C
GO TO 60
50 LOPTS = 0
ISCAL = 0
LSEQ = 0
C
60 CALL SEQOPT
C
IF (IPRUN.EQ.(-3)) GO TO 90
IF (IPRUN.EQ.3) GO TO 70
IF (IWLR.LT.100) CALL SEQSTA
GO TO 90
C
C JSTRU .EQ. 1 SECOND COEFFICIENT OPTIMIZATION
C
70 IF (JSTRU.EQ.2) GO TO 90
IPRUN = -3
GO TO 30
80 IF (NINP.EQ.0) GO TO 100
IF (NINP.GE.3) NGR(8) = 0
C
CALL ANALYS
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C HERE INSERT THE CALLS TO YOUR OWN DESIGN AND ANALYSIS PROCEDURES
C DEFINE NEW LABELS AND EXTEND THE STATEMENT '100+1'
C IPRUN HAS TO BE EQUIVALENT TO THE DEFINITIONS IN 'DORINP'
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
90 IF (LSPOUT.EQ.(-1)) CALL OUTSPE
IF (NINP.NE.2) GO TO 10
ENDFILE KA4
GO TO 10
C
100 CALL ERROR(26)
C
C THE PROGRAM WILL COME TO THE END IN THE SUBROUTINE DORINP
C OR IN THE SUBROUTINE ERROR
C
STOP
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: DEFIN0
C DEFINITION OF THE MACHINE PARAMETERS
C-----------------------------------------------------------------------
C
SUBROUTINE DEFIN0
C
DOUBLE PRECISION DPI, DOMI, D1MACH
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
COMMON /CONST2/ DPI, DOMI
COMMON /CIOFOR/ IOFO(3), ION
C
DIMENSION JOFO(3)
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
DATA JOFO /2H(4,2H0A,2H2)/, JON /40/
C
C THIS STATEMENT IS WRITTEN FOR THE PDP11/45 AND IS
C MACHINE DEPENDENT (SEE MANUAL 4.2.2.)
C E.G.: CDC CYBER 172
C DATA IOFO /10H(8A10) ,2*10H /,ION/8/
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
KA1 = I1MACH(1)
KA2 = I1MACH(2)
KA3 = 6
KA4 = 1
KA5 = 9
LINE = 65
C
PI = 4.*ATAN(1.0)
FLMA = 2.**(I1MACH(13)-2)
FLMI = 2.0*R1MACH(4)
FLER = 1.0E02*FLMI
C
DPI = 4.00*DATAN(1.0D00)
DOMI = 2.0D0*D1MACH(4)
C
MAXDEG = 32
IWLMI = 2
IWLMA = I1MACH(11)
MBL = 16
C
IOFO(1) = JOFO(1)
IOFO(2) = JOFO(2)
IOFO(3) = JOFO(3)
ION = JON
C
RETURN
C
END
C
C ======================================================================
C
C DOREDI - SUBROUTINES: ROOT ELEMENTS FOR OVERLAYS
C
C ======================================================================
C
C
C-----------------------------------------------------------------------
C SUBROUTINE: SYNTHE
C SUBROUTINE TO HANDLE THE NECESSARY COMMONS FOR SYNNOR AND SYNOPT
C FOR OVERLAYS
C-----------------------------------------------------------------------
C
SUBROUTINE SYNTHE
C
DOUBLE PRECISION DK, DKS, DCAP02, DCAP04
COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
COMMON /TOLNOR/ VSNN, NDEGN, NBN
COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
COMMON /RESIN2/ DK, DKS, DCAP02, DCAP04
COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
COMMON /OUTDAT/ IP, PRE, PIM, IZ, ZRE, ZIM
COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
C
CALL SYNNOR
IF (LOPTW.EQ.0) RETURN
CALL SYNOPT
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: SYNNOR
C SYNTHESIS OF RECURSIVE DIGITAL FILTERS
C CALCULATION OF THE ZEROS AND THE POLES
C-----------------------------------------------------------------------
C
SUBROUTINE SYNNOR
C
COMMON /CONTR/ IPRUN, IPCON, NINP, NOUT, NDOUT, LSPOUT, NSPOUT
COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
COMMON /OUTDAT/ IP, PRE, PIM, IZ, ZRE, ZIM
C
MOUT = NOUT
IF (IPRUN.EQ.(-3)) NOUT = 0
CALL DESIA
NOUT = MOUT
IF (LOPTW.NE.0) RETURN
CALL DESIB
IF (NDOUT.EQ.(-1)) CALL OUT012
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: SYNOPT
C SYNTHESIS OF RECURSIVE DIGITAL FILTERS WITH OPTIMIZED COEFFICIENTS
C CALCULATION OF THE POLES AND THE OPTIMIZATION PARAMETER
C-----------------------------------------------------------------------
C
SUBROUTINE SYNOPT
C
COMMON /CONTR/ IPRUN, IPCON, NINP, NOUT, NDOUT, LSPOUT, NSPOUT
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
* ISTAB, IDEPSL, IDEPSD
COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
COMMON /COPTCP/ IS1, IS2, IABO, IWLGS, IWLGN, IWLGP, ADEPSS,
* ADEPSN, ADEPSP, ACXS, ACXN, ACXP
COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
* LPOT2, JSTRU
COMMON /CREST1/ JSTRUS, JSTRUD
COMMON /SFILT/ NBS, SFACT, SB2(16), SB1(16), SB0(16), SC1(16),
* SC0(16)
COMMON /OUTDAT/ IP, PRE, PIM, IZ, ZRE, ZIM
COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
C
FLAB = 10.*FLMI
ACXN = ACXMI
ACXP = ACXMA
ADEPSS = FLMA
ADEPSP = FLMA
ADEPSN = -FLMA
ITER = 0
IWLGS = IWLMA
IWLGP = IWLMA
IWLGN = IWLMA
IWLL = IWLMA
IDEPSL = 0
IS1 = 0
IS2 = 0
IABO = 0
ITERM2 = ITERM
IF (IPRUN.EQ.3 .AND. JSTRU.EQ.1) ITERM = ITERM1
MOUT = NOUT
IF (NOUT.LE.3) NOUT = NOUT - 2
C
10 ITER = ITER + 1
IF (IPRUN.EQ.(-3)) CALL BLNUMZ
IF (NOUT.GE.2) CALL OUT023
CALL DESIB
IF (ITER.GT.1) GO TO 20
CALL GRID01
CALL GRID04
IF (NOUT.EQ.5) CALL OUT042
ACXS = ACX
20 IF (LOPTW.EQ.0) GO TO 60
CALL COEFW
CALL COPY02
CALL OPTPAR
IF (NOUT.GE.2) CALL OUT045
IF (ITER.GE.ITERM) GO TO 30
IF (ACX-ACXMI.LT.FLAB) GO TO 40
IF (ACXMA-ACX.LT.FLAB) GO TO 40
IF (IABO.LT.3) GO TO 10
IABO = 1
GO TO 50
30 IABO = 2
GO TO 50
40 IABO = 3
50 NOUT = MOUT
IF (NOUT.GE.1) CALL OUT022
ACX = ACXS
IF (IPRUN.EQ.3 .AND. JSTRU.EQ.1) GO TO 80
IF (IPRUN.EQ.(-3)) CALL BLNUMZ
CALL DESIB
IF (NOUT.EQ.1) CALL OUT011
IF (NOUT.EQ.2) CALL OUT038
60 IF (NOUT/2.EQ.1) CALL OUT039
IWL = IWLGS
IF (IPRUN.EQ.3 .AND. JSTRUS.EQ.1) GO TO 70
IF (NDOUT.EQ.(-1)) CALL OUT012
70 LOPTW = 1
MOUT = NOUT
IF (NOUT.NE.0) NOUT = NOUT + 2
CALL COEFW
IF (IPRUN.NE.(-3)) CALL COPY02
NOUT = MOUT
80 ITERM = ITERM2
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: SEQSTA
C START NOISE ANALYSIS AND OPTIMIZATION
C-----------------------------------------------------------------------
C
SUBROUTINE SEQSTA
C
COMMON /CONTR/ IPRUN, IPCON, NINP, NOUT, NDOUT, LSPOUT, NSPOUT
COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
* LPOT2, JSTRU
COMMON /CREST1/ JSTRUS, JSTRUD
COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
* ISTAB, IDEPSL, IDEPSD
COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
COMMON /COPTST/ LOPTS, ISTOR
C
IF (LSEQ.NE.(-1)) GO TO 10
IF (NOUT.GE.3) CALL OUT011
CALL COPY01
CALL ALLO02
GO TO 20
C
10 IF (NINP.GE.3 .AND. LOPTS.EQ.1) CALL FIXPAR
20 L = LSTAB
LSTAB = 0
CALL STABLE
IF (ISTAB.NE.(-1)) GO TO 50
LSTAB = L
IF (ISTRU.GE.30) GO TO 30
CALL TSTR01
GO TO 40
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C HERE INSERT TESTS FOR YOUR OWN STRUCTURE IMPLEMENTATIONS
C
30 CALL ERROR(8)
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
40 IF (IWLR.EQ.100) RETURN
IF (IWLR.GT.IWLMA) RETURN
L = IWL
IWL = IWLR
CALL RECO
CALL ROUND(1, 5)
CALL STABLE
IWL = L
RETURN
50 CALL ERROR(30)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: SEQOPT
C OPTIMIZATION OF THE PAIRING AND ORDERING OF THE SECOND-ORDER BLOCKS
C OPTIMIZATION PROCEDURE WITH LIMITED STORAGE
C-----------------------------------------------------------------------
C
SUBROUTINE SEQOPT
C
COMPLEX ZP, ZPS
COMMON /CONTR/ IPRUN, IPCON, NINP, NOUT, NDOUT, LSPOUT, NSPOUT
COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
* LPOT2, JSTRU
COMMON /CREST1/ JSTRUS, JSTRUD
COMMON /COPTST/ LOPTS, ISTOR
COMMON /CGRID/ GR(64), NGR(12)
COMMON /SFILT/ NBS, SFACT, SB2(16), SB1(16), SB0(16), SC1(16),
* SC0(16)
COMMON /CPOL/ ZP(16,2), ZPS(16,2)
COMMON /CNOISE/ RI, RIN, RE, REN, FAC
COMMON /COPST2/ PN(100,2), TF(100,2), TFA(100,2)
COMMON /CFFUNC/ PHI(5), BF2(5), BF1(5), BF0(5), IBB(5), ICOR
COMMON /CNFUNC/ AQC(5), BN2(5), BN1(5), BN0(5)
COMMON /CPOW/ PNU, PNC, AND, ITCORP
C
IF (ISCAL.EQ.0) GO TO 10
CALL DESCAL
10 CALL COPY01
IF (NOUT.GE.3) CALL OUT011
CALL POLLOC
IF (NGR(8).EQ.0) CALL GRID02
IF (NOUT.GE.5) CALL OUT042
IF (NOUT.GE.3) CALL OUT025
IF (LOPTS.NE.0) CALL OPTLST
CALL ANANOI
IF (NOUT.EQ.0) RETURN
CALL OUT027
IF (NOUT.EQ.1) RETURN
IF (JSTRUS.EQ.2) CALL OUT010
CALL OUT011
IF (JSTRUS.EQ.2) CALL DENORM
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: ANALYS
C ANALYSIS OF THE DESIGNED RECURSIVE DIGITAL FILTER
C-----------------------------------------------------------------------
C
SUBROUTINE ANALYS
C
COMPLEX ZP, ZPS
COMMON /CPLOT/ IPLOT, LNOR, IPAG, OMLO, OMUP, RMAX, RMIN
COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
* ISTAB, IDEPSL, IDEPSD
COMMON /COPTST/ LOPTS, ISTOR
COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
* LPOT2, JSTRU
COMMON /CREST1/ JSTRUS, JSTRUD
COMMON /CPOL/ ZP(16,2), ZPS(16,2)
COMMON /CNOISE/ RI, RIN, RE, REN, FAC
C
CALL ANAL01
IF (LOPTS.EQ.0) GO TO 10
CALL ANAL02
10 IF (IPLOT.EQ.0) GO TO 20
CALL PLOT00
C
20 RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: ANAL01
C FIRST PART OF ANALYSIS
C-----------------------------------------------------------------------
C
SUBROUTINE ANAL01
C
COMMON /CONTR/ IPRUN, IPCON, NINP, NOUT, NDOUT, LSPOUT, NSPOUT
COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
* ISTAB, IDEPSL, IDEPSD
COMMON /COPTST/ LOPTS, ISTOR
COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
* LPOT2, JSTRU
COMMON /CREST1/ JSTRUS, JSTRUD
COMMON /CGRID/ GR(64), NGR(12)
COMMON /CPLOT/ IPLOT, LNOR, IPAG, OMLO, OMUP, RMAX, RMIN
COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
C
IF (LSEQ.NE.(-1)) GO TO 10
IF (NOUT.GE.3) CALL OUT011
CALL COPY01
CALL ALLO02
C
10 L = LSTAB
LSTAB = 0
CALL STABLE
IF (ISTAB.NE.(-1)) GO TO 50
LSTAB = L
IF (LOPTW.EQ.0) GO TO 30
IF (NGR(1).NE.0) GO TO 20
IF (NINP.GE.3) CALL GRID03
CALL GRID01
IF (NINP.LE.2) CALL GRID04
20 IF (NOUT.EQ.5) CALL OUT042
CALL COEFW
CALL COPY02
LOPTW = 0
30 IF (LOPTS.NE.(-1) .AND. IPLOT.EQ.0) GO TO 60
IF (IWLR.EQ.100) GO TO 60
IWL = IWLR
IF (IWL.GT.IWLMA) GO TO 40
CALL RECO
CALL ROUND(1, 5)
CALL STABLE
GO TO 60
40 CALL ERROR(16)
GO TO 60
50 CALL ERROR(30)
60 RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: ANAL02
C SECOND PART OF THE ANALYSIS
C-----------------------------------------------------------------------
C
SUBROUTINE ANAL02
C
COMMON /CONTR/ IPRUN, IPCON, NINP, NOUT, NDOUT, LSPOUT, NSPOUT
COMMON /COPTST/ LOPTS, ISTOR
COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
* LPOT2, JSTRU
COMMON /CREST1/ JSTRUS, JSTRUD
COMMON /CNOISE/ RI, RIN, RE, REN, FAC
COMMON /COPST2/ PN(100,2), TF(100,2), TFA(100,2)
COMMON /CFFUNC/ PHI(5), BF2(5), BF1(5), BF0(5), IBB(5), ICOR
COMMON /CNFUNC/ AQC(5), BN2(5), BN1(5), BN0(5)
COMMON /CPOW/ PNU, PNC, AND, ITCORP
C
IF (ISTRU.GE.30) GO TO 10
CALL TSTR01
GO TO 20
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C HERE INSERT TESTS FOR YOUR OWN STRUCTURE IMPLEMENTATIONS
C
10 CALL ERROR(8)
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
20 IF (ISCAL.NE.0) CALL DESCAL
CALL COPY01
IF (NOUT.GE.3) CALL OUT011
LOPTS = 0
CALL POLLOC
CALL GRID02
IF (NOUT.GE.3) CALL OUT025
CALL ANANOI
CALL OUT027
IF (JSTRUS.EQ.2) CALL OUT010
CALL OUT011
IF (JSTRUS.EQ.2) CALL DENORM
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: PLOT00
C PRINTER PLOTS
C-----------------------------------------------------------------------
C
SUBROUTINE PLOT00
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /CPLOT/ IPLOT, LNOR, IPAG, OMLO, OMUP, RMAX, RMIN
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /CSTATV/ X1(16), X2(16)
COMMON /CLINE/ Y, X, IZ(111)
C
DIMENSION A(6)
EQUIVALENCE (A(1),IZ(1))
C
DATA ICHI, ICHP, ICHC, ICHM, ICHB /1HI,1H+,1H.,1H-,1H /
C
C NUMBER OF PRINTED LINES
C
IF (IPLOT.GT.4) GO TO 230
ML = IPAG*LINE - 6
ADOM = (OMUP-OMLO)/FLOAT(ML-1)
C
C SCALING
C
IF (LNOR.NE.(-1)) GO TO 60
C
RRMAX = -FLMA
RRMIN = FLMA
C
IF (IPLOT.EQ.4) Y = RESP(1.,1)
DO 30 L=1,ML
IF (IPLOT.EQ.4) GO TO 20
X = ADOM*FLOAT(L-1) + OMLO
IF (IPLOT.EQ.3) GO TO 10
Y = AMAGO(X)
IF (IPLOT.EQ.1) GO TO 20
IF (Y.LT.FLMI) Y = FLMI
Y = -20.*ALOG10(Y)
IF (Y.GT.RMAX) Y = RMAX
GO TO 20
10 Y = PHASE(X)
20 RRMAX = AMAX1(RRMAX,Y)
RRMIN = AMIN1(RRMIN,Y)
IF (IPLOT.NE.4) GO TO 30
Y = RESP(0.,0)
30 CONTINUE
C
C SCALING OF THE PLOTS
C
RR = ABS(RRMAX-RRMIN)
RN = 0.005*RR
RRMAX = RRMAX - RN
RRMIN = RRMIN + RN
RR = ABS(RRMAX-RRMIN)
IF (RR.LE.0.) RR = ABS(RRMAX)
IQ = INT(ALOG10(RR))
RN = 10.**FLOAT(IQ)
IF (RN.LT.RR) RN = RN*10.
AN = RR/RN
IF (AN.GT.0.0) AAN = 0.1
IF (AN.GT.0.1) AAN = 0.2
IF (AN.GT.0.2) AAN = 0.5
IF (AN.GT.0.5) AAN = 1.0
RR = AAN*RN
C
40 RN = RR*0.2
RMA = AINT(RRMAX/RN)*RN
IF (RMA.LT.RRMAX) RMA = RMA + RN
RMI = RMA - RR
IF (RMI.LE.RRMIN) GO TO 50
RMI = AINT(RRMIN/RN)*RN
IF (RMI.GT.RRMIN) RMI = RMI - RN
RMA = RMI + RR
IF (RMA.GE.RRMAX) GO TO 50
RR = RR*2.
IF (AAN.EQ.0.2) RR = RR*1.25
GO TO 40
C
50 RRMAX = RMA
RRMIN = RMI
GO TO 110
C
60 GO TO (70, 70, 80, 90), IPLOT
70 RRMAX = RMAX
RRMIN = 0.
GO TO 100
80 RRMAX = PI
RRMIN = -RRMAX
GO TO 100
90 RRMAX = RMAX
RRMIN = RMIN
100 RR = RRMAX - RRMIN
C
110 CALL OUT060
NULL = INT(-RRMIN*100./RR+1.5)
IF (NULL.LT.1 .OR. NULL.GT.111) NULL = 0
Q = RR/5.
DO 120 I=1,6
A(I) = Q*FLOAT(I-1) + RRMIN
120 CONTINUE
CALL OUT061
DO 130 I=1,111
IZ(I) = ICHM
130 CONTINUE
IF (IPLOT.EQ.4) Y = RESP(1.,1)
X = 1.
DO 220 L=1,ML
IF (IPLOT.EQ.4) GO TO 150
X = ADOM*FLOAT(L-1) + OMLO
IF (IPLOT.EQ.3) GO TO 140
Y = AMAGO(X)
IF (IPLOT.NE.2) GO TO 150
IF (Y.LT.FLMI) Y = FLMI
Y = -20.*ALOG10(Y)
GO TO 150
140 Y = PHASE(X)
150 DO 160 I=1,101,10
IZ(I) = ICHC
160 CONTINUE
IF (NULL.NE.0) IZ(NULL) = ICHC
Y1 = (Y-RRMIN)/RR
IY = INT(100.*Y1+1.5)
IF (IY.GT.0) GO TO 170
IZ(1) = ICHI
GO TO 190
170 IF (IY.LE.111) GO TO 180
IZ(111) = ICHI
GO TO 190
180 IZ(IY) = ICHP
190 CALL OUT062
IF (IPLOT.NE.4) GO TO 200
Y = RESP(0.,0)
X = FLOAT(L)
200 DO 210 I=1,111
IZ(I) = ICHB
210 CONTINUE
220 CONTINUE
GO TO 240
230 CALL ERROR(8)
240 RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: DORINP
C INPUT DATA COMPILER FOR THE DESIGN PROGRAM 'DOREDI'
C-----------------------------------------------------------------------
C
SUBROUTINE DORINP
C
COMMON /CONTR/ IPRUN, IPCON, NINP, NOUT, NDOUT, LSPOUT, NSPOUT
COMMON /CARD/ LCODE, IZ(40)
COMMON /INPDAT/ ICODE, JCODE, LINP, IINP, AINP(3), NPAR
COMMON /INTDAT/ JINP(16)
COMMON /CPLOT/ IPLOT, LNOR, IPAG, OMLO, OMUP, RMAX, RMIN
COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
COMMON /COPTST/ LOPTS, ISTOR
C
DIMENSION IPT1(5), IPT2(5)
DATA IPT1(1), IPT1(2), IPT1(3), IPT1(4), IPT1(5) /1HS,1HS,1HS,1HA,
* 1HF/
DATA IPT2(1), IPT2(2), IPT2(3), IPT2(4), IPT2(5) /1HY,1HE,1HS,1HN,
* 1HI/
DATA IPMAX /5/
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C HERE INSERT YOUR OWN PROGRAM DEFINITIONS
C CHANGE DIMENSION AND DATA STATEMENT IN FRONT
C 'FI' HAS TO BE LAST DEFINITION
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
LCODE = -1
IF (IPLOT.NE.0) GO TO 50
IF (IPCON.NE.0) GO TO 10
CALL DEFIN1
CALL DEFIN2
CALL DEFIN3
CALL DEFIN4
CALL DEFIN5
CALL DEFIN6
GO TO 60
C
10 IF (IPCON.EQ.99) STOP
IPRUN = IPCON
IF (IPCON.EQ.4) GO TO 40
CALL DEFIN2
IF (IPCON.EQ.2) GO TO 20
CALL DEFIN3
CALL DEFIN4
IF (IPCON.EQ.1) GO TO 60
20 CALL DEFIN5
30 IF (IPRUN.EQ.2 .OR. IPRUN.EQ.3) LOPTS = -1
IF (IPRUN.EQ.3) LOPTW = -1
GO TO 60
40 CALL DEFIN6
GO TO 60
C
50 IPLOT = 0
60 CALL INP001
IF (LCODE.EQ.0) GO TO 110
70 DO 80 IP=1,IPMAX
IF (ICODE.EQ.IPT1(IP) .AND. JCODE.EQ.IPT2(IP)) GO TO 90
80 CONTINUE
GO TO 190
90 IF (IPRUN.EQ.0) GO TO 100
IF (IP.EQ.IPMAX) IP = 99
IPCON = IP
IF (IPRUN.EQ.1 .OR. IPRUN.EQ.3) CALL INP031
RETURN
C
100 IPRUN = IP
GO TO 30
C
110 IF (IPRUN.EQ.0) GO TO 190
CALL INP002
IF (ICODE.EQ.0) GO TO 210
GO TO (120, 130, 140, 150, 160, 170), ICODE
120 CALL INP010
GO TO 180
130 IF (IPRUN.GE.4) GO TO 200
CALL INP020
GO TO 60
140 IF (IPRUN.EQ.2) GO TO 200
CALL INP030
GO TO 60
150 CALL INP040
GO TO 180
160 IF (IPRUN.LE.1) GO TO 200
CALL INP050
GO TO 180
170 IF (IPRUN.NE.4) GO TO 200
CALL INP060
IF (IPLOT.EQ.0) GO TO 60
RETURN
C
180 IF (LCODE.NE.0) GO TO 70
GO TO 60
C
190 I = 3
GO TO 220
200 I = 4
GO TO 220
210 I = 2
220 CALL ERROR(I)
GO TO 60
END
C
C ======================================================================
C
C DOREDI - SUBROUTINES: PART I
C
C ======================================================================
C
C
C-----------------------------------------------------------------------
C SUBROUTINE: DESIA
C FILTER DESIGN -- FIRST SECTION
C-----------------------------------------------------------------------
C
SUBROUTINE DESIA
C
COMMON /CONTR/ IPRUN, IPCON, NINP, NOUT, NDOUT, LSPOUT, NSPOUT
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /TOLSN/ VSN, VD, A
COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
COMMON /TOLNOR/ VSNN, NDEGN, NBN
COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
COMMON /OUTDAT/ IP, PRE, PIM, IZ, ZRE, ZIM
C
SFA = 0.
IF (NOUT.GE.3) CALL OUT030
CALL DESI00
IF (NOUT.GE.4) CALL OUT031
CALL DESI01
IF (NOUT.GE.3) CALL OUT033
GO TO (10, 20, 20, 30), IAPRO
10 CALL DESI11
GO TO 40
20 CALL DESI12
GO TO 40
30 CALL DESI14
C
40 VSNN = VSN
NDEGN = NDEG
NBN = NJ
IF (NOUT.GE.5) CALL OUT034
IF (NOUT.GE.3) CALL OUT035
IF (NOUT.LT.4) GO TO 50
IZ = 0
IP = 0
CALL OUT036
50 CALL TRANZE
IF (NOUT.LT.4) GO TO 60
IZ = 1
IP = 0
CALL OUT036
60 CALL TRBIZE
CALL BLNUMZ
CALL ROMEG
IF (NOUT.GE.3) CALL OUT032
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: DESIB
C FILTER DESIGN -- SECOND SECTION
C-----------------------------------------------------------------------
C
SUBROUTINE DESIB
C
COMMON /CONTR/ IPRUN, IPCON, NINP, NOUT, NDOUT, LSPOUT, NSPOUT
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /TOLSN/ VSN, VD, A
COMMON /TOLNOR/ VSNN, NDEGN, NBN
COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
COMMON /OUTDAT/ IP, PRE, PIM, IZ, ZRE, ZIM
COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
C
VSN = VSNN
NDEG = NDEGN
NJ = NBN
GO TO (10, 20, 20, 30), IAPRO
10 CALL DESI21
GO TO 40
20 CALL DESI22
GO TO 40
30 CALL DESI24
C
40 IF (NOUT.GE.3) CALL OUT037
IF (NOUT.LT.4) GO TO 50
IZ = 0
IP = 1
CALL OUT036
50 CALL TRANPO
IF (LSOUT.EQ.(-1)) GO TO 60
IF (NOUT.LT.4) GO TO 70
60 IZ = 1
IP = 1
CALL OUT036
IF (LSOUT.EQ.(-1)) CALL OUT016
70 CALL TRBIPO
IF (NOUT.GE.3) CALL OUT038
CALL BLDENZ
IF (NOUT.GE.2) CALL OUT011
IF (NOUT.GE.4) CALL OUT039
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: DESI00
C TRANSFORM TOLERANCE SCHEME
C-----------------------------------------------------------------------
C
SUBROUTINE DESI00
C
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
C
IF (ITYP.GE.3) NDEG = (NDEG+1)/2
IF (NDEG.NE.0) ADEG = FLOAT(NDEG)/(1.+EDEG)
IF (LVSN.LT.0) GO TO 20
IF (LVSN.EQ.0) GO TO 10
CALL PARCHA
CALL VSNMAX
10 CALL TRANSN
20 RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: DESI01
C DESIGN OF BUTTERWORTH, CHEBYSHEV (PASSBAND OR STOPBAND), AND
C ELLIPTIC FILTERS
C-----------------------------------------------------------------------
C
SUBROUTINE DESI01
C
COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /TOLSN/ VSN, VD, A
COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
C
IF (LVSN.NE.0) GO TO 20
CALL PARCHA
CALL DEGREE
Q = ADEG*(1.+EDEG) + 0.5
N = INT(Q)
M = INT(ADEG)
IF (FLOAT(M).LT.ADEG) M = M + 1
N = MAX0(M,N)
IF (NDEG.EQ.0) GO TO 10
IF (NDEG.GE.N) GO TO 20
CALL ERROR(15)
IF (NDEGF.EQ.(-1)) STOP
C
10 NDEG = N
C
20 IF (NDEG.LE.MAXDEG) RETURN
CALL ERROR(25)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: PARCHA
C COMPUTATION OF THE PARAMETERS OF THE CHARACTERISTIC FUNCTION
C-----------------------------------------------------------------------
C
SUBROUTINE PARCHA
C
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /TOLSN/ VSN, VD, A
COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
C
GD1 = 0.
GD2 = -1.
IF (ADELP.GT.0.) GD1 = SQRT((2.-ADELP)*ADELP)/(1.-ADELP)
IF (ADELS.GT.0.) GD2 = SQRT(1.-ADELS*ADELS)/ADELS
ACAP12 = GD1/GD2
IF (ACAP12.GT.0.) GO TO 60
GO TO (10, 20, 20, 30), IAPRO
10 ACAP12 = VSN**(-ADEG)
GO TO 40
20 Q = ARCOSH(VSN)*ADEG
ACAP12 = 1./COSH(Q)
GO TO 40
30 CALL BOUND
40 IF (GD2.EQ.(-1.)) GO TO 50
GD1 = ACAP12*GD2
ADELP = 1. - 1./SQRT(1.+GD1*GD1)
GO TO 60
50 GD2 = GD1/ACAP12
ADELS = 1./SQRT(1.+GD2*GD2)
60 RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: VSNMAX
C COMPUTATION OF THE NORMALIZED CUTOFF FREQUENCY FOR A FILTER
C DEGREE AND THE RIPPLES GIVEN
C-----------------------------------------------------------------------
C
SUBROUTINE VSNMAX
C
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /TOLSN/ VSN, VD, A
COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
C
AC = 1./ACAP12
E = 1./ADEG
GO TO (10, 20, 20, 30), IAPRO
10 VSN = AC**E
GO TO 40
20 VSN = COSH(E*ARCOSH(AC))
GO TO 40
30 CALL BOUND
40 RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: BOUND
C CALCULATION OF A BOUND FOR VSN OR ACAP12 FOR ELLIPTIC FILTERS
C-----------------------------------------------------------------------
C
SUBROUTINE BOUND
C
DOUBLE PRECISION DPI, DOMI
DOUBLE PRECISION DE, DKK, DDEG, DCAP12, DEG, DCAP14, DQ, DK1,
* DAB, DMAX, DDE
DOUBLE PRECISION DK, DF, DD
DOUBLE PRECISION DFF, DELLK
C
COMMON /CONST2/ DPI, DOMI
COMMON /TOLSN/ VSN, VD, A
COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
C
DIMENSION DK(3), DF(3)
C
DATA DE /1.D00/
DFF(DD) = (DELLK(DD)*DKK/DELLK(DSQRT(DE-DD*DD)))**II - DDEG
IF (ACAP12.LE.0.) GO TO 10
DCAP12 = DBLE(ACAP12)
DEG = DBLE(1./ADEG)
II = 1
GO TO 20
10 DCAP12 = DE/DBLE(VSN)
DEG = DBLE(ADEG)
II = -1
20 DCAP14 = DSQRT(DE-DCAP12*DCAP12)
DKK = DELLK(DCAP14)/DELLK(DCAP12)
DQ = DEXP(-DPI*DKK*DEG)
DK1 = 4.D00*DSQRT(DQ)
IF (DK1.LT.DE) GO TO 30
DQ = 2.D00*DQ
DQ = (DE-DQ)/(DE+DQ)
DQ = DQ*DQ
DK1 = DSQRT(DE-DQ*DQ)
30 DK(1) = DK1
DK(2) = (DE+DK(1))/2.D00
DDEG = DBLE(ADEG)
DF(1) = DFF(DK(1))
DF(2) = DFF(DK(2))
40 DK(3) = DK(1) - DF(1)*(DK(1)-DK(2))/(DF(1)-DF(2))
DF(3) = DFF(DK(3))
IF (DABS(DF(3)).LT.1.D-6) GO TO 60
DMAX = 0.D00
DO 50 J=1,3
DAB = DABS(DF(J))
IF (DMAX.GT.DAB) GO TO 50
JJ = J
DMAX = DAB
50 CONTINUE
IF (JJ.EQ.3) GO TO 40
DK(JJ) = DK(3)
DF(JJ) = DF(3)
GO TO 40
60 IF (ACAP12.LE.0.) GO TO 70
DDE = DE/DK(3)
VSN = SNGL(DDE)
RETURN
70 ACAP12 = SNGL(DK(3))
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: DEGREE
C COMPUTATION OF THE MINIMUM FILTER DEGREE (ADEG)
C-----------------------------------------------------------------------
C
SUBROUTINE DEGREE
C
DOUBLE PRECISION DE, DCAP02, DCAP04, DCAP12, DCAP14, DADEG
DOUBLE PRECISION DELLK
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /TOLSN/ VSN, VD, A
COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
C
GO TO (10, 20, 20, 30), IAPRO
C
10 ADEG = ALOG(1./ACAP12)/ALOG(VSN)
RETURN
C
20 ADEG = ARCOSH(1./ACAP12)/ARCOSH(VSN)
RETURN
C
30 DE = 1.D00
DCAP02 = DE/DBLE(VSN)
DCAP04 = DSQRT(DE-DCAP02*DCAP02)
DCAP12 = DBLE(ACAP12)
DCAP14 = DSQRT(DE-DCAP12*DCAP12)
DADEG = (DELLK(DCAP02)*DELLK(DCAP14))/(DELLK(DCAP04)*DELLK(DCAP12)
* )
ADEG = SNGL(DADEG)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: DESI11
C BUTTERWORTH FILTER
C COMPUTATION OF THE ZEROS AND LOCATIONS OF THE EXTREMA
C-----------------------------------------------------------------------
C
SUBROUTINE DESI11
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /TOLSN/ VSN, VD, A
COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
C
ADELTA = VSN**NDEG
C
NH = NDEG/2
NJ = (NDEG+1)/2
FDEG = FLOAT(NDEG)
FN = PI/2./FDEG
C
DO 10 I=1,NJ
NZERO(I) = 0
III = I + I - 1
Q = FN*FLOAT(III)
PREN(I) = SIN(Q)
PIMN(I) = COS(Q)
10 CONTINUE
C
FN = 2.*FN
NZERO(1) = NDEG
NZM(1) = 1
SM(1,1) = 0.
NZM(2) = 1
SM(1,2) = 1.
NZM(3) = 1
SM(1,3) = VSN
NZM(4) = 1
SM(1,4) = FLMA
C
UGC = GD2/ADELTA
OGC = GD1
SM(17,4) = 1.
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: DESI12
C CHEBYSHEV FILTER (PASSBAND OR STOPBAND)
C COMPUTATION OF THE ZEROS AND THE LOCATIONS OF THE EXTREMA
C-----------------------------------------------------------------------
C
SUBROUTINE DESI12
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /TOLSN/ VSN, VD, A
COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
C
ADELTA = COSH(FLOAT(NDEG)*ARCOSH(VSN))
C
FA = 1.
NH = NDEG/2
NJ = (NDEG+1)/2
FN = PI/(2.*FLOAT(NDEG))
C
DO 10 I=1,NJ
NZERO(I) = 0
INJ = I + I - 1
Q = FN*FLOAT(INJ)
PREN(I) = SIN(Q)
PIMN(I) = COS(Q)
10 CONTINUE
C
FN = 2.*FN
C
IF (IAPRO.EQ.3) GO TO 40
C
M = NJ + 1
DO 20 I=1,NJ
J = M - I
SM(I,1) = PIMN(J)
20 CONTINUE
NZM(1) = NJ
M = NH + 1
DO 30 I=1,M
MI = M - I
Q = FLOAT(MI)*FN
SM(I,2) = COS(Q)
30 CONTINUE
NZM(2) = M
SM(1,3) = VSN
NZM(3) = 1
SM(1,4) = FLMA
NZM(4) = 1
NZERO(1) = NDEG
C
UGC = GD2/ADELTA
OGC = GD1
GO TO 80
C
40 SM(1,1) = 0.
NZM(1) = 1
SM(1,2) = 1.
NZM(2) = 1
DO 50 I=1,NJ
INJ = NJ - I
Q = FLOAT(INJ)*FN
SM(INJ+1,3) = VSN/COS(Q)
50 CONTINUE
NZM(3) = NJ
DO 60 I=1,NH
NZERO(I) = 2
Q = PIMN(I)
FA = FA*Q*Q
SM(I,4) = VSN/Q
60 CONTINUE
IF (NH.EQ.NJ) GO TO 70
NZERO(NJ) = 1
SM(NJ,4) = FLMA
70 NZM(4) = NJ
C
UGC = GD2
OGC = GD1*ADELTA
80 ACK = FA
SM(17,4) = 1.
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: DESI14
C ELLIPTIC FILTER
C COMPUTATION OF THE ZEROS AND THE LOCATIONS OF THE EXTREMA
C-----------------------------------------------------------------------
C
SUBROUTINE DESI14
C
DOUBLE PRECISION DPI, DOMI
DOUBLE PRECISION DK, DKS, DCAP02, DCAP04, DE
DOUBLE PRECISION DSK(16)
DOUBLE PRECISION DCAP01, DQ, DN, DU, DM
DOUBLE PRECISION DELLK, DSN2, DEL1, DEL2, DDE, DDELTA, DDELT
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /CONST2/ DPI, DOMI
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /TOLSN/ VSN, VD, A
COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
COMMON /RESIN2/ DK, DKS, DCAP02, DCAP04
EQUIVALENCE (PREN(1),DSK(1))
C
DATA DE /1.D00/
C
DCAP02 = DE/DBLE(VSN)
DCAP01 = DSQRT(DCAP02)
DCAP04 = DSQRT(DE-DCAP02*DCAP02)
DK = DELLK(DCAP02)
DKS = DELLK(DCAP04)
C
DQ = DEXP(-DPI*DKS/DK)
C
NH = NDEG/2
M = NDEG + 1
NJ = M/2
MH = NH + 1
C
DN = DK/DBLE(FLOAT(NDEG))
C
DEL1 = DE
IF (NH.EQ.0) GO TO 20
DO 10 I=1,NH
INH = M - I - I
DU = DN*DBLE(FLOAT(INH))
DM = DSN2(DU,DK,DQ)
DEL1 = DEL1*DM*DCAP01
DSK(I) = DM
J = MH - I
SM(J,1) = SNGL(DM)
NZERO(I) = 2
DDE = DE/(DCAP02*DM)
SM(I,4) = SNGL(DDE)
10 CONTINUE
GO TO 30
20 SM(1,1) = 0.
30 KJ = NJ - 1
MJ = NJ + 1
DEL2 = DE
IF (KJ.EQ.0) GO TO 50
DO 40 I=1,KJ
NDEGI = NDEG - I - I
DU = DN*DBLE(FLOAT(NDEGI))
DM = DSN2(DU,DK,DQ)
J = MJ - I
SM(J,2) = SNGL(DM)
DDE = DE/(DCAP02*DM)
SM(I+1,3) = SNGL(DDE)
DEL2 = DEL2*DM*DCAP01
40 CONTINUE
GO TO 60
50 SM(NDEG,2) = 1.
SM(1,3) = VSN
C
60 DDELT = DEL1*DEL1
ADELTA = SNGL(DDELT)
ACK = 1./ADELTA
IF (NH.EQ.NJ) GO TO 80
ACK = ACK*SNGL(DCAP01)
DDELTA = DEL2*DEL2*DCAP01
ADELTA = SNGL(DDELTA)
DSK(NJ) = 0.D00
NZERO(NJ) = 1
SM(NJ,4) = FLMA
C
IF (NH.EQ.0) GO TO 90
DO 70 I=1,NH
J = MJ - I
SM(J,1) = SM(J-1,1)
SM(I,2) = SM(I+1,2)
70 CONTINUE
SM(1,1) = 0.
GO TO 90
C
80 SM(MH,3) = FLMA
SM(1,2) = 0.
C
90 NZM(1) = NJ
NZM(4) = NJ
NZM(2) = MH
NZM(3) = MH
SM(MH,2) = 1.
SM(1,3) = VSN
UGC = GD2*ADELTA
OGC = GD1/ADELTA
SM(17,4) = 1.
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: DESI21
C BUTTERWORTH FILTER
C COMPUTATION OF THE POLES
C-----------------------------------------------------------------------
C
SUBROUTINE DESI21
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
C
C COMPUTATION OF CONSTANT C AND REDUCED TOLERANCE SCHEME
C
IF (ACX.LT.999.) GO TO 20
IF ((OGC-UGC).LT.FLMI) GO TO 10
AC = (2.*ADELP/(ADELTA*ADELS))**0.33333
ACX = ALOG10(AC/UGC)/ALOG10(OGC/UGC)
IF (ACX.GE.0. .AND. ACX.LE.1.) GO TO 30
10 ACX = 0.5
20 AC = UGC*(OGC/UGC)**ACX
30 RDELP = 1. - SQRT(1./(1.+AC*AC))
Q = AC*ADELTA
RDELS = SQRT(1./(1.+Q*Q))
C
C COMPUTATION OF FACTOR SFA AND POLES
C
SFA = 1./AC
Q = AC**(-1./FLOAT(NDEG))
C
DO 40 I=1,NJ
SPR(I) = -Q*PREN(I)
SPI(I) = Q*PIMN(I)
40 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: DESI22
C CHEBYSHEV FILTER (PASSBAND OR STOPBAND)
C COMPUTATION OF THE POLES
C-----------------------------------------------------------------------
C
SUBROUTINE DESI22
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /TOLSN/ VSN, VD, A
COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
C
IF (ACX.LT.999.) GO TO 20
IF ((OGC-UGC).LT.FLMI) GO TO 10
IF (IAPRO.EQ.2) Q = 1./ADELTA
IF (IAPRO.EQ.3) Q = ADELTA*ADELTA
AC = (2.*ADELP*Q/ADELS)**0.33333
ACX = ALOG10(AC/UGC)/ALOG10(OGC/UGC)
IF (ACX.GE.0. .AND. ACX.LE.1.) GO TO 30
10 ACX = 0.5
20 AC = UGC*(OGC/UGC)**ACX
C
C COMPUTATION OF THE REDUCED TOLERANCE SCHEME
C
30 Q = AC
IF (IAPRO.EQ.3) Q = Q/ADELTA
Q = 1. + Q*Q
RDELP = 1. - SQRT(1./Q)
Q = AC
IF (IAPRO.EQ.2) Q = Q*ADELTA
Q = 1. + Q*Q
RDELS = SQRT(1./Q)
C
C COMPUTATION OF THE FACTOR SFA AND THE POLES
C
IF (IAPRO.EQ.3) GO TO 40
SFA = 2./(AC*2.**NDEG)
Q = -1./AC
GO TO 50
C
40 SFA = ACK
Q = AC
50 Q = ARSINH(Q)/FLOAT(NDEG)
QR = SINH(Q)
QI = COSH(Q)
IF (IAPRO.EQ.3) GO TO 70
DO 60 I=1,NJ
SPR(I) = QR*PREN(I)
SPI(I) = QI*PIMN(I)
60 CONTINUE
RETURN
C
70 DO 80 I=1,NH
Q = PIMN(I)*QI
QA = PREN(I)*QR
QQ = Q*Q
QQA = QA*QA
SFA = SFA/(QQ+QQA)
SPR(I) = -VSN/(QQ/QA+QA)
SPI(I) = VSN/(Q+QQA/Q)
80 CONTINUE
IF (NH.EQ.NJ) RETURN
SPI(NJ) = 0.
Q = VSN/QR
SFA = SFA*Q
SPR(NJ) = -Q
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: DESI24
C ELLIPTIC FILTER
C COMPUTATION OF THE REDUCED TOLERANCE SCHEME, THE FACTOR SFA, AND
C THE POLES
C-----------------------------------------------------------------------
C
SUBROUTINE DESI24
C
DOUBLE PRECISION DPI, DOMI
DOUBLE PRECISION DK, DKS, DCAP02, DCAP04, DE
DOUBLE PRECISION DSK(16)
DOUBLE PRECISION DU, DR, DQ, DUD, DUC, DRC, DRD, DM
DOUBLE PRECISION DELLK, DSN2
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /CONST2/ DPI, DOMI
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /TOLSN/ VSN, VD, A
COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
COMMON /RESIN2/ DK, DKS, DCAP02, DCAP04
EQUIVALENCE (PREN(1),DSK(1))
C
DATA DE /1.D00/
C
C IF ACX NOT DEFINED, COMPUTE A SYMMETRICAL USAGE OF THE TOLERANCE
C SCHEME
C
IF (ACX.LT.999.) GO TO 20
IF ((OGC-UGC).LT.FLMI) GO TO 10
AC = (2.*ADELP/(ADELTA*ADELS))**0.333333333
ACX = ALOG10(AC/UGC)/ALOG10(OGC/UGC)
IF (ACX.GE.0. .AND. ACX.LE.1.) GO TO 30
10 ACX = 0.5
C
20 AC = UGC*(OGC/UGC)**ACX
C
C COMPUTATION OF THE REDUCED TOLERANCE SCHEME
C
30 Q = AC*ADELTA
DU = DE/DBLE(Q)
RDELP = 1. - SQRT(1./(1.+Q*Q))
Q = 1. + AC*AC/(ADELTA*ADELTA)
RDELS = SQRT(1./Q)
C
C COMPUTATION OF THE FACTOR SFA AND THE POLES
C
Q = AC*ACK
IF (NH.EQ.NJ) Q = SQRT(1.+Q*Q)
SFA = 1./Q
C
DR = DBLE(ADELTA)
DR = DR*DR
DQ = DBLE(Q)
CALL DELI1(DQ, DU, DR)
DU = DQ
DQ = DSQRT(DE-DR*DR)
DQ = DELLK(DR)
DU = DK*DU/(DQ*DBLE(FLOAT(NDEG)))
DQ = DEXP(-DPI*DK/DKS)
DU = -DSN2(DU,DKS,DQ)
DQ = DU*DU
DUD = DE - DCAP04*DCAP04*DQ
DUD = DSQRT(DUD)
DUC = DSQRT(DE-DQ)
DO 40 I=1,NJ
DR = DSK(I)
DRC = DR*DR
DRD = DE - DCAP02*DCAP02*DRC
DRC = DSQRT(DE-DRC)
DM = DE - DQ*DRD
DRD = DSQRT(DRD)
DRD = DRD*DU*DUC*DRC/DM
SPR(I) = SNGL(DRD)
DR = DR*DUD/DM
SPI(I) = SNGL(DR)
40 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: TRANZE
C REACTANCE TRANSFORMATION OF THE ZEROS AND THE LOCATIONS OF THE
C EXTREMA
C-----------------------------------------------------------------------
C
SUBROUTINE TRANZE
C
DOUBLE PRECISION DR, DQI
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /TOLSN/ VSN, VD, A
COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
C
DIMENSION MSM(4)
FA = 1.
IF (ITYP.EQ.1) GO TO 190
IF (ITYP.EQ.3) GO TO 60
C
ME = NZM(4)
DO 10 I=1,ME
Q = SM(I,4)
IF (Q.LT.FLMA) FA = FA*Q
10 CONTINUE
C
FA = FA*FA
C
C LOWPASS - HIGHPASS
C
DO 50 J=1,4
ME = NZM(J)
DO 40 I=1,ME
QI = SM(I,J)
IF (ABS(QI).LT.FLMI) GO TO 20
QI = 1./QI
GO TO 30
20 QI = FLMA
30 SM(I,J) = QI
40 CONTINUE
50 CONTINUE
GO TO 90
60 DO 80 J=1,2
ME = NZM(J)
MA = ME + 1
ME = ME/2
DO 70 I=1,ME
QI = SM(I,J)
II = MA - I
SM(I,J) = SM(II,J)
SM(II,J) = QI
70 CONTINUE
80 CONTINUE
C
90 IF (ITYP.EQ.2) GO TO 190
C
C LOWPASS - BANDPASS TRANSFORMATION
C
QA = 2.*A
NN = NDEG + 1
IF (ITYP.EQ.4) GO TO 110
C
MSM(1) = 1
IF (NZM(1).NE.1) MSM(1) = NDEG
MSM(2) = 2
IF (NZM(2).NE.1) MSM(2) = NN
DO 100 J=3,4
MSM(J) = 2*NZM(J)
100 CONTINUE
GO TO 130
C
110 DO 120 J=1,2
MSM(J) = 2*NZM(J)
120 CONTINUE
MSM(3) = 2
IF (NZM(3).NE.1) MSM(3) = NN
MSM(4) = 1
IF (NZM(4).NE.1) MSM(4) = NDEG
C
130 S = 1.
DO 180 J=1,4
ME = NZM(J)
MA = MSM(J)
NZM(J) = MA
IF (J.EQ.3) S = -1.
DO 170 I=1,ME
QR = SM(I,J)
NU = NZERO(I)
IF (ABS(QR).LT.FLMA) GO TO 150
IF (J.NE.4) GO TO 140
FA = FA*(VD/A)**NU
140 QI = QR
GO TO 160
C
150 QR = QR/QA
DR = DBLE(QR)
DQI = DSQRT(DR*DR+1.D00)
QI = SNGL(DQI)
160 SM(I,J) = QI - S*QR
II = MA - I + 1
IF (ABS(QR).LT.FLMI) NU = 2*NU
IF (J.EQ.4) NZERO(II) = NU
SM(II,J) = QI + S*QR
170 CONTINUE
180 CONTINUE
C
190 DO 220 J=1,4
ME = NZM(J)
DO 210 I=1,ME
Q = SM(I,J)
IF (Q.LT.FLMA) GO TO 200
IF (J.NE.4 .OR. ITYP.GE.3) GO TO 210
NU = NZERO(I)
FA = FA*VD**NU
GO TO 210
200 SM(I,J) = Q*VD
210 CONTINUE
220 CONTINUE
SM(17,4) = SM(17,4)*FA
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: TRANPO
C REACTANCE TRANSFORMATION OF THE POLES
C-----------------------------------------------------------------------
C
SUBROUTINE TRANPO
C
DOUBLE PRECISION DR, DI, DQ
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /TOLSN/ VSN, VD, A
COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
C
IF (ITYP.EQ.1) GO TO 90
IF (ITYP.EQ.3) GO TO 40
C
DO 30 I=1,NJ
QR = SPR(I)
QI = SPI(I)
QH = QR*QR + QI*QI
IF (ABS(QI).GT.FLMI) GO TO 10
SFA = -SFA/QR
GO TO 20
10 SFA = SFA/QH
20 QI = QI/QH
IF (ABS(QI).LT.FLMI) QI = 0.
SPI(I) = QI
SPR(I) = QR/QH
30 CONTINUE
IF (ITYP.EQ.2) GO TO 90
C
40 QA = 2.*A
NN = NJ
NJ = NDEG
NDEG = 2*NDEG
C
ME = NJ
DO 80 I=1,NN
QR = SPR(I)/QA
QI = SPI(I)/QA
IF (ABS(QI).GE.FLMA) GO TO 70
DR = DBLE(QR)
DI = DBLE(QI)
DQ = DI*DI
DI = DI*DR*2.D00
DR = DR*DR - DQ - 1.D00
CALL DSQRTC(DR, DI, DR, DI)
QZ = SNGL(DR)
QN = SNGL(DI)
IF (ABS(QN).GT.FLMI) GO TO 60
JJ = NJ + ME
DO 50 II=ME,NJ
J = JJ - II
SPR(J+1) = SPR(J)
SPI(J+1) = SPI(J)
50 CONTINUE
NJ = NJ + 1
ME = ME + 1
60 SPR(I) = QR + QZ
SPI(I) = QI + QN
SPR(ME) = QR - QZ
SPI(ME) = QN - QI
ME = ME - 1
GO TO 80
C
70 SPR(I) = QR
SPI(I) = FLMA
NJ = NJ + 1
SPR(NJ) = QR
SPI(NJ) = 0.
80 CONTINUE
C
90 DO 100 I=1,NJ
SPR(I) = SPR(I)*VD
SPI(I) = SPI(I)*VD
100 CONTINUE
C
SFA = SFA*SM(17,4)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: TRBIZE
C BILINEAR TRANSFORMATION OF THE ZEROS AND THE LOCATIONS OF THE EXTREMA
C-----------------------------------------------------------------------
C
SUBROUTINE TRBIZE
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
COMMON /RESZ/ ZFA, ZM(18,4), ZZR(16), ZZI(16), ZPR(16), ZPI(16)
C
FA = 1.
DO 50 J=1,4
ME = NZM(J)
DO 40 I=1,ME
QI = SM(I,J)
ZM(I,J) = 2.*ATAN(QI)
IF (J.NE.4) GO TO 40
IF (QI.GE.FLMA) GO TO 10
IF (QI.LT.FLMI) GO TO 20
QQI = QI*QI
Q = 1. + QQI
ZZR(I) = (1.-QQI)/Q
ZZI(I) = 2.*QI/Q
NU = NZERO(I)/2
FA = FA*Q**NU
GO TO 40
C
10 ZZR(I) = -1.
GO TO 30
C
20 ZZR(I) = 1.
30 ZZI(I) = 0.
40 CONTINUE
50 CONTINUE
C
SM(17,1) = FA
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: TRBIPO
C BILINEAR TRANSFORMATION OF THE POLES
C-----------------------------------------------------------------------
C
SUBROUTINE TRBIPO
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
COMMON /RESZ/ ZFA, ZM(18,4), ZZR(16), ZZI(16), ZPR(16), ZPI(16)
COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
C
ZFA = SFA*SM(17,1)
C
DO 20 I=1,NJ
QR = SPR(I)
Q = 1. - QR
QI = SPI(I)
IF (ABS(QI).LT.FLMI) GO TO 10
QQR = QR*QR
QQI = QI*QI
ZFA = ZFA/(Q-QR+QQR+QQI)
Q = 1./(Q*Q+QQI)
ZPR(I) = (1.-QQR-QQI)*Q
ZPI(I) = 2.*QI*Q
GO TO 20
C
10 ZPR(I) = (1.+QR)/Q
ZPI(I) = 0.
ZFA = ZFA/Q
20 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: BLNUMZ
C BUILD NUMERATOR BLOCKS OF SECOND ORDER
C-----------------------------------------------------------------------
C
SUBROUTINE BLNUMZ
C
COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
COMMON /RESZ/ ZFA, ZM(18,4), ZZR(16), ZZI(16), ZPR(16), ZPI(16)
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
C
DIMENSION NZE(16)
COMMON /SCRAT/ ADUM(32)
EQUIVALENCE (ADUM(1),NZE(1))
C
N = 0
ME = NZM(4)
DO 10 I=1,ME
NZE(I) = NZERO(I)
10 CONTINUE
C
DO 70 I=1,ME
QR = ZZR(I)
NZ = NZE(I)
C
20 IF (NZ.EQ.0) GO TO 70
N = N + 1
B2(N) = 1.
IF (NZ.EQ.1) GO TO 30
B1(N) = -2.*QR
B0(N) = 1.
NZ = NZ - 2
IF (NZ.GT.0) GO TO 20
GO TO 70
C
30 IF (I.EQ.ME) GO TO 60
MA = I + 1
DO 40 II=MA,ME
IF (ZZI(II).EQ.0.) GO TO 50
40 CONTINUE
GO TO 60
C
50 QRR = ZZR(II)
B1(N) = -QR - QRR
B0(N) = QR*QRR
NZE(II) = NZE(II) - 1
GO TO 70
C
60 B1(N) = -QR
B0(N) = 0.
C
70 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: BLDENZ
C BUILD DENOMINATOR BLOCKS OF SECOND ORDER -- Z-DOMAIN
C-----------------------------------------------------------------------
C
SUBROUTINE BLDENZ
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
COMMON /RESZ/ ZFA, ZM(18,4), ZZR(16), ZZI(16), ZPR(16), ZPI(16)
COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
C
NB = (NDEG+1)/2
N = 0
FACT = ZFA
C
DO 40 I=1,NB
N = N + 1
QR = ZPR(N)
QI = ZPI(N)
IF (ABS(QI).LT.FLMI) GO TO 10
C1(I) = -2.*QR
C0(I) = QR*QR + QI*QI
GO TO 40
C
10 IF (N.GE.NJ) GO TO 20
IF (ABS(ZPI(N+1)).LT.FLMI) GO TO 30
20 C1(I) = -QR
C0(I) = 0.
GO TO 40
C
30 N = N + 1
QI = ZPR(N)
C1(I) = -QR - QI
C0(I) = QR*QI
40 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: ROMEG
C REALIZED FREQUENCIES OMEGA
C-----------------------------------------------------------------------
C
SUBROUTINE ROMEG
C
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
COMMON /RESZ/ ZFA, ZM(18,4), ZZR(16), ZZI(16), ZPR(16), ZPI(16)
C
N2 = NZM(2)
N3 = NZM(3)
GO TO (10, 20, 30, 40), ITYP
10 ROM(1) = ZM(N2,2)
ROM(2) = ZM(1,3)
GO TO 50
20 ROM(1) = ZM(1,3)
ROM(2) = ZM(N2,2)
GO TO 50
30 ROM(1) = ZM(N3,3)
ROM(2) = ZM(1,2)
ROM(3) = ZM(N2,2)
ROM(4) = ZM(1,3)
GO TO 50
40 N2 = N2/2
ROM(1) = ZM(N2,2)
ROM(4) = ZM(N2+1,2)
ROM(3) = ZM(1,3)
ROM(2) = ZM(N3,3)
50 RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: TRANSN
C COMPUTATION OF THE PARAMETERS OF THE NORMALIZED LOWPASS
C-----------------------------------------------------------------------
C
SUBROUTINE TRANSN
C
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /TOLSN/ VSN, VD, A
COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
C
TAN2(AA) = SIN(AA/2.)/COS(AA/2.)
C
V1 = TAN2(OM(1))
V2 = TAN2(OM(2))
IF (ITYP.LE.2) GO TO 210
V3 = TAN2(OM(3))
V4 = TAN2(OM(4))
IF (ITYP.EQ.3) GO TO 10
Q = V1
V1 = -V4
V4 = -Q
Q = V2
V2 = -V3
V3 = -Q
IF (LSYM.NE.0) LSYM = 5 - LSYM
IF (LVSN.EQ.0) GO TO 10
J = 5 - LVSN
LVSN = LSYM
LSYM = J
C
10 J = LSYM + 1
JJ = LVSN/2 + 1
GO TO (20, 80, 120, 160, 180), J
20 J = NORMA + 1
GO TO (30, 30, 40, 70), J
30 VDQ1 = V2*V3
VSN1 = VDQ1/V1 - V1
Q = V4 - VDQ1/V4
IF (Q.LT.VSN1) VSN1 = Q
A1 = 1./(V3-V2)
VSN1 = VSN1*A1
GO TO (40, 50, 40), J
40 VDQ = V1*V4
A = V2/(VDQ-V2*V2)
Q = V3/(V3*V3-VDQ)
IF (Q.LT.A) A = Q
VSN = A*(V4-V1)
IF (NORMA.EQ.2) GO TO 200
IF (VSN.GE.VSN1) GO TO 200
50 VDQ = VDQ1
60 VSN = VSN1
A = A1
GO TO 200
C
70 VDQ = SQRT(V1*V2*V3*V4)
A1 = V3/(V3*V3-VDQ)
VSN1 = (V4-VDQ/V4)*A1
A = V2/(VDQ-V2*V2)
VSN = (VDQ/V1-V1)*A
IF (VSN.GE.VSN1) GO TO 200
GO TO 60
C
80 GO TO (90, 100, 110), JJ
90 VDQ = V2*V3
VSN = V4 - VDQ/V4
GO TO 190
100 V3 = V4*(V4+VSN*V2)/(V2+VSN*V4)
110 VDQ = V2*V3
A = 1./(V3-V2)
GO TO 200
120 GO TO (130, 150, 140), JJ
130 VDQ = V1*V4
A = V3/(V3*V3-VDQ)
GO TO 170
140 V4 = V3*(V1+VSN*V3)/(V3+VSN*V1)
150 VDQ = V1*V4
A = VSN/(V4-V1)
GO TO 200
160 VDQ = V1*V4
A = V2/(VDQ-V2*V2)
170 VSN = (V4-V1)*A
GO TO 200
180 VDQ = V2*V3
VSN = VDQ/V1 - V1
190 A = 1./(V3-V2)
VSN = VSN*A
C
200 VD = SQRT(VDQ)
A = A*VD
IF (ITYP.LE.3) GO TO 270
A = A/VSN
GO TO 270
C
210 J = LVSN*2 + ITYP
GO TO (220, 220, 230, 240, 250, 260), J
220 VSN = V2/V1
GO TO (250, 240), J
230 VD = V2/VSN
GO TO 270
240 VD = V2
GO TO 270
250 VD = V1
GO TO 270
260 VD = V1*VSN
C
270 RETURN
END
C
C ======================================================================
C
C DOREDI - SUBROUTINES: PART II
C
C ======================================================================
C
C
C-----------------------------------------------------------------------
C SUBROUTINE: GRID01
C BUILD A START GRID FOR THE SEARCH OF THE GLOBAL EXTREMES OF THE
C TRANSFER FUNCTION WITH ROUNDED COEFFICIENTS
C-----------------------------------------------------------------------
C
SUBROUTINE GRID01
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
COMMON /RESZ/ ZFA, ZM(18,4), ZZR(16), ZZI(16), ZPR(16), ZPI(16)
COMMON /CGRID/ GR(64), NGR(12)
C
M = NDEG
IF (NDEG.GT.22) M = 22
IF (ITYP.GE.3) M = M/2
M1 = M/2
FM = PI/FLOAT(2*M)
M11 = 2*M1 + 1
FM1 = PI/FLOAT(M11)
M = M + 1
PH2 = 1.5*PI
C
NGR(8) = 1
C
C BANDPASS
C
N = 0
GO TO (70, 10, 20, 40), ITYP
C
C HIGHPASS
C
10 IA = 1
IB = NZM(2)
GO TO 30
C
C BANDPASS
C
20 IA = NZM(3)
IB = 1
30 N = N + 1
GR(N) = 0.
GO TO 50
C
C BANDSTOP
C
40 IA = 1
IB = NZM(2)/2 + 1
50 J = 1
PH1 = 0.
JA = 3
JB = 2
GO TO 350
60 NGR(1) = N + 1
70 K = 2
L = 1
MAK = 1
MAL = 1
MEK = NZM(K)
MEL = NZM(L)
J = 1
GO TO (80, 90, 100, 110), ITYP
C
C LOWPASS
C
80 ID = 1
NGR(1) = 1
PH = 0.
IF (MEK.GT.MEL) GO TO 370
IQ = MEK
MEK = MEL
MEL = IQ
K = 1
L = 2
GO TO 370
C
C HIGHPASS
C
90 ID = -1
GO TO 130
C
C BANDPASS
C
100 MEK = MEK/2
MEL = MEK
GO TO 120
C
C BANDSTOP
C
110 MAK = MEK/2 + 1
MAL = MEL/2 + 1
120 ID = 1
130 PH = PH2
GO TO 370
C
140 NGR(2) = N
GO TO (150, 160, 160, 160), ITYP
C
C LOWPASS
C
150 IA = NZM(2)
IB = 1
JA = 2
JB = 3
PH1 = PH2
J = 2
GO TO 350
160 NGR(3) = N
IF (ITYP.LE.2) GO TO 230
K = 2
L = 1
MEK = NZM(K)
MEL = NZM(L)
IF (ITYP.EQ.3) GO TO 180
IF (MEK.GT.MEL) GO TO 170
K = 1
L = 2
170 MEK = MEK/2
MEL = MEL/2
MAK = 1
MAL = 1
GO TO 190
180 MAK = MEK/2 + 1
MAL = MAK
190 ID = 1
PH = 0.
J = 2
GO TO 370
C
200 NGR(4) = N
IF (ITYP.EQ.3) GO TO 210
C
C BANDSTOP
C
IA = NZM(2)/2
IB = NZM(3)
GO TO 220
C
C BANDPASS
C
210 IA = NZM(2)
IB = 1
220 JA = 2
JB = 3
PH1 = PH2
J = 3
C
GO TO 350
230 IF (ABS(GR(N)-PI).LT.FLMI) GO TO 240
N = N + 1
GR(N) = PI
240 NGR(9) = N
NGR(5) = N
C
C BANDSTOP
C
K = 3
L = 4
MEK = NZM(K)
MEL = NZM(L)
MAK = 1
MAL = 1
J = 3
ID = 1
GO TO (280, 250, 270, 250), ITYP
C
C HIGHPASS
C
250 IF (MEK.GT.MEL) GO TO 260
IQ = MEK
MEK = MEL
MEL = IQ
K = 4
L = 3
260 PH = 0.
ID = -1
IF (ITYP.NE.4) GO TO 330
C
C BANDSTOP
C
MAK = MEK/2 + 1
MAL = MEL/2 + 1
GO TO 280
C
C BANDPASS
C
270 MEK = MEK/2
MEL = MEL/2
280 PH = PH2
GO TO 380
C
290 NGR(6) = N
J = 4
GO TO (500, 500, 300, 320), ITYP
C
C BANDPASS
C
300 IF (MEK.GT.MEL) GO TO 310
IQ = MEK
MEK = MEL
MEL = IQ
K = 4
L = 3
310 MAK = MEK + 1
MAL = MEL + 1
MEK = NZM(K)
MEL = NZM(L)
GO TO 330
C
C BANDSTOP
C
320 MEK = MEK/2
MEL = (MEL+1)/2
K = 4
L = 3
MAK = 1
MAL = 1
330 PH = 0.
GO TO 380
C
340 NGR(7) = N
GO TO 500
C
350 QA = ZM(IA,JA)
QD = ZM(IB,JB) - QA
IF (PH1.NE.0.) QA = QA + QD
DO 360 I=1,M1
N = N + 1
Q = FLOAT(I)*FM1 + PH1
GR(N) = QA + SIN(Q)*QD
360 CONTINUE
GO TO (60, 160, 230), J
C
370 GO TO (460, 390, 460, 390), IAPRO
380 GO TO (460, 460, 390, 390), IAPRO
C
C MULTIPLE EXTREMA
C
390 IF (ID.LT.0) GO TO 400
IK = MAK
IL = MAL
GO TO 410
C
400 IK = MEK
IL = MEL
410 DO 450 I=MAK,MEK
N = N + 1
GR(N) = ZM(IK,K)
IF (IL.GT.MEL .OR. IL.LT.MAL) GO TO 440
IF (NDEG.LE.22) GO TO 420
IF (IL.NE.MEL .AND. IL.NE.MAL) GO TO 430
420 N = N + 1
GR(N) = ZM(IL,L)
430 IL = IL + ID
440 IK = IK + ID
450 CONTINUE
GO TO 490
C
C SINGLE EXTREMA
C
460 QA = ZM(MEK,K)
QD = ZM(MEL,L)
IF (QD.GT.QA) GO TO 470
Q = QA
QA = QD
QD = Q
470 QD = QD - QA
IF (PH.NE.0.) QA = QA + QD
DO 480 I=1,M
N = N + 1
II = I - 1
Q = FLOAT(II)*FM + PH
GR(N) = QA + SIN(Q)*QD
480 CONTINUE
490 GO TO (140, 200, 290, 340), J
C
500 NGR(10) = NGR(3) + 1
NGR(11) = NGR(5) + 1
NGR(12) = NGR(6) + 1
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: GRID03
C START VALUES FOR THE CALCULATION OF A GRID
C-----------------------------------------------------------------------
C
SUBROUTINE GRID03
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /RESZ/ ZFA, ZM(18,4), ZZR(16), ZZI(16), ZPR(16), ZPI(16)
COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
C
IF (ITYP.LE.0 .OR. ITYP.GT.4) GO TO 70
IF (NDEG.LE.0) GO TO 70
IAPRO = 1
DO 10 I=1,4
NZM(I) = (ITYP+1)/2
10 CONTINUE
ME = 2
IF (ITYP.GE.3) ME = 4
DO 20 I=1,ME
IF (OM(I).LE.0.) GO TO 80
20 CONTINUE
C
GO TO (30, 40, 50, 60), ITYP
30 ZM(1,1) = 0.
ZM(1,2) = OM(1)
ZM(1,3) = OM(2)
ZM(1,4) = PI
GO TO 100
40 ZM(1,1) = PI
ZM(1,2) = OM(2)
ZM(1,3) = OM(1)
ZM(1,4) = 0.
GO TO 100
50 NZM(1) = 1
ZM(1,1) = (OM(2)+OM(3))/2.
ZM(1,2) = OM(2)
ZM(2,2) = OM(3)
ZM(1,3) = OM(4)
ZM(2,3) = OM(1)
ZM(1,4) = PI
ZM(2,4) = 0.
GO TO 100
60 ZM(1,1) = 0.
ZM(2,1) = PI
ZM(1,2) = OM(1)
ZM(2,2) = OM(4)
ZM(1,3) = OM(3)
ZM(2,3) = OM(2)
NZM(4) = 1
ZM(1,4) = (OM(2)+OM(3))/2.
GO TO 100
C
70 IE = 20
GO TO 90
80 IE = 23
90 CALL ERROR(IE)
100 RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: GRID04
C CHANGE GRID FOR UNSYMMETRICAL TOLERANCE SCHEME
C-----------------------------------------------------------------------
C
SUBROUTINE GRID04
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
COMMON /CGRID/ GR(64), NGR(12)
C
IF (ITYP.LE.2) GO TO 130
DO 10 I=1,4
IF (OM(I).EQ.0.) GO TO 130
10 CONTINUE
C
I = 1
MA = NGR(1)
ME = NGR(2)
IF (ITYP.EQ.3) OMM = OM(2)
IF (ITYP.EQ.4) OMM = OM(4)
N = 1
GO TO 50
C
20 MA = NGR(10)
ME = NGR(4)
IF (ITYP.EQ.3) OMM = OM(3)
IF (ITYP.EQ.4) OMM = OM(1)
N = 4
GO TO 90
C
30 I = 2
MA = NGR(11)
ME = NGR(6)
IF (ITYP.EQ.3) OMM = OM(4)
IF (ITYP.EQ.4) OMM = OM(2)
N = 11
GO TO 50
C
40 MA = NGR(12)
ME = NGR(7)
IF (ITYP.EQ.3) OMM = OM(1)
IF (ITYP.EQ.4) OMM = OM(3)
N = 7
GO TO 90
C
50 DO 60 II=MA,ME
J = ME + MA - II
IF ((GR(J)-OMM).LE.FLER) GO TO 70
60 CONTINUE
GO TO 80
70 GR(J) = OMM
NGR(N) = J
80 IF (I.EQ.1) GO TO 20
GO TO 40
C
90 DO 100 II=MA,ME
IF ((GR(II)-OMM).GE.(-FLER)) GO TO 110
100 CONTINUE
GO TO 120
110 GR(II) = OMM
NGR(N) = II
120 IF (I.EQ.1) GO TO 30
C
130 RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: COEFW
C SEARCH OF THE MINIMUM COEFFICIENT WORDLENGTH
C COMPUTATION OF THE OPTIMIZATION CRITERION
C-----------------------------------------------------------------------
C
SUBROUTINE COEFW
C
COMMON /CONTR/ IPRUN, IPCON, NINP, NOUT, NDOUT, LSPOUT, NSPOUT
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
* ISTAB, IDEPSL, IDEPSD
COMMON /CEPSIL/ EPS(2), PMAX
COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
* LPOT2, JSTRU
COMMON /CREST1/ JSTRUS, JSTRUD
COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
COMPLEX ZP, ZPS
COMMON /CPOL/ ZP(16,2), ZPS(16,2)
C
IWLRR = IWLR
IWLR = 100
IWLG = IWLMA
IWLD = 0
IWLM = IWL
II1 = 1
II2 = 0
II3 = 0
CALL COPY01
IF (IPRUN.GT.0) GO TO 20
IF (LSEQ.NE.(-1)) GO TO 10
CALL ALLO02
CALL COPY01
IF (NOUT.GE.2) CALL OUT011
10 CALL ALLO01
CALL COPY01
20 IF (NOUT.EQ.4) CALL OUT043
CALL RECO
30 CALL COPY02
IF (IPRUN.GT.0) GO TO 70
CALL ROUND(4, 5)
IF (ISCAL.NE.3) GO TO 40
CALL POLLOC
CALL COPY03
40 AMAX = 1.
FACT = 1.
DO 60 J=1,NB
IB = J
50 CALL SCAL01
CALL SCAL02
IF (JSTRUD.EQ.1) GO TO 60
IF (IB.NE.NB) GO TO 60
IB = IB + 1
GO TO 50
60 CONTINUE
CALL RECO
CALL ROUND(1, 3)
GO TO 80
70 CALL ROUND(1, 5)
80 CALL STABLE
IF (NOUT.GE.5) CALL OUT046
IF (ISTAB.NE.(-1)) GO TO 250
CALL EPSILO
IDEPS = 0
DO 90 K=1,2
IF (EPS(K).GT.1.) IDEPS = IDEPS + 1
90 CONTINUE
IF (NOUT.GE.5) CALL OUT043
IF (NOUT.GE.4) CALL OUT044
ADEPS = EPS(1) - EPS(2)
IF (LOPTW.NE.(-1)) GO TO 260
IF (IWL.NE.IWLL) GO TO 100
ADEPSL = ADEPS
IDEPSL = IDEPS
100 IF (IDEPS.EQ.0) GO TO 110
IF (IWL.LE.IWLD .OR. IWL.GT.IWLG) GO TO 120
IWLD = IWL
ADEPSD = ADEPS
IDEPSD = IDEPS
GO TO 120
110 IF (IWL.GE.IWLG) GO TO 120
IWLG = IWL
ADEPSG = ADEPS
IF (IWLD.LE.IWLG) GO TO 120
IWLD = 0
II3 = 0
120 GO TO (130, 140, 180), II1
130 II1 = 2
IF (IDEPS.EQ.2) II1 = 3
GO TO 120
140 IF (IDEPS.NE.0) GO TO 170
150 II2 = II2 + 1
160 IWL = IWL - 1
GO TO 220
170 II3 = II3 + 1
IF (II3.LT.2) GO TO 160
IF (II2.GE.2) GO TO 270
IWL = IWLM + 1
II1 = 3
GO TO 220
180 IF (IDEPS.EQ.0) GO TO 210
190 II3 = II3 + 1
200 IWL = IWL + 1
GO TO 220
210 II2 = II2 + 1
IF (II2.LT.2) GO TO 200
IF (II3.GE.2) GO TO 270
IWL = IWLM - 1
II1 = 2
220 IF (IWL.GE.IWLMI .AND. IWL.LE.IWLMA) GO TO 30
230 IF (II1.EQ.3) GO TO 240
IF (IDEPS.EQ.0) GO TO 270
II1 = 3
II2 = 2
GO TO 190
240 II1 = 2
IF (IDEPS.NE.0) GO TO 270
II3 = 2
GO TO 150
250 IF (LOPTW.NE.(-1)) GO TO 270
IF (II1.NE.1) GO TO 230
IWL = IWL + 1
IF (IWL.LE.IWLMA) GO TO 30
260 ADEPSD = ADEPS
IWLG = IWL
IWLD = IWL
270 RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: RECO
C TRANSFORM THE COEFFICIENTS FOR THE CHOSEN REALIZATION
C-----------------------------------------------------------------------
C
SUBROUTINE RECO
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /CRECO/ JRCO, JECO, JDCO, JJDCO(5), JMAXV, JTRB2, LREF
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
COMMON /FILTRE/ IRCO(5), IECO(16,5), IDCO(16,5), IECOM
C
DIMENSION ACO(16,5)
EQUIVALENCE (ACO(1,1),B2(1))
EQUIVALENCE (IRCO(1),IRB2), (IRCO(2),IRB1), (IRCO(3),IRB0)
EQUIVALENCE (IRCO(4),IRC1), (IRCO(5),IRC0)
C
ALOG2 = ALOG(2.)
C
C REPRESENTATION OF THE COEFFICIENTS WITH DIFFERENCES
C
IF (JDCO.EQ.0) GO TO 70
DO 60 J=1,5
ACOM = 0.
DO 10 I=1,NB
ACOM = AMAX1(ABS(ACO(I,J)),ACOM)
10 CONTINUE
IF (ACOM.GE.FLMI) GO TO 20
ID = -100
GO TO 30
20 ID = ALOG(ACOM)/ALOG2
Q = 2.**ID
IF (Q.GE.ACOM) GO TO 30
ID = ID + 1
Q = Q*2.
30 JJ = JJDCO(J)
IF (JJ.EQ.0) GO TO 60
BCOM = 0.5*Q
DO 50 I=1,NB
AC = ACO(I,J)
ACA = ABS(AC)
IF (ACA.GE.BCOM) GO TO 40
IF (JJ.NE.1) GO TO 50
40 IDCO(I,J) = ID
50 CONTINUE
60 CONTINUE
C
C RANGE OF COEFFICIENT NUMBERS
C
70 IF (JRCO.EQ.0) GO TO 180
DO 120 J=1,5
ACOM = 0.
BCOM = 0.
DO 80 I=1,NB
ACS = ACOEFS(ACO(I,J),0,IDCO(I,J),0)
ACOM = AMIN1(ACOM,ACS)
BCOM = AMAX1(BCOM,ACS)
80 CONTINUE
ACOM = AMAX1(ABS(ACOM),BCOM)
IF (ACOM.GE.FLMI) GO TO 90
IQ = -100
GO TO 110
90 IQ = ALOG(ACOM)/ALOG2
Q = 2.**IQ
IF (Q.LT.ACOM) GO TO 100
IF (JMAXV.LE.1) GO TO 110
IF (JMAXV.EQ.2 .AND. Q.GT.BCOM) GO TO 110
IF (JMAXV.EQ.3 .AND. Q.GT.ACOM) GO TO 110
100 IQ = IQ + 1
110 IF (JRCO.GT.0) IQ = MAX0(IQ,0)
IRCO(J) = IQ
120 CONTINUE
JJ = IABS(JRCO)
GO TO (180, 130, 130, 150, 160), JJ
130 IQ = MAX0(IRB2,IRB1,IRB0)
DO 140 J=1,3
IRCO(J) = IQ
140 CONTINUE
IF (JJ.EQ.2) GO TO 180
IQ = MAX0(IRC1,IRC0)
IRC1 = IQ
IRC0 = IQ
GO TO 180
150 IQ = MAX0(IRB1,IRC1)
IRB1 = IQ
IRC1 = IQ
GO TO 180
160 IQ = MAX0(IRB2,IRB1,IRB0,IRC1,IRC0)
DO 170 J=1,5
IRCO(J) = IQ
170 CONTINUE
C
C COMPUTE "PSEUDO" FLOATING-POINT EXPONENTS
C
180 IF (JECO.EQ.0) GO TO 230
IECOM = 0
DO 220 J=1,5
DO 210 I=1,NB
QA = ACOEFS(ACO(I,J),0,IDCO(I,J),IRCO(J))
QA = ABS(QA)
IF (QA.GE.FLMI) GO TO 190
IQ = JECO
GO TO 200
190 IQ = -ALOG(QA)/ALOG2
IF (IQ.GT.JECO) IQ = JECO
200 IECO(I,J) = IQ
IECOM = MAX0(IECOM,IQ)
210 CONTINUE
220 CONTINUE
230 RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: ROUND
C ROUNDING OF THE CHANGED COEFFICIENTS
C-----------------------------------------------------------------------
C
SUBROUTINE ROUND(IRA, IRB)
C
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
COMMON /CRECO/ JRCO, JECO, JDCO, JJDCO(5), JMAXV, JTRB2, LREF
COMMON /FILTRE/ IRCO(5), IECO(16,5), IDCO(16,5), IECOM
COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
* ISTAB, IDEPSL, IDEPSD
C
DIMENSION ACO(16,5)
EQUIVALENCE (ACO(1,1),B2(1))
C
ALSB = 2.**(1-IWL)
ELSB = 1. - ALSB
DO 120 I=1,NB
DO 110 J=IRA,IRB
AC = ACO(I,J)
IF (J.NE.2 .AND. J.NE.3) GO TO 10
IF (ACO(I,1).EQ.AC) GO TO 110
10 ACS = ACOEFS(AC,IECO(I,J),IDCO(I,J),IRCO(J))
ACSA = ABS(ACS)
IF (ACSA.LT.ALSB .AND. J.EQ.1) GO TO 110
IF (ACSA.GT.ELSB) GO TO 30
R = 0.5
IF (J.NE.1 .OR. JTRB2.LE.0) GO TO 20
R = 0.
IF (IDCO(I,J).GT.(-100)) R = 1.
20 ACSA = AINT(ACSA/ALSB+R)
ACSA = ACSA*ALSB
GO TO 50
30 IF (JMAXV.GE.1) GO TO 40
IF (JMAXV.EQ.(-2) .AND. AC.LT.0.) GO TO 40
ACSA = ELSB
GO TO 50
40 ACSA = 1.
50 BC = SIGN(ACSA,ACS)
BC = ACOEF(BC,IECO(I,J),IDCO(I,J),IRCO(J))
ACO(I,J) = BC
IF (J.NE.1) GO TO 110
F = BC/AC
JJ = IABS(JTRB2)
GO TO (100, 80, 60), JJ
60 DO 70 K=2,3
ACO(I,K) = ACO(I,K)*F
70 CONTINUE
80 II = I + 1
IF (II.GT.NB) GO TO 100
DO 90 K=1,3
ACO(II,K) = ACO(II,K)/F
90 CONTINUE
GO TO 110
100 FACT = FACT/F
110 CONTINUE
120 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: STABLE
C CHECK STABILITY OF THE ROUNDED SYSTEM
C-----------------------------------------------------------------------
C
SUBROUTINE STABLE
C
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
COMMON /SFILT/ NBS, SFACT, SB2(16), SB1(16), SB0(16), SC1(16),
* SC0(16)
COMMON /FILTRE/ IRCO(5), IECO(16,5), IDCO(16,5), IECOM
COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
* ISTAB, IDEPSL, IDEPSD
COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
C
DO 120 I=1,NB
IT = 0
CC0 = C0(I)
CC1 = C1(I)
IF (CC0.GE.1.) GO TO 20
AC1 = ABS(CC1)
10 IF (AC1.GE.1.+CC0) GO TO 30
IF (IT.EQ.0) GO TO 120
GO TO 110
20 IT = IT + 1
GO TO 40
30 IT = IT + 2
40 IF (LSTAB.NE.(-1)) GO TO 130
ALSB = 2.**(1-IWL)
Q0 = ACOEF(ALSB,IECO(I,5),-100,IRCO(5))
IF (IT.GE.2) GO TO 50
CC0 = CC0 - Q0
GO TO 10
50 IF (IT.GE.3) GO TO 80
AC0 = CC0
60 AC0 = AC0 + Q0
IF (AC0.GE.1.) GO TO 70
IF (AC1.GE.1.+AC0) GO TO 60
Q3 = AC1 - ABS(SC1(I))
Q2 = AC0 - SC0(I)
A = Q3*Q3 + Q2*Q2
A = SQRT(A)
GO TO 80
70 IT = 3
80 Q1 = ACOEF(ALSB,IECO(I,4),-100,IRCO(4))
BC1 = AC1
90 BC1 = BC1 - Q1
IF (BC1.GE.1.+AC0) GO TO 90
IF (IT.EQ.3) GO TO 100
Q3 = BC1 - ABS(SC1(I))
Q2 = CC0 - SC0(I)
B = Q3*Q3 + Q2*Q2
B = SQRT(B)
IF (A.GE.B) GO TO 100
BC1 = AC1
CC0 = AC0
100 C1(I) = SIGN(BC1,CC1)
110 C0(I) = CC0
120 CONTINUE
ISTAB = -1
RETURN
130 ISTAB = 0
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: EPSILO
C COMPUTATION OF THE GLOBAL EXTREMA OF THE TRANSFER FUNCTION
C WITH ROUNDED COEFFICIENTS
C COMPUTATION OF THE ERROR VALUES EPSILON
C-----------------------------------------------------------------------
C
SUBROUTINE EPSILO
C
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /CGRID/ GR(64), NGR(12)
COMMON /CEPSIL/ EPS(2), PMAX
C
C MAXIMUM IN THE PASSBAND(S)
C
ME = 3
IF (ITYP.EQ.3) ME = 5
CALL SMAX(PMAX, 8, ME, 1)
IF (ITYP.NE.4) GO TO 10
CALL SMAX(QMAX, 10, 5, 1)
PMAX = AMAX1(PMAX,QMAX)
C
C MINIMUM IN THE PASSBAND(S)
C
10 ME = 2
IF (ITYP.EQ.3) ME = 4
CALL SMAX(PMIN, 1, ME, -1)
IF (ITYP.NE.4) GO TO 20
CALL SMAX(QMIN, 10, 4, -1)
PMIN = AMIN1(PMIN,QMIN)
C
C MAXIMUM IN THE STOPBAND(S)
C
20 ME = 6
IF (ITYP.EQ.4) ME = 7
CALL SMAX(STMAX, 11, ME, 1)
IF (ITYP.NE.3) GO TO 30
CALL SMAX(QMAX, 12, 7, 1)
STMAX = AMAX1(STMAX,QMAX)
C
C COMPUTATION OF EPSILON
C
30 EPS(1) = (1.-PMIN/PMAX)/ADELP
EPS(2) = STMAX/PMAX/ADELS
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OPTPAR
C OPTIMIZATION OF THE COEFFICIENT WORDLENGTH BY VARIATION OF THE
C DESIGN PARAMETER ACX
C-----------------------------------------------------------------------
C
SUBROUTINE OPTPAR
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
* ISTAB, IDEPSL, IDEPSD
COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
COMMON /COPTCP/ IS1, IS2, IABO, IWLGS, IWLGN, IWLGP, ADEPSS,
* ADEPSN, ADEPSP, ACXS, ACXN, ACXP
C
IF (LOPTW.NE.(-1)) GO TO 40
IF (IWLD.GT.IWLL .AND. IDEPSL.LT.2) GO TO 20
IF (IDEPSD.NE.2) GO TO 10
ADEPSD = ADEPSG
IWLD = IWLG
10 IWLL = IWLD
GO TO 30
20 ADEPSD = ADEPSL
IWLD = IWLL
30 IWL = IWLD
40 IABO = IABO + 1
IF (ADEPSD.LT.0.) GO TO 80
IF (IS1.NE.(-1)) GO TO 50
IF (IS2.NE.(-1)) GO TO 70
GO TO 60
50 IF (IS2.NE.(-1)) GO TO 60
IS1 = -1
60 IS2 = -((IS2+2)/2)
ACXA = (ACX+ACXN)/2.
GO TO 120
70 ACXA = (ACX+ACXMI)/2.
GO TO 120
C
80 IF (IS1.EQ.(-1)) GO TO 90
IF (IS2.NE.(-1)) GO TO 110
GO TO 100
90 IF (IS2.NE.(-1)) GO TO 100
IS1 = 0
100 IS2 = -((IS2+2)/2)
ACXA = (ACX+ACXP)/2.
GO TO 120
110 ACXA = (ACX+ACXMA)/2.
C
120 IF (IWLG.LT.IWLGS) GO TO 130
IF (ADEPSS-ABS(ADEPSD).LE.FLMI .OR. IWLG.GT.IWLGS) GO TO 140
130 IWLGS = IWLG
ADEPSS = ABS(ADEPSD)
ACXS = ACX
140 IF (ADEPSD.LT.0.) GO TO 160
IF (IWLG.LT.IWLGP) GO TO 150
IF (ADEPSP-ADEPSD.LE.FLMI .OR. IWLG.GT.IWLGP) GO TO 190
150 IWLGP = IWLG
ADEPSP = ADEPSD
ACXP = ACX
GO TO 180
160 IF (IWLG.LT.IWLGN) GO TO 170
IF (ADEPSN-ADEPSD.GE.FLMI .OR. IWLG.GT.IWLGN) GO TO 190
170 IWLGN = IWLG
ADEPSN = ADEPSD
ACXN = ACX
180 IABO = 0
190 ACX = ACXA
RETURN
END
C
C ======================================================================
C
C DOREDI - SUBROUTINES: PART III
C
C ======================================================================
C
C
C-----------------------------------------------------------------------
C SUBROUTINE: GRID02
C BUILD A START GRID OUT OF THE POLE LOCATIONS FOR THE SEARCH OF THE
C GLOBAL EXTREMES OF THE TRANSFER FUNCTION
C-----------------------------------------------------------------------
C
SUBROUTINE GRID02
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
COMMON /CGRID/ GR(64), NGR(12)
C
NGR(1) = 0
J = 1
GR(1) = 0.
DO 10 I=1,NB
Q = -C1(I)/2.
QN = C0(I)
IF (ABS(QN).LT.FLMI) GO TO 10
QN = QN - Q*Q
IF (QN.LE.0.) GO TO 10
J = J + 1
GR(J) = ATAN2(SQRT(QN),Q)
10 CONTINUE
J = J + 1
GR(J) = PI
NGR(8) = 1
NGR(9) = J
C
C SORTING OF THE LOCATIONS
C
J = J - 1
20 JJ = 0
DO 30 I=1,J
Q1 = GR(I)
Q2 = GR(I+1)
IF (Q1.LT.Q2) GO TO 30
JJ = 1
GR(I) = Q2
GR(I+1) = Q1
30 CONTINUE
IF (JJ.NE.0) GO TO 20
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: ANANOI
C NOISE ANALYSIS OF A GIVEN CASCADE
C-----------------------------------------------------------------------
C
SUBROUTINE ANANOI
C
COMMON /CONTR/ IPRUN, IPCON, NINP, NOUT, NDOUT, LSPOUT, NSPOUT
COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
COMMON /COPTST/ LOPTS, ISTOR
COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
* LPOT2, JSTRU
COMMON /CREST1/ JSTRUS, JSTRUD
COMMON /CPOW/ PNU, PNC, AND, ITCORP
COMMON /CNFUNC/ AQC(5), BN2(5), BN1(5), BN0(5)
COMMON /CFFUNC/ PHI(5), BF2(5), BF1(5), BF0(5), IBB(5), ICOR
COMMON /CNOISE/ RI, RIN, RE, REN, FAC
COMMON /COPST2/ PN(100,2), TF(100,2), TFA(100,2)
C
IQ = IWLR
IF (IQ-1.GT.IWLMA) IQ = IWLMA + 1
IQ = IQ - 2
ALSBI = 2.**IQ
IB = NB
IF (LOPTS.NE.0) GO TO 20
DO 10 I=1,NB
JSEQN(I) = I
JSEQD(I) = I
10 CONTINUE
GO TO 40
20 DO 30 I=1,NB
JSEQN(I) = ISEQ(I,1)
JSEQD(I) = ISEQ(I,2)
30 CONTINUE
40 SAND = 0.
SPNU = 0.
SPNC = 0.
CALL ALLOND
CALL SMAX(Q, 8, 9, 1)
FAC1 = FACT/Q
IF (ISCAL.NE.0 .OR. JSTRUD.EQ.1) FACT = 1.
AMAX = 1.
ITCORP = -1
C
DO 80 I=1,NB
IB = I
IF (ISTRU.GE.30) GO TO 60
50 IF (ISCAL.EQ.0) SCA = 1.
IF (ISCAL.NE.0) CALL SCAL01
CALL ALNS01
CALL ALNS02
GO TO 70
60 CONTINUE
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C HERE INSERT THE CALLS TO YOUR OWN STRUCTURE IMPLEMENTATIONS
C SEE SUBROUTINE OPTLST
C
CALL ERROR(8)
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
70 Q = FAC1/FACT
CALL POWER
AND = AND*Q
PNC = Q
Q = Q*Q
PNU = PNU*Q
CALL PCORP
PNC = PNC*Q
SAND = SAND + AND
SPNU = SPNU + PNU
SPNC = SPNC + PNC
IF (NOUT.GE.3) CALL OUT026
IF (ISCAL.NE.0) CALL SCAL02
IF (JSTRUD.EQ.1 .OR. IB.NE.1) FAC1 = FAC1/SCA
IF (JSTRUD.EQ.1) GO TO 80
IF (IB.NE.NB) GO TO 80
IB = IB + 1
GO TO 50
80 CONTINUE
CALL TCORP
RIN = SPNU + SAND*SAND + PNC
FAC = FACT/FAC1
Q = FAC*FAC
RI = RIN*Q
IF (NOUT.LT.3) GO TO 90
IB = 0
AND = SAND
PNU = SPNU
CALL OUT026
90 CALL SPOW(RE)
REN = RE/Q
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OPTLST
C OPTIMIZATION OF THE PAIRING AND ORDERING BY THE PROCEDURE WITH
C LIMITED STORAGE
C DYNAMIC PROGRAMMING FOR SYSTEMS OF AN ORDER LESS THAN OR EQUAL TO TEN
C DIRECT PROCEDURE FOR THE PARAMETER ISTOR EQUAL TO ONE
C-----------------------------------------------------------------------
C
SUBROUTINE OPTLST
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
* LPOT2, JSTRU
COMMON /CREST1/ JSTRUS, JSTRUD
COMMON /COPTST/ LOPTS, ISTOR
COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
COMMON /COPST1/ LSEQN(16,100,2), LSEQD(16,100,2)
COMMON /COPST2/ PN(100,2), TF(100,2), TFA(100,2)
COMMON /CPOW/ PNU, PNC, AND, ITCORP
COMMON /CNFUNC/ AQC(5), BN2(5), BN1(5), BN0(5)
COMMON /CFFUNC/ PHI(5), BF2(5), BF1(5), BF0(5), IBB(5), ICOR
COMMON /CNOISE/ RI, RIN, RE, REN, FAC
C
DIMENSION ISEQN(16), ISEQD(16)
EQUIVALENCE (ISEQN(1),ISEQ(1,1)), (ISEQD(1),ISEQ(1,2))
C
ITCORP = 0
ALSBI = 2.**(IWLMA-1)
II = 1
PN(1,2) = 0.
TF(1,2) = 1.
TFA(1,2) = 1.
CALL SMAX(Q, 8, 9, 1)
FAC1 = FACT/Q
JS = 1
IS = 2
IN = 1
PNTM = FLMA
MA = NB
IF (LOPTS.NE.(-1)) MA = 1
C
DO 220 IK=1,NB
IIK = IK
MK = IK - 1
IQ = IS
IS = JS
JS = IQ
JN = IN
IN = 0
C
C EXTENSION OF THE ASSIGNED ORDERING
C
DO 210 JK=1,JN
IF (IK.EQ.1) GO TO 10
CALL CODE4(JSEQD, LSEQD(1,JK,JS), MK)
IF (LOPTS.EQ.(-1)) CALL CODE4(JSEQN, LSEQN(1,JK,JS), MK)
10 DO 200 I=1,NB
IF (IK.EQ.1) GO TO 30
DO 20 II=1,MK
IF (JSEQD(II).EQ.I) GO TO 200
20 CONTINUE
30 JSEQD(IK) = I
DO 190 J=1,MA
IF (LOPTS.NE.(-1)) GO TO 60
IF (IK.EQ.1) GO TO 50
DO 40 II=1,MK
IF (JSEQN(II).EQ.J) GO TO 190
40 CONTINUE
50 JSEQN(IK) = J
C
60 FACT = TF(JK,JS)
AMAX = TFA(JK,JS)
PNT = PN(JK,JS)
IB = IK
CALL ALLOND
70 IF (ISTRU.GE.30) GO TO 80
CALL SCAL01
CALL ALNS01
CALL ALNS02
GO TO 90
80 CONTINUE
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C HERE INSERT THE CALLS TO YOUR OWN STRUCTURE IMPLEMENTATIONS
C SCALING OF THE IB-TH BLOCK
C ALLOCATIONS OF THE NOISE SOURCES IN THE IB-TH BLOCK
C ALLOCATION OF THE INPUT TRANSFER FUNCTIONS
C FACT : PRODUCT OF THE IB-1 SCALING FACTORS IN FRONT
C AMAX : SCALING OF THE (IB-1)-TH BLOCK
C JSTRU = 1 : SECOND COEFFICIENT OPTIMIZATION
C JSTRU = 2 : NO SECOND OPTIMIZATION
C JSTRUS = 1 : NO SEPARATE SCALING FACTOR
C JSTRUS = 2 : SEPARATE SCALING FACTOR
C JSTRUD = 1 : NUMERATOR FIRST
C JSTRUD = 2 : DENOMINATOR FIRST
C BN2(I),BN1(I),BN0(I) : NOISE TRANSFER FUNCTION IN THE IB-TH BLOCK
C BF2(I),BF1(I),BF0(I) : TRANSFER FUNCTIONS IN THE IB-TH BLOCK
C FROM THE INPUT TO THE NONLINEARITIES
C IBB(I) .NE. IB : BLOCK IB ONLY NUMERATOR
C PHI(I) : VARIANCE OF THE INPUT TRANSFER FUNCTIONS
C
C REMOVE THE FOLLOWING STATEMENT
C
CALL ERROR(8)
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
90 Q = FAC1/FACT
CALL POWER
CALL PCORP
PNT = (PNU+PNC+AND*AND)*Q*Q + PNT
IF (JSTRUD.EQ.1) GO TO 100
IF (IB.NE.NB) GO TO 100
IB = IB + 1
GO TO 70
100 IF (IK.EQ.NB) GO TO 170
C
C STORE INTERMEDIATE RESULT
C CHECK IF COMBINATION IS IN STORAGE
C
IF (IK.EQ.1) GO TO 130
IF (IN.EQ.0) GO TO 150
DO 120 IL=1,IN
CALL CODE5(ISEQD, LSEQD(1,IL,IS), IIK, JSEQD, II)
IF (II.EQ.2) GO TO 120
IF (LOPTS.NE.(-1)) GO TO 110
CALL CODE5(ISEQN, LSEQN(1,IL,IS), IIK, JSEQN, II)
IF (II.EQ.2) GO TO 120
110 CONTINUE
IF (PN(IL,IS).LE.PNT) GO TO 190
CALL CODE3(JSEQD, LSEQD(1,IL,IS), IIK)
IF (LOPTS.EQ.(-1)) CALL CODE3(JSEQN, LSEQN(1,IL,IS),
* IIK)
PN(IL,IS) = PNT
TF(IL,IS) = FACT
TFA(IL,IS) = AMAX
GO TO 190
120 CONTINUE
130 IF (IN.LT.ISTOR) GO TO 150
C
C SEARCH IN THE STORAGE FOR THE ALLOCATION WITH GREATEST NOISE POWER
C
PNM = PNT
IQ = 0
DO 140 II=1,ISTOR
IF (PN(II,IS).LT.PNM) GO TO 140
IQ = II
PNM = PN(II,IS)
140 CONTINUE
IF (IQ.EQ.0) GO TO 190
GO TO 160
150 IN = IN + 1
IQ = IN
160 PN(IQ,IS) = PNT
TF(IQ,IS) = FACT
TFA(IQ,IS) = AMAX
CALL CODE3(JSEQD, LSEQD(1,IQ,IS), IIK)
IF (LOPTS.EQ.(-1)) CALL CODE3(JSEQN, LSEQN(1,IQ,IS), IIK)
GO TO 190
C
C CHECK RESULT
C
170 IF (PNT.GE.PNTM) GO TO 190
PNTM = PNT
FAC = FACT
DO 180 II=1,NB
ISEQD(II) = JSEQD(II)
ISEQN(II) = JSEQN(II)
180 CONTINUE
190 CONTINUE
200 CONTINUE
IF (IK.EQ.1) GO TO 220
210 CONTINUE
220 CONTINUE
FAC = FAC/FAC1
RIN = PNTM
RI = RIN*FAC*FAC
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: ALLOND
C ALLOCATION OF THE NUMERATOR AND DENOMINATOR BLOCKS
C-----------------------------------------------------------------------
C
SUBROUTINE ALLOND
C
COMPLEX ZP, ZPS
C
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
COMMON /SFILT/ NBS, SFACT, SB2(16), SB1(16), SB0(16), SC1(16),
* SC0(16)
COMMON /CPOL/ ZP(16,2), ZPS(16,2)
COMMON /COPTST/ LOPTS, ISTOR
COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
* LPOT2, JSTRU
COMMON /CREST1/ JSTRUS, JSTRUD
COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
C
DO 40 I=1,IB
J = JSEQD(I)
C1(I) = SC1(J)
C0(I) = SC0(J)
ZP(I,1) = ZPS(J,1)
ZP(I,2) = ZPS(J,2)
IF (LOPTS.NE.1) GO TO 20
IF (JSTRUD.EQ.1) GO TO 30
IF (I.NE.IB) GO TO 10
J = JSEQD(1)
GO TO 30
10 J = JSEQD(I+1)
GO TO 30
20 J = JSEQN(I)
30 B2(I) = SB2(J)
B1(I) = SB1(J)
B0(I) = SB0(J)
JSEQN(I) = J
40 CONTINUE
IF (IB.GE.NB) RETURN
II = IB
JJ = IB
DO 80 I=1,NB
DO 50 J=1,IB
IF (JSEQD(J).EQ.I) GO TO 60
50 CONTINUE
II = II + 1
C1(II) = SC1(I)
C0(II) = SC0(I)
ZP(II,1) = ZPS(I,1)
ZP(II,2) = ZPS(I,2)
60 DO 70 J=1,IB
IF (JSEQN(J).EQ.I) GO TO 80
70 CONTINUE
JJ = JJ + 1
B2(JJ) = SB2(I)
B1(JJ) = SB1(I)
B0(JJ) = SB0(I)
80 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: SCAL01
C SCALING OF THE IB-TH BLOCK
C-----------------------------------------------------------------------
C
SUBROUTINE SCAL01
C
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
* LPOT2, JSTRU
COMMON /CREST1/ JSTRUS, JSTRUD
COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
C
NNB = NB
BMAX = 0.
J = 1
IF (IB.GT.NB) GO TO 30
NB = IB
IF (JSTRUD.EQ.1) GO TO 30
J = 2
BB1 = B1(IB)
BB0 = B0(IB)
B1(IB) = 0.
B0(IB) = 0.
GO TO 30
10 IF (JSTRUS.EQ.1) GO TO 20
J = 3
B1(IB) = BB1
B0(IB) = BB0
GO TO 30
20 IF (ISTRU.LT.25) GO TO 80
J = 3
B1(IB) = 0.
B0(IB) = 1.
30 GO TO (40, 50, 60), ISCAL
40 CALL SMIMP(QMAX)
GO TO 70
50 CALL SMAX(QMAX, 8, 9, 1)
GO TO 70
60 CALL SPOW(QMAX)
QMAX = SQRT(QMAX)
70 BMAX = AMAX1(BMAX,QMAX)
GO TO (90, 10, 80), J
80 B1(IB) = BB1
B0(IB) = BB0
90 SCA = SCALM/BMAX
IF (JSTRUS.EQ.1) GO TO 100
IF (JSTRUD.EQ.2) GO TO 100
IF (AMAX*SCA.GT.SCALM) SCA = SCALM/AMAX
C
C SCALING WITH A POWER OF TWO
C
100 IF (LPOT2.NE.(-1)) GO TO 110
I = ALOG(SCA)/ALOG(2.)
Q = 2.**I
IF (Q.GT.SCA) Q = Q/2.
SCA = Q
GO TO 120
110 IF (IWLR.GE.100) GO TO 120
BB0 = 2.*ALSBI
Q = SCA*BB0
SCA = AINT(Q)/BB0
120 AMAX = BMAX*SCA
FACT = FACT*SCA
NB = NNB
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: SCAL02
C CALCULATE SCALING FACTOR INTO NUMERATOR
C-----------------------------------------------------------------------
C
SUBROUTINE SCAL02
C
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
* LPOT2, JSTRU
COMMON /CREST1/ JSTRUS, JSTRUD
C
I = IB
IF (JSTRUD.EQ.1) GO TO 10
IF (I.EQ.1) GO TO 20
I = I - 1
10 B2(I) = SCA
B1(I) = B1(I)*SCA
B0(I) = B0(I)*SCA
FACT = FACT/SCA
20 RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: ALNS01
C ALLOCATION OF THE NOISE SOURCES FOR THE MODIFIED STRUCTURES OF THE
C FIRST AND SECOND CANONIC FORMS
C-----------------------------------------------------------------------
C
SUBROUTINE ALNS01
C
COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
* LPOT2, JSTRU
COMMON /CREST1/ JSTRUS, JSTRUD
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
COMMON /CNFUNC/ AQC(5), BN2(5), BN1(5), BN0(5)
C
DO 10 I=1,5
AQC(I) = 1.
10 CONTINUE
C
BB2 = SCA
IF (IB.GT.NB) GO TO 120
CALL QUAN(C1(IB), AQC1)
CALL QUAN(C0(IB), AQC0)
BB1 = B1(IB)
BB0 = B0(IB)
C
J = ISTRU/10
IF (J.EQ.2) GO TO 130
BB2 = BB2*B2(IB)
CALL QUAN(BB2, AQB2)
IF (JSTRU.EQ.2) GO TO 20
BB1 = BB1*SCA
BB0 = BB0*SCA
20 CALL QUAN(BB1, AQB1)
CALL QUAN(BB0, AQB0)
J = ISTRU - 10
IIS = 1
GO TO (30, 30, 50, 50, 60, 70, 90, 100), J
30 AQC(5) = AQC0
AQC(4) = AQC1
AQC(3) = AQB0
AQC(2) = AQB1
40 AQC(1) = AQB2
GO TO 270
50 AQC(4) = AMIN1(AQB1,AQC1)
AQC(5) = AMIN1(AQB0,AQC0)
GO TO 40
60 AQC(3) = AQB0
GO TO 80
70 AQC(1) = AQB2
80 AQC(2) = AQB1
AQC(4) = AQC1
AQC(5) = AQC0
GO TO 270
90 AQC(5) = AMIN1(AQB0,AQC0)
GO TO 110
100 AQC(5) = AQC0
110 AQC(4) = AMIN1(AQB1,AQC1)
GO TO 40
C
120 AQC1 = 1.
AQC0 = 1.
AQB1 = 1.
AQB0 = 1.
GO TO 140
130 J = ISTRU - 20
CALL QUAN(BB1, AQB1)
CALL QUAN(BB0, AQB0)
140 IIS = 2
IF (IB.NE.1) BB2 = BB2*B2(IB-1)
CALL QUAN(BB2, AQB2)
GO TO (170, 150, 210, 200, 190, 160, 210, 220), J
150 AQC(3) = AQB0
160 AQC(2) = AQB1
170 AQC(1) = AQB2
180 AQC(4) = AQC1
AQC(5) = AQC0
GO TO 230
190 IF (IB.EQ.1) GO TO 170
GO TO 180
200 AQC(2) = AMIN1(AQB1,AQB0)
210 AQC(4) = AMIN1(AQB2,AQC1,AQC0)
IF (IB.LE.NB) GO TO 230
AQC(1) = AQC(4)
AQC(4) = 1.
GO TO 270
220 AQC(2) = AQB1
GO TO 210
C
230 IF (IB.EQ.1) GO TO 270
BB1 = B1(IB-1)
BB0 = B0(IB-1)
IF (JSTRUS.EQ.2) GO TO 240
BB1 = BB1*SCA
BB0 = BB0*SCA
240 CALL QUAN(BB1, AQB1)
CALL QUAN(BB0, AQB0)
GO TO (250, 270, 260, 270, 250, 270, 260, 270), J
250 AQC(3) = AQB0
AQC(2) = AQB1
GO TO 270
260 AQC(4) = AMIN1(AQC(4),AQB1,AQB0)
C
270 DO 460 I=1,5
BN2(I) = 0.
BN1(I) = 0.
BN0(I) = 0.
IF (AQC(I).EQ.1.) GO TO 460
IF (IB.GT.NB) GO TO 420
IF (IIS.EQ.2) GO TO 340
GO TO (280, 290, 300, 310, 320, 330, 300, 310), J
280 GO TO (420, 430, 440, 430, 440), I
290 GO TO (400, 430, 440, 430, 440), I
300 GO TO (420, 460, 460, 430, 440), I
310 GO TO (400, 460, 460, 430, 440), I
320 GO TO (460, 430, 410, 430, 440), I
330 GO TO (400, 430, 460, 430, 440), I
C
340 GO TO (390, 350, 360, 370, 390, 380, 360, 370), J
350 GO TO (390, 450, 450, 390, 390), I
360 GO TO (460, 460, 460, 390, 460), I
370 GO TO (460, 450, 460, 390, 460), I
380 GO TO (390, 450, 460, 390, 390), I
C
390 BN2(I) = B2(IB)
BN1(I) = B1(IB)
BN0(I) = B0(IB)
GO TO 460
400 BN2(I) = 1.
BN1(I) = B1(IB)/B2(IB)
BN0(I) = B0(IB)/B2(IB)
GO TO 460
410 BN0(I) = 1.
420 BN2(I) = 1.
GO TO 460
430 BN1(I) = 1.
GO TO 460
440 BN0(I) = 1.
GO TO 460
450 BN2(I) = B2(IB)
BN1(I) = B2(IB)*C1(IB)
BN0(I) = B2(IB)*C0(IB)
460 CONTINUE
C
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: ALNS02
C ALLOCATION OF THE TRANSFER FUNCTIONS FROM THE INPUT
C TO THE NONLINEARITIES
C-----------------------------------------------------------------------
C
SUBROUTINE ALNS02
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
* LPOT2, JSTRU
COMMON /CREST1/ JSTRUS, JSTRUD
COMMON /CRENO/ LCNO, ICNO(16,5)
COMMON /CNFUNC/ AQC(5), BN2(5), BN1(5), BN0(5)
COMMON /CFFUNC/ PHI(5), BF2(5), BF1(5), BF0(5), IBB(5), ICOR
C
NNB = NB
IF (IB.GT.NNB) GO TO 10
BB2 = B2(IB)
BB1 = B1(IB)
BB0 = B0(IB)
10 ICOR = 0
IS = ISTRU - 10
IIS = 1
IF (IS.LE.8) GO TO 20
IIS = 2
IS = IS - 10
C
20 DO 450 N=1,5
BF2(N) = 0.
BF1(N) = 0.
BF0(N) = 0.
IBB(N) = 0
AQ = AQC(N)
IF (AQ.EQ.1.) GO TO 400
IF (LCNO.EQ.(-1)) AQ = 0.
K = ICNO(IB,N)
GO TO (400, 30, 400, 400, 40, 50, 400), K
30 AQ = SQRT(3.)*AQ
GO TO 60
40 AQ = -SQRT(3.)
GO TO 60
50 AQ = SQRT(3.)*(AQ-1.)
60 IF (AQ.EQ.0.) GO TO 440
IBB(N) = IB
IF (IB.GT.NNB) GO TO 380
IF (IIS.EQ.2) GO TO 120
GO TO (70, 70, 80, 90, 100, 110, 90, 90), IS
70 GO TO (200, 210, 220, 240, 250), N
80 GO TO (260, 400, 400, 280, 290), N
90 GO TO (200, 400, 400, 280, 290), N
100 GO TO (400, 210, 220, 240, 250), N
110 GO TO (200, 210, 400, 240, 250), N
C
120 GO TO (130, 140, 150, 160, 170, 180, 150, 190), IS
130 GO TO (300, 310, 320, 340, 350), N
140 GO TO (330, 360, 370, 340, 350), N
150 GO TO (330, 400, 400, 380, 400), N
160 GO TO (330, 390, 400, 380, 400), N
170 GO TO (330, 310, 320, 340, 350), N
180 GO TO (330, 360, 400, 340, 350), N
190 GO TO (330, 360, 400, 380, 400), N
C
200 Q = BB2
GO TO 230
210 Q = BB1
GO TO 230
220 Q = BB0
230 BF2(N) = Q
BF1(N) = Q*C1(IB)
BF0(N) = Q*C0(IB)
GO TO 410
240 Q = -C1(IB)
GO TO 270
250 Q = -C0(IB)
GO TO 270
260 Q = 1.
270 BF2(N) = Q*BB2
BF1(N) = Q*BB1
BF0(N) = Q*BB0
GO TO 410
280 BF2(N) = BB1 - BB2*C1(IB)
BF1(N) = BB0 - BB2*C0(IB)
GO TO 410
290 BF2(N) = BB0 - BB2*C0(IB)
BF1(N) = BB0*C1(IB) - BB1*C0(IB)
GO TO 410
300 IF (IB.EQ.1) GO TO 330
IQ = IB - 1
IBB(N) = IQ
BF2(N) = B2(IQ)
GO TO 410
310 IQ = IB - 1
IBB(N) = IQ
BF1(N) = B1(IQ)
GO TO 410
320 IQ = IB - 1
IBB(N) = IQ
BF0(N) = B0(IQ)
GO TO 410
330 Q = 1.
GO TO 230
340 BF1(N) = -C1(IB)
GO TO 410
350 BF0(N) = -C0(IB)
GO TO 410
360 BF1(N) = BB1/BB2
GO TO 410
370 BF0(N) = BB0/BB2
GO TO 410
380 BF2(N) = 1.
GO TO 410
390 BF1(N) = BB1/BB2
GO TO 370
C
400 AQ = 0.
GO TO 440
C
410 NB = IBB(N)
IF (IB.GT.NNB) GO TO 430
IF (NB.EQ.IB) GO TO 420
BB22 = B2(NB)
BB12 = B1(NB)
BB02 = B0(NB)
420 B2(NB) = BF2(N)
B1(NB) = BF1(N)
B0(NB) = BF0(N)
430 CALL SPOW(Q)
PHI(N) = SQRT(Q)
ICOR = 1
IF (NB.EQ.IB) GO TO 440
B2(NB) = BB22
B1(NB) = BB12
B0(NB) = BB02
C
440 BF2(N) = BF2(N)*AQ
BF1(N) = BF1(N)*AQ
BF0(N) = BF0(N)*AQ
450 CONTINUE
C
IF (IB.GT.NNB) GO TO 460
B2(IB) = BB2
B1(IB) = BB1
B0(IB) = BB0
460 NB = NNB
C
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: TSTR01
C TEST ROUTINE FOR THE STRUCTURES NOS. 11-28
C-----------------------------------------------------------------------
C
SUBROUTINE TSTR01
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /COPTST/ LOPTS, ISTOR
COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
* LPOT2, JSTRU
COMMON /CREST1/ JSTRUS, JSTRUD
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
C
IF (LOPTS.NE.1) GO TO 10
IF (ISTRU.GT.20 .AND. JSTRUS.EQ.2) GO TO 30
10 L = ISTRU
IF (L.GT.20) L = L - 10
IF (L.LT.15) GO TO 60
DO 20 L=1,NB
Q = ABS(B0(L))
IF (Q.LT.FLMI) Q = ABS(B1(L))
IF ((B2(L)-Q).GT.FLMI) GO TO 40
20 CONTINUE
GO TO 60
30 L = 31
GO TO 50
40 L = 32
50 CALL ERROR(L)
60 RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: PCORP
C COMPUTATION OF THE CORRELATED NOISE OF THE BLOCK IB BY A
C LINEARIZED MODEL
C-----------------------------------------------------------------------
C
SUBROUTINE PCORP
C
COMPLEX ZBL1, ZQQ, ZOM, Z, ZA, ZAA
COMPLEX ZPHI(300)
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
COMMON /CPOW/ PNU, PNC, AND, ITCORP
COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
COMMON /CNFUNC/ AQC(5), BN2(5), BN1(5), BN0(5)
COMMON /CFFUNC/ PHI(5), BF2(5), BF1(5), BF0(5), IBB(5), ICOR
COMMON /COPST2/ PN(100,2), TF(100,2), TFA(100,2)
C
EQUIVALENCE (ZPHI(1),PN(1,1))
C
ZBL1(ZQQ,A2,A1,A0) = (ZQQ*A2+A1)*ZQQ + A0
C
FAC = PNC
IF (ITCORP.GE.0) GO TO 20
DO 10 I=1,300
ZPHI(I) = CMPLX(0.,0.)
10 CONTINUE
ITCORP = 1
C
20 PNC = 0.
IF (ICOR.EQ.0) GO TO 80
ADOM = PI/299.
DO 70 LL=1,300
OM = ADOM*FLOAT(LL-1)
ZOM = CMPLX(COS(OM),SIN(OM))
Z = CMPLX(FACT,0.)
C
DO 30 I=1,NB
IF (I.NE.IB) Z = Z*ZBL1(ZOM,B2(I),B1(I),B0(I))
Z = Z/ZBL1(ZOM,1.,C1(I),C0(I))
30 CONTINUE
C
IF (IB.LE.NB) Z = Z/ZBL1(ZOM,1.,C1(IB),C0(IB))
ZA = CMPLX(0.,0.)
DO 60 N=1,5
IF (IBB(N).EQ.0) GO TO 60
IF (IB.GT.NB) GO TO 50
ZAA = ZBL1(ZOM,BF2(N),BF1(N),BF0(N))*ZBL1(ZOM,BN2(N),BN1(N),
* BN0(N))
ZAA = ZAA/PHI(N)
IF (IBB(N).EQ.IB) GO TO 40
ZAA = ZAA*ZBL1(ZOM,1.,C1(IB),C0(IB))/ZBL1(ZOM,B2(IB),B1(IB),
* B0(IB))
40 ZA = ZA + ZAA
GO TO 60
50 ZA = ZA + BF2(N)
60 CONTINUE
C
Z = Z*ZA
IF (ITCORP.EQ.1) ZPHI(LL) = ZPHI(LL) + Z*FAC
Q = CABS(Z)
PNC = PNC + Q*Q
70 CONTINUE
C
PNC = PNC/PI/150.
80 RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: TCORP
C COMPUTATION OF THE TOTAL CORRELATED NOISE BY A LINEARIZED MODEL
C-----------------------------------------------------------------------
C
SUBROUTINE TCORP
C
COMPLEX ZPHI(300)
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /CPOW/ PNU, PNC, AND, ITCORP
COMMON /COPST2/ PN(100,2), TF(100,2), TFA(100,2)
COMMON /CRENO/ LCNO, ICNO(16,5)
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
C
EQUIVALENCE (ZPHI(1),PN(1,1))
C
PNC = 0.
C
DO 20 I=1,NB
DO 10 J=1,5
K = ICNO(I,J)
IF (K.EQ.2) GO TO 30
IF (K.EQ.5) GO TO 30
IF (K.EQ.6) GO TO 30
10 CONTINUE
20 CONTINUE
GO TO 50
C
30 DO 40 I=1,300
Q = CABS(ZPHI(I))
PNC = PNC + Q*Q
40 CONTINUE
C
PNC = PNC/PI/150.
50 RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: DESCAL
C REMOVE SCALING OF SECOND-ORDER BLOCKS
C-----------------------------------------------------------------------
C
SUBROUTINE DESCAL
C
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
C
DO 10 I=1,NB
Q = B2(I)
B2(I) = 1.
B1(I) = B1(I)/Q
B0(I) = B0(I)/Q
FACT = FACT*Q
10 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: DENORM
C TRANSFORM STRUCTURE TO STRUCTURE WITHOUT SEPARATE SCALING FACTOR
C-----------------------------------------------------------------------
C
SUBROUTINE DENORM
C
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
C
DO 10 I=1,NB
Q = B2(I)
B1(I) = B1(I)*Q
B0(I) = B0(I)*Q
10 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: ALLO01
C ALLOCATION OF NUMERATORS AND DENOMINATORS ACCORDING TO THE
C ALLOCATION FIELD ISEQ(I,J)
C-----------------------------------------------------------------------
C
SUBROUTINE ALLO01
C
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
COMMON /SFILT/ NBS, SFACT, SB2(16), SB1(16), SB0(16), SC1(16),
* SC0(16)
COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
* LPOT2, JSTRU
C
NB = NBS
FACT = SFACT
DO 10 I=1,NB
J = ISEQ(I,1)
B2(I) = SB2(J)
B1(I) = SB1(J)
B0(I) = SB0(J)
J = ISEQ(I,2)
C1(I) = SC1(J)
C0(I) = SC0(J)
10 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: ALLO02
C ALLOCATION OF NUMERATORS AND DENOMINATORS ACCORDING TO THE
C ALLOCATION FIELD KSEQ(I,J), FOR A GIVEN START ORDERING
C-----------------------------------------------------------------------
C
SUBROUTINE ALLO02
C
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
COMMON /SFILT/ NBS, SFACT, SB2(16), SB1(16), SB0(16), SC1(16),
* SC0(16)
COMMON /CREST2/ KSEQ(16,2)
C
NB = NBS
FACT = SFACT
DO 10 I=1,NB
J = KSEQ(I,1)
B2(I) = SB2(J)
B1(I) = SB1(J)
B0(I) = SB0(J)
J = KSEQ(I,2)
C1(I) = SC1(J)
C0(I) = SC0(J)
10 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: FIXPAR
C FIXED PAIRING OF POLES AND ZEROS
C USED TOGETHER WITH A FREE INPUT OF THE COEFFICIENTS
C-----------------------------------------------------------------------
C
SUBROUTINE FIXPAR
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
COMMON /SFILT/ NBS, SFACT, SB2(16), SB1(16), SB0(16), SC1(16),
* SC0(16)
COMMON /SCRAT/ ADUM(32)
C
DIMENSION IPOL(16,2)
EQUIVALENCE (ADUM(1),IPOL(1,1))
C
DO 10 I=1,NB
IPOL(I,1) = 0
IPOL(I,2) = 0
10 CONTINUE
C
DO 30 I=1,NB
PM = 0.
DO 20 J=1,NB
IF (IPOL(J,1).NE.0) GO TO 20
P2O = POW2O(1.0,0.,0.,C1(J),C0(J))
IF (P2O.LE.PM) GO TO 20
PM = P2O
IPM = J
20 CONTINUE
IPOL(IPM,1) = I
30 CONTINUE
C
DO 70 I=1,NB
DO 40 K=1,NB
IF (IPOL(K,1).NE.I) GO TO 40
KK = K
GO TO 50
40 CONTINUE
C
50 CC0 = C0(KK)
CC1 = C1(KK)
DO 60 J=1,NB
IF (IPOL(J,2).NE.0) GO TO 60
P2O = POW2O(1.0,BB1,BB0,CC1,CC0)
PM = P2O
IZM = J
60 CONTINUE
IPOL(IZM,2) = KK
70 CONTINUE
C
CALL COPY01
DO 80 I=1,NB
IPM = IPOL(I,2)
B2(IPM) = SB2(I)
B1(IPM) = SB1(I)
B0(IPM) = SB0(I)
80 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: POLLOC
C CALCULATE THE POLES IN THE CLOSED UPPER Z-DOMAIN
C-----------------------------------------------------------------------
C
SUBROUTINE POLLOC
C
COMPLEX ZP, ZPS
C
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
COMMON /CPOL/ ZP(16,2), ZPS(16,2)
C
DO 20 I=1,NB
Q = -C1(I)/2.
QN = C0(I)
QN = Q*Q - QN
IF (QN.LT.0.) GO TO 10
QN = SQRT(QN)
ZPS(I,1) = CMPLX(Q+QN,0.)
ZPS(I,2) = CMPLX(Q-QN,0.)
GO TO 20
10 QN = SQRT(-QN)
ZPS(I,1) = CMPLX(Q,QN)
ZPS(I,2) = CMPLX(Q,-QN)
20 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C FUNCTION: POW2O
C NOISE POWER OF A SECOND-ORDER BLOCK
C-----------------------------------------------------------------------
C
FUNCTION POW2O(B2, B1, B0, C1, C0)
C
BB0 = B2*B2 + B1*B1 + B0*B0
BB1 = 2.*B1*(B2+B0)
BB2 = 2.*B2*B0
E1 = 1. + C0
Q = BB0*E1 - BB1*C1 + BB2*(C1*C1-C0*E1)
Q = Q/((1.-C0*C0)*E1-(C1-C1*C0)*C1)
POW2O = Q
RETURN
END
C
C ======================================================================
C
C DOREDI - SUBROUTINES: PART IV
C
C ======================================================================
C
C
C-----------------------------------------------------------------------
C SUBROUTINE: DEFIN1
C DEFAULT VALUES FOR CLASS 01 INPUT DATA
C-----------------------------------------------------------------------
C
SUBROUTINE DEFIN1
C
COMMON /CONTR/ IPRUN, IPCON, NINP, NOUT, NDOUT, LSPOUT, NSPOUT
COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
COMMON /CGRID/ GR(64), NGR(12)
C
NINP = 0
NOUT = 3
NDOUT = 0
LSPOUT = 0
NSPOUT = 0
DO 10 I=1,12
NGR(I) = 0
10 CONTINUE
DO 20 I=1,4
ROM(I) = 0.
20 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: INP010
C COMMAND CARDS FOR INPUT AND OUTPUT DIRECTIVES
C PROCESSING OF THE COMMAND CARDS CLASS 01
C-----------------------------------------------------------------------
C
SUBROUTINE INP010
C
COMMON /CONTR/ IPRUN, IPCON, NINP, NOUT, NDOUT, LSPOUT, NSPOUT
COMMON /INPDAT/ ICODE, JCODE, LINP, IINP, AINP(3), NPAR
COMMON /CGRID/ GR(64), NGR(12)
C
J = JCODE/10
IF (J.EQ.0) GO TO 110
IF (J.GT.4) GO TO 110
GO TO (10, 80, 90, 100), J
C
C *0110 B,NINP
C
10 IF (IINP.GT.4) GO TO 120
GO TO (20, 30, 40, 50), IINP
20 IF (NINP.EQ.0) GO TO 130
GO TO 70
30 CALL INP012
GO TO 70
40 NINP = IINP
CALL INP013
GO TO 60
50 NINP = IINP
CALL INP014
60 NGR(1) = 0
RETURN
70 NINP = IINP
RETURN
C
C *0120 B,NOUT
C
80 IF (IINP.LT.0 .OR. IINP.GT.5) GO TO 120
NOUT = IINP
RETURN
C
C *0130 NDOUT
C
90 NDOUT = LINP
RETURN
C
C *0140 LSPOUT,ISPOUT
C
100 LSPOUT = LINP
NSPOUT = IINP
RETURN
C
110 I = 1
GO TO 140
120 I = 4
GO TO 140
130 I = 12
140 CALL ERROR(I)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: INP012
C INPUT OF THE FILTER DESCRIPTION AND DATA FROM DISK: (CHANNEL KA4)
C-----------------------------------------------------------------------
C
SUBROUTINE INP012
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
COMMON /RESZ/ ZFA, ZM(18,4), ZZR(16), ZZI(16), ZPR(16), ZPI(16)
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
COMMON /FILTRE/ IRCO(5), IECO(16,5), IDCO(16,5), IECOM
COMMON /CRECO/ JRCO, JECO, JDCO, JJDCO(5), JMAXV, JTRB2, LREF
COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
* ISTAB, IDEPSL, IDEPSD
C
9999 FORMAT (10I5)
9998 FORMAT (4E15.7)
C
DO 10 I=1,3
READ (KA4,9999)
10 CONTINUE
READ (KA4,9999) ITYP, IAPRO, NDEG
READ (KA4,9998) SF, OM, ADELP, ADELS
READ (KA4,9998) AC, ROM, RDELP, RDELS
READ (KA4,9999) NZM, M
DO 20 I=1,M
READ (KA4,9998) (ZM(I,J),J=1,4)
20 CONTINUE
READ (KA4,9999)
READ (KA4,9999) NB
READ (KA4,9998) FACT
DO 30 I=1,NB
READ (KA4,9998) B2(I), B1(I), B0(I), C1(I), C0(I)
30 CONTINUE
READ (KA4,9999) IRCO
DO 40 I=1,NB
READ (KA4,9999) (IECO(I,J),J=1,5)
40 CONTINUE
DO 50 I=1,NB
READ (KA4,9999) (IDCO(I,J),J=1,5)
50 CONTINUE
READ (KA4,9999) IWL, IECOM, JMAXV, JTRB2
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: INP013
C INPUT OF SECOND-ORDER BLOCKS
C-----------------------------------------------------------------------
C
SUBROUTINE INP013
C
COMMON /CARD/ LCODE, IZ(40)
COMMON /INPDAT/ ICODE, JCODE, LINP, IINP, AINP(3), NPAR
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
C
CALL INP001
IF (LCODE.NE.0) RETURN
CALL INP002
IF (ICODE.NE.0 .OR. JCODE.NE.10) CALL ERROR(2)
NB = IINP
FACT = AINP(1)
C
MA = 2*NB
DO 20 I=1,MA
CALL INP001
IF (LCODE.NE.0) RETURN
CALL INP002
IF (ICODE.NE.0 .OR. JCODE.NE.10) CALL ERROR(2)
IF (I.GT.NB) GO TO 10
B2(I) = AINP(1)
B1(I) = AINP(2)
B0(I) = AINP(3)
GO TO 20
10 J = I - NB
C1(J) = AINP(2)
C0(J) = AINP(3)
20 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: INP014
C INPUT OF POLES AND ZEROS
C-----------------------------------------------------------------------
C
SUBROUTINE INP014
C
COMMON /CARD/ LCODE, IZ(40)
COMMON /INPDAT/ ICODE, JCODE, LINP, IINP, AINP(3), NPAR
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
C
CALL INP001
IF (LCODE.NE.0) RETURN
CALL INP002
IF (ICODE.NE.0 .OR. JCODE.NE.10) CALL ERROR(2)
NB = IINP
FACT = AINP(1)
NN = 0
JJ = NB
C
10 CALL INP001
IF (LCODE.NE.0) RETURN
CALL INP002
IF (ICODE.NE.0 .OR. JCODE.NE.10) CALL ERROR(2)
QR = AINP(1)
QI = AINP(2)
IF (QI.EQ.0.) GO TO 20
A0 = QR*QR + QI*QI
A1 = -2.*QR
GO TO 40
20 IF (JJ.NE.NB) GO TO 30
JJ = NB - 1
R = QR
GO TO 10
30 JJ = NB
A0 = R*QR
A1 = -R - QR
40 NN = NN + 1
IF (K.EQ.2) GO TO 50
B2(NN) = 1.
B1(NN) = A1
B0(NN) = A0
GO TO 60
50 C1(NN) = A1
C0(NN) = A0
60 IF (NN.LT.JJ) GO TO 10
IF (JJ.EQ.NB) GO TO 70
A1 = -R
A0 = 0.
JJ = JJ + 1
GO TO 40
70 IF (K.EQ.2) RETURN
K = 2
NN = 0
GO TO 10
C
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: DEFIN2
C DEFAULT VALUES FOR THE COMMAND CARDS OF CLASS 02
C-----------------------------------------------------------------------
C
SUBROUTINE DEFIN2
C
COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
* LPOT2, JSTRU
COMMON /CREST2/ KSEQ(16,2)
COMMON /COPTST/ LOPTS, ISTOR
C
LOPTW = 0
LSTAB = -1
IWLR = 100
ACXMI = 0.
ACXMA = 1.
LOPTS = 0
ISTOR = 100
LSEQ = 0
ISCAL = 2
SCALM = 1.
ITERM = 10
ITERM1 = 4
DO 20 J=1,2
DO 10 I=1,MBL
KSEQ(I,J) = I
ISEQ(I,J) = I
10 CONTINUE
20 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: INP020
C COMMAND CARDS FOR CONTROLLING THE OPTIMIZATION
C PROCESSING OF THE COMMAND CARDS OF CLASS 02
C-----------------------------------------------------------------------
C
SUBROUTINE INP020
C
COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
COMMON /INPDAT/ ICODE, JCODE, LINP, IINP, AINP(3), NPAR
COMMON /INTDAT/ JINP(16)
COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
* ISTAB, IDEPSL, IDEPSD
COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
* LPOT2, JSTRU
COMMON /CREST2/ KSEQ(16,2)
COMMON /COPTST/ LOPTS, ISTOR
C
IF (JCODE.LT.10) GO TO 110
J = JCODE/10
IF (J.GT.4) GO TO 110
GO TO (10, 50, 60, 70), J
C
C *0210 LWLF,ITERM,ACX,ACXMI,ACXMA
C
10 LOPTW = -1
IF (LINP.EQ.(-1)) LOPTW = 1
IF (IINP.EQ.0) GO TO 20
IF (IINP.LT.0) GO TO 120
ITERM = IINP
20 IF (NPAR.EQ.0) GO TO 30
Q = AINP(1)
IF (Q.LT.0. .OR. Q.GT.1.) GO TO 120
ACX = Q
30 Q = AINP(2)
IF (Q.EQ.0.) GO TO 40
IF (Q.LT.0. .OR. Q.GE.1.) GO TO 120
ACXMI = Q
40 Q = AINP(3)
IF (Q.EQ.0.) RETURN
IF (Q.LE.0. .OR. Q.GT.1.) GO TO 120
ACXMA = Q
RETURN
C
C *0220 B,ITERM1
C
50 IF (IINP.LT.0) GO TO 120
ITERM1 = IINP
RETURN
C
C *0230 PAIRF,ISTOR
C
60 LOPTS = -1
IF (LINP.EQ.(-1)) LOPTS = 1
IF (IINP.GT.100) GO TO 120
IF (IINP.LT.0) GO TO 120
IF (IINP.NE.0) ISTOR = IINP
RETURN
C
C *0240 LSEQ,ISCAL,SCALM
C
70 IF (LINP.NE.(-1)) GO TO 100
LSEQ = -1
DO 90 J=1,2
CALL INP001
IF (ICODE.NE.0) GO TO 130
IF (JCODE.NE.40) GO TO 130
CALL INP003
DO 80 I=1,MBL
KSEQ(I,J) = JINP(I)
80 CONTINUE
90 CONTINUE
C
100 IF (IINP.NE.0) ISCAL = IINP
IF (ISCAL.LT.0) ISCAL = 0
IF (AINP(1).NE.0.) SCALM = AINP(1)
RETURN
C
110 I = 4
GO TO 140
120 I = 5
GO TO 140
130 I = 2
140 CALL ERROR(I)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: DEFIN3
C SET DEFAULT VALUES
C-----------------------------------------------------------------------
C
SUBROUTINE DEFIN3
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /TOLSN/ VSN, VD, A
COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
C
ITYP = 1
IAPRO = 4
NDEG = 0
NDEGF = 0
MDEG = 0
NORMA = 0
LSOUT = 0
C
DO 10 I=1,4
OM(I) = 0.
10 CONTINUE
SF = 0.
A = 0.
ADELP = 0.
ADELS = 1./FLMA
VSN = 0.
VD = 0.
ACX = 1000.
EDEG = ACX
C
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: INP030
C INPUT ROUTINE FOR THE FILTER TYPE AND FOR THE TOLERANCE SCHEME
C PROCESSING OF THE COMMAND CARDS OF CLASS 03
C-----------------------------------------------------------------------
C
SUBROUTINE INP030
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /INPDAT/ ICODE, JCODE, LINP, IINP, AINP(3), NPAR
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /TOLSN/ VSN, VD, A
COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
C
SOM(AA) = 2.*ATAN(AA)
C
J = JCODE/10
IF (J.GT.7) GO TO 240
GO TO (10, 20, 30, 40, 170, 220, 230), J
C
C *0310 B,FILTER TYPE,SAMPLING FREQUENCY
C
10 IF (IINP.NE.0) ITYP = IINP
IF (AINP(1).NE.0.) SF = AINP(1)
RETURN
C
C *0320 B,IAPRO,,ACX
C
20 IF (IINP.NE.0) IAPRO = IINP
IF (NPAR.NE.0) ACX = AINP(1)
RETURN
C
C *0330 NDEGF,EDEG
C
30 IF (LINP.EQ.(-1)) NDEGF = -1
IF (IINP.NE.0) NDEG = IINP
IF (NPAR.NE.0) EDEG = AINP(1)
RETURN
C
C CUTOFF FREQUENCIES OF THE TOLERANCE SCHEME
C
40 K = JCODE - J*10
IF (K.GT.7) GO TO 250
GO TO (50, 50, 100, 110, 120, 140, 160), K
C
C *0341 B,B,FR(1),FR(2), :: FREQUENCY DOMAIN
C *0342 B,B,FR(3),FR(4),
C
50 IF (AINP(3).EQ.0.) GO TO 70
IF (SF.NE.0.) GO TO 260
60 SF = AINP(3)
70 IF (SF.EQ.0.) GO TO 270
Q = 2.*PI/SF
GO TO (80, 90), K
80 OM(1) = AINP(1)*Q
OM(2) = AINP(2)*Q
RETURN
C
90 OM(3) = AINP(1)*Q
OM(4) = AINP(2)*Q
RETURN
C
C *0343 B,B,OM(1),OM(2) :: Z-DOMAIN
C
100 OM(1) = AINP(1)
OM(2) = AINP(2)
RETURN
C
C *0344 B,B,OM(3),OM(4) :: Z-DOMAIN SEC. CARD
C
110 OM(3) = AINP(1)
OM(4) = AINP(2)
RETURN
C
C *0345 B,B,S(1),S(2) :: S-DOMAIN
C
120 DO 130 I=1,2
OM(I) = SOM(AINP(I))
130 CONTINUE
RETURN
C
C *0346 B,B,S(3),S(4) :: S-DOMAIN SEC. CARD
C
140 DO 150 I=1,2
OM(I+2) = SOM(AINP(I))
150 CONTINUE
RETURN
C
C *0347 B,B,VSN,VD,A :: NORMALIZED PARAMETER
C
160 VSN = AINP(1)
VD = AINP(2)
A = AINP(3)
RETURN
C
C DELTA P, DELTA S
C
170 K = JCODE - J*10
IF (K.GT.3) GO TO 250
GO TO (180, 190, 210), K
C
C *0351 B,B,DELP,DELS
C
180 ADELP = AINP(1)
ADELS = AINP(2)
RETURN
C
C *0352 B,B,AP,AS (IN DB)
C
190 ADELP = 1. - 10.**(-0.05*AINP(1))
200 ADELS = 10.**(-0.05*AINP(2))
RETURN
C
C *0353 B,B,P,AS (AS IN DB)
C
210 Q = AINP(1)
ADELP = 1. - SQRT(1.-Q*Q)
GO TO 200
C
C *0360 B,NORMA
C
220 IF (IINP.LT.0 .OR. IINP.GT.3) GO TO 250
NORMA = IINP
RETURN
C
C *0370 LSOUT
C
230 LSOUT = LINP
RETURN
C
C OUTPUT OF ERROR MESSAGES
C
240 I = 4
GO TO 280
250 I = 5
GO TO 280
260 I = 7
GO TO 280
270 I = 6
280 CALL ERROR(I)
IF (I.EQ.7) GO TO 60
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: INP031
C CHECK ROUTINE FOR THE FILTER INPUT DATA
C CHECK DATA CARDS OF CLASS 03
C-----------------------------------------------------------------------
C
SUBROUTINE INP031
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /TOLSN/ VSN, VD, A
COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
C
IF (ITYP.GT.4) GO TO 90
IF (IAPRO.GT.4) GO TO 100
C
IF (EDEG.GT.999.) EDEG = 0.2
C
C CHECK OM(1) ... OM(4)
C
LVSN = 0
LSYM = 0
J = 0
Q = 0.
ME = 2
IF (ITYP.GE.3) ME = 4
DO 20 I=1,ME
QQ = OM(I)
IF (QQ.EQ.0.) GO TO 10
IF (QQ.GT.PI) GO TO 120
IF (QQ.LT.Q) GO TO 110
Q = QQ
GO TO 20
10 LSYM = LVSN
LVSN = I
J = J + 1
20 CONTINUE
IF (J.EQ.0) GO TO 80
IF (ITYP.GT.2) GO TO 30
IF (J.EQ.1) GO TO 70
GO TO 50
30 GO TO (40, 60, 50, 50), J
40 LSYM = LVSN
LVSN = 0
J = 0
GO TO 80
C
50 IF (VSN.EQ.0.) GO TO 120
IF (VD.EQ.0.) GO TO 120
IF (ITYP.GT.2 .AND. A.EQ.0.) GO TO 110
LVSN = -1
GO TO 80
C
60 J = 1
IQ = LSYM + LVSN
IF (IQ.LE.3 .OR. IQ.GE.7) GO TO 120
70 IF (NDEG.EQ.0) GO TO 120
C
C CHECK FOR PASSBAND AND STOPBAND RIPPLE
C
80 IF (ADELP.EQ.0.) J = J + 1
IF (ADELS.LE.(1./FLMA)) J = J + 1
IF (J.EQ.0) RETURN
IF (J.NE.1 .OR. NDEG.EQ.0) GO TO 130
RETURN
C
C OUTPUT OF ERROR MESSAGE
C
90 I = 20
GO TO 140
100 I = 21
GO TO 140
110 I = 22
GO TO 140
120 I = 23
GO TO 140
130 I = 24
140 CALL ERROR(I)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: DEFIN4
C SET DEFAULT VALUES FOR COMMAND CARDS OF CLASS 04
C-----------------------------------------------------------------------
C
SUBROUTINE DEFIN4
C
COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
COMMON /CRECO/ JRCO, JECO, JDCO, JJDCO(5), JMAXV, JTRB2, LREF
COMMON /FILTRE/ IRCO(5), IECO(16,5), IDCO(16,5), IECOM
COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
* ISTAB, IDEPSL, IDEPSD
C
LSTAB = 0
LREF = 0
IWL = 16
JRCO = 2
JECO = 0
JDCO = 0
JMAXV = -3
JTRB2 = 3
DO 20 J=1,5
IRCO(J) = 0
DO 10 I=1,MBL
IECO(I,J) = 0
IDCO(I,J) = -100
10 CONTINUE
20 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: INP040
C INPUT ROUTINE FOR THE SPECIFICATION OF THE COEFFICIENT REALIZATION
C-----------------------------------------------------------------------
C
SUBROUTINE INP040
C
COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
COMMON /CARD/ LCODE, IZ(40)
COMMON /INPDAT/ ICODE, JCODE, LINP, IINP, AINP(3), NPAR
COMMON /INTDAT/ JINP(16)
COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
COMMON /CRECO/ JRCO, JECO, JDCO, JJDCO(5), JMAXV, JTRB2, LREF
COMMON /FILTRE/ IRCO(5), IECO(16,5), IDCO(16,5), IECOM
COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
* ISTAB, IDEPSL, IDEPSD
C
J = JCODE/10
IF (J.EQ.0) GO TO 170
IF (J.GT.6) GO TO 170
GO TO (10, 20, 50, 90, 150, 160), J
C
C *0410 LSTAB,IWL
C
10 LSTAB = LINP
IF (IINP.NE.0) IWL = IINP
RETURN
C
C *0420 RCOINP,JRCO
C
20 IF (LINP.NE.(-1)) GO TO 40
CALL INP001
IF (LCODE.EQ.1) RETURN
CALL INP003
IF (ICODE.NE.0 .OR. JCODE.NE.20) CALL ERROR(2)
DO 30 I=1,5
IRCO(I) = JINP(I)
30 CONTINUE
JRCO = 0
IF (IINP.NE.0) GO TO 180
40 JRCO = IINP
RETURN
C
C *0430 ECOINP,JECO
C
50 IF (LINP.NE.(-1)) GO TO 80
JECO = 0
IECOM = 0
DO 70 J=1,5
CALL INP001
IF (LCODE.EQ.1) RETURN
CALL INP003
IF (ICODE.NE.0 .OR. JCODE.NE.30) CALL ERROR(2)
DO 60 I=1,MBL
IECOM = MAX0(IECOM,JINP(I))
IECO(I,J) = JINP(I)
60 CONTINUE
70 CONTINUE
RETURN
C
80 JECO = IINP
RETURN
C
C *0440 DCOINP,JDCO
C
90 IF (LINP.EQ.(-1)) GO TO 120
DO 110 J=1,5
II = IINP/10
M = 6 - J
JJ = IINP - 10*II
IINP = II
JJDCO(M) = JJ
IF (II.NE.0) GO TO 110
DO 100 I=1,16
IDCO(I,J) = -100
100 CONTINUE
110 CONTINUE
JDCO = -1
RETURN
C
120 JDCO = 0
DO 140 J=1,5
CALL INP001
IF (LCODE.EQ.1) RETURN
CALL INP003
IF (ICODE.NE.0 .OR. JCODE.NE.40) CALL ERROR(2)
DO 130 I=1,16
IDCO(I,J) = JINP(I)
130 CONTINUE
140 CONTINUE
RETURN
C
C *0450 B,JMAXV
C
150 JMAXV = IINP
RETURN
C
C *0460 B,JTRB2
C
160 JTRB2 = IINP
RETURN
C
C OUTPUT OF ERROR MESSAGES
C
170 I = 1
GO TO 190
180 I = 9
190 CALL ERROR(I)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: DEFIN5
C DEFAULT VALUES FOR COMMAND CARDS OF CLASS 05
C-----------------------------------------------------------------------
C
SUBROUTINE DEFIN5
C
COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
* LPOT2, JSTRU
COMMON /CREST1/ JSTRUS, JSTRUD
COMMON /CRENO/ LCNO, ICNO(16,5)
C
LCNO = 0
LPOT2 = 0
ISTRU = 13
JSTRU = 1
JSTRUS = 1
JSTRUD = 1
IWLR = 100
DO 20 J=1,5
DO 10 I=1,MBL
ICNO(I,J) = 3
10 CONTINUE
20 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: INP050
C INPUT DATA FOR THE STRUCTURE REALIZATION
C COMMAND CARDS OF CLASS 05
C-----------------------------------------------------------------------
C
SUBROUTINE INP050
C
COMMON /CONST1/ MAXDEG, IWLMI, IWLMA, MBL
COMMON /CARD/ LCODE, IZ(40)
COMMON /INPDAT/ ICODE, JCODE, LINP, IINP, AINP(3), NPAR
COMMON /INTDAT/ JINP(16)
COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
* LPOT2, JSTRU
COMMON /CREST1/ JSTRUS, JSTRUD
COMMON /CRENO/ LCNO, ICNO(16,5)
C
J = JCODE/10
IF (J.EQ.0) GO TO 100
IF (J.GT.3) GO TO 100
GO TO (10, 20, 90), J
C
C *0510 LPOT2,ISTRU
C
10 LPOT2 = LINP
IF (IINP.EQ.0) RETURN
IF (IINP.GE.30) GO TO 110
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C HERE INSERT THE PARAMETERS JSTRU,JSTRUS,JSTRUD FOR YOUR STRUCTURES
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
ISTRU = IINP
JSTRU = 1
IF (ISTRU/2.EQ.(ISTRU+1)/2) JSTRU = 2
JSTRUS = JSTRU
JSTRUD = 1
IF (ISTRU.GE.20) JSTRUD = 2
RETURN
C
C *0520 CNOINP,JCNO
C
20 IF (LINP.NE.(-1)) GO TO 50
DO 40 J=1,5
CALL INP001
IF (LCODE.NE.0) RETURN
CALL INP003
IF (ICODE.NE.0 .OR. JCODE.NE.20) CALL ERROR(2)
DO 30 I=1,MBL
ICNO(I,J) = JINP(I)
30 CONTINUE
40 CONTINUE
IF (IINP.NE.0) GO TO 130
RETURN
C
50 IF (IINP.EQ.0) GO TO 120
DO 80 J=1,5
II = IINP/10
M = 6 - J
JJ = IINP - 10*II
IINP = II
IF (JJ.NE.0) GO TO 60
IF (J.EQ.1) GO TO 130
JJ = ICNO(1,M+1)
60 DO 70 I=1,16
ICNO(I,M) = JJ
70 CONTINUE
80 CONTINUE
RETURN
C
C *0530 LCNO,IWLR
C
90 LCNO = LINP
IF (IINP.LT.0) GO TO 110
IF (IINP.NE.0) IWLR = IINP
RETURN
C
100 I = 1
GO TO 140
110 I = 5
GO TO 140
120 I = 9
GO TO 140
130 I = 4
140 CALL ERROR(I)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: DEFIN6
C DEFAULT VALUES FOR THE COMMAND CARDS OF CLASS 06
C-----------------------------------------------------------------------
C
SUBROUTINE DEFIN6
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
COMMON /COPTST/ LOPTS, ISTOR
COMMON /CPLOT/ IPLOT, LNOR, IPAG, OMLO, OMUP, RMAX, RMIN
C
LOPTS = 0
LOPTW = 0
IPOT = 0
LNOR = 0
IPAG = 2
OMLO = 0.
OMUP = PI
RMAX = 1.
RMIN = -1.
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: INP060
C COMMAND CARDS FOR THE ANALYSIS
C PROCESSING OF THE COMMAND CARDS OF CLASS 06
C-----------------------------------------------------------------------
C
SUBROUTINE INP060
C
COMMON /INPDAT/ ICODE, JCODE, LINP, IINP, AINP(3), NPAR
COMMON /COPTST/ LOPTS, ISTOR
COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
* ISTAB, IDEPSL, IDEPSD
COMMON /CPLOT/ IPLOT, LNOR, IPAG, OMLO, OMUP, RMAX, RMIN
C
J = JCODE/10
IF (J.EQ.0) GO TO 90
IF (J.GT.6) GO TO 90
GO TO (10, 20, 30, 60, 70, 80), J
C
C *0610 LWLM,IWL
C
10 LOPTW = 1
IF (LINP.EQ.(-1)) LOPTW = -1
IF (IINP.LT.0) GO TO 100
IF (IINP.NE.0) IWL = IINP
RETURN
C
C *0620 LSEQ,LSCAL,SCALM
C
20 LOPTS = -1
JCODE = 40
CALL INP020
RETURN
C
C *0630 LNOR,IPAG,OMLO,OMUP,RMAX
C
30 IPLOT = 1
40 LNOR = LINP
IF (IINP.GT.0) IPAG = IINP
IF (IPLOT.EQ.4) RETURN
IF (NPAR.NE.0) OMLO = AINP(1)
Q = AINP(2)
IF (Q.EQ.0.) GO TO 50
IF (Q.LE.OMLO) GO TO 110
OMUP = Q
50 IF (AINP(3).NE.0.) RMAX = AINP(3)
RETURN
C
C *0640 LNOR,IPAG,OMLO,OMUP,RMAX
C
60 IPLOT = 2
RMAX = 100.
GO TO 40
C
C *0650 LNOR,IPAG,OMLO,OMUP
C
70 IPLOT = 3
GO TO 40
C
C *0660 LNOR,IPAG,RMIN,RMAX
C
80 IPLOT = 4
RMAX = 1.
IF (NPAR.NE.0) RMIN = AINP(1)
Q = AINP(2)
IF (Q.EQ.0.) GO TO 40
IF (Q.LE.RMIN) GO TO 110
RMAX = Q
GO TO 40
C
90 I = 1
GO TO 120
100 I = 5
GO TO 120
110 I = 10
IPLOT = 0
120 CALL ERROR(I)
C
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: INP001
C READ COMMAND CARD
C
C IF 'DECODE' DOES NOT CONFORM TO YOUR COMPILER,
C SET THE ALTERNATIVE INPUT VERSION
C
C ELIMINATE THE STATEMENTS ENCLOSED BY THE COMMENT LINES C111
C
C AND REPLACE THE COMMENT MARKS 'C222' WITH ' '
C
C-----------------------------------------------------------------------
C
SUBROUTINE INP001
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /CARD/ LCODE, IZ(40)
COMMON /INPDAT/ ICODE, JCODE, LINP, IINP, AINP(3), NPAR
C111
COMMON /CIOFOR/ IOFO(3), ION
C111
C
DIMENSION IZ1(2)
DATA IPR, ISL /1H*,1H//
C
9999 FORMAT (1X, 4A1)
9998 FORMAT (1X//11H DATA INPUT)
C
IF (LCODE.EQ.(-1)) WRITE (KA2,9998)
C
10 CONTINUE
C111
READ (KA1,IOFO) (IZ(I),I=1,ION)
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C DECODE IS A MACHINE-DEPENDENT FUNCTION
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
DECODE(2,21,IZ) IZ1(1),IZ1(2)
C111
C
C222 READ (KA1,21) (IZ1(I),I=1,2)
21 FORMAT (2A1)
IF (IZ1(1).EQ.IPR) GO TO 30
IF (IZ1(1).EQ.ISL .AND. IZ1(2).EQ.ISL) GO TO 20
CALL ERROR(1)
GO TO 10
20 CONTINUE
C111
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C DECODE IS A MACHINE-DEPENDENT FUNCTION
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
DECODE (4,22,IZ) ICODE,JCODE
C111
C
C222 READ (KA1,22) ICODE,JCODE
22 FORMAT (2X, 2A1)
WRITE (KA2,9999) IZ1(1), IZ1(2), ICODE, JCODE
LCODE = 1
RETURN
30 LCODE = 0
ICODE = 0
RETURN
C
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: INP002
C DECODE COMMAND CARD
C
C IF 'DECODE' DOES NOT CONFORM TO YOUR COMPILER, SEE INP001
C
C-----------------------------------------------------------------------
C
SUBROUTINE INP002
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /CARD/ LCODE, IZ(40)
COMMON /INPDAT/ ICODE, JCODE, LINP, IINP, AINP(3), NPAR
COMMON /SCRAT/ ADUM(32)
DIMENSION IZZ(15)
EQUIVALENCE (ADUM(1),IZZ(1))
C
DATA IBL, IT /1H ,1HT/
C
NPAR = 0
C111
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C DECODE IS A MACHINE-DEPENDENT FUNCTION
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
DECODE (30,11,IZ) IZZ
11 FORMAT (15X, 15A1)
DO 10 I=1,15
IF (IZZ(I).NE.IBL) NPAR = 1
10 CONTINUE
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C DECODE IS A MACHINE-DEPENDENT FUNCTION
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
DECODE (60,12,IZ) ICODE,JCODE,JINP,IINP,AINP
C111
C
C222 READ (KA1,12) ICODE,JCODE,JINP,IINP,(AINP(I),I=1,3)
C222 NPAR=1
C
WRITE (KA2,9999) ICODE, JCODE, JINP, IINP, AINP
C
12 FORMAT (1X, 2I2, 2X, 1A1, I6, 1X, 3E15.0)
9999 FORMAT (2H *, 2I4, 2X, 1A1, I6, 1X, 3E15.6)
C
LINP = 0
IF (JINP.EQ.IT) LINP = -1
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: INP003
C INTEGER DATA INPUT
C
C
C IF 'DECODE' DOES NOT CONFORM TO YOUR COMPILER, SEE INP001
C-----------------------------------------------------------------------
C
SUBROUTINE INP003
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /CARD/ LCODE, IZ(40)
COMMON /INTDAT/ JINP(16)
C
12 FORMAT (1X, 2I2, 3X, 16I4)
9999 FORMAT (1H*, 2I3, 1X, 16I4)
C111
C
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C DECODE IS A MACHINE-DEPENDENT FUNCTION
C >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C
DECODE (72,12,IZ) ICODE,JCODE,(JINP(I),I=1,16)
C111
C
C222 READ (KA1,12) ICODE,JCODE,(JINP(I),I=1,16)
C
WRITE (KA2,9999) ICODE, JCODE, JINP
RETURN
END
C
C ======================================================================
C
C DOREDI - SUBROUTINES: PART V
C
C ======================================================================
C
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUT010
C LEADER TO OUT011 FOR STRUCTURES WITH SEPARATE SCALING
C-----------------------------------------------------------------------
C
SUBROUTINE OUT010
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
C
WRITE (KA2,9999)
9999 FORMAT (1X//46H STRUCTURE WITH SEPARATE SCALING FACTORS B2(L)//
* 45H FOR STRUCTURES 21 TO 28 B2(0) = GAIN FACTOR)
C
DO 10 I=1,NB
Q = B2(I)
B1(I) = B1(I)/Q
B0(I) = B0(I)/Q
10 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUT011
C OUTPUT OF BLOCKS OF SECOND ORDER
C-----------------------------------------------------------------------
C
SUBROUTINE OUT011
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
C
WRITE (KA2,9999) FACT
9999 FORMAT (//23H BLOCKS OF SECOND ORDER//24H CONSTANT GAIN FACTOR = ,
* E15.7//3H L, 5X, 5HB2(L), 8X, 5HB1(L), 8X, 5HB0(L), 10X,
* 5HC1(L), 8X, 5HC0(L)/)
WRITE (KA2,9998) (I,B2(I),B1(I),B0(I),C1(I),C0(I),I=1,NB)
9998 FORMAT (1X, I3, 3F13.8, 2X, 2F13.8)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUT012
C OUTPUT OF THE RESULTS TO THE DISK: (CHANNEL KA4)
C-----------------------------------------------------------------------
C
SUBROUTINE OUT012
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
COMMON /RESZ/ ZFA, ZM(18,4), ZZR(16), ZZI(16), ZPR(16), ZPI(16)
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
COMMON /FILTRE/ IRCO(5), IECO(16,5), IDCO(16,5), IECOM
COMMON /CRECO/ JRCO, JECO, JDCO, JJDCO(5), JMAXV, JTRB2, LREF
COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
* ISTAB, IDEPSL, IDEPSD
C
9999 FORMAT (10I5)
9998 FORMAT (4E15.7)
C
WRITE (KA4,9997)
9997 FORMAT (19H FILTER DESCRIPTION)
WRITE (KA4,9999) ITYP, IAPRO, NDEG
WRITE (KA4,9998) SF, OM, ADELP, ADELS
WRITE (KA4,9998) AC, ROM, RDELP, RDELS
M = 0
DO 10 I=1,4
M = MAX0(M,NZM(I))
10 CONTINUE
WRITE (KA4,9999) NZM, M
DO 20 I=1,M
WRITE (KA4,9998) (ZM(I,J),J=1,4)
20 CONTINUE
WRITE (KA4,9996)
9996 FORMAT (5H DATA)
WRITE (KA4,9999) NB
WRITE (KA4,9998) FACT
DO 30 I=1,NB
WRITE (KA4,9998) B2(I), B1(I), B0(I), C1(I), C0(I)
30 CONTINUE
WRITE (KA4,9999) IRCO
DO 40 I=1,NB
WRITE (KA4,9999) (IECO(I,J),J=1,5)
40 CONTINUE
DO 50 I=1,NB
WRITE (KA4,9999) (IDCO(I,J),J=1,5)
50 CONTINUE
WRITE (KA4,9999) IWL, IECOM, JMAXV, JTRB2
WRITE (KA4,9995)
9995 FORMAT (12H END OF DATA)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUT016
C OUTPUT OF BLOCKS OF SECOND ORDER IN THE S-DOMAIN
C-----------------------------------------------------------------------
C
SUBROUTINE OUT016
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
C
COMMON /SCRAT/ ADUM(32)
DIMENSION NZE(16)
EQUIVALENCE (ADUM(1),NZE(1))
C
WRITE (KA2,9999) SFA
9999 FORMAT (//39H BLOCKS OF SECOND ORDER IN THE S-DOMAIN//9H CONSTANT,
* 15H GAIN FACTOR = , E15.7//3H L, 5X, 5HA2(L), 8X, 5HA1(L),
* 8X, 5HA0(L), 10X, 5HD1(L), 8X, 5HD0(L)/)
C
N = 0
ME = NZM(4)
DO 10 I=1,ME
NZE(I) = NZERO(I)
10 CONTINUE
NB = (NDEG+1)/2
II = 1
NZ = NZE(II)
C
DO 180 I=1,NB
20 IF (NZ.GT.0) GO TO 30
II = II + 1
NZ = NZE(II)
GO TO 20
30 QI = SM(II,4)
IF (NZ.EQ.1) GO TO 70
IF (QI.GE.FLMA) GO TO 50
40 B2 = 1.
B1 = 0.
B0 = QI*QI
GO TO 60
50 B2 = 0.
B1 = 0.
B0 = 1.
60 NZ = NZ - 2
GO TO 130
70 IF (II.EQ.ME) GO TO 120
MA = II + 1
DO 80 J=MA,ME
IF (SM(J,4).GE.FLMA) GO TO 90
IF (SM(J,4).LE.FLMI) GO TO 100
80 CONTINUE
GO TO 120
90 NZE(J) = NZE(J) - 1
IF (QI.GE.FLMA) GO TO 50
GO TO 110
100 NZE(J) = NZE(J) - 1
IF (QI.LE.FLMI) GO TO 40
110 B2 = 0.
B1 = 1.
B0 = 0.
NZ = NZ - 1
GO TO 130
120 IF (QI.GE.FLMA) GO TO 50
GO TO 110
C
130 N = N + 1
QR = SPR(N)
QI = SPI(N)
IF (ABS(QI).LT.FLMI) GO TO 140
C1 = -2.*QR
C0 = QR*QR + QI*QI
GO TO 170
140 IF (N.GE.NJ) GO TO 150
IF (ABS(SPI(N+1)).LT.FLMI) GO TO 160
150 C1 = -QR
C0 = 0.
GO TO 170
160 N = N + 1
QI = SPR(N)
C1 = -QR - QI
C0 = QR*QI
C
170 WRITE (KA2,9998) I, B2, B1, B0, C1, C0
9998 FORMAT (1X, I3, 3F13.8, 2X, 2F13.8)
180 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUT022
C OUTPUT OF THE TERMINATION KIND
C-----------------------------------------------------------------------
C
SUBROUTINE OUT022
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
COMMON /COPTCP/ IS1, IS2, IABO, IWLGS, IWLGN, IWLGP, ADEPSS,
* ADEPSN, ADEPSP, ACXS, ACXN, ACXP
C
WRITE (KA2,9999) ITER
9999 FORMAT (///38H OPTIMIZATION IS TERMINATED AFTER THE , I4,
* 6H. STEP)
GO TO (10, 20, 30), IABO
10 WRITE (KA2,9998)
9998 FORMAT (25H THREE UNSUCCESSFUL STEPS)
GO TO 40
20 WRITE (KA2,9997)
9997 FORMAT (24H MAXIMUM NUMBER OF STEPS)
GO TO 40
30 WRITE (KA2,9996)
9996 FORMAT (54H OPTIMIZATION PARAMETER IS COMING TO THE END OF THE RE,
* 4HGION)
40 RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUT023
C OUTPUT OF THE ITERATION NUMBER
C-----------------------------------------------------------------------
C
SUBROUTINE OUT023
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /COPTCO/ LOPTW, LSTAB, ACXMI, ACXMA, ITER, ITERM, ITERM1
C
WRITE (KA2,9999) ITER
9999 FORMAT (1X////1X, I3, 13H-TH ITERATION)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUT025
C OUTPUT OF THE OPTIMIZATION OPTIONS
C-----------------------------------------------------------------------
C
SUBROUTINE OUT025
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /COPTST/ LOPTS, ISTOR
COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
* LPOT2, JSTRU
COMMON /CRENO/ LCNO, ICNO(16,5)
C
DIMENSION IES(4), NO(4)
DATA IES(1), IES(2), IES(3), IES(4) /1H ,1HY,1HE,1HS/
DATA NO(1), NO(2), NO(3), NO(4) /1H ,1H ,1HN,1HO/
C
IF (LOPTS.EQ.0) GO TO 10
WRITE (KA2,9999)
9999 FORMAT (1X//41H OPTIMIZATION OF THE PAIRING AND ORDERING//)
WRITE (KA2,9998) ISTOR
9998 FORMAT (46H STORAGE FOR THE INTERMEDIATE RESULTS ISTOR =, I4/)
IF (LOPTS.EQ.(-1)) WRITE (KA2,9997) NO
IF (LOPTS.NE.(-1)) WRITE (KA2,9997) IES
9997 FORMAT (/14H FIXED PAIRING, 32X, 4A1)
C
10 WRITE (KA2,9996) ISTRU
9996 FORMAT (1X/20H REALIZED STRUCTURE , 19X, 7HISTRU =, I4/)
WRITE (KA2,9995) ISCAL
9995 FORMAT (/16H SCALING OPTION, 23X, 7HISCAL =, I4)
IF (LPOT2.EQ.(-1)) WRITE (KA2,9994) IES
IF (LPOT2.NE.(-1)) WRITE (KA2,9994) NO
9994 FORMAT (10X, 24HBY A FACTOR OF POWER TWO, 12X, 4A1)
WRITE (KA2,9993) SCALM
9993 FORMAT (10X, 21HCHOSEN MAXIMUM OF THE/10X, 15HOVERFLOW POINTS,
* 14X, 7HSCALM =, F6.3/)
WRITE (KA2,9992)
9992 FORMAT (/33H REALIZATION OF THE COEFFICIENTS )
IF (LCNO.EQ.(-1)) WRITE (KA2,9991) NO, IWLR
IF (LCNO.NE.(-1)) WRITE (KA2,9991) IES, IWLR
9991 FORMAT (5X, 33HCONSIDERATION OF THE QUANTIZATION, 8X, 4A1/5X,
* 10HWORDLENGTH, 24X, 7HIWLR =, I4/24X, 19H(FOR IWLR=100 NO RO,
* 7HUNDING)/)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUT026
C OUTPUT OF THE NOISE VALUES OF THE SECOND-ORDER BLOCKS
C-----------------------------------------------------------------------
C
SUBROUTINE OUT026
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /CREST/ ISTRU, ISCAL, SCALM, ISEQ(16,2), LSEQ, IWLR,
* LPOT2, JSTRU
COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
COMMON /CNFUNC/ AQC(5), BN2(5), BN1(5), BN0(5)
COMMON /CRENO/ LCNO, ICNO(16,5)
COMMON /CPOW/ PNU, PNC, AND, ITCORP
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
C
DIMENSION NCHA(5), NQUAN(5)
C
IF (IB.LE.0) GO TO 50
IF (IB.EQ.1) WRITE (KA2,9999)
9999 FORMAT (1X//47H NUM DEN SCALING UNCOR. NOISE COR. NOISE,
* 14H DC-OFFSET, 12H CHA CO-QUAN)
ALOG2 = ALOG(2.)
IF (IB.GT.NB) GO TO 10
II = IB
IN = ISEQ(IB,1)
ID = ISEQ(IB,2)
GO TO 20
10 IN = 0
ID = 0
II = NB
20 DO 40 J=1,5
Q = AQC(J)
IF (Q.EQ.1.) GO TO 30
Q = ALOG(Q)/ALOG2
NQUAN(J) = -IFIX(Q-0.5)
NCHA(J) = ICNO(II,J)
GO TO 40
30 NQUAN(J) = 0
NCHA(J) = 0
40 CONTINUE
WRITE (KA2,9998) IN, ID, SCA, PNU, PNC, AND, (NCHA(J),NQUAN(J),J=
* 1,5)
9998 FORMAT (/1X, 2(I3, 2X), F9.6, 1X, 3E13.5, 2I4, 4H BIT/(60X, 2I4,
* 4H BIT))
RETURN
C
50 WRITE (KA2,9997) AND
9997 FORMAT (46X, 14(1H-)/46X, E14.5/)
Q = AND*AND
PNT = PNU + PNC + Q
WRITE (KA2,9996) PNU, PNC, Q, PNT
9996 FORMAT (12H TOTAL NOISE, 9X, 3E13.5, 1H=, E11.4)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUT027
C OUTPUT OF THE DIFFERENT NOISE VALUES
C-----------------------------------------------------------------------
C
SUBROUTINE OUT027
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /CNOISE/ RI, RIN, RE, REN, FAC
C
ANP = RI/12.
ALANP = 10.*ALOG10(ANP)
WRITE (KA2,9999) ANP, ALANP
9999 FORMAT (1H //22H ABSOLUTE OUTPUT NOISE, 7X, 3HANP, 3X, 1H=,
* E13.5, 12H (Q*Q) =, F6.2, 3H DB/)
WRITE (KA2,9998) FAC
9998 FORMAT (/22H SCALING AT THE OUTPUT, 7X, 3HSCA, 3X, 1H=, E13.5/)
RNP = RIN/12.
ALRNP = 10.*ALOG10(RNP)
WRITE (KA2,9997) RNP, ALRNP
9997 FORMAT (/22H RELATIVE OUTPUT NOISE, 7X, 3HRNP, 3X, 1H=, E13.5,
* 12H (Q*Q) ,F6.2, 3H DB)
WRITE (KA2,9995)
ALRIN = 10.*ALOG10(RIN)
WRITE (KA2,9996) RIN, ALRIN
9996 FORMAT (/19H INNER NOISE FIGURE, 10X, 3HRIN, 3X, 1H=, E13.5,
* 12H (Q*Q)/12 = , F6.2, 3H DB)
WRITE (KA2,9995)
9995 FORMAT (/23H RELATED TO MAX /H/ = 1/)
WRITE (KA2,9994) REN
9994 FORMAT (/22H ENTRANCE NOISE FIGURE, 7X, 3HREN, 3X, 1H=, E13.5/)
WRITE (KA2,9995)
ADW = RIN/(1.+REN)
ADW = ALOG(ADW)/(2.*ALOG(2.))
IDW = IFIX(ADW)
Q = FLOAT(IDW)
IF (Q.LT.ADW) IDW = IDW + 1
WRITE (KA2,9993) ADW, IDW
9993 FORMAT (/22H ADDITIONAL WORDLENGTH, 6X, 11HDELTA W = /, F7.3,
* 4H/ = , I4/)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUT030
C OUTPUT OF THE TOLERANCE SCHEME
C-----------------------------------------------------------------------
C
SUBROUTINE OUT030
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
C
DIMENSION FE(4)
C
AL(ADEL) = -20.*ALOG10(ADEL)
C
ME = 2
IF (ITYP.GT.2) ME = 4
WRITE (KA2,9999)
9999 FORMAT (//17H TOLERANCE SCHEME)
WRITE (KA2,9998)
9998 FORMAT (//13H FILTER-TYPE )
GO TO (10, 20, 30, 40), ITYP
10 WRITE (KA2,9997)
9997 FORMAT (25X, 7HLOWPASS)
GO TO 50
20 WRITE (KA2,9996)
9996 FORMAT (25X, 8HHIGHPASS)
GO TO 50
30 WRITE (KA2,9995)
9995 FORMAT (25X, 8HBANDPASS)
GO TO 50
40 WRITE (KA2,9994)
9994 FORMAT (25X, 8HBANDSTOP)
50 WRITE (KA2,9993)
9993 FORMAT (/15H APPROXIMATION )
GO TO (60, 70, 80, 90), IAPRO
60 WRITE (KA2,9992)
9992 FORMAT (25X, 11HBUTTERWORTH)
GO TO 100
70 WRITE (KA2,9991)
9991 FORMAT (25X, 11HCHEBYSHEV I)
GO TO 100
80 WRITE (KA2,9990)
9990 FORMAT (25X, 12HCHEBYSHEV II)
GO TO 100
90 WRITE (KA2,9989)
9989 FORMAT (25X, 8HELLIPTIC)
C
100 IF (SF.EQ.0.) GO TO 120
WRITE (KA2,9988) SF
9988 FORMAT (/18H SAMPLING FREQ. , F15.6)
Q = SF*0.5/PI
DO 110 I=1,ME
FE(I) = OM(I)*Q
110 CONTINUE
WRITE (KA2,9987) (FE(I),I=1,ME)
9987 FORMAT (20H CUTOFF FREQUENCIES , 4F13.6)
C
120 WRITE (KA2,9986) (OM(I),I=1,ME)
9986 FORMAT (/20H NORM. CUTOFF FREQ. , 4F13.6)
DO 130 I=1,ME
Q = OM(I)*0.5
FE(I) = SIN(Q)/COS(Q)
130 CONTINUE
WRITE (KA2,9985) (FE(I),I=1,ME)
9985 FORMAT (/20H CUTOFF FREQ. S-DOM., 4F13.6)
P = 1. - ADELP
Q = AL(P)
P = SQRT(1.-P*P)
WRITE (KA2,9984) ADELP, Q, P
9984 FORMAT (//20H PASSBAND RIPPLE(S) , F15.6, F12.4, 9H DB P =,
* F12.4)
Q = FLMA
IF (ADELS.GT.0.) Q = AL(ADELS)
WRITE (KA2,9983) ADELS, Q
9983 FORMAT (20H STOPBAND RIPPLE(S) , F15.6, F12.4, 3H DB)
IF (NDEG.EQ.0) RETURN
WRITE (KA2,9982) NDEG
9982 FORMAT (//32H FILTER DEGREE ASSIGNED TO ORDER, I3)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUT031
C OUTPUT OF THE NORMALIZED PARAMETER IN THE S-DOMAIN
C-----------------------------------------------------------------------
C
SUBROUTINE OUT031
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /TOLSN/ VSN, VD, A
C
WRITE (KA2,9999)
9999 FORMAT (//37H NORMALIZED PARAMETER IN THE S-DOMAIN)
WRITE (KA2,9998) VD, VSN
9998 FORMAT (/3H VD, 6X, 1H=, F13.6/4H VSN, 5X, 1H=, F13.6)
IF (ITYP.LE.2) RETURN
C
WRITE (KA2,9997) A
9997 FORMAT (2H A, 7X, 1H=, F13.6)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUT032
C OUTPUT OF THE REALIZED CUTOFF FREQUENCIES
C-----------------------------------------------------------------------
C
SUBROUTINE OUT032
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
C
ME = 2
IF (ITYP.GE.3) ME = 4
WRITE (KA2,9999) (ROM(I),I=1,ME)
9999 FORMAT (/9H REALIZED/20H NORM. CUTOFF FREQ. , 4F13.6)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUT033
C OUTPUT OF THE FILTER DEGREE
C-----------------------------------------------------------------------
C
SUBROUTINE OUT033
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /TOLSN/ VSN, VD, A
COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
C
Q = ADEG*EDEG
WRITE (KA2,9999) ADEG, Q, NDEG
9999 FORMAT (//20H MIN. FILTER DEGREE , F12.4//18H DEGREE EXTENSION ,
* F14.4//20H CHOSEN FILTER DEG. , I7)
IF (ITYP.LE.2) RETURN
IQ = 2*NDEG
WRITE (KA2,9998) IQ
9998 FORMAT (26H FOR THE REFERENCE LOWPASS//21H ACTUAL FILTER DEGREE,
* I6)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUT034
C OUTPUT OF THE CHARACTERISTIC FUNCTION /K(JV)/**2
C-----------------------------------------------------------------------
C
SUBROUTINE OUT034
C
DOUBLE PRECISION DSK(16), DK, DKS, DCAP02, DCAP04
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /OUTDAT/ IP, PRE, PIM, IZ, ZRE, ZIM
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /TOLCHA/ GD1, GD2, ACAP12, ADELTA, ADEG
COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
COMMON /RESIN2/ DK, DKS, DCAP02, DCAP04
EQUIVALENCE (PREN(1),DSK(1))
C
WRITE (KA2,9999) ADELTA
9999 FORMAT (/11H CAP. DELTA, 10X, F15.6)
IF (IAPRO.LT.4) RETURN
WRITE (KA2,9998)
9998 FORMAT (//50H ZEROS OF THE CHARACTERISTIC FUNCTION /K(J*V)/**2/)
IP = 0
IZ = 2
ZIM = 0.
DO 10 I=1,NJ
ZRE = SNGL(DSK(I))
CALL OUT002
10 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUT035
C OUTPUT OF THE TOLERANCE
C OUTPUT OF THE BOUND PAIR OF THE DESIGN PARAMETER AC
C-----------------------------------------------------------------------
C
SUBROUTINE OUT035
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
C
WRITE (KA2,9999) UGC, OGC
9999 FORMAT (//38H BOUND PAIR OF THE DESIGN PARAMETER C , F9.4,
* 16H .LE. C .LE. , F9.4)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUT036
C OUTPUT OF THE ZEROS AND POLES OF THE NORMALIZED REFERENCE LOWPASS
C-----------------------------------------------------------------------
C
SUBROUTINE OUT036
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /OUTDAT/ IP, PRE, PIM, IZ, ZRE, ZIM
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
C
NZ = NZM(4)
IF (IZ.NE.0) GO TO 10
WRITE (KA2,9999)
9999 FORMAT (////50H POLES AND ZEROS OF THE NORMALIZED REFERENCE LOWPA,
* 18HSS IN THE S-DOMAIN)
IF (IP.NE.0) NZ = 0
GO TO 20
10 WRITE (KA2,9998)
9998 FORMAT (////50H POLES AND ZEROS OF THE REFERENCE FILTER IN THE S-,
* 6HDOMAIN)
20 CONTINUE
PRE = SFA
CALL OUT001
ME = MAX0(NZ,NJ)
DO 60 I=1,ME
IF (IP.EQ.0) GO TO 30
IP = 2
Q = SPI(I)
IF (ABS(Q).LT.FLMI) IP = 1
PRE = SPR(I)
PIM = Q
IF (NZ.LT.I) GO TO 40
IF (NJ.LT.I) IP = 0
30 IZ = NZERO(I)
ZRE = 0.
ZIM = SM(I,4)
GO TO 50
40 IZ = 0
50 CALL OUT002
60 CONTINUE
C
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUT037
C OUTPUT OF THE CHOSEN DESIGN PARAMETER
C-----------------------------------------------------------------------
C
SUBROUTINE OUT037
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /TOL/ ITYP, IAPRO, NDEG, OM(4), SF, ADELP, ADELS
COMMON /DESIGN/ NDEGF, EDEG, ACX, NORMA, LSOUT, LVSN, LSYM
COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
C
WRITE (KA2,9999) ACX, AC
9999 FORMAT (/26H CHOSEN DESIGN PAR. CX =, F10.7, 10H (=) C =,
* F15.7)
QP = 100.*RDELP/ADELP
QS = 100.*RDELS/ADELS
WRITE (KA2,9998) RDELP, QP, RDELS, QS
9998 FORMAT (/40H UTILIZATION OF THE PASSBAND DELTA P =, F10.7,
* 3H =, F10.5, 8H PERCENT/20X, 20HSTOPBAND DELTA S =, F10.7,
* 3H =, F10.5, 8H PERCENT)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUT038
C OUTPUT OF THE POLES AND ZEROS IN THE Z-DOMAIN
C-----------------------------------------------------------------------
C
SUBROUTINE OUT038
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /OUTDAT/ IP, PRE, PIM, IZ, ZRE, ZIM
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
COMMON /RESZ/ ZFA, ZM(18,4), ZZR(16), ZZI(16), ZPR(16), ZPI(16)
COMMON /RESIN1/ PREN(16), PIMN(16), UGC, OGC, ACK, NJ, NH
C
WRITE (KA2,9999)
9999 FORMAT (////32H POLES AND ZEROS IN THE Z-DOMAIN)
NZ = NZM(4)
PRE = ZFA
CALL OUT001
ME = MAX0(NJ,NZ)
C
DO 30 I=1,ME
IP = 0
IF (NJ.LT.I) GO TO 10
IP = 2
Q = ZPI(I)
IF (ABS(Q).LT.FLMI) IP = 1
PRE = ZPR(I)
PIM = Q
10 IZ = 0
IF (NZ.LT.I) GO TO 20
IZ = NZERO(I)
ZRE = ZZR(I)
ZIM = ZZI(I)
20 CALL OUT002
30 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUT039
C OUTPUT OF THE EXTREMA OF THE MAGNITUDE OF THE TRANSFER FUNCTION
C-----------------------------------------------------------------------
C
SUBROUTINE OUT039
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /RES/ AC, ROM(4), RDELP, RDELS, NZM(4)
COMMON /RESS/ SFA, SM(18,4), NZERO(16), SPR(16), SPI(16)
COMMON /RESZ/ ZFA, ZM(18,4), ZZR(16), ZZI(16), ZPR(16), ZPI(16)
C
WRITE (KA2,9999)
9999 FORMAT (///51H EXTREMES OF THE MAGNITUDE OF THE TRANSFER FUNCTION,
* 1H , 20H(COEFS NOT ROUNDED) )
FN = 180./PI
WRITE (KA2,9998)
9998 FORMAT (/1X, 26HMAXIMA IN THE PASSBAND (P)/1X, 15HMAXIMA IN THE S,
* 11HTOPBAND (S)/1X, 14HBAND S-DOMAIN, 7X, 8HZ-DOMAIN, 10X,
* 9HMAGNITUDE/18X, 6HIN RAD, 3X, 9HIN DEGREE/)
C
II = 1
JJ = 3
10 I = NZM(II)
J = NZM(JJ)
K = MAX0(I,J)
DO 40 L=1,K
IF (I.LT.L) GO TO 20
OMEG1 = FN*ZM(L,II)
AMAG1 = AMAGO(ZM(L,II))
IF (J.LT.L) GO TO 30
OMEG2 = FN*ZM(L,JJ)
AMAG2 = AMAGO(ZM(L,JJ))
WRITE (KA2,9997) SM(L,II), ZM(L,II), OMEG1, AMAG1, SM(L,JJ),
* ZM(L,JJ), OMEG2, AMAG2
9997 FORMAT (2H P, 3(F10.4, 2X), F12.6/2H S, 3(F10.4, 2X), F12.6)
9996 FORMAT (2H P, 3(F10.4, 2X), F12.6)
GO TO 40
20 OMEG2 = FN*ZM(L,JJ)
AMAG2 = AMAGO(ZM(L,JJ))
WRITE (KA2,9995) SM(L,JJ), ZM(L,JJ), OMEG2, AMAG2
9995 FORMAT (2H S, 3(F10.4, 2X), F12.6)
GO TO 40
30 WRITE (KA2,9996) SM(L,II), ZM(L,II), OMEG1, AMAG1
40 CONTINUE
IF (II.EQ.2) RETURN
WRITE (KA2,9994)
9994 FORMAT (/1X, 26HMINIMA IN THE PASSBAND (P)/1X, 15HMINIMA IN THE S,
* 11HTOPBAND (S)/1X, 14HBAND S-DOMAIN, 7X, 8HZ-DOMAIN, 10X,
* 9HMAGNITUDE/18X, 6HIN RAD, 3X, 9HIN DEGREE/)
II = 2
JJ = 4
GO TO 10
C
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUT042
C OUTPUT OF THE START GRID FOR THE SEARCH OF THE GLOBAL EXTREMA
C OF THE TRANSFER FUNCTION
C-----------------------------------------------------------------------
C
SUBROUTINE OUT042
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /CGRID/ GR(64), NGR(12)
C
WRITE (KA2,9999)
9999 FORMAT (///36H GRID FOR THE SEARCH OF THE EXTREMES)
J = 1
IF (NGR(1).NE.0) GO TO 10
MA = NGR(8)
ME = NGR(9)
J = 4
10 GO TO (20, 30, 40, 50, 80), J
20 WRITE (KA2,9998)
9998 FORMAT (/10H PASS BAND/)
MA = 1
ME = NGR(3)
GO TO 60
30 MA = ME + 1
ME = NGR(5)
IF (ME.LT.MA) GO TO 70
GO TO 60
40 WRITE (KA2,9997)
9997 FORMAT (/10H STOP BAND/)
MA = MAX0(NGR(3),NGR(5)) + 1
ME = NGR(6)
GO TO 60
50 MA = ME + 1
ME = NGR(7)
IF (ME.LT.MA) GO TO 80
60 WRITE (KA2,9996) (GR(I),I=MA,ME)
9996 FORMAT (1X, 7F10.4)
70 J = J + 1
GO TO 10
80 RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUT043
C HEADLINE FOR THE SEARCH OF THE MINIMUM WORDLENGTH
C-----------------------------------------------------------------------
C
SUBROUTINE OUT043
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
C
WRITE (KA2,9999)
9999 FORMAT (//29H SEARCH OF MINIMUM WORDLENGTH//4H IWL, 6X, 5HEPS-P,
* 7X, 5HEPS-S, 7X, 4HPMAX)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUT044
C OUTPUT OF ERROR EPSILON AND MAXIMUM MAGNITUDE
C-----------------------------------------------------------------------
C
SUBROUTINE OUT044
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
* ISTAB, IDEPSL, IDEPSD
COMMON /CEPSIL/ EPS(2), PMAX
C
9999 FORMAT (1X, I4, 3F12.6)
C
WRITE (KA2,9999) IWL, EPS, PMAX
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUT045
C OUTPUT OF MINIMUM WORDLENGTH
C-----------------------------------------------------------------------
C
SUBROUTINE OUT045
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
* ISTAB, IDEPSL, IDEPSD
C
WRITE (KA2,9999) IWLG, ADEPSG, IWLD, ADEPSD
9999 FORMAT (/19H MINIMUM WORDLENGTH, I4, 13H DELTA EPS =,
* F12.6/19H TEST WORDLENGTH, I4, 13H DELTA EPS =, F12.6)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUT046
C OUTPUT OF THE LAYOUT OF THE ROUNDED COEFFICIENTS
C-----------------------------------------------------------------------
C
SUBROUTINE OUT046
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /CCOEFW/ IWL, IWLG, IWLD, IWLL, ADEPSG, ADEPSD, ADEPSL,
* ISTAB, IDEPSL, IDEPSD
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
COMMON /FILTRE/ IRCO(5), IECO(16,5), IDCO(16,5), IECOM
COMMON /SCRAT/ ADUM(32)
C
DIMENSION ACO(16,5)
DIMENSION ITEXT(5), JTEXT(5)
DIMENSION IOCT(12)
EQUIVALENCE (B2(1),ACO(1,1))
EQUIVALENCE (ADUM(1),IOCT(1))
C
DATA ITEXT(1), ITEXT(2), ITEXT(3), ITEXT(4), ITEXT(5) /1HB,1HB,
* 1HB,1HC,1HC/
DATA JTEXT(1), JTEXT(2), JTEXT(3), JTEXT(4), JTEXT(5) /1H2,1H1,
* 1H0,1H1,1H0/
C
WRITE (KA2,9999) IWL
9999 FORMAT (///35H LAYOUT OF THE ROUNDED COEFFICIENTS//11H WORDLENGTH,
* 2H , I4///11X, 4HCOEF, 10X, 5HCOEFS, 10X, 2HIR, 3X, 2HIE,
* 10X, 2HID, 5X, 5HOCTAL/)
DO 20 J=1,5
WRITE (KA2,9998)
9998 FORMAT (/)
DO 10 I=1,NB
AC = ACO(I,J)
ACS = ACOEFS(AC,IECO(I,J),IDCO(I,J),IRCO(J))
ID = SIGN(2.,AC)
Q = ACS
CALL OCTAL(Q, IOCT, 12)
WRITE (KA2,9997) ITEXT(J), JTEXT(J), I, AC, ACS, IRCO(J),
* IECO(I,J), ID, IDCO(I,J), (IOCT(II),II=1,12)
9997 FORMAT (1X, 2A1, 1H(, I2, 2H) , F13.8, 1H=, F13.8, 5H*2**(, I3,
* 2H -, I3, 1H), 2X, I3, 2H**, I4, 1X, 12I1)
10 CONTINUE
20 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUT060
C NEW PAGE
C-----------------------------------------------------------------------
C
SUBROUTINE OUT060
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
WRITE (KA2,9999)
9999 FORMAT (1H1)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUT061
C HEAD OF PRINTER PLOT
C-----------------------------------------------------------------------
C
SUBROUTINE OUT061
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /CPLOT/ IPLOT, LNOR, IPAG, OMLO, OMUP, RMAX, RMIN
COMMON /CLINE/ Y, X, IZ(111)
DIMENSION A(6)
EQUIVALENCE (A(1),IZ(1))
GO TO (10, 20, 30, 40), IPLOT
10 WRITE (KA2,9999)
9999 FORMAT (1H /30X, 35HMAGNITUDE OF THE FREQUENCY RESPONSE/)
GO TO 50
20 WRITE (KA2,9998)
9998 FORMAT (1H /28X, 37HATTENUATION OF THE FREQUENCY RESPONSE/)
GO TO 50
30 WRITE (KA2,9997)
9997 FORMAT (1H /31X, 14HPHASE RESPONSE/)
GO TO 50
40 WRITE (KA2,9996)
9996 FORMAT (1H /31X, 16HIMPULSE RESPONSE/)
50 WRITE (KA2,9995) A(1), A(3), A(5)
9995 FORMAT (10H MAGNITUDE, 2X, 1HX, F8.3, 2(11X, F9.3)/)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUT062
C OUTPUT OF THE PLOTTER LINE
C-----------------------------------------------------------------------
C
SUBROUTINE OUT062
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /CPLOT/ IPLOT, LNOR, IPAG, OMLO, OMUP, RMAX, RMIN
COMMON /CLINE/ Y, X, IZ(111)
DIMENSION IZDUM(56)
DATA ICHP, ICHB /1H+,1H /
DO 10 I=1,110,2
J = (I+1)/2
IZDUM(J) = ICHB
IF (IZ(I).NE.ICHB) IZDUM(J) = IZ(I)
IF (IZ(I+1).NE.ICHB) IZDUM(J) = IZ(I+1)
IF ((IZ(I).EQ.ICHP) .OR. (IZ(I+1).EQ.ICHP)) IZDUM(J) = ICHP
10 CONTINUE
IZDUM(56) = IZ(111)
IF (IPLOT.EQ.4) GO TO 20
WRITE (KA2,9999) Y, X, IZDUM
9999 FORMAT (1X, E8.2, F6.3, 1X, 56A1)
GO TO 30
20 IX = INT(X)
WRITE (KA2,9998) Y, IX, IZDUM
9998 FORMAT (1X, E8.2, I6, 1X, 56A1)
30 RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: HEAD
C OUTPUT OF A HEAD LINE
C-----------------------------------------------------------------------
C
SUBROUTINE HEAD(KAN)
C
WRITE (KAN,9999)
9999 FORMAT (41H === DOREDI === VERSION V005 === )
C
C HERE ON LINE FOR THE OUTPUT OF DATE AND TIME;
C REMOVE THE FOLLOWING STATEMENT
WRITE (KAN,9998)
9998 FORMAT (1H )
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUT001
C HEADLINE OF THE POLES AND ZERO OUTPUT
C-----------------------------------------------------------------------
C
SUBROUTINE OUT001
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /OUTDAT/ IP, PRE, PIM, IZ, ZRE, ZIM
C
IF (PRE.NE.0.) WRITE (KA2,9999) PRE
9999 FORMAT (/24H CONSTANT GAIN FACTOR = , E15.7)
WRITE (KA2,9998)
9998 FORMAT (//3X, 4HNUM., 11X, 5HPOLES, 15X, 4HNUM., 11X, 5HZEROS/)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUT002
C OUTPUT OF POLE AND ZERO
C-----------------------------------------------------------------------
C
SUBROUTINE OUT002
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
COMMON /OUTDAT/ IP, PRE, PIM, IZ, ZRE, ZIM
C
IF (IP.EQ.0) GO TO 20
IF (IZ.EQ.0) GO TO 10
WRITE (KA2,9999) IP, PRE, PIM, IZ, ZRE, ZIM
9999 FORMAT (2(1X, I5, 1X, F11.6, 5H +-J*, F11.6, 1X))
RETURN
10 WRITE (KA2,9999) IP, PRE, PIM
RETURN
20 IF (IZ.EQ.0) RETURN
WRITE (KA2,9998) IZ, ZRE, ZIM
9998 FORMAT (36X, I5, 1X, F11.6, 5H +-J*, F11.6, 1X)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OCTAL
C OCTAL CONVERSION
C-----------------------------------------------------------------------
C
SUBROUTINE OCTAL(AC, IOCT, N)
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
C
DIMENSION IOCT(1)
C
BC = AC
IF (BC.LT.(-1.0)) BC = -1.0
IF (BC.LT.0.) BC = -BC
IF (BC.GE.1.) BC = 0.
C
IOCT(1) = 0
DO 10 J=2,N
Q = 8.*BC
IQ = INT(Q)
IOCT(J) = IQ
BC = Q - FLOAT(IQ)
10 CONTINUE
C
IF (AC.GE.FLMI) RETURN
IA = 1
IOCT(1) = 1
IQ = N - 1
DO 20 J=1,IQ
JJ = N - J + 1
IF (IA.EQ.1 .AND. IOCT(JJ).EQ.0) GO TO 20
IOCT(JJ) = 7 - IOCT(JJ) + IA
IA = 0
20 CONTINUE
C
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: OUTSPE
C DUMMY SUBROUTINE FOR 'SPECIAL OUTPUT'
C-----------------------------------------------------------------------
C
SUBROUTINE OUTSPE
C
COMMON /CONTR/ IPRUN, IPCON, NINP, NOUT, NDOUT, LSPOUT, ISPOUT
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
C
C DATA TRANSFER BY NAMED COMMON BLOCKS:
C
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
C
C ISPMAX: MAXIMUM NUMBER OF IMPLEMENTED OPTIONS
ISPMAX=0
C
IF(ISPOUT.LE.0) RETURN
IF(ISPOUT.GT.ISPMAX) GO TO 900
C
C EXAMPLE FOR ISPMAX=3
C
C GO TO (100,200,300),ISPOUT
C 100 CONTINUE
C .
C OPTION NO. 1
C .
C RETURN
C
C 200 CONTINUE
C .
C OPTION NO. 2
C .
C RETURN
C
C 300 CONTINUE
C .
C OPTION NO. 3
C .
C RETURN
C
900 CALL ERROR (8)
RETURN
END
C
C ======================================================================
C
C DOREDI - SUBROUTINES: PART VI
C
C ======================================================================
C
C
C-----------------------------------------------------------------------
C SUBROUTINE: SPOW
C CALCULATION OF THE L-2 NORM
C-----------------------------------------------------------------------
C
SUBROUTINE SPOW(BMAX)
C
COMMON /CNFUNC/ AQC(5), BN2(5), BN1(5), BN0(5)
COMMON /CPOW/ PNU, PNC, AND, ITCORP
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
COMMON /CRENO/ LCNO, ICNO(16,5)
COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
C
BB2 = BN2(1)
BB1 = BN1(1)
BB0 = BN0(1)
BN2(1) = B2(1)
BN1(1) = B1(1)
BN0(1) = B0(1)
LQ = LCNO
LCNO = 2
IQ = IB
IB = 1
Q = PNU
CALL POWER
BMAX = PNU*FACT*FACT
PNU = Q
IB = IQ
LCNO = LQ
BN2(1) = BB2
BN1(1) = BB1
BN0(1) = BB0
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: POWER
C CALCULATION OF THE NOISE POWER GENERATED IN THE IB-TH BLOCK
C-----------------------------------------------------------------------
C
SUBROUTINE POWER
C
COMPLEX ZA, ZB, ZP, ZPS, ZQ, ZQI, Z0, Z1, ZBL, ZBL1, ZZ, ZZ1
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /CPOL/ ZP(16,2), ZPS(16,2)
COMMON /CRENO/ LCNO, ICNO(16,5)
COMMON /CNFUNC/ AQC(5), BN2(5), BN1(5), BN0(5)
COMMON /CPOW/ PNU, PNC, AND, ITCORP
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
C
DIMENSION EU(5)
C
DATA Z0,Z1 /(0.,0.),(1.,0.)/
C
ZBL(ZZ,A1,A0) = (ZZ+A1)*ZZ + A0
ZBL1(ZZ1,A21,A11,A01) = (A21*ZZ1+A11)*ZZ1 + A01
C
ZA = Z0
PNU = 0.
IF (LCNO.NE.2) GO TO 20
EU(1) = 1.
DO 10 I=2,5
EU(I) = 0.
10 CONTINUE
GO TO 160
20 DO 130 N=1,5
AQ = AQC(N)
IF (AQ.EQ.1.) GO TO 120
IF (LCNO.EQ.(-1)) AQ = 0.
KK = IB
IF (KK.GT.NB) KK = NB
K = ICNO(KK,N)
QS = 1.
QC1 = 0.
GO TO (30, 40, 50, 60, 70, 80, 90), K
30 QC2 = -1.
Q = AQ/2.
GO TO 100
40 QC2 = 2.*(1.-3./PI)
GO TO 110
50 QC2 = 2.
GO TO 110
60 QC2 = -1.
Q = (AQ-1.)/2.
GO TO 100
70 QS = 4. - 6./PI
GO TO 30
80 QC1 = -6.*(1.-2./PI)
GO TO 40
90 QS = 4.
GO TO 30
100 ZA = ZA + Q*ZBL1(Z1,BN2(N),BN1(N),BN0(N))
110 EU(N) = QS + (QC2*AQ+QC1)*AQ
GO TO 130
120 EU(N) = 0.
130 CONTINUE
IF (IB.GT.NB) GO TO 340
IF (CABS(ZA).EQ.0.) GO TO 150
ZA = ZA/ZBL(Z1,C1(IB),C0(IB))
IF (IB.EQ.NB) GO TO 150
MA = IB + 1
DO 140 I=MA,NB
ZA = ZA*ZBL1(Z1,B2(I),B1(I),B0(I))/ZBL(Z1,C1(I),C0(I))
140 CONTINUE
150 AND = REAL(ZA)
C
160 MARK = 0
R1 = 0.
DO 170 N=1,5
IF (EU(N).EQ.0.) GO TO 170
R1 = R1 + BN2(N)*BN0(N)*EU(N)
170 CONTINUE
R2 = 1.
DO 190 J=IB,NB
IF (J.NE.IB) R2 = R2*B2(J)*B0(J)
Q = C0(J)
IF (ABS(Q).GT.FLMI) GO TO 180
IF (MARK.NE.0) GO TO 360
MARK = J
Q = C1(J)
180 R2 = R2/Q
190 CONTINUE
IF (MARK.NE.0) GO TO 200
PNU = PNU + R1*R2
GO TO 280
C
200 DO 270 I=IB,NB
IF (I.EQ.MARK) GO TO 210
Q1 = C0(I)
Q2 = C1(I)
GO TO 220
210 Q1 = C1(I)
Q2 = 1.
220 IF (I.NE.IB) GO TO 240
R = 0.
DO 230 N=1,5
IF (EU(N).EQ.0.) GO TO 230
QQ1 = BN2(N)*(BN1(N)*Q1-BN0(N)*Q2)/Q1
QQ2 = BN0(N)*(BN1(N)-BN2(N)*C1(IB))
R = R + (QQ1+QQ2)*EU(N)
230 CONTINUE
R = R*R2
GO TO 260
240 R = R1*(B2(I)*(B1(I)*Q1-B0(I)*Q2)/Q1+B0(I)*(B1(I)-C1(I)*B2(I)))
DO 250 J=IB,NB
IF (J.NE.I .AND. J.NE.IB) R = R*B2(J)*B0(J)
Q = C0(J)
IF (J.EQ.MARK) Q = C1(J)
R = R/Q
250 CONTINUE
260 PNU = PNU + R
270 CONTINUE
C
280 DO 330 J=IB,NB
DO 320 K=1,2
ZQ = ZP(J,K)
IF (CABS(ZQ).LT.FLMI) GO TO 320
KK = 3 - K
ZQI = ZP(J,KK)
ZB = Z1/(ZQ-ZQI)
ZQI = Z1/ZQ
ZB = ZB*ZQI
ZA = Z0
DO 290 N=1,5
IF (EU(N).EQ.0.) GO TO 290
ZA = ZA + EU(N)*ZBL1(ZQ,BN2(N),BN1(N),BN0(N))*ZBL1(ZQI,
* BN2(N),BN1(N),BN0(N))
290 CONTINUE
ZB = ZB*ZA
DO 310 JJ=IB,NB
IF (JJ.EQ.IB) GO TO 300
ZB = ZB*ZBL1(ZQ,B2(JJ),B1(JJ),B0(JJ))*ZBL1(ZQI,B2(JJ),B1(JJ)
* ,B0(JJ))
C MAY 15, 1980 Q CHANGED TO ZQ IN FOLLOWING STATEMENT 300
300 IF (JJ.NE.J) ZB = ZB/ZBL(ZQ,C1(JJ),C0(JJ))
ZB = ZB/ZBL(ZQI,C1(JJ),C0(JJ))
310 CONTINUE
PNU = PNU + REAL(ZB)
320 CONTINUE
330 CONTINUE
RETURN
C
340 AND = REAL(ZA)
DO 350 N=1,5
PNU = PNU + EU(N)
350 CONTINUE
RETURN
C
360 CALL ERROR(33)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: SMAX
C COMPUTE THE MAXIMUM (IC=1) OR MINIMUM (IC=-1) OF THE TRANSFER FUNCTION
C-----------------------------------------------------------------------
C
SUBROUTINE SMAX(BMAX, MMA, MME, IC)
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /CGRID/ GR(64), NGR(12)
C
DIMENSION O(5), B(5)
C
FLM = 1./FLMA
FLER1 = FLER*0.05
MA = NGR(MMA)
ME = NGR(MME)
BMAX = 0.
MAX = 0
DO 30 I=MA,ME
QR = AMAGO(GR(I))
IF (IC.GT.0) GO TO 10
IF (QR.LT.FLM) GO TO 220
QR = 1./QR
10 IF (QR.LT.BMAX) GO TO 20
BMAX1 = BMAX
MAX1 = MAX
BMAX = QR
MAX = I
GO TO 30
20 IF (QR.LT.BMAX1) GO TO 30
BMAX1 = QR
MAX1 = I
30 CONTINUE
JJ = 1
IF ((ME-MA).GT.1) GO TO 40
O(1) = GR(MA)
O(3) = (GR(MA)+GR(ME))/2.
O(5) = GR(ME)
JJ = 3
GO TO 60
40 IF (IABS(MAX-MAX1).LT.2) JJ = 3
50 IF (MAX.EQ.MA) MAX = MA + 1
IF (MAX.EQ.ME) MAX = MAX - 1
O(1) = GR(MAX-1)
O(3) = GR(MAX)
O(5) = GR(MAX+1)
60 DO 80 I=1,5,2
Q = AMAGO(O(I))
IF (IC.GT.0) GO TO 70
IF (Q.LT.FLM) GO TO 220
Q = 1./Q
70 B(I) = Q
80 CONTINUE
GO TO 150
C
90 EMAX = 0.
DO 100 I=1,5
Q = B(I)
IF (Q.LT.EMAX) GO TO 100
MAX = I
EMAX = Q
100 CONTINUE
IF (EMAX.EQ.BMAX) GO TO 110
IF (EMAX/BMAX-1..LT.FLER) GO TO 180
110 IF (MAX.EQ.3) GO TO 130
IF (MAX.LT.3) GO TO 120
O(1) = O(3)
B(1) = B(3)
O(3) = O(4)
B(3) = B(4)
GO TO 140
120 O(5) = O(3)
B(5) = B(3)
O(3) = O(2)
B(3) = B(2)
GO TO 140
130 O(1) = O(2)
B(1) = B(2)
O(5) = O(4)
B(5) = B(4)
140 BMAX = EMAX
150 O(2) = (O(3)+O(1))/2.
Q = AMAGO(O(2))
IF (IC.GT.0) GO TO 160
IF (Q.LT.FLM) GO TO 220
Q = 1./Q
160 B(2) = Q
O(4) = (O(3)+O(5))/2.
Q = AMAGO(O(4))
IF (IC.GT.0) GO TO 170
IF (Q.LT.FLM) GO TO 220
Q = 1./Q
170 B(4) = Q
IF (O(5)-O(1).GE.FLER1) GO TO 90
180 GO TO (190, 200, 210), JJ
190 EMAX = BMAX
BMAX = BMAX1
BMAX1 = EMAX
MAX = MAX1
JJ = 2
GO TO 50
200 BMAX = AMAX1(BMAX,BMAX1)
210 IF (IC.LT.0) BMAX = 1./BMAX
RETURN
220 BMAX = 0.
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: SMIMP
C CALCULATION OF THE SUM OF THE IMPULSE RESPONSE MAGNITUDES
C ABSOLUTE SCALING CRITERION
C-----------------------------------------------------------------------
C
SUBROUTINE SMIMP(SUM)
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
C
DIMENSION X1(16), X2(16)
COMMON /SCRAT/ ADUM(32)
EQUIVALENCE (ADUM(1),X1(1)), (ADUM(17),X2(1))
C
DO 10 I=1,NB
X1(I) = 0.
X2(I) = 0.
10 CONTINUE
C
BMAX = 0.
SUM = 0.
IN = 0
JN = 0
Y = FACT
IM = 0
GO TO 30
20 Y = 0.
30 DO 40 I=1,NB
U = Y
Y = U*B2(I) + X1(I)
X1(I) = U*B1(I) - Y*C1(I) + X2(I)
X2(I) = U*B0(I) - Y*C0(I)
40 CONTINUE
JM = SIGN(1.,Y)
Y = ABS(Y)
SUM = SUM + Y
BMAX = AMAX1(BMAX,Y)
IF (IM.NE.0) GO TO 50
EMAX = BMAX
IM = JM
GO TO 20
50 IF (IM.EQ.JM) GO TO 70
IF (JN.GT.3) RETURN
IN = 0
Q = EMAX/BMAX
EMAX = Y
IM = JM
IF (Q.LT.FLER) GO TO 60
JN = 0
GO TO 20
60 JN = JN + 1
GO TO 20
70 IF (IN.GT.3) RETURN
EMAX = AMAX1(EMAX,Y)
IF (Y/EMAX.LT.FLER) IN = IN + 1
GO TO 20
END
C
C-----------------------------------------------------------------------
C FUNCTION: AMAGO
C COMPUTATION OF THE MAGNITUDE OF THE TRANSFER FUNCTION
C PARAMETER OMEGA IN RADIANS
C-----------------------------------------------------------------------
C
FUNCTION AMAGO(OMEG)
C
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
C
COMPLEX C, CH
C
C = CMPLX(COS(OMEG),SIN(OMEG))
CH = CMPLX(FACT,0.)
C
DO 10 I=1,NB
CH = CH*((B2(I)*C+B1(I))*C+B0(I))/((C+C1(I))*C+C0(I))
10 CONTINUE
AMAGO = CABS(CH)
RETURN
END
C
C-----------------------------------------------------------------------
C FUNCTION: PHASE
C PHASE OF OMEGA
C-----------------------------------------------------------------------
C
FUNCTION PHASE(OM)
C
COMPLEX CQ, CF
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
CF = CMPLX(FACT,0.)
CQ = CMPLX(COS(OM),SIN(OM))
DO 10 I=1,NB
CF = CF*((B2(I)*CQ+B1(I))*CQ+B0(I))/((CQ+C1(I))*CQ+C0(I))
10 CONTINUE
Y = AIMAG(CF)
X = REAL(CF)
PHASE = ATAN2(Y,X)
RETURN
END
C
C-----------------------------------------------------------------------
C FUNCTION: RESP
C RESPONSE OF AN INPUT SIGNAL
C-----------------------------------------------------------------------
C
FUNCTION RESP(X, IS)
C
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
COMMON /CSTATV/ X1(16), X2(16)
C
IF (IS.EQ.0) GO TO 20
C CLEAR STATE VARIABLES
DO 10 I=1,NB
X1(I) = 0.
X2(I) = 0.
10 CONTINUE
C
20 Y = X*FACT
DO 30 I=1,NB
U = Y
Y = U*B2(I) + X1(I)
X1(I) = U*B1(I) - Y*C1(I) + X2(I)
X2(I) = U*B0(I) - Y*C0(I)
30 CONTINUE
RESP = Y
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: CODE3
C STORE PAIRING AND ORDERING
C-----------------------------------------------------------------------
C
SUBROUTINE CODE3(L, J, N)
C
DIMENSION J(16), L(16)
DO 10 I=1,N
J(I) = L(I)
10 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: CODE4
C RESTORE PAIRING AND ORDERING
C-----------------------------------------------------------------------
C
SUBROUTINE CODE4(L, J, N)
C
DIMENSION J(16), L(16)
DO 10 I=1,N
L(I) = J(I)
10 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: CODE5
C SEARCH FOR A GIVEN CODE ARRAY IN LL
C RESTORE PAIRING AND ORDERING TO L
C-----------------------------------------------------------------------
C
SUBROUTINE CODE5(L, J, N, LL, M)
C
DIMENSION L(16), J(16), LL(16)
M = 2
DO 30 I=1,N
IQ = J(I)
DO 10 II=1,N
IF (IQ.EQ.LL(II)) GO TO 20
10 CONTINUE
RETURN
20 L(I) = IQ
30 CONTINUE
M = 1
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: QUAN
C CHECK THE ACTUAL COEFFICIENT QUANTIZATION
C-----------------------------------------------------------------------
C
SUBROUTINE QUAN(ACO, ACQ)
C
COMMON /COPTSP/ IB, JSEQN(16), JSEQD(16), AMAX, SCA, ALSBI
C
Q = ABS(ACO)
IF (Q.NE.AINT(Q)) GO TO 10
ACQ = 1.
RETURN
10 Q = Q*ALSBI
ACQ = 0.5/ALSBI
GO TO 30
20 Q = Q*0.5
ACQ = ACQ*2.
30 IF (Q.EQ.AINT(Q)) GO TO 20
RETURN
END
C
C-----------------------------------------------------------------------
C FUNCTION: ACOEF
C RETRANSFORM ONE COEFFICIENT
C-----------------------------------------------------------------------
C
FUNCTION ACOEF(AS, IE, ID, IR)
C
Q = AS*2.**(IR-IE)
ACOEF = Q
IF (ID.LE.(-100)) RETURN
AD = 2.**ID
ACOEF = Q - SIGN(AD,AS)
RETURN
END
C
C-----------------------------------------------------------------------
C FUNCTION: ACOEFS
C TRANSFORM ONE COEFFICIENT
C-----------------------------------------------------------------------
C
FUNCTION ACOEFS(AC, IE, ID, IR)
C
AD = 0.
IF (ID.GT.(-100)) AD = 2.**ID
ACOEFS = (AC-SIGN(AD,AC))*2.**(IE-IR)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: COPY01
C COPY COEFFICIENT FIELD /FILT/ INTO THE FIELD /SFILT/
C-----------------------------------------------------------------------
C
SUBROUTINE COPY01
C
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
COMMON /SFILT/ NBS, SFACT, SB2(16), SB1(16), SB0(16), SC1(16),
* SC0(16)
C
DIMENSION ACO(16,5), SACO(16,5)
EQUIVALENCE (ACO(1,1),B2(1)), (SACO(1,1),SB2(1))
C
NBS = NB
SFACT = FACT
DO 20 J=1,5
DO 10 I=1,NB
SACO(I,J) = ACO(I,J)
10 CONTINUE
20 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: COPY02
C COPY COEFFICIENT FIELD /SFILT/ INTO THE FIELD /FILT/
C-----------------------------------------------------------------------
C
SUBROUTINE COPY02
C
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
COMMON /SFILT/ NBS, SFACT, SB2(16), SB1(16), SB0(16), SC1(16),
* SC0(16)
C
DIMENSION ACO(16,5), SACO(16,5)
EQUIVALENCE (ACO(1,1),B2(1)), (SACO(1,1),SB2(1))
C
NB = NBS
FACT = SFACT
DO 20 J=1,5
DO 10 I=1,NB
ACO(I,J) = SACO(I,J)
10 CONTINUE
20 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: COPY03
C COPY POLE LOCATIONS
C-----------------------------------------------------------------------
C
SUBROUTINE COPY03
C
COMPLEX ZP, ZPS
C
COMMON /FILT/ NB, FACT, B2(16), B1(16), B0(16), C1(16), C0(16)
COMMON /CPOL/ ZP(16,2), ZPS(16,2)
C
DO 20 I=1,NB
DO 10 J=1,2
ZP(I,J) = ZPS(I,J)
10 CONTINUE
20 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C FUNCTION: DELLK
C CALCULATE COMPLETE ELLIPTIC INTEGRAL OF FIRST KIND
C-----------------------------------------------------------------------
C
DOUBLE PRECISION FUNCTION DELLK(DK)
C
DOUBLE PRECISION DPI, DOMI
DOUBLE PRECISION DE, DGEO, DK, DRI, DARI, DTEST
C
COMMON /CONST/ PI, FLMA, FLMI, FLER
COMMON /CONST2/ DPI, DOMI
C
DATA DE /1.D00/
C
DGEO = DE - DK*DK
IF (DGEO) 10, 10, 20
10 DELLK = DBLE(FLMA)
RETURN
C
20 DGEO = DSQRT(DGEO)
DRI = DE
30 DARI = DRI
DTEST = DARI*DOMI
DRI = DGEO + DRI
C
IF (DARI-DGEO-DTEST) 50, 50, 40
40 DGEO = DSQRT(DARI*DGEO)
DRI = 0.5D00*DRI
GO TO 30
50 DELLK = DPI/DRI
RETURN
END
C
C-----------------------------------------------------------------------
C FUNCTION: DSN2
C CALCULATION OF THE JACOBI'S ELLIPTIC FUNCTION SN(U,K)
C
C EXTERNAL CALCULATION OF THE PARAMETER NECESSARY
C DK = K($K)
C DQ = EXP(-PI*K'/K) ... (JACOBI'S NOME)
C-----------------------------------------------------------------------
C
DOUBLE PRECISION FUNCTION DSN2(DU, DK, DQ)
C
DOUBLE PRECISION DPI, DOMI
DOUBLE PRECISION DE, DZ, DPI2, DQ, DM, DU, DK, DC, DQQ, DH, DQ1,
* DQ2
C
COMMON /CONST2/ DPI, DOMI
C
DATA DE, DZ /1.D00,2.D00/
C
DPI2 = DPI/DZ
IF (DABS(DQ).GE.DE) GO TO 30
C
DM = DPI2*DU/DK
DC = DZ*DM
DC = DCOS(DC)
C
DM = DSIN(DM)*DK/DPI2
DQQ = DQ*DQ
DQ1 = DQ
DQ2 = DQQ
C
DO 10 I=1,100
DH = (DE-DQ1)/(DE-DQ2)
DH = DH*DH
DH = DH*(DE-DZ*DQ2*DC+DQ2*DQ2)
DH = DH/(DE-DZ*DQ1*DC+DQ1*DQ1)
DM = DM*DH
C
DH = DABS(DE-DH)
IF (DH.LT.DOMI) GO TO 20
C
DQ1 = DQ1*DQQ
DQ2 = DQ2*DQQ
10 CONTINUE
C
GO TO 30
C
20 DSN2 = DM
RETURN
C
30 DSN2 = 0.D00
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: DELI1
C ELLIPTIC FUNCTION
C-----------------------------------------------------------------------
C
SUBROUTINE DELI1(RES, X, CK)
C
DOUBLE PRECISION RES, X, CK, ANGLE, GEO, ARI, PIM, SQGEO, AARI,
* TEST, DPI
DOUBLE PRECISION DOMI
C
COMMON /CONST2/ DPI, DOMI
IF (X) 20, 10, 20
10 RES = 0.D0
RETURN
C
20 IF (CK) 40, 30, 40
30 RES = DLOG(DABS(X)+DSQRT(1.D0+X*X))
GO TO 130
C
40 ANGLE = DABS(1.D0/X)
GEO = DABS(CK)
ARI = 1.D0
PIM = 0.D0
50 SQGEO = ARI*GEO
AARI = ARI
ARI = GEO + ARI
ANGLE = -SQGEO/ANGLE + ANGLE
SQGEO = DSQRT(SQGEO)
IF (ANGLE) 70, 60, 70
C
C REPLACE 0 BY A SMALL VALUE, TEST
C
60 ANGLE = SQGEO*DOMI
70 TEST = AARI*DOMI*1.D+05
IF (DABS(AARI-GEO)-TEST) 100, 100, 80
80 GEO = SQGEO + SQGEO
PIM = PIM + PIM
IF (ANGLE) 90, 50, 50
90 PIM = PIM + DPI
GO TO 50
100 IF (ANGLE) 110, 120, 120
110 PIM = PIM + DPI
120 RES = (DATAN(ARI/ANGLE)+PIM)/ARI
130 IF (X) 140, 150, 150
140 RES = -RES
150 RETURN
END
C
C-----------------------------------------------------------------------
C FUNCTION: SINH
C-----------------------------------------------------------------------
C
FUNCTION SINH(X)
C
SINH = (EXP(X)-EXP(-X))/2.
RETURN
END
C
C-----------------------------------------------------------------------
C FUNCTION: COSH
C-----------------------------------------------------------------------
C
FUNCTION COSH(X)
C
COSH = (EXP(X)+EXP(-X))/2.
RETURN
END
C
C-----------------------------------------------------------------------
C FUNCTION: ARSINH
C-----------------------------------------------------------------------
C
FUNCTION ARSINH(X)
C
ARSINH = ALOG(X+SQRT(X*X+1.))
RETURN
END
C
C-----------------------------------------------------------------------
C FUNCTION: ARCOSH
C-----------------------------------------------------------------------
C
FUNCTION ARCOSH(X)
C
IF (X.LT.1.) GO TO 10
ARCOSH = ALOG(X+SQRT(X*X-1.))
RETURN
10 ARCOSH = 0.
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: DSQRTC
C COMPUTATION OF THE COMPLEX SQUARE ROOT IN DOUBLE PRECISION
C DU + J*DV = SQRT ( DX + J*DY )
C-----------------------------------------------------------------------
C
SUBROUTINE DSQRTC(DX, DY, DU, DV)
C
DOUBLE PRECISION DPI, DOMI, D1MACH
DOUBLE PRECISION DX, DU, DY, DV, DQ, DP
COMMON /CONST2/ DPI, DOMI
C
DQ = DX
DP = DY
C
DV = 0.5D00*DQ
DU = DQ*DQ + DP*DP
DU = DSQRT(DU)
DU = 0.5D00*DU
DV = DU - DV
DU = DV + DQ
IF (DABS(DU).LE.3.D0*D1MACH(3)) DU = 0.D0
DU = DSQRT(DU)
IF (DABS(DV).LE.3.D0*D1MACH(3)) DV = 0.D0
DV = DSQRT(DV)
IF (DP.LT.(-DOMI)) DU = -DU
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: ERROR
C OUTPUT OF AN ERROR MESSAGE
C-------------------------------------------------------------------
C
SUBROUTINE ERROR(IERCOD)
C
COMMON /CANPAR/ KA1, KA2, KA3, KA4, KA5, LINE
C
WRITE (KA2,9999) IERCOD
9999 FORMAT (22H *** ERROR *** NUMBER , I3)
IF (IERCOD.GE.20) STOP
RETURN
END
Return to Main Software Page