PROGRAM BESSEL_LOWPASS_FILTER
C
C  BESSEL_LOWPASS_FILTER.F   ver 1.0  November 28, 2012
C
C  By Tom Irvine
C  Email:  tom@vibrationdata.com
C
C  Fourier transform of a time history
C
C     Variables:
C
C        T = input time
C        Y = input amplitude
C
C********************************************************
C23456789
C
      DOUBLE PRECISION, ALLOCATABLE :: TT(:),Y(:),YF(:)
      DOUBLE PRECISION AMP,DT,DF,DUR
      DOUBLE PRECISION MS,SR
      DOUBLE PRECISION TTA
      DOUBLE PRECISION NY
      DOUBLE PRECISION SCALE,FC
C
      INTEGER I
      INTEGER ios
      INTEGER stat_alloc
C
      CHARACTER :: A_file*20,B_file*20
C
      AVE=0
C
      PRINT '("")'
      PRINT '("  The input file must have two columns: ")'
      PRINT '("  time(sec) & amplitude (engineering unit) ")'
      PRINT '("")'
C
      ios=1
      DO WHILE(ios.NE.0)
        PRINT '( " Enter the input filename:")'
        READ(*,*)A_file
C
      OPEN(UNIT=10,FILE=A_file,STATUS='OLD',ACTION='READ',IOSTAT=ios)
C
        IF(ios.EQ.0)THEN
            WRITE(*,720)
720         FORMAT(/,'File opened.',/)
        ELSE
            WRITE(*,721)
721         FORMAT(/,' Error: File does not exist.',/)
        ENDIF
C
      ENDDO
C
      I=1
      READ(10,*,IOSTAT=IOstatus) TT1,AMP
      I=2
      IFLAG=0
      DO WHILE(IFLAG.EQ.0)
            READ(10,*,IOSTAT=IOstatus) TTA,AMP
            IF (IOstatus.EQ.0)THEN
                NUM=I
                I=I+1
            ELSE
                NUM=I-1
                IFLAG=1
            ENDIF
      ENDDO
C
      DUR=TTA-TT1
      DT=DUR/FLOAT(NUM-1)
      SR=1/DT
C
      WRITE(*,701)DUR
      WRITE(*,702)DT
      WRITE(*,703)SR
      WRITE(*,704)NUM
C
701   FORMAT(/,'  DUR = ',G9.3,' sec ')
702   FORMAT('   DT = ',G9.3,' sec ')
703   FORMAT('   SR = ',G11.4,' samples/sec ')
704   FORMAT('  NUM = ',I8)
C
      REWIND(10)
C
      ALLOCATE(Y(NUM))
      ALLOCATE(YF(NUM))
      ALLOCATE(TT(NUM))
C
      DO I=1,NUM
        READ(10,*,IOSTAT=IOstatus) TT(I),Y(I)
      ENDDO
C
      CLOSE(10)
C
      WRITE(*,*)' '
      WRITE(*,*)' Input Time History Statistics '
      CALL STATISTICS(Y,NUM,MEAN,STD_DEV,RMS)
C
      SCALE=SQRT(3.*(-1.+SQRT(5.))/2.)
C
      NY=SR*0.5
C
      WRITE(*,*)
      DO WHILE(1==1)
         WRITE(*,*)' Enter lowpass frequency (Hz) '
         READ(*,*)FC
         IF(FC.LT.NY)THEN
            EXIT
         ELSE
            WRITE(*,202)NY
         ENDIF
      ENDDO
202   FORMAT(' Error: the maximum lowpass frequency must be < ',G9.4,/)
C
      CALL BESSEL_TRANSFER_FUNCTIONS(FC,SCALE)
C
      CALL BESSEL_FILTER_COEFFICIENTS(FC,DT,SCALE,A1,A2,B0,B1,B2)
C
      CALL BESSEL_APPLY_FILTER(NUM,A1,A2,B0,B1,B2,Y,YF)
C
      PRINT '(" ")'
      PRINT '( " Enter the filtered time history output filename. ")'
      READ(*,*)B_file
      OPEN(UNIT=11,FILE=B_file,ACTION='WRITE')
C
      DO I=1,NUM
          WRITE(11,200)TT(I),YF(I)
200       FORMAT(2(1X,E12.6))
      ENDDO
C
      WRITE(*,*)' '
      WRITE(*,*)' Filtered Time History Statistics '
      CALL STATISTICS(YF,NUM,MEAN,STD_DEV,RMS)
C
      DEALLOCATE(TT,stat=stat_alloc)
      DEALLOCATE(Y,stat=stat_alloc)
      DEALLOCATE(YF,stat=stat_alloc)
C
      CLOSE(11)
