Make your own free website on Tripod.com
FASTFILT - An FFT Based Filtering Program

Return to Main Software Page

C
C-----------------------------------------------------------------------
C MAIN PROGRAM: FASTFILT - 8/1/75
C AUTHOR:       JONT B. ALLEN
C               BELL LABS,MURRAY HILL N.J., 07974
C
C INPUT:
C               A)0=FILTERING;1=CORRELATION
C               B)LENGTH OF FILTER SPECIFIED IN DISK BLOCKS
C               C)IMPULSE RESPONSE INPUT FILENAME
C               THIS MAY BE 'TTI' OR ANY DISK FILENAME
C               D)IF INPUT IS FROM 'TTI',THE IMPULSE RESP
C               IS NOW TYPED IN AS DELAY AND GAIN.
C               GAIN MUST BE LESS THAN OR EQUAL TO 1. IN
C               MAGNITUDE.
C               D)IF INPUT IS FROM DISK FILE, AN OVERALL GAIN
C               VALUE IS REQUESTED; THIS WILL BE USED TO SCALE
C               THE IMPULSE RESPONSE AFTER IT IS READ FROM DISK.
C               E)THE INPUT DISK FILENAME WHERE THE SIGNAL TO BE
C               FILTERED IS STORED IN INTEGER (BINARY) FORMAT.
C               F)OUTPUT DISK FILENAME WHERE FILTERED OUTPUT IS TO
C               G)NUMBER OF DISK BLOCKS TO BE FILTERED
C
C  IF THE USER OPTS THE CORRELATION OPTION, THE IMPULSE RESP. IS
C      TIME REVERSED.
C
C
C      USES BERGLAND REAL FFT TO INCREASE SIZE OF MAXIMUM FILTER LENGTH
C      PGM TO FILTER DATA OFF OF DISK BY FFT OVERLAP-SAVE FILTERING
C      THE FILTER IMPULSE RESPONSE DATA IS EITHER SPECIFIED IN A DISK
C      FILE, OR FROM THE TTY
C
C SUBROUTINES:  RFILT,RBLK,WBLK,FOPEN,CLOSE,DKIO,GNAME,I1MACH,FAST,FSST
C
C RFILT:        DOES OVERLAP-ADD FILTERING;MAY BE USED ALONE
C RBLK,WBLK,FOPEN,CLOSE,DKIO:   DISK FILE HANDELING ROUTINES.
C
C*** THERE ARE 2 MACHINE DEPENDENT SUBROUTINES WHICH MUST BE
C*** SUPPLIED BY THE USER - DKIO AND GNAME - COMPLETE INSTRUCTIONS
C*** FOR WRITING THEM ARE GIVEN IN THE ROUTINES. FORTRAN 5
C*** CODE IS SUPPLIED FOR A DATA GENERAL COMPUTER.
C***    DKIO - OPENS, CLOSES, READS AND WRITES A DISK FILE
C***    GNAME - READS AN ASCII FILENAME FROM THE INPUT DEVICE(TTY)
C      I1MACH: MACHINE DEPENDENT; SETS MACHINE CONSTANTS
C      FAST,FSST: BERGLAND RADIX-4 REAL TO COMPLEX FFT'S
C-----------------------------------------------------------------------
C
      DIMENSION IS(4096), F(4098), R(4098), S(2049)
      DIMENSION IFIN(10), IFOUT(10)
      LOGICAL COR
      EQUIVALENCE (S(1),IS(1))
      DATA IDSK0 /0/, IDSK1 /1/, IB0 /0/
      DATA NBLKL /256/, NMAX /4096/, COR /.FALSE./
      ITTO = I1MACH(2)
      ITTI = I1MACH(1)
      SCALE = I1MACH(9)
      NMAXD2 = NMAX/2
C
C INPUT INFO FROM TTY
C
      WRITE (ITTO,9999)
9999  FORMAT (33H DO YOU WISH TO DO FILTERING (=0)/17H OR CORRELATION (,
     *    25HMATCHED FILTERING) (=1) =)
      READ (ITTI,9998) I
9998  FORMAT (I1)
      IF (I.EQ.1) COR = .TRUE.
      IF (COR) WRITE (ITTO,9997)
9997  FORMAT (28H PROGRAM WILL DO CORRELATION)
      IF (.NOT.COR) WRITE (ITTO,9996)
9996  FORMAT (26H PROGRAM WILL DO FILTERING)
  10  CONTINUE
