Make your own free website on Tripod.com
A Program for Multistage Decimation, Interpolation, and Narrow Band Filtering

Return to Main Software Page

C
C-----------------------------------------------------------------------
C MAIN PROGRAM: TEST PROGRAM FOR DECIMATION,INTERPOLATION AND FILTERING
C AUTHOR:       R E CROCHIERE, L R RABINER
C               BELL LABORATORIES, MURRAY HILL, NEW JERSEY, 07974
C INPUT:        COEF1 = ARRAY OF FILTER COEFFICIENTS FOR FIRST STAGE
C                       FOR TEST PROBLEM 1.
C               COEF2 = ARRAY OF FILTER COEFFICIENTS FOR SECOND STAGE
C                       FOR TEST PROBLEM 1.
C               COEF11 = ARRAY OF FILTER COEFFICIENTS FOR FIRST STAGE
C                        FOR 1-STAGE EXAMPLE IN TEST PROBLEM 3.
C               COEF1A = ARRAY OF FILTER COEFFICIENTS FOR FIRST STAGE
C                        FOR 3-STAGE EXAMPLE IN TEST PROBLEM 3.
C               COEF2A = ARRAY OF FILTER COEFFICIENTS FOR SECOND STAGE
C                        FOR 3-STAGE EXAMPLE IN TEST PROBLEM 3.
C               COEF3A = ARRAY OF FILTER COEFFICIENTS FOR THIRD STAGE
C                        FOR 3-STAGE EXAMPLE IN TEST PROBLEM 3.
C-----------------------------------------------------------------------
C
      COMMON /DIFC1/ K, ID1, ID2, ID3, N1, N2, N3, ITYPE, ID, ISB, N1P,
     *    N2P, N3P
      COMMON /DIFC2/ IDD, J1, J2, J3, J4, J5, J6, J7, J8, J9, J10, J11,
     *    J12, J1S, J2S, J3S, NCF1, NCF2, NCF3, NF1, NF2, NF3, IQ1,
     *    IQ2, IQ3, K1S, K2S, K3S
      DIMENSION X(60), COEF1(13), COEF2(14), SBUFF(106), BUFF(20),
     *    COEF3(1)
      DIMENSION COEF11(60)
      DIMENSION COEF1A(25), COEF2A(20), COEF3A(75)
      INTEGER TTI, TTO
C
C DEFINE I/O DEVICE CODES
C INPUT: INPUT TO THIS PROGRAM IS USER INTERACTIVE
C        THAT IS - A QUESTION IS WRITTEN ON THE USER
C        TERMINAL (TTO) AND THE USER TYPES IN THE ANSWER.
C
C OUTPUT: ALL OUTPUT IS WRITTEN ON THE STANDARD
C         OUTPUT UNIT (LPT).
C
      TTI = I1MACH(1)
      TTO = I1MACH(4)
      LPT = I1MACH(2)
C
C TEST EXAMPLE 1
C
      A = 0.9
      X(1) = -A
      X(2) = 1. - A*A
      DO 10 I=3,60
        X(I) = X(I-1)*A
  10  CONTINUE
      READ (TTI,9999) (COEF1(I),I=1,13)
      READ (TTI,9999) (COEF2(I),I=1,14)
      WRITE (TTO,9998)
      WRITE (TTO,9999) (COEF1(I),I=1,13)
      WRITE (TTO,9997)
      WRITE (TTO,9999) (COEF2(I),I=1,14)
9999  FORMAT (7F10.6)
9998  FORMAT (25H COEFFICIENTS FOR STAGE 1)
9997  FORMAT (25H COEFFICIENTS FOR STAGE 2)
      CALL DIINIT(2, 5, 2, 1, 25, 28, 0, COEF1, COEF2, COEF3, 1, BUFF,
     *    20, SBUFF, 106, IER)
      IF (IER.NE.0) STOP
      DO 30 I=1,3
        DO 20 J=1,20
          JX = (I-1)*20 + J
          BUFF(J) = X(JX)
  20    CONTINUE
        WRITE (TTO,9996)
9996    FORMAT (15H INPUT SEQUENCE)
        WRITE (TTO,9995) (BUFF(J),J=1,20)
9995    FORMAT (1H , 10F7.4)
        CALL DIFILT(COEF1, COEF2, COEF3, BUFF, SBUFF)
        WRITE (TTO,9994)
