add fftfreq function

Author: cval26

src/Makefile

@@ -60,7 +60,8 @@ SRCF90 = \
 	fftpack_qct.f90\
 	fftpack_iqct.f90\
 	fftpack_dct.f90\
-	rk.f90
+	rk.f90\
+	fftpack_utils.f90
 
 OBJF := $(SRCF:.f90=.o)
 OBJF90 := $(SRCF90:.f90=.o)
@@ -86,6 +87,7 @@ fftpack_iqct.o: fftpack.o rk.o
 fftpack_dct.o: fftpack.o rk.o
 fftpack_fftshift.o: fftpack.o rk.o
 fftpack_ifftshift.o: fftpack.o rk.o
+fftpack_utils.o: fftpack.o
 
 zfftb.f90: rk.o
 cfftb1.f90: rk.o

src/fftpack.f90

@@ -7,6 +7,7 @@ module fftpack
     public :: zffti, zfftf, zfftb
     public :: fft, ifft
     public :: fftshift, ifftshift
+    public :: fftfreq
 
     public :: dffti, dfftf, dfftb
     public :: rfft, irfft
@@ -175,6 +176,11 @@ pure subroutine dcost(n, x, wsave)
             real(kind=rk), intent(in) :: wsave(*)
         end subroutine dcost
 
+        pure module function fftfreq(n) result(out)
+            integer, intent(in) :: n
+            integer, dimension(n) :: out
+        end function fftfreq
+
     end interface
 
     !> Version: experimental

src/fftpack_utils.f90

@@ -0,0 +1,35 @@
+submodule(fftpack) fftpack_utils
+
+contains
+
+    !> Returns an integer array with the frequency values involved in the
+    !> performed FFT, ordered in the standard way (zero first, then positive
+    !> frequencies, then the negative ones).
+    pure module function fftfreq(n) result(out)
+        integer, intent(in) :: n
+        integer, dimension(n) :: out
+        integer :: i
+
+        out(1) = 0
+        if (n == 1) return
+
+        if (mod(n, 2) == 0) then  !> n even, smallest n = 2
+            do i = 2, n/2
+                out(i) = i-1
+            end do
+            out(n/2+1) = -n/2
+            do i = n/2+2, n  !> only enters if n/2+2 <= n
+                out(i) = out(i-1) + 1
+            end do
+        else  !> n odd, smallest n = 3
+            do i = 2, n/2+1
+                out(i) = i-1
+            end do
+            out(n/2+2) = -out(n/2+1)
+            do i = n/2+3, n  !> only enters if n/2+3 <= n
+                out(i) = out(i-1) + 1
+            end do
+        end if
+    end function fftfreq
+
+end submodule fftpack_utils