program lowpassf

implicit none

integer,parameter :: num=200000,  filternum=100  ! 2*fiulternum+!
integer :: i 
real :: x(num), xl(num), t(num)
real :: filt(-filternum:filternum)
real :: cutoff, dt
cutoff = 25.0		! Hz
dt = 0.005


 do i=1, num
  read(*,*) t(i), x(i)
 end do


 call setfilter( filt, filternum,  cutoff, dt )

 call convolution( x, xl, num, filt, filternum )

 do i=1, num
  write(*,*) t(i), x(i), xl(i)
 end do

end


subroutine setfilter( filt, filternum,  cutoff, dt )
implicit none

 integer :: filternum
 real :: pi = 3.14159265358979323846
 real :: filt(-filternum:filternum), cutoff, dt
 real :: wc, sum
 integer :: n, i

  wc= 2*pi*cutoff


 do i=1, filternum
	filt(i) = 2.0/(i*dt) * sin( wc*dt*i )/(2*pi)
	filt(-i) = 2.0/(i*dt) * sin( wc*dt*i )/(2*pi)
 end do
   filt(0) = 2*wc/(2*pi)

!    check sum

sum=0
 do i=-filternum, filternum
	sum = sum + filt(i)
 end do

 write(0,*) 'sum=', sum

 do i=-filternum, filternum
   filt(i) = filt(i)/sum
	write(0, *) i, filt(i)
 
 end do

end subroutine 


subroutine  convolution( x, xl, num, filt, filternum )
implicit none
integer :: num, filternum
real :: x(num), xl(num), filt(-filternum:filternum)
real :: xtmp
integer :: i, j

 do i=1, num
	xl(i) = 0
   do j=-filternum, filternum
	if( i+j <= 0 .OR. i+j > num ) then
		xtmp = 0
        else 
                xtmp = x(i+j)
        end if 
       xl(i) = xl(i) + filt(j)*xtmp
   end do
 end do

end 


