Power of 2 FFT in place

Author: dannys4

Project.toml

@@ -4,6 +4,7 @@ authors = ["Danny Sharp <[email protected]> and contributors"]
 version = "0.1.0"
 
 [deps]
+LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 

src/NSFFT.jl

@@ -1,6 +1,6 @@
 module NSFFT
 
-using Primes, StaticArrays
+using Primes, StaticArrays, LoopVectorization
 export pow2FFT!, myfft_sv
 
 include("algos.jl")

src/algos.jl

@@ -1,5 +1,9 @@
 @enum Direction FFT_FORWARD FFT_BACKWARD
 
+"""
+Power of 2 FFT in place, forward
+
+"""
 function pow2FFT!(out::AbstractVector{T}, in::AbstractVector{T}, ::Val{FFT_FORWARD}) where {T<:Complex}
     N = length(out)
     if N == 1
@@ -19,4 +23,29 @@ function pow2FFT!(out::AbstractVector{T}, in::AbstractVector{T}, ::Val{FFT_FORWA
         out[j+m] = out_j - wj*out[j+m]
         wj *= w1
     end
+end
+
+"""
+Power of 2 FFT in place, backward
+
+"""
+function pow2FFT!(out::AbstractVector{T}, in::AbstractVector{T}, ::Val{FFT_BACKWARD}) where {T<:Complex}
+    N = length(out)
+    if N == 1
+        out[1] = in[1]
+        return
+    end
+    pow2FFT!(@view(out[1:(end÷2)]), @view(in[1:2:end]), Val(FFT_BACKWARD))
+    pow2FFT!(@view(out[(end÷2+1):end]), @view(in[2:2:end]), Val(FFT_BACKWARD))
+
+    inc = 2*π/N
+    w1 = T(cos(inc), sin(inc));
+    wj = T(1,0)
+    m = N ÷ 2
+    for j in 1:m
+        out_j    = out[j]
+        out[j]   = out_j + wj*out[j+m]
+        out[j+m] = out_j - wj*out[j+m]
+        wj *= w1
+    end
 end