Skip to content

Commit

Permalink
Resolved third round of requested changes
Browse files Browse the repository at this point in the history
  • Loading branch information
dpaulzc committed Nov 13, 2024
1 parent 09009f7 commit b9caf40
Show file tree
Hide file tree
Showing 2 changed files with 7 additions and 11 deletions.
16 changes: 6 additions & 10 deletions src/misc.f90
Original file line number Diff line number Diff line change
Expand Up @@ -1657,39 +1657,35 @@ end subroutine subtitle
subroutine Hilbert_transform(fx, Hfx)
!! Does Hilbert tranform for a given function
!! Ref - EQ (4)3, R. Balito et. al.
!! "An algorithm for fast Hilbert transform of real
!! functions"
!! "An algorithm for fast Hilbert transform of real functions"

real(r64), intent(in) :: fx(:)
real(r64), allocatable, intent(out) :: Hfx(:)

! Local variables
integer :: n, k, nfx
real(r64) :: term2, term3, b

nfx = size(fx)
allocate(Hfx(nfx))

! Hilbert function is zero at the edges
Hfx(1) = 0.0_r64
Hfx(nfx) = 0.0_r64
! Note: In the reference, we have N + 1 points and indexing is 0-based
! whereas here we have N points and indexing is 1-based
do k = 1, nfx - 2 ! Run over the internal points
term2 = 0.0_r64 ! 2nd term in Bilato Eq. 4
term3 = 0.0_r64 ! 3rd term in Bilato Eq. 4
do k = 1, nfx - 2 ! Run over the internal points
term2 = 0.0_r64 ! 2nd term in Bilato Eq. 4
term3 = 0.0_r64 ! 3rd term in Bilato Eq. 4

do n = 1, nfx - 2 - k ! Partial sum over internal points
do n = 1, nfx - 2 - k ! Partial sum over internal points
b = log((n + 1.0_r64)/n)
term2 = term2 - (1.0_r64 - (n + 1.0_r64)*b)*fx(k + n + 1) + &
(1.0_r64 - n*b)*fx(k + n + 2)
end do
do n = 1, k - 1 ! Partial sum over internal points
do n = 1, k - 1 ! Partial sum over internal points
b = log((n + 1.0_r64)/n)
term3 = term3 + (1.0_r64 - (n + 1.0_r64)*b)*fx(k - n + 1) - &
(1.0_r64 - n*b)*fx(k - n)
end do

Hfx(k + 1) = -(fx(k + 2) - fx(k) + term2 + term3)/pi
end do
end subroutine Hilbert_transform
Expand Down
2 changes: 1 addition & 1 deletion test/test_misc.f90
Original file line number Diff line number Diff line change
Expand Up @@ -409,6 +409,6 @@ end function fx2
! Hfx2 is actual hilbert transform of fx2, is an even function
pure elemental real(r64) function hfx2(x)
real(r64), intent(in) :: x
hfx2 = (1/exp(1.0_r64) - cos(x))/(1.0_r64 + x**2)
hfx2 = (exp(-1.0_r64) - cos(x))/(1.0_r64 + x**2)
end function hfx2
end program test_misc

0 comments on commit b9caf40

Please sign in to comment.