From 673da9a755b70f929869f182eaadc4d902a8470f Mon Sep 17 00:00:00 2001
From: Dwaipayan Paul <dwaipayanpaul6@gmail.com>
Date: Tue, 12 Nov 2024 15:36:26 +0100
Subject: [PATCH] Added unit test for hilbert transform

---
 src/misc.f90       |  2 +-
 test/test_misc.f90 | 86 ++++++++++++++++++++++++++++++++++++++++++++--
 2 files changed, 85 insertions(+), 3 deletions(-)

diff --git a/src/misc.f90 b/src/misc.f90
index 2c85447..156d22e 100644
--- a/src/misc.f90
+++ b/src/misc.f90
@@ -1658,7 +1658,7 @@ 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
-     !!   function
+     !!   functions
      !!
      !! fx is the function 
      !! Hfx is the Hilbert transform of the function
diff --git a/test/test_misc.f90 b/test/test_misc.f90
index 5cab0e4..e363557 100644
--- a/test/test_misc.f90
+++ b/test/test_misc.f90
@@ -7,12 +7,12 @@ program test_misc
        twonorm, binsearch, mux_vector, demux_vector, interpolate, coarse_grained, &
        unique, linspace, compsimps, mux_state, demux_state, demux_mesh, expm1, &
        Fermi, Bose, Pade_continued, precompute_interpolation_corners_and_weights, &
-       interpolate_using_precomputed, operator(.umklapp.), shrink
+       interpolate_using_precomputed, operator(.umklapp.), shrink, hilbert_transform
   
   implicit none
 
   integer :: itest
-  integer, parameter :: num_tests = 28
+  integer, parameter :: num_tests = 32
   type(testify) :: test_array(num_tests), tests_all
   integer(i64) :: index, quotient, remainder, int_array(5), v1(3), v2(3), &
        v1_muxed, v2_muxed, ik, ik1, ik2, ik3, ib1, ib2, ib3, wvmesh(3), &
@@ -23,6 +23,9 @@ program test_misc
        real_array(5), result, q1(3, 4), q2(3, 4), q3(3, 4)
   real(r64), allocatable :: integrand(:), domain(:), im_axis(:), real_func(:), &
        widc(:, :), f_coarse(:), f_interp(:), array_of_reals(:)
+  real(r64), allocatable :: hfx1_even(:), hfx1_odd(:), hfx2_even(:), hfx2_odd(:), &
+       ind_even(:), ind_odd(:), x_even(:), x_odd(:), xmin, xmax
+  integer(i64) :: n_even, n_odd  
 
   print*, '<<module misc unit tests>>'
   
@@ -333,9 +336,88 @@ program test_misc
   array_of_reals = [1, 2, 3, 4, 5]*1.0_r64
   call shrink(array_of_reals, 2_i64)
   call test_array(itest)%assert(array_of_reals, [1, 2]*1.0_r64)
+
+  ! Hilbert transform tests (H)
+  ! fx1 -> function 1, fx2 -> function 2
+  ! - hfx1_even stores hilbert transform calculated for fx1, and for even number
+  ! - of points 
+  xmin = -30.0
+  xmax = 30.0
+  n_even = 4000
+  n_odd = 4001
+  ! ind_even are indices to compare in case of even number of points
+  ! ind_odd are indices to compare in case of odd number of points 
+  allocate(ind_even(6),ind_odd(5))
+  ind_even = [801, 1201, 1601, 2001, 2401, 2801]
+  ind_odd = [889, 1333, 1777, 2221, 2665]
+   
+  itest = itest + 1
+  test_array(itest) = testify("Hilbert transform: even function, 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)), &
+       tol = 2e-4_r64)
+
+  itest = itest + 1
+  test_array(itest) = testify("Hilbert transform: odd function, 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)), &
+       tol = 1e-4_r64)
+  
+  itest = itest + 1
+  test_array(itest) = testify("Hilbert transform: even function, 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)), &
+       tol = 4e-4_r64)
+  
+  itest = itest + 1
+  test_array(itest) = testify("Hilbert transform: odd function, 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)), &
+       tol = 1e-5_r64)
   
   tests_all = testify(test_array)
   call tests_all%report
   
   if(tests_all%get_status() .eqv. .false.) error stop -1
+
+  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
+     
+     ! 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
+     
+     ! 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
+     
+     ! 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
 end program test_misc