9994    FORMAT (26H DECIMATED OUTPUT SEQUENCE)
        WRITE (TTO,9995) (BUFF(J),J=1,2)
  30  CONTINUE
C
C TEST EXAMPLE 2
C RESTORE COEF1(1)
C
      COEF1(1) = COEF1(1)*2.
      CALL DIINIT(2, 5, 2, 1, 25, 28, 0, COEF1, COEF2, COEF3, 2, BUFF,
     *    20, SBUFF, 106, IER)
      IF (IER.NE.0) STOP
      WRITE (TTO,9993)
9993  FORMAT (1H /1H )
      DO 50 I=1,3
        DO 40 J=1,2
          JX = (I-1)*2 + J
          BUFF(J) = X(JX)
  40    CONTINUE
        WRITE (TTO,9996)
        WRITE (TTO,9995) (BUFF(J),J=1,2)
        CALL DIFILT(COEF1, COEF2, COEF3, BUFF, SBUFF)
        WRITE (TTO,9992)
9992    FORMAT (29H INTERPOLATED OUTPUT SEQUENCE)
        WRITE (TTO,9995) (BUFF(J),J=1,20)
  50  CONTINUE
C
C TEST EXAMPLE 3
C
      WRITE (TTO,9993)
      READ (TTI,9999) (COEF11(I),I=1,60)
      WRITE (TTO,9991)
9991  FORMAT (34H COEFFICIENTS FOR A 1-STAGE DESIGN)
      WRITE (TTO,9999) (COEF11(I),I=1,60)
      WRITE (TTO,9993)
      READ (TTI,9999) (COEF1A(I),I=1,25)
      READ (TTI,9999) (COEF2A(I),I=1,20)
      READ (TTI,9999) (COEF3A(I),I=1,75)
      WRITE (TTO,9990)
9990  FORMAT (34H COEFFICIENTS FOR A 3-STAGE DESIGN)
      WRITE (TTO,9989)
      WRITE (TTO,9999) (COEF1A(I),I=1,25)
      WRITE (TTO,9988)
      WRITE (TTO,9999) (COEF2A(I),I=1,20)
      WRITE (TTO,9987)
      WRITE (TTO,9999) (COEF3A(I),I=1,75)
9989  FORMAT (8H STAGE 1)
9988  FORMAT (8H STAGE 2)
9987  FORMAT (8H STAGE 3)
      STOP
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: DIINIT
C INITIALIZATION FOR DIFILT WHICH DECIMATES, INTERPOLATES OR
C FILTERS A SIGNAL
C-----------------------------------------------------------------------
C
      SUBROUTINE DIINIT(KD, ID1D, ID2D, ID3D, N1D, N2D, N3D, COEF1,
     *    COEF2, COEF3, ITYPED, BUFF, IDJD, SBUFF, ISBD, IERR)
      COMMON /DIFC1/ K, ID1, ID2, ID3, N1, N2, N3, ITYPE, ID, ISB, N1P,
     *    N2P, N3P
      COMMON /DIFC2/ IDD, J1, J2, J3, J4, J5, J6, J7, J8, J9, J10, J11,
     *    J12, J1S, J2S, J3S, NCF1, NCF2, NCF3, NF1, NF2, NF3, IQ1,
     *    IQ2, IQ3, K1S, K2S, K3S
      DIMENSION COEF1(1), COEF2(1), COEF3(1), SBUFF(1), BUFF(1)