C
      WRITE(*,*)' Output Files: '
      WRITE(*,*)
      WRITE(*,*)' ',B_file,' - time(sec) & amplitude'
      WRITE(*,*)
      WRITE(*,*)' Bessel_transfer_mag.txt   - freq(Hz) & magnitude'
      WRITE(*,*)' Bessel_transfer_phase.txt - freq(Hz) & phase(deg)'
C
100   STOP
      END
C
      SUBROUTINE STATISTICS(Y,NUM,MEAN,STD_DEV,RMS)
      INTEGER NUM
      DOUBLE PRECISION Y(NUM)
      DOUBLE PRECISION MEAN,MS,STD_DEV,RMS
C
      MEAN = SUM(Y)/FLOAT(NUM)
C
      MS=0.
      DO I=1,NUM
          MS=MS+Y(I)**2.
      ENDDO
C
      STD_DEV=SQRT(MS/FLOAT(NUM))
      RMS=SQRT(STD_DEV**2.+MEAN**2.)
C
      WRITE(*,111)MEAN,STD_DEV,RMS
111   FORMAT(/,' Mean= ',G10.4,'  Std Dev= ',G10.4,'  RMS= ',G10.4,/)
C
      RETURN
      END
C
      SUBROUTINE BESSEL_TRANSFER_FUNCTIONS(FC,SCALE)
      PARAMETER(PI=3.141592653589793D0)
      DOUBLE PRECISION F(1000)
      DOUBLE PRECISION FMIN,FMAX,FC,FFF
      DOUBLE PRECISION SCALE
      DOUBLE PRECISION H_MAG,H_PHASE
      DOUBLE COMPLEX S,S2,H,TERM
      INTEGER I,NF
C
      OPEN(UNIT=12,FILE='Bessel_transfer_mag.txt',ACTION='WRITE')
      OPEN(UNIT=13,FILE='Bessel_transfer_phase.txt',ACTION='WRITE')
C
      FMIN=FC/100.
C
      IF(FMIN>1.)FMIN=1.
      FMAX=5*FC
C
      NF=NINT(48.*LOG(FMAX/FMIN)/LOG(2.))
      IF(NF>1000)NF=1000
C
      F(1)=FMIN
      DO I=2,NF
         F(I)=F(I-1)*2.**(1./48.)
      ENDDO
C
      DO I=1,NF
         FFF=SCALE*F(I)/FC
         S=CMPLX(0.,FFF)
         S2=S**2.
         TERM=S2+3.*S+3.
         H=3./TERM
         H_PHASE=(180./PI)*ATAN2(IMAG(H),REAL(H))
         H_MAG=ABS(H)
         WRITE(12,401)F(I),H_MAG
         WRITE(13,401)F(I),H_PHASE
401      FORMAT(2(1X,E12.5))
      ENDDO
C
      CLOSE(12)
      CLOSE(13)
C
      RETURN
      END
C
      SUBROUTINE BESSEL_FILTER_COEFFICIENTS(FC,DT,SCALE,A1,A2,B0,B1,B2)
      PARAMETER(PI=3.141592653589793D0)
      DOUBLE PRECISION FC,DT,SCALE
      DOUBLE PRECISION OM,OM2,DEN
      DOUBLE PRECISION A1,A2
      DOUBLE PRECISION B0,B1,B2
C
      OM=TAN(PI*FC*DT/SCALE)
      OM2=OM**2.
C
      DEN=1.+3.*OM+3.*OM2
C
      B0=3.*OM2/DEN
      B1=2.*B0
      B2=B0
C
      A1=2.*(-1.+3.*OM2)/DEN
      A2=(1.-3.*OM+3.*OM2)/DEN
C
      RETURN
      END
C
      SUBROUTINE BESSEL_APPLY_FILTER(NUM,A1,A2,B0,B1,B2,Y,YF)
      INTEGER I,NUM
      DOUBLE PRECISION A1,A2,B0,B1,B2
      DOUBLE PRECISION Y(NUM),YF(NUM)
C
      WRITE(*,*)' '
      WRITE(*,*)' Filter Coefficients '
      WRITE(*,601)A1,A2
      WRITE(*,602)B0,B1,B2
601   FORMAT('  A1= ',G11.5,' A2= ',G11.5)
602   FORMAT('  B0= ',G11.5,' B1= ',G11.5,' B2= ',G11.5)
C
C  Digital Recursive Filtering Relationship
C
      YF(1)=B0*Y(1)
      YF(2)=B0*Y(2)+B1*Y(1)-A1*YF(1)
C
      DO I=1,NUM
C
         YF(i)= B0*Y(I) + B1*Y(I-1) +B2*Y(I-2)
     $                 -A1*YF(I-1) -A2*YF(I-2)
	  ENDDO
C
      RETURN
      END