program konvolucija
      implicit none
      integer(4),parameter :: nmax=4096
      integer(4) :: i
      real(8)    :: data(nmax),filter(nmax),resitev(2*nmax),r

      open(1,file='funkcija.dat')
      open(2,file='filter.dat')
      do i=1,nmax
         read(1,*)r,data(i)
         read(2,*)r,filter(i)
      enddo
      close(1)
      close(2)

      call convlv(data,nmax,filter,nmax,1,resitev)

      open(3,file='konvol.dat')
      do i=1,2*nmax
         write(3,*)i,resitev(i) !resitev(2*i-1)
      end do

      close(3)
    
      end program konvolucija


      SUBROUTINE convlv(data,n,respns,m,isign,ans)
      INTEGER(4),PARAMETER :: NMAX=4096
      INTEGER(4) :: isign,m,n
      REAL(8)    :: data(n),respns(n)
      COMPLEX(8) :: ans(n)
!CU    USES realft,twofft

!Convolves or deconvolves a real data set data(1:n) (including any user-supplied zero
!padding) with a response function respns, stored in wrap-around order in a real array of
!length m ≤ n. (m should be an odd integer.) Wrap-around order means that the first half
!of the array respns contains the impulse response function at positive times, while the
!second half of the array contains the impulse response function at negative times, counting
!down from the highest element respns(m). On input isign is +1 for convolution, −1
!for deconvolution. The answer is returned in the first n components of ans. However, ans
!must be supplied in the calling program with length at least 2*n, for consistency with
!twofft. n MUST be an integer power of two.

      INTEGER(4) :: i,no2
      COMPLEX(8) :: fft(NMAX)
      do 11 i=1,(m-1)/2
        respns(n+1-i)=respns(m+1-i)
11    continue
      do 12 i=(m+3)/2,n-(m-1)/2
        respns(i)=0.0
12    continue
      call twofft(data,respns,fft,ans,n)
      no2=n/2
      do 13 i=1,no2+1
        if (isign.eq.1) then
          ans(i)=fft(i)*ans(i)/no2
        else if (isign.eq.-1) then
          if (abs(ans(i)).eq.0.0) write(*,*)'deconvolving at response zero in convlv'
          ans(i)=fft(i)/ans(i)/no2
        else
          pause 'no meaning for isign in convlv'
        endif
13    continue
      ans(1)=cmplx(real(ans(1)),real(ans(no2+1)))
      call realft(ans,n,-1)
      return
      END


      SUBROUTINE twofft(data1,data2,fft1,fft2,n)
      INTEGER(4) :: n
      REAL(8)    :: data1(n),data2(n)
      COMPLEX(8) :: fft1(n),fft2(n)
!CU    USES four1
      INTEGER(4) :: j,n2
      COMPLEX(8) :: h1,h2,c1,c2
      c1=cmplx(0.5,0.0)
      c2=cmplx(0.0,-0.5)
      do 11 j=1,n
        fft1(j)=cmplx(data1(j),data2(j))
11    continue
      call four1(fft1,n,1)
      fft2(1)=cmplx(aimag(fft1(1)),0.0)
      fft1(1)=cmplx(real(fft1(1)),0.0)
      n2=n+2
      do 12 j=2,n/2+1
        h1=c1*(fft1(j)+conjg(fft1(n2-j)))
        h2=c2*(fft1(j)-conjg(fft1(n2-j)))
        fft1(j)=h1
        fft1(n2-j)=conjg(h1)
        fft2(j)=h2
        fft2(n2-j)=conjg(h2)
12    continue
      return
      END

      SUBROUTINE FOUR1(DATA,NN,ISIGN)
      REAL(8) :: WR,WI,WPR,WPI,WTEMP,THETA
      real(8) :: DATA(NN)
!      real(8) :: DATA(4096)
!      DIMENSION DATA(*)
      N=2*NN
      J=1
      DO 11 I=1,N,2
        IF(J.GT.I)THEN
          TEMPR=DATA(J)
          TEMPI=DATA(J+1)
          DATA(J)=DATA(I)
          DATA(J+1)=DATA(I+1)
          DATA(I)=TEMPR
          DATA(I+1)=TEMPI
        ENDIF
        M=N/2
1       IF ((M.GE.2).AND.(J.GT.M)) THEN
          J=J-M
          M=M/2
        GO TO 1
        ENDIF
        J=J+M
11    CONTINUE
      MMAX=2
2     IF (N.GT.MMAX) THEN
        ISTEP=2*MMAX
        THETA=6.28318530717959D0/(ISIGN*MMAX)
        WPR=-2.D0*DSIN(0.5D0*THETA)**2
        WPI=DSIN(THETA)
        WR=1.D0
        WI=0.D0
        DO 13 M=1,MMAX,2
          DO 12 I=M,N,ISTEP
            J=I+MMAX
            TEMPR=SNGL(WR)*DATA(J)-SNGL(WI)*DATA(J+1)
            TEMPI=SNGL(WR)*DATA(J+1)+SNGL(WI)*DATA(J)
            DATA(J)=DATA(I)-TEMPR
            DATA(J+1)=DATA(I+1)-TEMPI
            DATA(I)=DATA(I)+TEMPR
            DATA(I+1)=DATA(I+1)+TEMPI
12        CONTINUE
          WTEMP=WR
          WR=WR*WPR-WI*WPI+WR
          WI=WI*WPR+WTEMP*WPI+WI
13      CONTINUE
        MMAX=ISTEP
      GO TO 2
      ENDIF
      RETURN
      END


      SUBROUTINE realft(data,n,isign)
!      INTEGER(4),parameter :: n=4096
      INTEGER(4) :: isign,n
      REAL(8)    :: data(n)
! CU    USES four1
      INTEGER(4) :: i,i1,i2,i3,i4,n2p3
      REAL(8)    :: c1,c2,h1i,h1r,h2i,h2r,wis,wrs
      DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp
      theta=3.141592653589793d0/dble(n/2)
      c1=0.5
      if (isign.eq.1) then
        c2=-0.5
        call four1(data,n/2,+1)
      else
        c2=0.5
        theta=-theta
      endif
      wpr=-2.0d0*sin(0.5d0*theta)**2
      wpi=sin(theta)
      wr=1.0d0+wpr
      wi=wpi
      n2p3=n+3
      do 11 i=2,n/4
        i1=2*i-1
        i2=i1+1
        i3=n2p3-i2
        i4=i3+1
        wrs=sngl(wr)
        wis=sngl(wi)
        h1r=c1*(data(i1)+data(i3))
        h1i=c1*(data(i2)-data(i4))
        h2r=-c2*(data(i2)+data(i4))
        h2i=c2*(data(i1)-data(i3))
        data(i1)=h1r+wrs*h2r-wis*h2i
        data(i2)=h1i+wrs*h2i+wis*h2r
        data(i3)=h1r-wrs*h2r+wis*h2i
        data(i4)=-h1i+wrs*h2i+wis*h2r
        wtemp=wr
        wr=wr*wpr-wi*wpi+wr
        wi=wi*wpr+wtemp*wpi+wi
11    continue
      if (isign.eq.1) then
        h1r=data(1)
        data(1)=h1r+data(2)
        data(2)=h1r-data(2)
      else
        h1r=data(1)
        data(1)=c1*(h1r+data(2))
        data(2)=c1*(h1r-data(2))
        call four1(data,n/2,-1)
      endif
      return
      END