C
C     KD = NO. OF STAGES OF DECIMATION AND /OR INTERPOLATION
C   ID1D = DECIMATION RATIO FOR STAGE 1
C   ID2D = DECIMATION RATIO FOR STAGE 2 = 1 IF KD=1
C   ID3D = DECIMATION RATIO FOR STAGE 3 = 1 IF KD=1 OR 2
C    N1D = FILTER LENGTH FOR STAGE 1
C    N2D = FILTER LENGTH FOR STAGE 2 = 0 IF KD = 1
C    N3D = FILTER LENGTH FOR STAGE 3 = 0 IF KD = 1 OR 2
C  COEF1 = COEF. ARRAY FOR STAGE 1     (SIZE = ?(N1D+1)/2! )?X!=INTEGER
C  COEF2 = COEF. ARRAY FOR STAGE 2     (SIZE = ?(N2D+1)/2! )    PART OF
C  COEF3 = COEF. ARRAY FOR STAGE 3     (SIZE = ?(N3D+1)/2! )
C ITYPED = 1 FOR DECIMATOR
C        = 2 FOR INTERPOLATOR
C        = 3 FOR FILTER
C   BUFF = STORAGE BUFFER FOR INPUT/OUTPUT
C   IDJD = SIZE OF BUFF = INTEGER*ID1D*ID2D*ID3D
C  SBUFF = SCRATCH STORAGE BUFFER FOR INTERNAL VARIABLES
C   ISBD = SIZE OF SBUFF
C        = 2*(N1P+N2P+N3P)                FOR DECIMATOR
C        = (N1P+N2P+N3P)+ID1D+ID2D        FOR INTERPOLATOR
C        = 3*(N1P+N2P+N3P)+ID1D+ID2D      FOR FILTER
C WHERE  N1P = (?N1D/ID1D!+1)*ID1D       ?X!=INTEGER PART OF X
C        N2P = (?N2D/ID2D!+1)*ID2D           OR X-1 IF X=INTEGER
C        N3P = (?N3D/ID3D!+1)*ID3D
C   IERR = ERROR CODE FOR DEBUGGING PURPOSES
C        = 0    NO ERRORS ENCOUNTERED IN INITIALIZATION
C        = 1    KD NOT EQUAL TO 1, 2, OR 3
C        = 2    IDJD NOT AN INTEGER MULTIPLE OF ID1D*ID2D*ID3D
C
C INITIALIZATION AND SETUP
C TRANSFER OF DUMMY ARGUMENTS TO COMMON
C
      K = KD
      ID1 = ID1D
      ID2 = ID2D
      ID3 = ID3D
      N1 = N1D
      N2 = N2D
      N3 = N3D
      ITYPE = ITYPED
      ID = IDJD
      ISB = ISBD
      IERR = 0
C
C PRELIMINARY CHECKS
C
      IF ((K.GT.3) .OR. (K.LT.1)) IERR = 1
      M = 1
      IF (K.EQ.3) GO TO 10
      N3 = 0
      ID3 = 0
      IQ3 = 0
      IF (K.EQ.2) GO TO 20
      N2 = 0
      ID2 = 0
      IQ2 = 0
      GO TO 30
  10  M = M*ID3
      IQ3 = N3/ID3
      IF (N3.NE.(IQ3*ID3)) IQ3 = IQ3 + 1
  20  M = M*ID2
      IQ2 = N2/ID2
      IF (N2.NE.(IQ2*ID2)) IQ2 = IQ2 + 1
  30  M = M*ID1
      IQ1 = N1/ID1
      IF (N1.NE.(IQ1*ID1)) IQ1 = IQ1 + 1
      IDD = ID/M
      IF (ID.NE.(M*IDD)) IERR = 2
      N1P = IQ1*ID1
      N2P = IQ2*ID2
      N3P = IQ3*ID3
C
C SETUP OF ADDRESS LOCATIONS IN SBUFF FOR INTERNAL STORAGE BUFFERS
C AND ZERO OUT SBUFF FOR INITIALIZATION
C
      J1 = 0
      J2 = J1 + 2*N1P
      J3 = J2 + 2*N2P
      J4 = J3 + 2*N3P
      IF (ITYPE.EQ.2) J4 = 0
      J5 = J4 + IQ1*ID1
      J6 = J5 + IQ2*ID2
      J7 = J6 + IQ3*ID3
      J8 = J7 + 2*IQ1
      J9 = J8 + 2*IQ2
      J10 = J9 + 2*IQ3
      J11 = J10 + ID2
      J12 = J11 + ID3
      IF (ITYPE.EQ.1) J12 = J4
      IF (ISB.LT.J12) IERR = 3
      DO 40 M=1,J12
        SBUFF(M) = 0.0
  40  CONTINUE
