Skip to content

Commit

Permalink
Resolved requested changes
Browse files Browse the repository at this point in the history
  • Loading branch information
@dpaulzc
dpaulzc committed Nov 13, 2024
1 parent 4d3d4ed commit eda1e54
Show file tree
Hide file tree
Showing 2 changed files with 40 additions and 49 deletions.
16 changes: 6 additions & 10 deletions src/misc.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
73 changes: 34 additions & 39 deletions test/test_misc.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)
Expand All @@ -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

0 comments on commit eda1e54

Please sign in to comment.