C
C INPUT FILTER LENGTH
C
      WRITE (ITTO,9995) NBLKL
9995  FORMAT (18H FILTER LENGTH IN , I3/28H WORD DISK BLOCKS (1,2,4 OR ,
     *    3H8)=)
      READ (ITTI,9994) NB
9994  FORMAT (I1)
      IF (NB.LE.0) GO TO 10
      NP = NB*NBLKL
      IF (NP.GT.NMAXD2) WRITE (ITTO,9993)
9993  FORMAT (23H FILTER LENGTH TOO LONG)
      IF (NP.GT.NMAXD2) GO TO 10
      WRITE (ITTO,9992) NP
9992  FORMAT (16H FILTER WILL BE , I4, 12H POINTS LONG)
      NPT2 = NP*2
      NPP1 = NP + 1
C
C ZERO FILTER
C
      DO 20 I=1,NPT2
        F(I) = 0
  20  CONTINUE
      WRITE (ITTO,9991)
9991  FORMAT (37H INPUT FILTER FROM DISK(1) OR TTY(2)=)
      READ (ITTI,9990) IRESP
9990  FORMAT (I1)
      IF (IRESP.EQ.1) GO TO 60
C
C TTY INPUT OF IMPULSE RESPONSE
C
  30  CONTINUE
      WRITE (ITTO,9989)
9989  FORMAT (32H DELAY=0 TO TERMINATE INPUT MODE/12H DELAY= (I4))
      READ (ITTI,9988) IDEL
9988  FORMAT (I4)
C
C CHECK FOR INPUT TERMINATION
C
      IF (IDEL.LE.0) GO TO 40
C
C CHECK FOR TO LARGE A DELAY
C
      IF (IDEL.GT.NP) GO TO 30
C
C READ IN GAIN VALUE OF IMPULSE RESP. AT DELAY IDEL
C
      WRITE (ITTO,9987)
9987  FORMAT (36H GAIN=(F10.1), ABS(GAIN) LESS THAN 1)
      READ (ITTI,9986) F(IDEL)
9986  FORMAT (F10.1)
      IF (ABS(F(IDEL)).GT.1.) WRITE (ITTO,9985)
9985  FORMAT (15H GAIN TOO LARGE)
      IF (ABS(F(IDEL)).GT.1.) F(IDEL) = 0
C
C GET A NEW RESP VALUE
C
      GO TO 30
  40  CONTINUE
C
C WRITE IMPULSE RESPONSE DATA ON DISK IN FILE "FILTDATA"
C OPEN A DISK FILE CALLED "FILTDATA" ON CHANNEL IDSK0
C AND SCALE DATA
C
      CALL FOPEN(IDSK0, 8HFILTDATA)
      DO 50 I=1,NP
        IS(I) = SCALE*F(I)
  50  CONTINUE
C
C PUT DATA ON DISK IN BINARY FORMAT
C WRITE NB BLOCKS STARTING AT BLOCK IB0 ON CHANNEL IDSK0
C FROM ARRAY IS; IER IS UNUSED ERROR FLAG
C
      CALL WBLK(IDSK0, IB0, IS, NB, IER)
      CALL CLOSE(IDSK0)
      WRITE (ITTO,9984)
9984  FORMAT (47H IMPULSE RESPONSE IS ON DISK IN FILE - FILTDATA)
      GO TO 80
C
C READ FILTER H(T) FROM DISK AND PUT IN ARRAY F
C
  60  CONTINUE
      WRITE (ITTO,9983)
9983  FORMAT (17H FILTER FILENAME=)
      CALL GNAME(IFIN)
      WRITE (ITTO,9982)
9982  FORMAT (20H SYSTEM GAIN=(F10.1))
      READ (ITTI,9981) GAIN
9981  FORMAT (F10.1)
C
C OPEN DISK ON CHANNEL IDSK0 ; NAME OF FILE IS IN ARRAY IFIN
C AND READ IN IMPULSE RESPONSE
C
      CALL FOPEN(IDSK0, IFIN)
      CALL RBLK(IDSK0, IB0, IS, NB, IER)
      CALL CLOSE(IDSK0)
C
C SCALE H(T) TO 'GAIN' MAXIMUM
C
      DO 70 I=1,NP
        F(I) = GAIN*FLOAT(IS(I))/SCALE
  70  CONTINUE
C
C FILTER RESP HAS BEEN DEFINED
C
  80  CONTINUE
      WRITE (ITTO,9980)