C
C SETUP OF SCRAMBLED COEFFICIENT SETS FOR INTERPOLATION,
C INITIALIZATION OF MOVING ADDRESS POINTERS AND
C HALVING OF FIRST COEFFICIENT IN SETS COEF1 TO COEF3
C
C STAGE 1
C
      IDX = J4 + 1
      NCF1 = (N1+1)/2
      SUM = ID1
      IF (ITYPE.EQ.1) GO TO 70
      DO 60 MD=1,ID1
        DO 50 MQ=1,IQ1
          M = (MD-1) + (MQ-1)*ID1
          IF (M.LT.NCF1) MM = NCF1 - M
          IF (M.GE.NCF1) MM = M - (N1-NCF1-1)
          IF (MM.LE.NCF1) SBUFF(IDX) = COEF1(MM)*SUM
          IF (MM.GT.NCF1) SBUFF(IDX) = 0.
          IDX = IDX + 1
  50    CONTINUE
  60  CONTINUE
  70  IF (N1.EQ.(2*NCF1-1)) COEF1(1) = COEF1(1)/2.0
      NF1 = (N1+2)/2
      J1S = J1 + N1P
      K1S = J7 + IQ1
      IF (K.LT.2) GO TO 140
C
C STAGE 2
C
      IDX = J5 + 1
      NCF2 = (N2+1)/2
      SUM = ID2
      IF (ITYPE.EQ.1) GO TO 100
      DO 90 MD=1,ID2
        DO 80 MQ=1,IQ2
          M = (MD-1) + (MQ-1)*ID2
          IF (M.LT.NCF2) MM = NCF2 - M
          IF (M.GE.NCF2) MM = M - (N2-NCF2-1)
          IF (MM.LE.NCF2) SBUFF(IDX) = COEF2(MM)*SUM
          IF (MM.GT.NCF2) SBUFF(IDX) = 0.
          IDX = IDX + 1
  80    CONTINUE
  90  CONTINUE
 100  IF (N2.EQ.(2*NCF2-1)) COEF2(1) = COEF2(1)/2.0
      NF2 = (N2+2)/2
      J2S = J2 + N2P
      K2S = J8 + IQ2
      IF (K.LT.3) GO TO 140
C
C STAGE 3
C
      IDX = J6 + 1
      NCF3 = (N3+1)/2
      SUM = ID3
      IF (ITYPE.EQ.1) GO TO 130
      DO 120 MD=1,ID3
        DO 110 MQ=1,IQ3
          M = (MD-1) + (MQ-1)*ID3
          IF (M.LT.NCF3) MM = NCF3 - M
          IF (M.GE.NCF3) MM = M - (N3-NCF3-1)
          IF (MM.LE.NCF3) SBUFF(IDX) = COEF3(MM)*SUM
          IF (MM.GT.NCF3) SBUFF(IDX) = 0.
          IDX = IDX + 1
 110    CONTINUE
 120  CONTINUE
 130  IF (N3.EQ.(2*NCF3-1)) COEF3(1) = COEF3(1)/2.0
      NF3 = (N3+2)/2
      J3S = J3 + N3P
      K3S = J9 + IQ3
 140  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: DIFILT
C DIFILT DECIMATES, INTERPOLATES OR FILTERS A SIGNAL,
C DIINIT MUST BE CALLED PRIOR TO CALLING THIS ROUTINE
C-----------------------------------------------------------------------
C
C FOR DECIMATOR, FILL BUFF WITH IDJD SAMPLES, CALL DIFILT, AND
C     RECEIVE (IDJD/ID1D*ID2D*ID3D) SAMPLES
C FOR INTERPOLATOR, SUPPLY (IDJD/ID1D*ID2D*ID3D) INPUT SAMPLES
C     TO BUFF, CALL DIFILT, AND RECEIVE IDJD OUTPUT SAMPLES
C FOR FILTER, SUPPLY AND RECEIVE IDJD SAMPLES FROM BUFF
C
      SUBROUTINE DIFILT(COEF1, COEF2, COEF3, BUFF, SBUFF)
      COMMON /DIFC1/ K, ID1, ID2, ID3, N1, N2, N3, ITYPE, ID, ISB, N1P,
     *    N2P, N3P
      COMMON /DIFC2/ IDD, J1, J2, J3, J4, J5, J6, J7, J8, J9, J10, J11,
     *    J12, J1S, J2S, J3S, NCF1, NCF2, NCF3, NF1, NF2, NF3, IQ1,
     *    IQ2, IQ3, K1S, K2S, K3S
      DIMENSION COEF1(1), COEF2(1), COEF3(1), SBUFF(1), BUFF(1)
