Comply with plan_inv function of AbstractFFTs (#85)

Author: mfherbst

Project.toml

@@ -6,19 +6,20 @@ version = "0.7.0"
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
 BinaryProvider = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
-FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 ToeplitzMatrices = "c751599d-da0a-543b-9d20-d0a503d91d24"
 
 [compat]
 AbstractFFTs = "0.4"
 DSP = "0.6"
-FastGaussQuadrature = "0.4"
 FFTW = "1"
+FastGaussQuadrature = "0.4"
 SpecialFunctions = "0.8"
 ToeplitzMatrices = "0.6"
 julia = "1"

src/FastTransforms.jl

@@ -1,24 +1,31 @@
 module FastTransforms
 
-using AbstractFFTs, DSP, FastGaussQuadrature, FFTW, Libdl, LinearAlgebra, SpecialFunctions, ToeplitzMatrices
+using DSP, FastGaussQuadrature, Libdl, LinearAlgebra, SpecialFunctions, ToeplitzMatrices
+using Reexport
+@reexport using AbstractFFTs
+@reexport using FFTW
 
 import Base: unsafe_convert, eltype, ndims, adjoint, transpose, show, *, \,
              inv, size, view
 
 import Base.GMP: Limb
 import Base.MPFR: BigFloat, _BigFloat
 
-import AbstractFFTs: Plan
+import AbstractFFTs: Plan, ScaledPlan,
+                     fft, ifft, bfft, fft!, ifft!, bfft!,
+                     plan_fft, plan_ifft, plan_bfft, plan_fft!, plan_ifft!, plan_bfft!,
+                     rfft, irfft, brfft, plan_rfft, plan_irfft, plan_brfft,
+                     fftshift, ifftshift,
+                     rfft_output_size, brfft_output_size,
+                     plan_inv, normalization
+
+import FFTW: dct, dct!, idct, idct!, plan_dct!, plan_idct!,
+             plan_dct, plan_idct, fftwNumber
 
 import DSP: conv
 
 import FastGaussQuadrature: unweightedgausshermite
 
-import FFTW: dct, dct!, idct, idct!,
-             plan_fft!, plan_ifft!, plan_dct!, plan_idct!,
-             plan_fft, plan_ifft, plan_rfft, plan_irfft, plan_dct, plan_idct,
-			 plan_bfft, plan_bfft!, plan_brfft, fftwNumber
-
 import LinearAlgebra: mul!, lmul!, ldiv!
 
 export leg2cheb, cheb2leg, ultra2ultra, jac2jac,

src/fftBigFloat.jl

@@ -9,7 +9,7 @@ const ComplexFloats = Complex{T} where T<:AbstractFloat
 # To add more types, add them in the union of the function's signature.
 
 function generic_fft(x::Vector{T}) where T<:AbstractFloats
-	T <: FFTW.fftwNumber && (@warn("Using generic fft for FFTW number type."))
+    T <: FFTW.fftwNumber && (@warn("Using generic fft for FFTW number type."))
     n = length(x)
     ispow2(n) && return generic_fft_pow2(x)
     ks = range(zero(real(T)),stop=n-one(real(T)),length=n)
@@ -37,8 +37,8 @@ end
 
 generic_bfft(x::Vector{T}) where {T <: AbstractFloats} = conj!(generic_fft(conj(x)))
 function generic_bfft!(x::Vector{T}) where {T <: AbstractFloats}
-	x[:] = generic_bfft(x)
-	return x
+    x[:] = generic_bfft(x)
+    return x
 end
 
 generic_brfft(v::Vector, n::Integer) = generic_irfft(v, n)*n
@@ -113,9 +113,9 @@ function generic_ifft_pow2(x::Vector{Complex{T}}) where T<:AbstractFloat
 end
 
 function generic_dct(a::AbstractVector{Complex{T}}) where {T <: AbstractFloat}
-	T <: FFTW.fftwNumber && (@warn("Using generic dct for FFTW number type."))
-	N = length(a)
-	twoN = convert(T,2) * N
+    T <: FFTW.fftwNumber && (@warn("Using generic dct for FFTW number type."))
+    N = length(a)
+    twoN = convert(T,2) * N
     c = generic_fft([a; reverse(a, dims=1)]) # c = generic_fft([a; flipdim(a,1)])
     d = c[1:N]
     d .*= exp.((-im*convert(T, pi)).*(0:N-1)./twoN)
@@ -126,9 +126,9 @@ end
 generic_dct(a::AbstractArray{T}) where {T <: AbstractFloat} = real(generic_dct(complex(a)))
 
 function generic_idct(a::AbstractVector{Complex{T}}) where {T <: AbstractFloat}
-	T <: FFTW.fftwNumber && (@warn("Using generic idct for FFTW number type."))
-	N = length(a)
-	twoN = convert(T,2)*N
+    T <: FFTW.fftwNumber && (@warn("Using generic idct for FFTW number type."))
+    N = length(a)
+    twoN = convert(T,2)*N
     b = a * sqrt(twoN)
     b[1] = b[1] * sqrt(convert(T,2))
     shift = exp.(-im * 2 * convert(T, pi) * (N - convert(T,1)/2) * (0:(2N-1)) / twoN)
@@ -154,32 +154,37 @@ end
 
 # dummy plans
 abstract type DummyPlan{T} <: Plan{T} end
-struct DummyFFTPlan{T,inplace} <: DummyPlan{T} end
-struct DummyiFFTPlan{T,inplace} <: DummyPlan{T} end
-struct DummybFFTPlan{T,inplace} <: DummyPlan{T} end
-struct DummyDCTPlan{T,inplace} <: DummyPlan{T} end
-struct DummyiDCTPlan{T,inplace} <: DummyPlan{T} end
-struct DummyrFFTPlan{T,inplace} <: DummyPlan{T}
-	n::Integer
-end
-struct DummyirFFTPlan{T,inplace} <: DummyPlan{T}
-	n::Integer
+for P in (:DummyFFTPlan, :DummyiFFTPlan, :DummybFFTPlan, :DummyDCTPlan, :DummyiDCTPlan)
+    # All plans need an initially undefined pinv field
+    @eval begin
+        mutable struct $P{T,inplace} <: DummyPlan{T}
+            pinv::DummyPlan{T}
+            $P{T,inplace}() where {T<:AbstractFloats, inplace} = new()
+        end
+    end
 end
-struct DummybrFFTPlan{T,inplace} <: DummyPlan{T}
-	n::Integer
+for P in (:DummyrFFTPlan, :DummyirFFTPlan, :DummybrFFTPlan)
+    @eval begin
+        mutable struct $P{T,inplace} <: DummyPlan{T}
+            n::Integer
+            pinv::DummyPlan{T}
+            $P{T,inplace}(n::Integer) where {T<:AbstractFloats, inplace} = new(n)
+        end
+    end
 end
 
 for (Plan,iPlan) in ((:DummyFFTPlan,:DummyiFFTPlan),
                      (:DummyDCTPlan,:DummyiDCTPlan))
    @eval begin
-       Base.inv(::$Plan{T,inplace}) where {T,inplace} = $iPlan{T,inplace}()
-       Base.inv(::$iPlan{T,inplace}) where {T,inplace} = $Plan{T,inplace}()
+       plan_inv(::$Plan{T,inplace}) where {T,inplace} = $iPlan{T,inplace}()
+       plan_inv(::$iPlan{T,inplace}) where {T,inplace} = $Plan{T,inplace}()
     end
 end
 
 # Specific for rfft, irfft and brfft:
-Base.inv(::DummyirFFTPlan{T,inplace}) where {T,inplace} = DummyrFFTPlan{T,Inplace}(p.n)
-Base.inv(::DummyrFFTPlan{T,inplace}) where {T,inplace} = DummyirFFTPlan{T,Inplace}(p.n)
+plan_inv(p::DummyirFFTPlan{T,inplace}) where {T,inplace} = DummyrFFTPlan{T,Inplace}(p.n)
+plan_inv(p::DummyrFFTPlan{T,inplace}) where {T,inplace} = DummyirFFTPlan{T,Inplace}(p.n)
+
 
 
 for (Plan,ff,ff!) in ((:DummyFFTPlan,:generic_fft,:generic_fft!),
@@ -202,14 +207,14 @@ end
 *(p::DummyirFFTPlan{T,true}, x::StridedArray{T,N}) where {T<:AbstractFloats,N} = generic_irfft!(x, p.n)
 *(p::DummyirFFTPlan{T,false}, x::StridedArray{T,N}) where {T<:AbstractFloats,N} = generic_irfft(x, p.n)
 function mul!(C::StridedVector, p::DummyirFFTPlan, x::StridedVector)
-	C[:] = generic_irfft(x, p.n)
-	C
+    C[:] = generic_irfft(x, p.n)
+    C
 end
 *(p::DummybrFFTPlan{T,true}, x::StridedArray{T,N}) where {T<:AbstractFloats,N} = generic_brfft!(x, p.n)
 *(p::DummybrFFTPlan{T,false}, x::StridedArray{T,N}) where {T<:AbstractFloats,N} = generic_brfft(x, p.n)
 function mul!(C::StridedVector, p::DummybrFFTPlan, x::StridedVector)
-	C[:] = generic_brfft(x, p.n)
-	C
+    C[:] = generic_brfft(x, p.n)
+    C
 end
 
 

test/fftBigFloattests.jl

@@ -26,6 +26,14 @@ using FastTransforms, FFTW, Test
     @test norm(idct(dct(c))-c,Inf) < 1000eps(BigFloat)
     @test norm(dct(idct(c))-c,Inf) < 1000eps(BigFloat)
 
+    c = randn(ComplexF16, 20)
+    p = plan_fft(c)
+    @test inv(p) * (p * c) ≈ c
+
+    c = randn(ComplexF16, 20)
+    pinpl = plan_fft!(c)
+    @test inv(pinpl) * (pinpl * c) ≈ c
+
     # Make sure we don't accidentally hijack any FFTW plans
     for T in (Float32, Float64)
         @test plan_fft(rand(BigFloat,10)) isa FastTransforms.DummyPlan