9980  FORMAT (21H INPUT DATA FILENAME=)
      CALL GNAME(IFIN)
      WRITE (ITTO,9979)
9979  FORMAT (22H OUTPUT DATA FILENAME=)
      CALL GNAME(IFOUT)
      CALL FOPEN(IDSK0, IFIN)
      CALL FOPEN(IDSK1, IFOUT)
      WRITE (ITTO,9978)
9978  FORMAT (43H NUMBER OF DATA BLOCKS TO BE FILTERED= (I5))
      READ (ITTI,9977) NBLKS
9977  FORMAT (I5)
C
C FILTER DATA FROM FILE IFIN AND PUT IN FILE IFOUT
C IF COR IS TRUE, REVERSE IMPULSE RESP
C
      IF (.NOT.COR) GO TO 100
      NPD2 = NP/2
      DO 90 I=1,NPD2
        J = NP + 1 - I
        F0 = F(J)
        F(J) = F(I)
        F(I) = F0
  90  CONTINUE
 100  CONTINUE
C
C FFT FILTER H(T) TO GET FILTER SPECTRUM
C FORCE SECOND HALF OF FILTER TO ZERO; ZERO R SECOND HALF
C
      DO 110 I=NPP1,NPT2
        F(I) = 0
        R(I) = 0
 110  CONTINUE
      CALL FAST(F, NPT2)
C
C START LOOP OVER FRAMES OF DATA
C
      IFRAME = 0
      DO 140 IB=1,NBLKS,NB
C
C READ IN FRAME OF DATA FROM DISK CHANNEL IDSK0
C AND STORE FLOATED DATA IN R ARRAY
C
        CALL RBLK(IDSK0, IFRAME, IS, NB, IER)
        DO 120 I=1,NP
          R(I) = IS(I)
 120    CONTINUE
C
C FILTER FRAME OF DATA
C
        CALL RFILT(R, F, S, NP)
C
C OUTPUT PRESENT FRAME OF DATA TO OUTPUT DISK FILE IDSK1
C
        DO 130 I=1,NP
          IF (ABS(R(I)).GT.SCALE) WRITE (ITTO,9976) IB, I
9976      FORMAT (32H OUTPUT OVERFLOW BLOCK AND WORD , 2I5)
          IS(I) = R(I)
 130    CONTINUE
C
C WRITE DATA ON DISK IN BINARY FORMAT
C
        CALL WBLK(IDSK1, IFRAME, IS, NB, IER)
        WRITE (ITTO,9975) IB
9975    FORMAT (7H BLOCK , I6, 12H IS FINISHED)
        IFRAME = IFRAME + NB
 140  CONTINUE
C
C WRITE OUT LAST FRAME
C
      DO 150 I=NPP1,NPT2
        IF (ABS(R(I)).GT.SCALE) WRITE (ITTO,9976) IB, I
        J = I - NP
        IS(J) = R(I)
 150  CONTINUE
      CALL WBLK(IDSK1, IFRAME, IS, NB, IER)
      CALL CLOSE(IDSK0)
      CALL CLOSE(IDSK1)
      STOP
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: RFILT
C FILTER ONE FRAME (I.E., NP POINTS) OF DATA
C
C PROGRAM ASSUMES:
C      1. R(NP+1) TO R(2*NP) HAS NOT BEEN USED SINCE LAST CALL TO
C      'RFILT' AND WAS ZERO BEFORE FIRST CALL.
C      2. INPUT AND OUTPUT DATA (TIME SERIES) ARE IN R(1) TO R(NP)
C      3. F(1) TO F(2*NP+2) IS THE FILTER FREQUENCY RESP.
C      4. DATA ARE STORED IN F AS: DC,0.,F1REAL,F1IMAG,
C      F2REAL,F2IMAG,...,FSAMPLING/2.,0.
C      5. S(1) TO S(NP) IS A SCRATCH ARRAY; IT MAY BE USED IN 'MAIN
C      6. IMPULSE RESP. OF FILTER F IS ZERO FOR NP+1 THRU 2*NP
C      (I.E. SECOND HALF OF H(T)=0, WHERE F(W)=FFT(H) )
C-----------------------------------------------------------------------
C
      SUBROUTINE RFILT(R, F, S, NP)
      DIMENSION S(1), R(1), F(1)
      NPT2 = NP*2