C
C  COEF1 = COEF. ARRAY FOR STAGE 1     (SIZE = ?(N1D+1)/2! )?X!=INTEGER
C  COEF2 = COEF. ARRAY FOR STAGE 2     (SIZE = ?(N2D+1)/2! )    PART OF
C  COEF3 = COEF. ARRAY FOR STAGE 3     (SIZE = ?(N3D+1)/2! )
C   BUFF = STORAGE BUFFER FOR INPUT/OUTPUT
C  SBUFF = SCRATCH STORAGE BUFFER FOR INTERNAL VARIABLES
C
      MM = 1
      NN = 1
      IF (ITYPE.NE.2) GO TO 20
      MM = ID - IDD + 1
      DO 10 M=1,IDD
        BUFF(MM) = BUFF(M)
        MM = MM + 1
  10  CONTINUE
      MM = ID - IDD + 1
  20  IF (K.EQ.2) GO TO 100
      IF (K.EQ.3) GO TO 220
C
C ONE STAGE DECIMATION AND/OR INTERPOLATION
C THIS SECTION OF CODE IS REMOVABLE IF 1 STAGE (KD=1) IMPLEMENTATION
C ARE NOT USED        (SUBSTITUTE       200 GO TO 2000)
C ONE STAGE DECIMATION
C READ ID1 SAMPLES FROM BUFF INTO S1D BUFFER
C
C
  30  IF (ITYPE.NE.2) GO TO 40
      SUM = BUFF(MM)
      MM = MM + 1
      GO TO 70
  40  M1 = J1S - ID1
      MD = J1S + N1P
  50  SBUFF(J1S) = BUFF(MM)
      SBUFF(MD) = BUFF(MM)
      MM = MM + 1
      J1S = J1S - 1
      MD = MD - 1
      IF (M1.LT.J1S) GO TO 50
C
C COMPUTE ONE FILTER OUTPUT FOR STAGE 1D
C
      JD = J1S + NCF1
      JU = J1S + NF1
      SUM = 0.0
      M1 = 0
  60  M1 = M1 + 1
      SUM = SUM + COEF1(M1)*(SBUFF(JU)+SBUFF(JD))
      JD = JD - 1
      JU = JU + 1
      IF (M1.LT.NCF1) GO TO 60
      IF (J1S.LE.J1) J1S = J1 + N1P
      IF (ITYPE.NE.1) GO TO 70
      BUFF(NN) = SUM
      NN = NN + 1
      IF (MM.LE.ID) GO TO 40
      GO TO 390
C
C ONE STAGE INTERPOLATION
C STORE DATA INTO S1I
C
  70  SBUFF(K1S) = SUM
      MD = K1S + IQ1
      SBUFF(MD) = SUM
      MQ = J4
C
C COMPUTE ID1 SAMPLES FROM STAGE 1I
C
  80  M1 = K1S
      SUM = 0.0
  90  MQ = MQ + 1
      SUM = SUM + SBUFF(MQ)*SBUFF(M1)
      M1 = M1 + 1
      IF (M1.LT.MD) GO TO 90
C
C STORE OUTPUT INTO BUFF
C
      BUFF(NN) = SUM
      NN = NN + 1
      IF (MQ.LT.J5) GO TO 80
      K1S = K1S - 1
      IF (K1S.LE.J7) K1S = J7 + IQ1
      IF (MM.LE.ID) GO TO 30
      GO TO 390
C
C TWO STAGE DECIMATION AND/OR INTERPOLATION:
C THIS SECTION OF CODE IS REMOVABLE IF 2 STAGE (KD=2) IMPLEMENTATION
C ARE NOT USED        (SUBSTITUTE       300 GO TO 2000)
C TWO STAGE DECIMATION
C
 100  IF (ITYPE.NE.2) GO TO 110
      SUM = BUFF(MM)
      MM = MM + 1
      GO TO 160
 110  M2 = J2S - ID2
      MQ = J2S + N2P
C
C READ ID1 SAMPLES FROM BUFF INTO S1D BUFFER
C
 120  M1 = J1S - ID1
      MD = J1S + N1P
 130  SBUFF(J1S) = BUFF(MM)
      SBUFF(MD) = BUFF(MM)
      MM = MM + 1
      J1S = J1S - 1
      MD = MD - 1
      IF (M1.LT.J1S) GO TO 130
