Mistake in real FFT

Author: dannys4

src/FFTA.jl

@@ -45,6 +45,13 @@ function bfft(x::AbstractVector{T}) where {T}
     y
 end
 
+function bfft(x::AbstractVector{T}) where {T <: Real}
+    y = similar(x, Complex{T})
+    g = CallGraph{Complex{T}}(length(x))
+    fft!(y, x, 1, 1, FFT_BACKWARD(), g[1].type, g, 1)
+    y
+end
+
 function bfft(x::AbstractMatrix{T}) where {T}
     M,N = size(x)
     y1 = similar(x)

src/algos.jl

@@ -62,8 +62,7 @@ end
 Power of 2 FFT in place, complex
 
 """
-function fft_pow2!(out::AbstractVector{T}, in::AbstractVector{T}, N::Int, start_out::Int, stride_out::Int, start_in::Int, stride_in::Int, d::Direction) where {T}
-
+function fft_pow2!(out::AbstractVector{T}, in::AbstractVector{U}, N::Int, start_out::Int, stride_out::Int, start_in::Int, stride_in::Int, d::Direction) where {T, U}
     if N == 2
         out[start_out]              = in[start_in] + in[start_in + stride_in]
         out[start_out + stride_out] = in[start_in] - in[start_in + stride_in]
@@ -86,35 +85,6 @@ function fft_pow2!(out::AbstractVector{T}, in::AbstractVector{T}, N::Int, start_
     end
 end
 
-"""
-Power of 2 FFT in place, real
-
-"""
-function fft_pow2!(out::AbstractVector{Complex{T}}, in::AbstractVector{T}, N::Int, start_out::Int, stride_out::Int, start_in::Int, stride_in::Int, d::Direction) where {T<:Real}
-    if N == 2
-        out[1] = in[1] + in[2]
-        out[2] = in[1] - in[2]
-        return
-    end
-    m = N ÷ 2
-    fft_pow2!(out, in, m, start_out    , stride_out, start_in    , stride_in*2, d)
-    fft_pow2!(out, in, m, start_out + m, stride_out, start_in + 1, stride_in*2, d)
-
-    w1 = convert(Complex{T}, cispi(direction_sign(d)*2/N))
-    wj = one(Complex{T})
-    @inbounds @turbo for j in 1:m-1
-        j1_out = start_out + j*stride_out
-        j2_out = start_out + (j+m)*stride_out
-        out[j1_out] = out[j1_out] + wj*out[j2_out]
-        wj *= w1
-    end
-    @inbounds @turbo for j in 1:m-1
-        j1_out = start_out + (j+m)*stride_out
-        j2_out = start_out + (m-j+1)*stride_out
-        out[j1_out] = conj(out[j2_out])
-    end
-end
-
 function fft_dft!(out::AbstractVector{T}, in::AbstractVector{T}, N::Int, start_out::Int, stride_out::Int, start_in::Int, stride_in::Int, d::Direction) where {T}
     tmp = in[start_in]
     @inbounds for j in 1:N-1
@@ -139,22 +109,26 @@ function fft_dft!(out::AbstractVector{Complex{T}}, in::AbstractVector{T}, N::Int
     halfN = N÷2
     wk = wkn = w = convert(Complex{T}, cispi(direction_sign(d)*2/N))
 
-    out[start_out + 1:stride_out:start_out+stride_out*N] .= in[1]
-    out[1] = sum(@view in[start_in:stride_in:start_int+stride_out*N])
-    iseven(N) && (out[start_out + stride_out*halfN] = alternatingSum(@view in[start_in:stride_in:start_int+stride_out*N]))
+    tmpBegin = tmpHalf = in[start_in]
+    @inbounds for j in 1:N-1
+        tmpBegin += in[start_in + stride_in*j]
+        iseven(j) ? tmpHalf += in[start_in + stride_in*j] : tmpHalf -= in[start_in + stride_in*j]
+    end
+    out[start_out] = convert(Complex{T}, tmpBegin)
+    iseven(N) && (out[start_out + stride_out*halfN] = convert(Complex{T}, tmpHalf))
 
-    @inbounds for d in 2:halfN+1
+    @inbounds for d in 1:halfN
         tmp = in[start_in]
-        @inbounds for k in 2:N
+        @inbounds for k in 1:N-1
             tmp += wkn*in[start_in + k*stride_in]
             wkn *= wk
         end
         out[start_out + d*stride_out] = tmp
         wk *= w
         wkn = wk
     end
-    @inbounds @turbo for i in 0:halfN-1
-        out[start_out + stride_out*(N-i)] = conj(out[start_out + stride_out*(halfN-i)])
+    @inbounds @turbo for i in 1:halfN-1
+        out[start_out + stride_out*(N-i)] = conj(out[start_out + stride_out*i])
     end
 end
 

test/complex_backward.jl

@@ -4,8 +4,8 @@ test_nums = [8, 11, 15, 100]
     for N in test_nums
         x = ones(ComplexF64, N)
         y = bfft(x)
-        b1 = isapprox(y[1], N, atol=1e-12)
-        b2 = isapprox(y[2:end], 0*x[2:end], atol=1e-12)
-        @test b1 && b2
+        y_ref = 0*y
+        y_ref[1] = N
+        @test y ≈ y_ref atol=1e-12
     end
 end

test/complex_forward.jl

@@ -2,9 +2,10 @@ using FFTA, Test
 test_nums = [8, 11, 15, 100]
 @testset verbose = true " forward" begin
     for N in test_nums
-        x = zeros(ComplexF64, N)
-        x[1] = 1
+        x = ones(ComplexF64, N)
         y = fft(x)
-        @test y ≈ ones(size(x))
+        y_ref = 0*y
+        y_ref[1] = N
+        @test y ≈ y_ref atol=1e-12
     end
 end

test/real_backward.jl

@@ -0,0 +1,11 @@
+using FFTA, Test
+test_nums = [8, 11, 15, 100]
+@testset "backward" begin
+    for N in test_nums
+        x = ones(Float64, N)
+        y = bfft(x)
+        y_ref = 0*y
+        y_ref[1] = N
+        @test y_ref ≈ y atol=1e-12
+    end
+end

test/real_forward.jl

@@ -1,8 +1,11 @@
 using FFTA, Test
 test_nums = [8, 11, 15, 100]
-@testset verbose = true "fft 1D real, size $(padnum(maximum(test_nums),N))" for N in test_nums
-    x = zeros(Float64, N)
-    x[1] = 1
-    y = fft(x)
-    @test y ≈ ones(size(x))
+@testset verbose = true "forward" begin 
+    for N in test_nums
+        x = ones(Float64, N)
+        y = fft(x)
+        y_ref = 0*y
+        y_ref[1] = N
+        @test y ≈ y_ref atol=1e-12
+    end
 end

test/runtests.jl

@@ -14,7 +14,7 @@ end
         include("complex_backward.jl")
     end
     @testset verbose = true "Real" begin
-        # include("real_forward.jl")
+        include("real_forward.jl")
         include("real_backward.jl")
     end
 end