C
C STORE PREVIOUS TAIL IN SCRATCH ARRAY S(.);
C ZERO SECOND HALF OF R--FIRST HALF OF R IS NEW DATA TO BE FILTERED
C
      DO 10 K=1,NP
        NPPK = NP + K
        S(K) = R(NPPK)
        R(NPPK) = 0.
  10  CONTINUE
C
C TAKE FFT OF DATA
C
      CALL FAST(R, NPT2)
C
C HANDLE VALUES AS COMPLEX VALUES
C
      KMAX = NPT2 + 1
      DO 20 K=1,KMAX,2
        X = F(K)*R(K) - F(K+1)*R(K+1)
        Y = F(K)*R(K+1) + F(K+1)*R(K)
        R(K) = X
        R(K+1) = Y
  20  CONTINUE
C
C INVERSE TRANSFORM PRODUCT
C
      CALL FSST(R, NPT2)
C
C ADD IN TAIL FROM PREVIOUS FRAME WHICH WAS STORED IN S(.)
C
      DO 30 K=1,NP
        R(K) = R(K) + S(K)
  30  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C PROGRAM: RBLK
C PGM TO READ DATA FROM DISK
C-----------------------------------------------------------------------
C
      SUBROUTINE RBLK(ICN, IB0, IARRAY, NBLK, IER)
      DIMENSION NAME(1), IARRAY(1)
      CALL DKIO(NAME, ICN, IB0, IARRAY, NBLK, IER, 2)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: WBLK
C PGM TO WRITE DATA ON DISK
C-----------------------------------------------------------------------
C
      SUBROUTINE WBLK(ICN, IB0, IARRAY, NBLK, IER)
      DIMENSION NAME(1), IARRAY(1)
      CALL DKIO(NAME, ICN, IB0, IARRAY, NBLK, IER, 3)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  FOPEN
C PGM TO OPEN A DISC FILE ON CHANNEL ICN WITH NAME IN ARRAY "NAME"
C-----------------------------------------------------------------------
C
      SUBROUTINE FOPEN(ICN, NAME)
      DIMENSION NAME(1), IARRAY(1)
      CALL DKIO(NAME, ICN, 0, IARRAY, 0, IER, 0)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: CLOSE
C PGM TO CLOSE A DISK FILE ON CHANNEL ICN
C-----------------------------------------------------------------------
C
      SUBROUTINE CLOSE(ICN)
      DIMENSION NAME(1), IARRAY(1)
      CALL DKIO(NAME, ICN, 0, IARRAY, 0, IER, 1)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: DKIO
C *** THIS PROGRAM IS MACHINE DEPENDENT
C *** FORTRAN CODE HAS BEEN SUPPLIED FOR A DATA GENERAL COMPUTER
C *** WITH A FORTRAN 5 COMPILER
C ***
      SUBROUTINE DKIO(NAME, ICN, IB0, IARRAY, NBLK, IER, IOPER)
      DIMENSION NAME(10), IARRAY(1)
C PGM TO OPEN, CLOSE, READ, AND WRITE DATA FROM BULK STORAGE
C DEVICE.
C
C     IOPER    OPERATION
C      0       OPEN
C      1       CLOSE DISC
C      2       READ FROM DISK
C      3       WRITE ON DISC
C
C      ICN=FORTRAN CHANNEL NUMBER FOR BULK STORAGE
C      NAME=FORTRAN ARRAY CONTAINING ASCII NAME OF DISC FILE
C      IB0=BLOCK TO BE READ OR WRITTEN FROM DISK FILE
C      IARRAY=ARRAY WHERE DATA IS TO BE TRANSFERED TO OR FROM
C      NBLK=NUMBER OF BLOCKS TO BE TRANSFERED EACH CALL
C      IER=ERROR FLAG (PARAMETER NOT USED IN THIS PROGRAM)
C-----------------------------------------------------------------------
C
      LENBLK = 256
      IOPER1 = IOPER + 1
      GO TO (10, 20, 30, 40), IOPER1
C
C *** OPEN - INTRODUCE A FILE ON A SPECIFIC CHANNEL TO THE
C ***        OPERATING SYSTEM PRIOR TO READING OR WRITING DATA IN IT
C
  10  CONTINUE
       OPEN ICN,NAME
      RETURN
C
C *** CLOSE - RELEASE THE FILE AND CHANNEL NUMBER FROM THE
C ***         OPERATING SYSTEM
C
  20  CLOSE ICN
      RETURN
