Merge pull request #5 from fortran-fans/add_shift

Add `fftshift/ifftshift` utility function.

Author: zoziha

doc/specs/fftpack.md

@@ -32,7 +32,7 @@ Pure function.
 This argument is `intent(in)`.  
 The length of the sequence to be transformed.
 
-`wsave`: Shall be a `real` array.
+`wsave`: Shall be a `real` and rank-1 array.
 This argument is `intent(out)`.  
 A work array which must be dimensioned at least `4*n+15`.
 The same work array can be used for both `zfftf` and `zfftb`
@@ -88,7 +88,7 @@ Pure function.
 This argument is `intent(in)`.  
 The length of the `complex` sequence `c`. The method is more efficient when `n` is the product of small primes.
 
-`c`: Shall be a `complex` array.
+`c`: Shall be a `complex` and rank-1 array.
 This argument is `intent(inout)`.  
 A `complex` array of length `n` which contains the sequence.
 ```
@@ -101,7 +101,7 @@ for j=1,...,n
                        where i=sqrt(-1)
 ```
 
-`wsave`: Shall be a `real` array.
+`wsave`: Shall be a `real` and rank-1 array.
 This argument is `intent(inout)`.  
 A `real` work array which must be dimensioned at least `4n+15` in the program that calls `zfftf`. 
 The wsave array must be initialized by calling subroutine `zffti(n,wsave)` and a different `wsave` array must be used for each different value of `n`.  
@@ -162,7 +162,7 @@ Pure function.
 This argument is `intent(in)`.  
 The length of the `complex` sequence `c`. The method is more efficient when `n` is the product of small primes.
 
-`c`: Shall be a `complex` array.
+`c`: Shall be a `complex` and rank-1 array.
 This argument is `intent(inout)`.  
 A `complex` array of length `n` which contains the sequence.
 ```
@@ -175,7 +175,7 @@ for j=1,...,n
                        where i=sqrt(-1)
 ```
 
-`wsave`: Shall be a `real` array.
+`wsave`: Shall be a `real` and rank-1 array.
 This argument is `intent(inout)`.  
 A `real` work array which must be dimensioned at least `4n+15` in the program that calls `zfftf`. The `wsave` array must be initialized by calling subroutine `zffti(n,wsave)` and a different `wsave` array must be used for each different value of `n`. This initialization does not have to be repeated so long as `n` remains unchanged thus subsequent transforms can be obtained faster than the first. The same `wsave` array can be used by `zfftf` and `zfftb`.  
 Contains initialization calculations which must not be destroyed between calls of subroutine `zfftf` or `zfftb`.
@@ -218,7 +218,7 @@ Pure function.
 
 ### Argument
 
-`x`: Shall be a `complex` array.
+`x`: Shall be a `complex` and rank-1 array.
 This argument is `intent(in)`. 
 
 `n`: Shall be an `integer` scalar.
@@ -270,7 +270,7 @@ Pure function.
 
 ### Argument
 
-`x`: Shall be a `complex` array.
+`x`: Shall be a `complex` and rank-1 array.
 This argument is `intent(in)`. 
 
 `n`: Shall be an `integer` scalar.
