Add easy-to-use interface: fft & ifft for complex sequences.

Author: zoziha

doc/specs/fftpack.md

@@ -196,3 +196,98 @@ program demo_zfftb
     call zfftb(4,x,w)   !! `x` returns [(4.0,0.0), (8.0,0.0), (12.0,0.0), (16.0,0.0)].
 end program demo_zfftb
 ```
+
+## `fft`
+
+### Description
+
+Computes the forward complex discrete fourier
+transform (the fourier analysis).   
+
+### Status
+
+Experimental.
+
+### Class
+
+Pure function.
+
+### Snytax
+
+`call [[fftpack(module):fft(interface)]](x [, n])`
+
+### Argument
+
+`x`: Shall be a `complex` array.
+This argument is `intent(in)`. 
+
+`n`: Shall be an `integer` scalar.
+This argument is `intent(in)` and `optional`.  
+The needed length of the `complex` sequence `c`.
+
+#### Warning
+
+if `n <= size(x)`, the first `n` elements of `x` will be included in the calculation.
+
+if `n > size(x)`, the all elements of `x` and `n-size(x)` elements filled with zeros will be included in the calculation.
+
+### Example
+
+```fortran
+program demo_fft
+    use fftpack, only: fft
+    complex(kind=8) :: x(4)
+    x = [real(kind=8) :: 1.0, 2.0, 3.0, 4.0]
+    print *, fft(x)     !! [(10.0,0.0), (-2.0,2.0), (-2.0,0.0), (-2.0,-2.0)].
+    print *, fft(x,3)   !! [(6.0,0.0), (-1.5,0.86), (-1.5,0.86)].
+    print *, fft(x,5)   !! [(10.0,0.0), (-4.0,1.3), (1.5,-2.1), (1.5,2.1), (-4.0,1.3)].
+end program demo_fft
+```
+
+## `ifft`
+
+### Description
+
+Unnormalized inverse of `fft`.
+
+### Status
+
+Experimental.
+
+### Class
+
+Pure function.
+
+### Snytax
+
+`call [[fftpack(module):ifft(interface)]](x [, n])`
+
+### Argument
+
+`x`: Shall be a `complex` array.
+This argument is `intent(in)`. 
+
+`n`: Shall be an `integer` scalar.
+This argument is `intent(in)` and `optional`.  
+The needed length of the `complex` sequence `c`.
+
+#### Warning
+
+if `n <= size(x)`, the first `n` elements of `x` will be included in the calculation.
+
+if `n > size(x)`, the all elements of `x` and `n-size(x)` elements filled with zeros will be included in the calculation.
+
+### Example
+
+```fortran
+program demo_ifft
+    use fftpack, only: fft, ifft
+    complex(kind=8) :: x(4)
+    x = [real(kind=8) :: 1.0, 2.0, 3.0, 4.0]
+    print *, fft(x)             !! [(10.0,0.0), (-2.0,2.0), (-2.0,0.0), (-2.0,-2.0)].
+    print *, fft(x,3)           !! [(6.0,0.0), (-1.5,0.86), (-1.5,0.86)].
+    print *, fft(x,5)           !! [(10.0,0.0), (-4.0,1.3), (1.5,-2.1), (1.5,2.1), (-4.0,1.3)].
+    print *, ifft(fft(x))/4.0   !! [(1.0,0.0), (2.0,0.0), (3.0,0.0), (4.0,0.0)]
+    print *, ifft(fft(x), 3)    !! [(6.0,2.0), (10.3,-1.0), (13.73,-1.0)]
+end program demo_ifft
+```

fpm.toml

@@ -23,4 +23,14 @@ main = "tstfft.f"
 [[test]]
 name = "fftpack_zfft"
 source-dir = "test"
-main = "test_fftpack_zfft.f90"
+main = "test_fftpack_zfft.f90"
+
+[[test]]
+name = "fftpack_fft"
+source-dir = "test"
+main = "test_fftpack_fft.f90"
+
+[[test]]
+name = "fftpack_ifft"
+source-dir = "test"
+main = "test_fftpack_ifft.f90"

src/fftpack.f90

@@ -5,6 +5,7 @@ module fftpack
     private
 
     public :: zffti, zfftf, zfftb
+    public :: fft, ifft
 
     interface
 
@@ -42,4 +43,28 @@ end subroutine zfftb
 
     end interface
 
-end module fftpack
+    !> Version: experimental
+    !>
+    !> Forward transform of a double complex periodic sequence.
+    !> ([Specifiction](../page/specs/fftpack.html#fft))
+    interface fft
+        pure module function fft_cdp(x, n) result(result)
+            complex(kind=dp), intent(in) :: x(:)
+            integer, intent(in), optional :: n
+            complex(kind=dp), allocatable :: result(:)
+        end function fft_cdp
+    end interface fft
+
+    !> Version: experimental
+    !>
+    !> Backward transform of a double complex periodic sequence.
+    !> ([Specifiction](../page/specs/fftpack.html#ifft))
+    interface ifft
+        pure module function ifft_cdp(x, n) result(result)
+            complex(kind=dp), intent(in) :: x(:)
+            integer, intent(in), optional :: n
+            complex(kind=dp), allocatable :: result(:)
+        end function ifft_cdp
+    end interface ifft
+
+end module fftpack

src/fftpack_fft.f90

@@ -0,0 +1,36 @@
+submodule(fftpack) fftpack_fft
+
+contains
+
+    !> Forward transform of a double complex periodic sequence.
+    pure module function fft_cdp(x, n) result(result)
+        complex(kind=dp), intent(in) :: x(:)
+        integer, intent(in), optional :: n
+        complex(kind=dp), allocatable :: result(:)
+
+        integer :: lenseq, lensav, lenwrk
+        real(kind=dp), allocatable :: wsave(:)
+
+        if (present(n)) then
+            lenseq = n
+            if (lenseq <= size(x)) then
+                result = x(:lenseq)
+            else if (lenseq > size(x)) then
+                result = [x, ((0.0_dp, 0.0_dp), i=1, lenseq - size(x))]
+            end if
+        else
+            lenseq = size(x)
+            result = x
+        end if
+
+        !> Initialize FFT
+        lensav = 4*lenseq + 15
+        allocate (wsave(lensav))
+        call zffti(lenseq, wsave)
+
+        !> Forward transformation
+        call zfftf(lenseq, result, wsave)
+
+    end function fft_cdp
+
+end submodule fftpack_fft

src/fftpack_ifft.f90

@@ -0,0 +1,36 @@
+submodule(fftpack) fftpack_ifft
+
+contains
+
+    !> Backward transform of a double complex periodic sequence.
+    pure module function ifft_cdp(x, n) result(result)
+        complex(kind=dp), intent(in) :: x(:)
+        integer, intent(in), optional :: n
+        complex(kind=dp), allocatable :: result(:)
+
+        integer :: lenseq, lensav, lenwrk
+        real(kind=dp), allocatable :: wsave(:)
+
+        if (present(n)) then
+            lenseq = n
+            if (lenseq <= size(x)) then
+                result = x(:lenseq)
+            else if (lenseq > size(x)) then
+                result = [x, ((0.0_dp, 0.0_dp), i=1, lenseq - size(x))]
+            end if
+        else
+            lenseq = size(x)
+            result = x
+        end if
+
+        !> Initialize FFT
+        lensav = 4*lenseq + 15
+        allocate (wsave(lensav))
+        call zffti(lenseq, wsave)
+
+        !> Backward transformation
+        call zfftb(lenseq, result, wsave)
+
+    end function ifft_cdp
+
+end submodule fftpack_ifft

test/test_fftpack_fft.f90

@@ -0,0 +1,32 @@
+program tester
+
+    call test_fftpack_fft()
+    print *, "All tests in `test_fftpack_fft` passed."
+
+contains
+
+    subroutine check(condition, msg)
+        logical, intent(in) :: condition
+        character(*), intent(in) :: msg
+        if (.not. condition) error stop msg
+    end subroutine check
+
+    subroutine test_fftpack_fft
+        use fftpack, only: fft
+        use iso_fortran_env, only: dp => real64
+        real(kind=dp) :: eps = 1.0e-10
+
+        complex(kind=dp) :: x(3)
+
+        x = [real(kind=dp) :: 1.0, 2.0, 3.0]
+
+        call check(sum(abs(fft(x, 2) - [complex(kind=dp) ::(3.0, 0.0), (-1.0, 0.0)])) < eps, &
+                   msg="sum(abs(fft(x,2) - [complex(kind=dp) ::(3.0, 0.0), (-1.0, 0.0)])) < eps failed.")
+        call check(sum(abs(fft(x, 3) - fft(x))) < eps, &
+                   msg="sum(abs(fft(x,3) - fft(3))) < eps failed.")
+        call check(sum(abs(fft(x, 4) - [complex(kind=dp) ::(6.0, 0.0), (-2.0, -2.0), (2.0, 0.0), (-2.0, 2.0)])) < eps, &
+                   msg="sum(abs(fft(x,4) - [complex(kind=dp) ::(6.0, 0.0), (-2.0, -2.0), (2.0,0.0), (-2.0,2.0)])) < eps failed.")
+
+    end subroutine test_fftpack_fft
+
+end program tester

test/test_fftpack_ifft.f90

@@ -0,0 +1,32 @@
+program tester
+
+    call test_fftpack_ifft()
+    print *, "All tests in `test_fftpack_fft` passed."
+
+contains
+
+    subroutine check(condition, msg)
+        logical, intent(in) :: condition
+        character(*), intent(in) :: msg
+        if (.not. condition) error stop msg
+    end subroutine check
+
+    subroutine test_fftpack_ifft
+        use fftpack, only: fft, ifft
+        use iso_fortran_env, only: dp => real64
+        real(kind=dp) :: eps = 1.0e-10
+
+        complex(kind=dp) :: x(4)
+
+        x = [real(kind=dp) :: 1.0, 2.0, 3.0, 4.0]
+        
+        call check(sum(abs(ifft(fft(x))/4.0 - [complex(kind=dp) ::(1.0, 0.0), (2.0, 0.0), (3.0, 0.0), (4.0, 0.0)])) < eps, &
+        msg="abs(sum(ifft(fft(x))/4.0 - [complex(kind=dp) ::(1.0, 0.0), (2.0, 0.0), (3.0, 0.0), (4.0, 0.0)])) < eps failed.")
+        call check(sum(abs(ifft(fft(x),2) - [complex(kind=dp) ::(8.0, 2.0), (12.0, -2.0)])) < eps, &
+                   msg="abs(sum(ifft(fft(x),2) - [complex(kind=dp) ::(8.0, 2.0), (12.0, -2.0)])) < eps failed.")
+        call check(sum(abs(ifft(fft(x,2),4) - [complex(kind=dp) ::(2.0, 0.0), (3.0, -1.0), (4.0,0.0), (3.0,1.0)])) < eps, &
+        msg="abs(sum(ifft(fft(x,2),4) - [complex(kind=dp) ::(2.0, 0.0), (3.0, -1.0), (4.0,0.0), (3.0,1.0)])) < eps failed.")
+
+    end subroutine test_fftpack_ifft
+
+end program tester