From eda1e54d5cde18d7e1f59309add982434f7c3ce3 Mon Sep 17 00:00:00 2001 From: Dwaipayan Paul Date: Wed, 13 Nov 2024 15:01:43 +0100 Subject: [PATCH] Resolved requested changes --- src/misc.f90 | 16 ++++------ test/test_misc.f90 | 73 +++++++++++++++++++++------------------------- 2 files changed, 40 insertions(+), 49 deletions(-) diff --git a/src/misc.f90 b/src/misc.f90 index 563a527..d6cab81 100644 --- a/src/misc.f90 +++ b/src/misc.f90 @@ -1657,15 +1657,12 @@ 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 - !! - !! fx is the function - !! Hfx is the Hilbert transform of the function + !! "An algorithm for fast Hilbert transform of real + !! functions" !! - real(r64), allocatable, intent(out) :: Hfx(:) real(r64), intent(in) :: fx(:) + real(r64), allocatable, intent(out) :: Hfx(:) ! Local variables integer :: n, k, nfx @@ -1677,8 +1674,8 @@ subroutine Hilbert_transform(fx, Hfx) ! Hilbert function is zero at the edges Hfx(1) = 0.0_r64 Hfx(nfx) = 0.0_r64 - - ! Both Hfx and fx follow 1-indexing system unlike the paper + ! 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 @@ -1688,14 +1685,13 @@ subroutine Hilbert_transform(fx, Hfx) 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 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) = (-1.0_r64/pi)*(fx(k + 2) - fx(k) + term2 + term3) + Hfx(k + 1) = -(fx(k + 2) - fx(k) + term2 + term3)/pi end do end subroutine Hilbert_transform end module misc diff --git a/test/test_misc.f90 b/test/test_misc.f90 index e77acfb..788c168 100644 --- a/test/test_misc.f90 +++ b/test/test_misc.f90 @@ -339,8 +339,8 @@ program test_misc ! Hilbert transform tests (H) ! fx1 -> function 1, fx2 -> function 2 - ! - hfx1_even stores hilbert transform calculated for fx1, and for even number - ! - of points + ! hfx1_even stores hilbert transform calculated for fx1, and for even number + ! of points xmin = -30.0 xmax = 30.0 n_even = 4000 @@ -352,33 +352,33 @@ program test_misc ind_odd = [889, 1333, 1777, 2221, 2665] itest = itest + 1 - test_array(itest) = testify("Hilbert transform: even function, even points") + test_array(itest) = testify("Hilbert transform: f(x) = 1/(1 + x^2), even points") allocate(x_even(n_even), hfx1_even(n_even)) call linspace(x_even, xmin, xmax, n_even) - call Hilbert_transform(fx1_array(x_even), hfx1_even) - call test_array(itest)%assert(hfx1_even(ind_even), hfx1_array(x_even(ind_even)), & + call Hilbert_transform(fx1_el(x_even), hfx1_even) + call test_array(itest)%assert(hfx1_even(ind_even), hfx1_el(x_even(ind_even)), & tol = 2e-4_r64) itest = itest + 1 - test_array(itest) = testify("Hilbert transform: odd function, even points") + test_array(itest) = testify("Hilbert transform: f(x) = sin(x)/(1 + x^2), even points") allocate(hfx2_even(n_even)) - call Hilbert_transform(fx2_array(x_even), hfx2_even) - call test_array(itest)%assert(hfx2_even(ind_even), hfx2_array(x_even(ind_even)), & + call Hilbert_transform(fx2_el(x_even), hfx2_even) + call test_array(itest)%assert(hfx2_even(ind_even), hfx2_el(x_even(ind_even)), & tol = 1e-4_r64) itest = itest + 1 - test_array(itest) = testify("Hilbert transform: even function, odd points") + test_array(itest) = testify("Hilbert transform: f(x) = 1/(1 + x^2), odd points") allocate(x_odd(n_odd), hfx1_odd(n_odd)) call linspace(x_odd, xmin, xmax, n_odd) - call Hilbert_transform(fx1_array(x_odd), hfx1_odd) - call test_array(itest)%assert(hfx1_odd(ind_odd), hfx1_array(x_odd(ind_odd)), & + call Hilbert_transform(fx1_el(x_odd), hfx1_odd) + call test_array(itest)%assert(hfx1_odd(ind_odd), hfx1_el(x_odd(ind_odd)), & tol = 4e-4_r64) itest = itest + 1 - test_array(itest) = testify("Hilbert transform: odd function, odd points") + test_array(itest) = testify("Hilbert transform: f(x) = sin(x)/(1 + x^2), odd points") allocate(hfx2_odd(n_odd)) - call Hilbert_transform(fx2_array(x_odd), hfx2_odd) - call test_array(itest)%assert(hfx2_odd(ind_odd), hfx2_array(x_odd(ind_odd)), & + call Hilbert_transform(fx2_el(x_odd), hfx2_odd) + call test_array(itest)%assert(hfx2_odd(ind_odd), hfx2_el(x_odd(ind_odd)), & tol = 1e-5_r64) tests_all = testify(test_array) @@ -389,35 +389,30 @@ program test_misc contains ! reference functions for the Hilbert transform test ! fx1 is an even function - function fx1_array(x) result(fx1) - real(r64), intent(in) :: x(:) - real(r64), allocatable :: fx1(:) - allocate(fx1(size(x))) - fx1 = 1/(1 + x**2) - end function fx1_array + elemental function fx1_el(x) result(fx1) + real(r64), intent(in) :: x + real(r64) :: fx1 + fx1 = 1/(1.0_r64 + x**2) + end function fx1_el ! Hfx1 is actual hilbert transform of fx1, is an odd function - function hfx1_array(x) result(hfx1) - real(r64), intent(in) :: x(:) - real(r64), allocatable :: hfx1(:) - allocate(hfx1(size(x))) - hfx1 = x/(1 + x**2) - end function hfx1_array + elemental function hfx1_el(x) result(hfx1) + real(r64), intent(in) :: x + real(r64) :: hfx1 + hfx1 = x/(1.0_r64 + x**2) + end function hfx1_el ! fx2 is an odd function - function fx2_array(x) result(fx2) - real(r64), intent(in) :: x(:) - real(r64), allocatable :: fx2(:) - allocate(fx2(size(x))) - fx2 = sin(x)/(1 + x**2) - end function fx2_array + elemental function fx2_el(x) result(fx2) + real(r64), intent(in) :: x + real(r64) :: fx2 + fx2 = sin(x)/(1.0_r64 + x**2) + end function fx2_el ! Hfx2 is actual hilbert transform of fx2, is an even function - function hfx2_array(x) result(hfx2) - real(r64), intent(in) :: x(:) - real(r64), allocatable :: hfx2(:) - real(r64), parameter :: e = 2.718281 - allocate(hfx2(size(x))) - hfx2 = (1/e - cos(x))/(1 + x**2) - end function hfx2_array + elemental function hfx2_el(x) result(hfx2) + real(r64), intent(in) :: x + real(r64) :: hfx2 + hfx2 = (1/exp(1.0_r64) - cos(x))/(1.0_r64 + x**2) + end function hfx2_el end program test_misc