From b9caf4091dd375cca32651a01ae0e877fc92825f Mon Sep 17 00:00:00 2001 From: Dwaipayan Paul Date: Wed, 13 Nov 2024 16:12:37 +0100 Subject: [PATCH] Resolved third round of requested changes --- src/misc.f90 | 16 ++++++---------- test/test_misc.f90 | 2 +- 2 files changed, 7 insertions(+), 11 deletions(-) diff --git a/src/misc.f90 b/src/misc.f90 index 8957ba0..efc6fd8 100644 --- a/src/misc.f90 +++ b/src/misc.f90 @@ -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 diff --git a/test/test_misc.f90 b/test/test_misc.f90 index cd389c1..ac254e0 100644 --- a/test/test_misc.f90 +++ b/test/test_misc.f90 @@ -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