@@ -303,3 +303,85 @@ program demo_ifft
     print *, ifft(fft(x), 3)    !! [(6.0,2.0), (10.3,-1.0), (13.73,-1.0)]
 end program demo_ifft
 ```
+
+## `fftshift`
+
+### Description
+
+Rearranges the Fourier transform by moving the zero-frequency component to the center of the array.  
+
+### Status
+
+Experimental.
+
+### Class
+
+Pure function.
+
+### Snytax
+
+`result = [[fftpack(module):fftshift(interface)]](x)`
+
+### Argument
+
+`x`: Shall be a `complex/real` and rank-1 array.
+This argument is `intent(in)`. 
+
+### Return value
+
+Returns the `complex/real` and rank-1 Fourier transform by moving the zero-frequency component to the center of the array.
+
+### Example
+
+```fortran
+program demo_fftshift
+    use fftpack, only: fftshift
+    complex(kind=8) :: c(5) = [1, 2, 3, 4, 5]
+    real(kind=8) :: x(5) = [1, 2, 3, 4, 5]
+    print *, fftshift(c(1:4))   !! [(3.0,0.0), (4.0,0.0), (1.0,0.0), (2.0,0.0)]
+    print *, fftshift(c)        !! [(4.0,0.0), (5.0,0.0), (1.0,0.0), (2.0,0.0), (3.0,0.0)]
+    print *, fftshift(x(1:4))   !! [3.0, 4.0, 1.0, 2.0]
+    print *, fftshift(x)        !! [4.0, 5.0, 1.0, 2.0, 3.0]
+end program demo_fftshift
+```
+
+## `ifftshift`
+
+### Description
+
+Rearranges the Fourier transform with zero frequency shifting back to the original transform output. In other words, `ifftshift` is the result of undoing `fftshift`.
+
+### Status
+
+Experimental.
+
+### Class
+
+Pure function.
+
+### Snytax
+
+`result = [[fftpack(module):ifftshift(interface)]](x)`
+
+### Argument
+
+`x`: Shall be a `complex/real` and rank-1 array.
+This argument is `intent(in)`. 
+
+### Return value
+
+Returns the `complex/real` and rank-1 Fourier transform with zero frequency shifting back to the original transform output.
+
+### Example
+
+```fortran
+program demo_ifftshift
+    use fftpack, only: fftshift, ifftshift
+    complex(kind=8) :: c(5) = [1, 2, 3, 4, 5]
+    real(kind=8) :: x(5) = [1, 2, 3, 4, 5]
+    print *, ifftshift(fftshift(c(1:4)))   !! [(1.0,0.0), (2.0,0.0), (3.0,0.0), (4.0,0.0)]
+    print *, ifftshift(fftshift(c) )       !! [(1.0,0.0), (2.0,0.0), (3.0,0.0), (4.0,0.0), (5.0,0.0)]
+    print *, ifftshift(fftshift(x(1:4)))   !! [1.0, 2.0, 3.0, 4.0]
+    print *, ifftshift(fftshift(x))        !! [1.0, 2.0, 3.0, 4.0, 5.0]
+end program demo_ifftshift
+```

fpm.toml

@@ -33,4 +33,14 @@ main = "test_fftpack_fft.f90"
 [[test]]
 name = "fftpack_ifft"
 source-dir = "test"
-main = "test_fftpack_ifft.f90"
+main = "test_fftpack_ifft.f90"
+
+[[test]]
+name = "fftpack_fftshift"
+source-dir = "test"
+main = "test_fftpack_fftshift.f90"
+
+[[test]]
+name = "fftpack_ifftshift"
+source-dir = "test"
+main = "test_fftpack_ifftshift.f90"

src/fftpack.f90

@@ -1,11 +1,12 @@
 module fftpack
 
-    use iso_fortran_env, only: dp => real64
+    use, intrinsic :: iso_fortran_env, only: dp => real64
     implicit none
     private
 
     public :: zffti, zfftf, zfftb
     public :: fft, ifft
+    public :: fftshift, ifftshift
 
     interface
 
@@ -67,4 +68,34 @@ pure module function ifft_cdp(x, n) result(result)
         end function ifft_cdp
     end interface ifft
 
+    !> Version: experimental
+    !>
+    !> Shifts zero-frequency component to center of spectrum.
+    !> ([Specifiction](../page/specs/fftpack.html#fftshift))
+    interface fftshift
+        pure module function fftshift_cdp(x) result(result)
+            complex(kind=dp), intent(in) :: x(:)
+            complex(kind=dp), allocatable :: result(:)
+        end function fftshift_cdp
+        pure module function fftshift_rdp(x) result(result)
+            real(kind=dp), intent(in) :: x(:)
+            real(kind=dp), allocatable :: result(:)
+        end function fftshift_rdp
+    end interface fftshift
+
+    !> Version: experimental
+    !>
+    !> Shifts zero-frequency component to beginning of spectrum.
+    !> ([Specifiction](../page/specs/fftpack.html#ifftshift))
+    interface ifftshift
+        pure module function ifftshift_cdp(x) result(result)
+            complex(kind=dp), intent(in) :: x(:)
+            complex(kind=dp), allocatable :: result(:)
+        end function ifftshift_cdp
+        pure module function ifftshift_rdp(x) result(result)
+            real(kind=dp), intent(in) :: x(:)
+            real(kind=dp), allocatable :: result(:)
+        end function ifftshift_rdp
+    end interface ifftshift
+
 end module fftpack

src/fftpack_fftshift.f90

@@ -0,0 +1,23 @@
+submodule(fftpack) fftpack_fftshift
+
+contains
+
+    !> Shifts zero-frequency component to center of spectrum for `complex` type.
+    pure module function fftshift_cdp(x) result(result)
+        complex(kind=dp), intent(in) :: x(:)
+        complex(kind=dp), allocatable :: result(:)
+
+        result = cshift(x, shift=-floor(0.5_dp*size(x)))
+
+    end function fftshift_cdp
+
+    !> Shifts zero-frequency component to center of spectrum for `real` type.
+    pure module function fftshift_rdp(x) result(result)
+        real(kind=dp), intent(in) :: x(:)
+        real(kind=dp), allocatable :: result(:)
+
+        result = cshift(x, shift=-floor(0.5_dp*size(x)))
+
+    end function fftshift_rdp
+
+end submodule fftpack_fftshift

src/fftpack_ifftshift.f90

@@ -0,0 +1,23 @@
+submodule(fftpack) fftpack_ifftshift
+
+contains
+
+    !> Shifts zero-frequency component to beginning of spectrum for `complex` type.
+    pure module function ifftshift_cdp(x) result(result)
+        complex(kind=dp), intent(in) :: x(:)
+        complex(kind=dp), allocatable :: result(:)
+
+        result = cshift(x, shift=-ceiling(0.5_dp*size(x)))
+
+    end function ifftshift_cdp
+
+    !> Shifts zero-frequency component to beginning of spectrum for `real` type.
+    pure module function ifftshift_rdp(x) result(result)
+        real(kind=dp), intent(in) :: x(:)
+        real(kind=dp), allocatable :: result(:)
+
+        result = cshift(x, shift=-ceiling(0.5_dp*size(x)))
+
+    end function ifftshift_rdp
+
+end submodule fftpack_ifftshift

test/test_fftpack_fftshift.f90

@@ -0,0 +1,43 @@
+program tester
+
+    call test_fftpack_fftshift_complex
+    call test_fftpack_fftshift_real
+    print *, "All tests in `test_fftpack_fftshift` 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_fftshift_complex
+        use fftpack, only: fftshift
+        use iso_fortran_env, only: dp => real64
+
+        complex(kind=dp) :: xeven(4) = [1, 2, 3, 4]
+        complex(kind=dp) :: xodd(5) = [1, 2, 3, 4, 5]
+
+        call check(all(fftshift(xeven) == [complex(kind=dp) :: 3, 4, 1, 2]), &
+                   msg="all(fftshift(xeven) == [complex(kind=dp) :: 3, 4, 1, 2]) failed.")
+        call check(all(fftshift(xodd) == [complex(kind=dp) :: 4, 5, 1, 2, 3]), &
+                   msg="all(fftshift(xodd) == [complex(kind=dp) :: 4, 5, 1, 2, 3]) failed.")
+
+    end subroutine test_fftpack_fftshift_complex
+
+    subroutine test_fftpack_fftshift_real
+        use fftpack, only: fftshift
+        use iso_fortran_env, only: dp => real64
+
+        real(kind=dp) :: xeven(4) = [1, 2, 3, 4]
+        real(kind=dp) :: xodd(5) = [1, 2, 3, 4, 5]
+
+        call check(all(fftshift(xeven) == [real(kind=dp) :: 3, 4, 1, 2]), &
+                   msg="all(fftshift(xeven) == [real(kind=dp) :: 3, 4, 1, 2]) failed.")
+        call check(all(fftshift(xodd) == [real(kind=dp) :: 4, 5, 1, 2, 3]), &
+                   msg="all(fftshift(xodd) == [real(kind=dp) :: 4, 5, 1, 2, 3]) failed.")
+
+    end subroutine test_fftpack_fftshift_real
+
+end program tester

test/test_fftpack_ifft.f90

@@ -1,7 +1,7 @@
 program tester
 
     call test_fftpack_ifft()
-    print *, "All tests in `test_fftpack_fft` passed."
+    print *, "All tests in `test_fftpack_ifft` passed."
 
 contains
 

test/test_fftpack_ifftshift.f90

@@ -0,0 +1,45 @@
+program tester
+
+    call test_fftpack_ifftshift_complex
+    call test_fftpack_ifftshift_real
+    print *, "All tests in `test_fftpack_ifftshift` 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_ifftshift_complex
+        use fftpack, only: ifftshift
+        use iso_fortran_env, only: dp => real64
+        integer :: i
+
+        complex(kind=dp) :: xeven(4) = [3, 4, 1, 2]
+        complex(kind=dp) :: xodd(5) = [4, 5, 1, 2, 3]
+
+        call check(all(ifftshift(xeven) == [complex(kind=dp) ::(i, i=1, 4)]), &
+                   msg="all(ifftshift(xeven) == [complex(kind=dp) ::(i, i=1, 4)]) failed.")
+        call check(all(ifftshift(xodd) == [complex(kind=dp) ::(i, i=1, 5)]), &
+                   msg="all(ifftshift(xodd) == [complex(kind=dp) ::(i, i=1, 5)]) failed.")
+
+    end subroutine test_fftpack_ifftshift_complex
+
+    subroutine test_fftpack_ifftshift_real
+        use fftpack, only: ifftshift
+        use iso_fortran_env, only: dp => real64
+        integer :: i
+
+        real(kind=dp) :: xeven(4) = [3, 4, 1, 2]
+        real(kind=dp) :: xodd(5) = [4, 5, 1, 2, 3]
+
+        call check(all(ifftshift(xeven) == [real(kind=dp) ::(i, i=1, 4)]), &
+                   msg="all(ifftshift(xeven) == [real(kind=dp) ::(i, i=1, 4)]) failed.")
+        call check(all(ifftshift(xodd) == [real(kind=dp) ::(i, i=1, 5)]), &
+                   msg="all(ifftshift(xodd) == [real(kind=dp) ::(i, i=1, 5)]) failed.")
+
+    end subroutine test_fftpack_ifftshift_real
+
+end program tester