A Computer Program for Digital Interpolator Design
Return to Main Software Page
C
C-----------------------------------------------------------------------
C MAIN PROGRAM: TEST PROGRAM FOR OPTIMAL DIGITAL INTERPOLATING FILTER
C AUTHORS: GEERD OETKEN, H. W. SCHUSSLER AND T. W. PARKS
C T. W. PARKS
C DEPT. OF ELECTRICAL ENGINEERING
C RICE UNIV.
C HOUSTON, TEXAS 77001
C
C G. OETKEN AND H. W. SCHUSSLER
C INSTITUT FUER NACHRICHTENTECHNIK
C UNIVERSITAET ERLANGEN-NUERNBERG
C 8520 ERLANGEN
C CAUERSTRASSE 7
C WEST GERMANY
C
C INPUT: NUM IS THE NUMBER OF TIMES TO LOOP CALLING DODIF
C L DETERMINES THE DEGREE OF THE FILTER
C THE DEGREE IS 2*L*R
C ALPHA IS THE CUTOFF FREQUENCY (RADIANS) OF THE FILTER
C INPUT DATA. RANGE: 0.D0 < ALPHA <= 1.D0
C R IS THE WANTED INCREASE IN THE SAMPLING RATE
C RANGE: AN INTEGER > 1
C-----------------------------------------------------------------------
C
DOUBLE PRECISION ALPHA, DIF(100)
INTEGER R
C
C DEFINE I/O DEVICE CODES
C
IN = I1MACH(1)
IOUT = I1MACH(2)
C
READ (IN,9999) NUM
9999 FORMAT (I5)
DO 20 II=1,NUM
READ (IN,9998) L, ALPHA, R
9998 FORMAT (I5, F5.3, I5)
LR = L*R + 1
WRITE (IOUT,9997)
9997 FORMAT (1H0, 35H OPTIMAL INTERPOLATING FILTER DATA)
WRITE (IOUT,9996)
9996 FORMAT (1H0)
WRITE (IOUT,9995) L, R, ALPHA
9995 FORMAT (1H0, 5H L= , I2, 10H R= , I2, 14H ALPHA= ,
* F4.2)
WRITE (IOUT,9996)
CALL DODIF(L, R, ALPHA, DIF)
WRITE (IOUT,9994) LR
9994 FORMAT (1H0, 10H THERE ARE, I5, 28H VALUES IN HALF OF THE FILTE,
* 1HR)
WRITE (IOUT,9996)
J = -LR
DO 10 I=1,LR
J = J + 1
WRITE (IOUT,9993) J, DIF(I)
10 CONTINUE
9993 FORMAT (1H0, 4H H(, I3, 2H)=, E16.7)
20 CONTINUE
STOP
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: DODIF
C SUBROUTINE TO CALCULATE THE OPTIMAL DIGITAL INTERPOLATING
C FILTER
C-----------------------------------------------------------------------
C
SUBROUTINE DODIF(L, R, ALPHA, DIF)
DOUBLE PRECISION AP(100), AM(100), D(200)
DOUBLE PRECISION DBLE, DFLOAT, S1, S2, DPHI, ALPHA, H, X, Y,
* DIF(1)
INTEGER R
C
C INPUT:
C R -WANTED INCREASE IN SAMPLING RATE, I.E. THE FILTER
C CALCULATES R-1 SAMPLES BETWEEN EACH PAIR OF
C INPUT SAMPLES. RANGE FOR R/ ANY INTEGER > 1
C L -INTEGER WHICH DETERMINES THE DEGREE OF THE FILTER.
C THE DEGREE IS 2*L*R.
C RANGE FOR L/ 1 TO 10, MAY BE INCREASED AT THE
C RISK OF NUMERICAL INSTABILITY (CHANGE LMAX AND
C ARRAY BOUNDS OF AP,AM,D FOR L > 10).
C ALPHA -CUTOFF FREQUENCY (RADIANS) OF THE FILTER INPUT DATA
C DIVIDED BY PI/R.
C RANGE FOR ALPHA/ 0.D0< ALPHA <= 1.D0 (DOUBLE
C PRECISION).
C DIF -DOUBLE PRECISION ARRAY IN WHICH THE FIRST HALF (L*R+1
C VALUES) OF THE UNIT SAMPLE RESPONSE IS RETURNED.
C (THE OPTIMAL FILTER HAS LINEAR PHASE).
C
C FORTRAN 4 USAGE/
C CALL DODIF(L,R,ALPHA,DIF)
C
C REMARKS/
C BIG VALUES OF L AND SMALL VALUES OF ALPHA INCREASE THE
C POSSIBILITY OF PROBLEMS WITH MATRIX INVERSION (DEPENDING UPON
C THE WORDLENGTH OF THE COMPUTER). A VERY GOOD POSSIBILITY TO
C CONTROL THE QUALITY OF THE MATRIX INVERSION IS TO CHECK THE
C MAGNITUDE OF THE UNIT SAMPLE RESPONSE DIF(K) AT K=M*R+1,
C M=0 TO L-1, WHICH IDEALLY HAVE TO BE ZERO.
C
C SUBROUTINES AND FUNCTIONS REQUIRED/
C DPHI, SMINVD
C
DATA LMAX /10/
C
C (MAXIMUM VALUE OF L SPECIFIED BY ARRAY BOUNDS)
C
DATA EPS /1.E-6/
C
C (RELATIVE TOLERANCE FOR LOSS OF SIGNIFICANCE TEST IN MATRIX)
C
DFLOAT(IF) = DBLE(FLOAT(IF))
IODEV = I1MACH(2)
C
C (DEVICE NO. ON WHICH ERROR MESSAGES ARE TO HAPPEN)
C
C
C TEST IF L WITHIN CURRENT BOUNDS
C
IF (L.GT.LMAX) GO TO 100
C
C CALCULATE THE TWO L*L MATRICES
C
LR = 0
DO 20 M=1,L
DO 10 K=1,M
K1 = K + LR
S1 = K - M
S1 = DABS(S1)
S2 = DFLOAT(2*L+1-M-K)
S1 = DPHI(S1,ALPHA)
S2 = DPHI(S2,ALPHA)
AP(K1) = S1 + S2
AM(K1) = S1 - S2
10 CONTINUE
LR = LR + M
20 CONTINUE
C
C INVERT BOTH MATRICES USING IBM/SSP SUBROUTINES
C
CALL SMINVD(AP, L, EPS, IER1)
CALL SMINVD(AM, L, EPS, IER2)
IF (IER1.EQ.(-1) .OR. IER2.EQ.(-1)) GO TO 110
C
C CALCULATE THE L*2L MATRIX D
C
DO 40 M=1,L
M1 = (M-1)*L
DO 30 K=1,L
K1 = 2*(M1+K) - 1
K2 = (K-1)*K/2 + M
IF (M.GT.K) K2 = (M-1)*M/2 + K
S1 = AP(K2) + AM(K2)
S2 = AP(K2) - AM(K2)
D(K1) = S1
D(K1+1) = S2
30 CONTINUE
40 CONTINUE
C
C CALCULATE UNIT SAMPLE RESPONSE FROM D AND DPHI
C
J = 0
DO 90 I=1,L
I1 = (I-1)*2*L
DO 80 M=1,R
M1 = M - 1
GO TO 60
50 M1 = R
60 X = DFLOAT(M1)/DFLOAT(R)
H = 0.D0
DO 70 K=1,L
Y = DFLOAT(L-(K-1)) - X
IT = I1 + 2*K - 1
H = H + D(IT)*DPHI(Y,ALPHA)
Y = DFLOAT(L-K) + X
IT = I1 + 2*K
H = H + D(IT)*DPHI(Y,ALPHA)
70 CONTINUE
J = J + 1
C
C FIRST HALF OF UNIT SAMPLE RESPONSE IS STORED IN DIF(J)
C
DIF(J) = .5D0*H
IF (I.EQ.L .AND. M1.EQ.(R-1)) GO TO 50
80 CONTINUE
90 CONTINUE
RETURN
C
C ERROR MESSAGES
C
100 WRITE (IODEV,9999) LMAX
9999 FORMAT (18H L IS GREATER THAN, I4)
GO TO 120
110 WRITE (IODEV,9998)
9998 FORMAT (54H ERROR DURING MATRIX INVERSION, ILL CONDITIONED SYSTEM)
120 RETURN
END
C
C-----------------------------------------------------------------------
C FUNCTION: DPHI
C FUNCTION REQUIRED BY DODIF
C SPECIAL VERSION FOR LOWPASS INTERPOLATORS
C-----------------------------------------------------------------------
C
DOUBLE PRECISION FUNCTION DPHI(X, ALPHA)
DOUBLE PRECISION X, ALPHA, PI
PI = 4.D0*DATAN(1.D0)
IF (DABS(X).LT.1.D-10) GO TO 10
X = ALPHA*PI*X
DPHI = DSIN(X)/X
RETURN
10 CONTINUE
DPHI = 1.D0
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: SMINVD
C AUTHOR: DEAN KOLBA
C DEPT. OF ELEC. ENG.
C RICE UNIVERSITY
C HOUSTON, TX. 77001
C
C SUBROUTINE TO INVERT A POSITIVE-DEFINITE SYMMETRIC MATRIX
C-----------------------------------------------------------------------
C
SUBROUTINE SMINVD(A, N, EPS, IER)
DOUBLE PRECISION A(1), SUM
C
C CALCULATE THE CHOLESKY DECOMPOSITION
C
IER = 0
ID = 0
DO 70 K=1,N
ID = ID + K
IV = ID
KM1 = K - 1
TEST = EPS*SNGL(A(ID))
DO 60 I=K,N
SUM = 0.D0
IF (KM1.EQ.0) GO TO 20
DO 10 J=1,KM1
J1 = ID - J
J2 = IV - J
SUM = SUM + A(J1)*A(J2)
10 CONTINUE
20 SUM = A(IV) - SUM
IF (I.GT.K) GO TO 40
IF (SNGL(SUM).GT.ABS(TEST)) GO TO 30
IF (SUM.LE.0.D0) GO TO 150
IF (IER.GT.0) GO TO 30
IER = KM1
30 A(ID) = DSQRT(SUM)
GO TO 50
40 A(IV) = SUM/A(ID)
50 IV = IV + I
60 CONTINUE
70 CONTINUE
C
C INVERT UPPER TRIANGULAR MATRIX RESULTING FROM CHOLESKY DECOMPOSITION
C
NV = N*(N+1)/2
ID = NV
DO 110 K=1,N
A(ID) = 1.D0/A(ID)
JH = N
KM1 = K - 1
IF (KM1.EQ.0) GO TO 100
JL = N - KM1
IV = NV - K + 1
DO 90 I=1,KM1
SUM = 0.D0
JH = JH - 1
J1 = ID
J2 = IV
DO 80 J=JL,JH
J1 = J1 + J
J2 = J2 + 1
SUM = SUM + A(J1)*A(J2)
80 CONTINUE
A(IV) = -SUM*A(ID)
IV = IV - JH
90 CONTINUE
100 ID = ID - JH
110 CONTINUE
C
C CALCULATE THE INVERSE OF A USING THE INVERTED TRIANGULAR MATRIX
C
ID = 0
DO 140 K=1,N
ID = ID + K
IV = ID
DO 130 I=K,N
SUM = 0.D0
J1 = IV
JL = I - K
DO 120 J=I,N
J2 = J1 + JL
SUM = SUM + A(J1)*A(J2)
J1 = J1 + J
120 CONTINUE
A(IV) = SUM
IV = IV + I
130 CONTINUE
140 CONTINUE
RETURN
150 IER = -1
RETURN
C
C IER=-1 INVERSION FAILED
C IER=0 GOOD INVERSION
C IER=I LOSS OF SIGNIFICANCE AT THE I+1 STAGE
C
END
Return to Main Software Page