C
C COMPUTE ONE FILTER OUTPUT FOR STAGE 1D
C
      JD = J1S + NCF1
      JU = J1S + NF1
      SUM = 0.0
      M1 = 0
 140  M1 = M1 + 1
      SUM = SUM + COEF1(M1)*(SBUFF(JD)+SBUFF(JU))
      JD = JD - 1
      JU = JU + 1
      IF (M1.LT.NCF1) GO TO 140
      IF (J1S.LE.J1) J1S = J1 + N1P
C
C STORE DATA INTO S2D
C
      SBUFF(J2S) = SUM
      SBUFF(MQ) = SUM
      J2S = J2S - 1
      MQ = MQ - 1
      IF (M2.LT.J2S) GO TO 120
C
C COMPUTE ONE FILTER OUTPUT FOR STAGE 2D
C
      JD = J2S + NCF2
      JU = J2S + NF2
      SUM = 0.0
      M2 = 0
 150  M2 = M2 + 1
      SUM = SUM + COEF2(M2)*(SBUFF(JD)+SBUFF(JU))
      JD = JD - 1
      JU = JU + 1
      IF (M2.LT.NCF2) GO TO 150
      IF (J2S.LE.J2) J2S = J2 + N2P
      IF (ITYPE.NE.1) GO TO 160
      BUFF(NN) = SUM
      NN = NN + 1
      IF (MM.LE.ID) GO TO 110
      GO TO 390
C
C TWO STAGE INTERPOLATION
C STORE DATA INTO S2I
C
 160  SBUFF(K2S) = SUM
      MD = K2S + IQ2
      SBUFF(MD) = SUM
      IDX = J10 + 1
      MQ = J5
C
C COMPUTE OUTPUT FOR STAGE 2I
C
 170  M2 = K2S
      SUM = 0.0
 180  MQ = MQ + 1
      SUM = SUM + SBUFF(MQ)*SBUFF(M2)
      M2 = M2 + 1
      IF (M2.LT.MD) GO TO 180
C
C STORE OUTPUT IN T2I
C
      SBUFF(IDX) = SUM
      IDX = IDX + 1
      IF (IDX.LE.J11) GO TO 170
      K2S = K2S - 1
      IF (K2S.LE.J8) K2S = J8 + IQ2
      IDX = J10 + 1
C
C STORE DATA INTO S1I
C
 190  SBUFF(K1S) = SBUFF(IDX)
      MD = K1S + IQ1
      SBUFF(MD) = SBUFF(IDX)
      IDX = IDX + 1
      MQ = J4
 200  M1 = K1S
      SUM = 0.0
 210  MQ = MQ + 1
      SUM = SUM + SBUFF(MQ)*SBUFF(M1)
      M1 = M1 + 1
      IF (M1.LT.MD) GO TO 210
      BUFF(NN) = SUM
      NN = NN + 1
      IF (MQ.LT.J5) GO TO 200
      K1S = K1S - 1
      IF (K1S.LE.J7) K1S = J7 + IQ1
      IF (IDX.LE.J11) GO TO 190
      IF (MM.LE.ID) GO TO 100
      GO TO 390
C
C THREE STAGE DECIMATION AND/OR INTERPOLATION
C THIS SECTION OF CODE IS REMOVABLE IF 3 STAGE (KD=3) IMPLEMENTATION
C THREE STAGE DECIMATION
C
 220  IF (ITYPE.NE.2) GO TO 230
      SUM = BUFF(MM)
      MM = MM + 1
      GO TO 300
 230  M3 = J3S - ID3
      M = J3S + N3P
 240  M2 = J2S - ID2
      MQ = J2S + N2P
C
C READ ID1 SAMPLES FROM BUFF INTO S1D BUFFER
C
 250  M1 = J1S - ID1
      MD = J1S + N1P
 260  SBUFF(J1S) = BUFF(MM)
      SBUFF(MD) = BUFF(MM)
      MM = MM + 1
      J1S = J1S - 1
      MD = MD - 1
      IF (M1.LT.J1S) GO TO 260
C
C COMPUTE ONE FILTER OUTPUT FOR STAGE 1D
C
      JD = J1S + NCF1
      JU = J1S + NF1
      SUM = 0.0
      M1 = 0
 270  M1 = M1 + 1
      SUM = SUM + COEF1(M1)*(SBUFF(JD)+SBUFF(JU))
      JD = JD - 1
      JU = JU + 1
      IF (M1.LT.NCF1) GO TO 270
      IF (J1S.LE.J1) J1S = J1 + N1P
