PROGRAM AUTOCORRELATION C C AUTOCORRELATION.F ver 1.0 November 24, 2012 C C By Tom Irvine C Email: tom@vibrationdata.com C C Autocorrelation of a time history C C Variables: C C T = input time C Y = input amplitude C C******************************************************** C23456789 C DOUBLE PRECISION, ALLOCATABLE :: Y(:) DOUBLE PRECISION, ALLOCATABLE :: XC(:,:) DOUBLE PRECISION AMP,DT,DF,DUR DOUBLE PRECISION MS,SR DOUBLE PRECISION TT,TT1 DOUBLE PRECISION MAXA,MAXT C INTEGER I,IJK,J,K INTEGER N,N2,NUM,M INTEGER IFLAG INTEGER ios INTEGER stat_alloc C CHARACTER :: A_file*20 C PARAMETER(PI=3.141592653589793D0,TPI=2D0*PI) C 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 OPEN(UNIT=11,FILE='autocorrelation.out',ACTION='WRITE') C I=1 READ(10,*,IOSTAT=IOstatus) TT1,AMP I=2 IFLAG=0 DO WHILE(IFLAG.EQ.0) READ(10,*,IOSTAT=IOstatus) TT,AMP IF (IOstatus.EQ.0)THEN NUM=I I=I+1 ELSE NUM=I-1 IFLAG=1 ENDIF ENDDO C DUR=TT-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 NH=NINT(NUM/2.) C REWIND(10) C DF=1/(NUM*DT) C ALLOCATE(Y(NUM)) C DO I=1,NUM READ(10,*,IOSTAT=IOstatus) TT,Y(I) ENDDO C CLOSE(10) C C**************************************************************** C C23456789 C N=NUM M=N N2=2*N C ALLOCATE(XC(N2-1,2)) C DO I=1,N2-1 XC(I,1)=-DUR+(I-1)*DT XC(I,2)=0. ENDDO C DO I=1,N C IJK=I+N-1 C DO J=1,M C IF((I+J-1).GT.N)THEN EXIT ELSE XC(IJK,2)=XC(IJK,2)+Y(J)*Y(I+J-1) ENDIF C ENDDO C ENDDO C C23456789 DO I=1,N-1 C IJK=N-I C DO J=1,M C IF((I+J)>M)THEN EXIT ELSE XC(IJK,2)=XC(IJK,2)+Y(J)*Y(I+J) ENDIF C ENDDO C ENDDO C MAXA=0. MAXT=0. C DO I=1,N2-1 XC(I,2)=XC(I,2)/FLOAT(N-1) WRITE(11,121)XC(I,1),XC(I,2) 121 FORMAT(1X,E12.6,2X,E12.6) C IF(XC(I,2).GT.MAXA)THEN MAXA=XC(I,2) MAXT=XC(I,1) ENDIF C ENDDO C CLOSE(11) C DEALLOCATE(XC,stat=stat_alloc) DEALLOCATE(Y,stat=stat_alloc) C WRITE(*,333)MAXA,MAXT 333 FORMAT(/,' Max Autocorrelation = ',G11.5,' at ',G11.5,' sec') C PRINT'(" ")' PRINT'(" Calculation complete.")' PRINT'(" ")' PRINT'(" The output file has two columns: ")' PRINT'(" Delay(sec) & Autocorrleation ")' PRINT'(" ")' PRINT'(" The output file is: autocorrelation.out")' PRINT'(" ")' C C**************************************************************** 100 STOP END