Make your own free website on Tripod.com
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