add hilbert (transform) function.

doc/specs/fftpack.md

@@ -1170,3 +1170,47 @@ program demo_ifftshift
     print *, ifftshift(fftshift(x))        !! [1.0, 2.0, 3.0, 4.0, 5.0]
 end program demo_ifftshift
 ```
+
+### `hilbert`
+
+#### Description
+
+Computes the discrete-time analytic signal using Hilbert.
+
+#### Status
+
+Experimental.
+
+#### Class
+
+Pure function.
+
+#### Syntax
+
+`result = [[fftpack(module):hilbert(interface)]](x [, n])`
+
+#### Arguments
+
+`x`: Shall be a `complex` and rank-1 array.
+This argument is `intent(in)`. 
+
+`n`: Shall be an `integer` scalar.
+This argument is `intent(in)` and `optional`.  
+Defines the length of the Hilbert transform. If `n` is not specified (the default) then `n = size(x)`. If `n <= size(x)`, `x` is truncated, if `n > size(x)`, `x` is zero-padded.
+
+#### Return value
+
+Returns the `complex` and rank-1 Hilbert transform.
+
+#### Example
+
+```fortran
+program demo_fftpack_hilbert
+
+    use fftpack, only: hilbert, dp
+    complex(kind=dp) :: x(4) = [1, 2, 3, 4]
+
+    print *, hilbert(x)     !! [(1.000,1.000), (2.000,-1.000), (3.000,-1.000), (4.000,1.000)]
+
+end program demo_fftpack_hilbert
+```

example/demo_fftpack_hilbert.f90

@@ -0,0 +1,11 @@
+program demo_fftpack_hilbert
+
+    use fftpack, only: hilbert, dp
+
+    print *, hilbert([complex(kind=dp) :: 1, 2, 3, 4   ])
+    print *, hilbert([complex(kind=dp) :: 1, 2, 3, 4, 5])
+
+    !! [(1.000,1.000), (2.000,-1.000),  (3.000,-1.000),  (4.000,1.000)]
+    !! [(1.000,1.701), (2.000,-1.376), (3.000,-0.6498), (4.000,-1.376), (5.000,1.701)]
+
+end program demo_fftpack_hilbert

fpm.toml

@@ -15,6 +15,12 @@ auto-executables = false
 auto-tests = false
 auto-examples = false
 
+# Examples
+[[example]]
+name       = "fftpack_hilbert"
+source-dir = "example"
+main       = "demo_fftpack_hilbert.f90"
+
 # Original test
 [[test]]
 name = "tstfft"
@@ -82,3 +88,8 @@ main = "test_fftpack_fftshift.f90"
 name = "fftpack_ifftshift"
 source-dir = "test"
 main = "test_fftpack_ifftshift.f90"
+
+[[test]]
+name       = "fftpack_hilbert"
+source-dir = "test"
+main       = "test_fftpack_hilbert.f90"

src/Makefile

@@ -58,7 +58,8 @@ SRCF90 = \
 	fftpack_fftshift.f90\
 	fftpack_ifftshift.f90\
 	fftpack_qct.f90\
-	fftpack_iqct.f90
+	fftpack_iqct.f90\
+	fftpack_hilbert.f90
 
 OBJF := $(SRCF:.f=.o)
 OBJF90 := $(SRCF90:.f90=.o)

src/fftpack.f90

@@ -5,6 +5,7 @@ module fftpack
     integer, parameter :: dp = kind(1.0d0)
 
     public :: dp
+
     public :: zffti, zfftf, zfftb
     public :: fft, ifft
     public :: fftshift, ifftshift
@@ -17,6 +18,8 @@ module fftpack
     public :: dcosqi, dcosqf, dcosqb
     public :: qct, iqct
 
+    public :: hilbert
+
     interface
 
         !> Version: experimental
@@ -255,4 +258,16 @@ pure module function ifftshift_rdp(x) result(result)
         end function ifftshift_rdp
     end interface ifftshift
 
+    !> Version: experimental
+    !> 
+    !> Computes the discrete-time analytic signal using Hilbert.
+    !> ([Specification](../page/specs/fftpack.html#hilbert))
+    interface hilbert
+        pure module function hilbert_dp(x, n) result(result)
+            complex(kind=dp), intent(in)  :: x(:)
+            integer, intent(in), optional :: n
+            complex(kind=dp), allocatable :: result(:)
+        end function hilbert_dp
+    end interface hilbert
+
 end module fftpack

src/fftpack_hilbert.f90

@@ -0,0 +1,59 @@
+!> Reference: https://ww2.mathworks.cn/help/signal/ref/hilbert.html
+submodule(fftpack) fftpack_hilbert
+
+    use iso_fortran_env, only: int8
+
+contains
+
+    !> Computes the discrete-time analytic signal using Hilbert.
+    pure module function hilbert_dp(x, n) result(result)
+        complex(kind=dp), intent(in)  :: x(:)
+        integer, intent(in), optional :: n
+        complex(kind=dp), allocatable :: result(:)
+
+        integer                         :: n_, lensav, i
+        real(kind=dp), allocatable      :: wsave(:)
+        integer(kind=int8), allocatable :: H(:)
+
+        if (present(n)) then
+            n_ = n
+            if (n_ <= size(x)) then
+                result = x(:n_)
+            else if (n_ > size(x)) then
+                result = [x, ((0.0_dp, 0.0_dp), i=1, n_ - size(x))]
+            end if
+        else
+            n_ = size(x)
+            result = x
+        end if
+
+        allocate(H(n_))
+
+        !> Create a vector H whose elements H(i) have the values:
+        !> 1 for i = 1
+        !> 2 for i = 2, 3, ... , (n/2)
+        !> 1(2) for i = n/2 + 1, mod(n, 2) ==(/=) 0
+        !> 0 for i = (n/2)+2, ... , n
+        H(1)           = 1_int8
+        H(2:n_/2)      = 2_int8
+        H(n_/2 + 1)    = merge(1_int8, 2_int8, mod(n_, 2) == 0)
+        H(n_/2 + 2:n_) = 0_int8
+
+        !> Initialize FFT
+        lensav = 4*n_ + 15
+        allocate(wsave(lensav))
+        call zffti(n_, wsave)
+
+        !> Forward transformation
+        call zfftf(n_, result, wsave)
+
+        result = result*H
+
+        !> Backward transformation
+        call zfftb(n_, result, wsave)
+
+        result = result/n_
+
+    end function hilbert_dp
+
+end submodule fftpack_hilbert

test/Makefile

@@ -8,7 +8,8 @@ all: tstfft \
 	fftpack_dzfft \
 	fftpack_dcosq \
 	fftpack_qct \
-	fftpack_iqct
+	fftpack_iqct \
+	fftpack_hilbert
 
 # Orginal test
 tstfft: tstfft.o
@@ -57,5 +58,9 @@ fftpack_ifftshift: test_fftpack_ifftshift.f90
 	$(FC) $(FFLAGS) $< -L../src -l$(LIB) -I../src -o [email protected]
 	./fftpack_ifftshift.x
 
+fftpack_hilbert: test_fftpack_hilbert.f90
+	$(FC) $(FFLAGS) $< -L../src -l$(LIB) -I../src -o [email protected]
+	./fftpack_hilbert.x
+
 clean:
 	rm -f -r *.o *.x

test/test_fftpack_hilbert.f90

@@ -0,0 +1,31 @@
+program tester
+
+    call test_fftpack_hilbert
+    print *, "All tests in `test_fftpack_hilbert` passed."
+
+contains
+
+    subroutine check(condition, msg)
+        logical, intent(in)          :: condition
+        character(len=*), intent(in) :: msg
+        if (.not. condition) error stop msg
+    end subroutine check
+
+    subroutine test_fftpack_hilbert()
+
+        use fftpack, only: hilbert, dp
+        real(kind=dp) :: eps = 1.0e-10_dp
+
+        call check(sum(abs(hilbert([complex(kind=dp) :: 1, 2, 3, 4   ]) - &
+                                   [complex(kind=dp) :: (1.0_dp, 1.0_dp),  (2.0_dp, -1.0_dp), &
+                                                        (3.0_dp, -1.0_dp), (4.0_dp, 1.0_dp)])) < eps, &
+                   msg="hilbert([complex(kind=dp) :: 1, 2, 3, 4]) failed.")
+        call check(sum(abs(hilbert([complex(kind=dp) :: 1, 2, 3, 4, 5]) - &
+                                   [complex(kind=dp) :: (1.0_dp, 1.7013016167040795_dp),   (2.0_dp, -1.3763819204711736_dp), &
+                                                        (3.0_dp, -0.64983939246581257_dp), (4.0_dp, -1.3763819204711731_dp), &
+                                                        (5.0_dp, 1.7013016167040800_dp)])) < eps, &
+                   msg="hilbert([complex(kind=dp) :: 1, 2, 3, 4, 5]) failed.")
+
+    end subroutine test_fftpack_hilbert
+
+end program tester