C
C *** RDBLK - READS NBLK BLOCKS (1 BLOCK=256 WORDS; 1WORD=16 BITS) OF
C ***         DATA STARTING AT BLOCK IB0 INTO ARRAY IARRAY (FIRST
C ***         BLOCK IN FILE IS NUMBER 0)
C
  30  CALL RDBLK(ICN, IB0, IARRAY, NBLK, IER)
      RETURN
C
C *** WRBLK - WRITES NBLK BLOCKS OF DATA STARTING AT BLOCK IB0 FROM
C ***         ARRAY IARRAY INTO DISK FILE OPENED ON CHANNEL ICN
C
  40  CALL WRBLK(ICN, IB0, IARRAY, NBLK, IER)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE: GNAME
C ***
C ***  THIS PROGRAM IS MACHINE DEPENDENT
C ***  FORTRAN CODE HAS BEEN SUPPLIED FOR A DATA GENERAL COMPUTER
C ***  WITH A FORTRAN 5 COMPILER
C ***
C THIS PROGRAM READS ASCII DATA INTO AN ARRAY 'NAME(I)'
C IN A FORMAT THAT MAY BE USED BY : CALL FOPEN(ICN,NAME)
C WHICH OPENS DISK FILE 'NAME' ON FORTRAN CHANNEL ICN
C-----------------------------------------------------------------------
C
      SUBROUTINE GNAME(NAME)
      DIMENSION NAME(10)
      ITTI = I1MACH(1)
C
C READ UP TO 10 CHARACTERS FROM DEVICE ITTI IN S (STRING) FORMAT
C THE CHARACTERS ARE PACKED 2 PER 16 BIT WORD AND ARE LEFT
C JUSTIFIED IN THE ARRAY NAME.
C
      READ (ITTI,9999) NAME(1)
9999  FORMAT (S10)
      RETURN
      END
C
C-----------------------------------------------------------------------
C MAIN PROGRAM: TESTFILT - TEST PROGRAM FOR RFILT SUBROUTINE - 1/8/79
C SUBROUTINES:  RFILT,FAST,FSST,R1MACH,I1MACH
C-----------------------------------------------------------------------
C
      DIMENSION Y(4096), F(258), R(258), S(128)
C
C DEFINE FUNCTION STATEMENT FOR INPUT "FILE" X(*)
C
      X(K) = (AMOD(FLOAT(K-1),10.) - 5.) / 5.
C
C DEFINE CONSTANTS
C
      ITTO = I1MACH(2)
      NPTS = 4096
      NFILT = 128
      NFFT = 2*NFILT
      NBLKL = 64
      NBLKS = NPTS/NBLKL
      NB = NFILT/NBLKL
C
C INITIALIZE FILTER: DELAY OF 1 SAMPLE,CHANGE OF SIGN
C
      DO 10 I=1,NFFT
        F(I) = 0.
        R(I) = 0.
  10  CONTINUE
      F(2) = -1.
C
C FFT FILTER
C
      CALL FAST(F, NFFT)
C
C START LOOP ON DATA
C
      DO 40 IB=1,NBLKS,NB
C
C
C GET BLK OF DATA FROM INPUT "FILE" FUNCTION X(I)
C
        DO 20 I=1,NFILT
          J = I + (IB-1)*NBLKL
          R(I) = X(J)
  20    CONTINUE
C
C FILTER DATA
C
        CALL RFILT(R, F, S, NFILT)
C
C WRITE OUTPUT INTO OUTPUT "FILE" ARRAY Y(I)
C
        DO 30 I=1,NFILT
          J = I + (IB-1)*NBLKL
          Y(J) = R(I)
  30    CONTINUE
  40  CONTINUE
C
C CHECK RESULTS
C
      RMACH4 = R1MACH(4)
      ERMAX = 0.
      FNFFT = NFFT
      RMAX = SQRT(FNFFT)*RMACH4
      NPTS1 = NPTS - 1
      DO 50 I=1,NPTS1
        ERR = X(I) + Y(I+1)
        ERMAX = AMAX1(ERMAX,ABS(ERR))
  50  CONTINUE
      WRITE (ITTO,9999) RMACH4, ERMAX
9999  FORMAT (18H MACHINE ACCURACY=, E12.5, 11H MAX ERROR=, E12.5)
      IF (RMAX.GT.ERMAX) WRITE (ITTO,9998)
9998  FORMAT (12H TEST PASSES)
      IF (RMAX.LE.ERMAX) WRITE (ITTO,9997)
9997  FORMAT (11H TEST FAILS)
      STOP
      END
Return to Main Software Page