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