C
C STORE DATA INTO S2D
C
      SBUFF(J2S) = SUM
      SBUFF(MQ) = SUM
      J2S = J2S - 1
      MQ = MQ - 1
      IF (M2.LT.J2S) GO TO 250
C
C COMPUTE ONE FILTER OUTPUT FOR STAGE 2D
C
      JD = J2S + NCF2
      JU = J2S + NF2
      SUM = 0.0
      M2 = 0
 280  M2 = M2 + 1
      SUM = SUM + COEF2(M2)*(SBUFF(JD)+SBUFF(JU))
      JD = JD - 1
      JU = JU + 1
      IF (M2.LT.NCF2) GO TO 280
      IF (J2S.LE.J2) J2S = J2 + N2P
C
C STORE DATA INTO S3D
C
      SBUFF(J3S) = SUM
      SBUFF(M) = SUM
      J3S = J3S - 1
      M = M - 1
      IF (M3.LT.J3S) GO TO 240
C
C COMPUTE ONE FILTER OUTPUT FOR STAGE 3D
C
      JD = J3S + NCF3
      JU = J3S + NF3
      SUM = 0.0
      M3 = 0
 290  M3 = M3 + 1
      SUM = SUM + COEF3(M3)*(SBUFF(JD)+SBUFF(JU))
      JD = JD - 1
      JU = JU + 1
      IF (M3.LT.NCF3) GO TO 290
      IF (J3S.LE.J3) J3S = J3 + N3P
      IF (ITYPE.NE.1) GO TO 300
      BUFF(NN) = SUM
      NN = NN + 1
      IF (MM.LE.ID) GO TO 230
      GO TO 390
C
C THREE STAGE INTERPOLATION
C STORE DATA INTO S3I
C
 300  SBUFF(K3S) = SUM
      MD = K3S + IQ3
      SBUFF(MD) = SUM
      M = J11 + 1
      MQ = J6
C
C COMPUTE OUTPUT FOR STAGE 3I
C
 310  M3 = K3S
      SUM = 0.0
 320  MQ = MQ + 1
      SUM = SUM + SBUFF(MQ)*SBUFF(M3)
      M3 = M3 + 1
      IF (M3.LT.MD) GO TO 320
C
C STORE OUTPUT IN T3I
C
      SBUFF(M) = SUM
      M = M + 1
      IF (M.LE.J12) GO TO 310
      K3S = K3S - 1
      IF (K3S.LE.J9) K3S = J9 + IQ3
      M = J11 + 1
C
C STORE DATA INTO S2I
C
 330  SBUFF(K2S) = SBUFF(M)
      MD = K2S + IQ2
      SBUFF(MD) = SBUFF(M)
      M = M + 1
      IDX = J10 + 1
      MQ = J5
C
C COMPUTE OUTPUT FOR STAGE 2I
C
 340  M2 = K2S
      SUM = 0.0
 350  MQ = MQ + 1
      SUM = SUM + SBUFF(MQ)*SBUFF(M2)
      M2 = M2 + 1
      IF (M2.LT.MD) GO TO 350
C
C STORE OUTPUT IN T2I
C
      SBUFF(IDX) = SUM
      IDX = IDX + 1
      IF (IDX.LE.J11) GO TO 340
      K2S = K2S - 1
      IF (K2S.LE.J8) K2S = J8 + IQ2
      IDX = J10 + 1
C
C STORE DATA INTO S1I
C
 360  SBUFF(K1S) = SBUFF(IDX)
      MD = K1S + IQ1
      SBUFF(MD) = SBUFF(IDX)
      IDX = IDX + 1
      MQ = J4
C
C COMPUTE OUTPUT FOR STAGE 1I
C
 370  M1 = K1S
      SUM = 0.0
 380  MQ = MQ + 1
      SUM = SUM + SBUFF(MQ)*SBUFF(M1)
      M1 = M1 + 1
      IF (M1.LT.MD) GO TO 380
C
C STORE OUTPUT IN BUFF
C
      BUFF(NN) = SUM
      NN = NN + 1
      IF (MQ.LT.J5) GO TO 370
      K1S = K1S - 1
      IF (K1S.LE.J7) K1S = J7 + IQ1
      IF (IDX.LE.J11) GO TO 360
      IF (M.LE.J12) GO TO 330
      IF (MM.LE.ID) GO TO 220
 390  RETURN
      END
Return to Main Software Page