Skip to content

Commit

Permalink
Resolved second 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 eda1e54 commit 09009f7
Show file tree
Hide file tree
Showing 2 changed files with 16 additions and 21 deletions.
1 change: 0 additions & 1 deletion src/misc.f90
Original file line number Diff line number Diff line change
Expand Up @@ -1659,7 +1659,6 @@ subroutine Hilbert_transform(fx, Hfx)
!! Ref - EQ (4)3, R. Balito et. al.
!! "An algorithm for fast Hilbert transform of real
!! functions"
!!

real(r64), intent(in) :: fx(:)
real(r64), allocatable, intent(out) :: Hfx(:)
Expand Down
36 changes: 16 additions & 20 deletions test/test_misc.f90
Original file line number Diff line number Diff line change
Expand Up @@ -355,30 +355,30 @@ program test_misc
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_el(x_even), hfx1_even)
call test_array(itest)%assert(hfx1_even(ind_even), hfx1_el(x_even(ind_even)), &
call Hilbert_transform(fx1(x_even), hfx1_even)
call test_array(itest)%assert(hfx1_even(ind_even), hfx1(x_even(ind_even)), &
tol = 2e-4_r64)

itest = itest + 1
test_array(itest) = testify("Hilbert transform: f(x) = sin(x)/(1 + x^2), even points")
allocate(hfx2_even(n_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)), &
call Hilbert_transform(fx2(x_even), hfx2_even)
call test_array(itest)%assert(hfx2_even(ind_even), hfx2(x_even(ind_even)), &
tol = 1e-4_r64)

itest = itest + 1
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_el(x_odd), hfx1_odd)
call test_array(itest)%assert(hfx1_odd(ind_odd), hfx1_el(x_odd(ind_odd)), &
call Hilbert_transform(fx1(x_odd), hfx1_odd)
call test_array(itest)%assert(hfx1_odd(ind_odd), hfx1(x_odd(ind_odd)), &
tol = 4e-4_r64)

itest = itest + 1
test_array(itest) = testify("Hilbert transform: f(x) = sin(x)/(1 + x^2), odd points")
allocate(hfx2_odd(n_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)), &
call Hilbert_transform(fx2(x_odd), hfx2_odd)
call test_array(itest)%assert(hfx2_odd(ind_odd), hfx2(x_odd(ind_odd)), &
tol = 1e-5_r64)

tests_all = testify(test_array)
Expand All @@ -389,30 +389,26 @@ program test_misc
contains
! reference functions for the Hilbert transform test
! fx1 is an even function
elemental function fx1_el(x) result(fx1)
pure elemental real(r64) function fx1(x)
real(r64), intent(in) :: x
real(r64) :: fx1
fx1 = 1/(1.0_r64 + x**2)
end function fx1_el
end function fx1

! Hfx1 is actual hilbert transform of fx1, is an odd function
elemental function hfx1_el(x) result(hfx1)
pure elemental real(r64) function hfx1(x)
real(r64), intent(in) :: x
real(r64) :: hfx1
hfx1 = x/(1.0_r64 + x**2)
end function hfx1_el
end function hfx1

! fx2 is an odd function
elemental function fx2_el(x) result(fx2)
pure elemental real(r64) function fx2(x)
real(r64), intent(in) :: x
real(r64) :: fx2
fx2 = sin(x)/(1.0_r64 + x**2)
end function fx2_el
end function fx2

! Hfx2 is actual hilbert transform of fx2, is an even function
elemental function hfx2_el(x) result(hfx2)
pure elemental real(r64) function hfx2(x)
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 function hfx2
end program test_misc

0 comments on commit 09009f7

Please sign in to comment.