Merge pull request #64 from daanhb/genericfft

Genericfft

Author: MikaelSlevinsky

src/fftBigFloat.jl

@@ -1,57 +1,59 @@
-const BigFloats = Union{BigFloat,Complex{BigFloat}}
+const AbstractFloats = Union{AbstractFloat,Complex{T} where T<:AbstractFloat}
 
 if VERSION < v"0.7-"
-    import Base.FFTW: fft, fft!, rfft, irfft, ifft, conv, dct, idct, dct!, idct!,
+    import Base.FFTW: fft, fft!, rfft, irfft, ifft, ifft!, conv, dct, idct, dct!, idct!,
                         plan_fft!, plan_ifft!, plan_dct!, plan_idct!,
                         plan_fft, plan_ifft, plan_rfft, plan_irfft, plan_dct, plan_idct
 else
-    import FFTW: fft, fft!, rfft, irfft, ifft, dct, idct, dct!, idct!,
-                        plan_fft!, plan_ifft!, plan_dct!, plan_idct!,
-                        plan_fft, plan_ifft, plan_rfft, plan_irfft, plan_dct, plan_idct
+    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
     import DSP: conv
 end
 
 # The following implements Bluestein's algorithm, following http://www.dsprelated.com/dspbooks/mdft/Bluestein_s_FFT_Algorithm.html
 # To add more types, add them in the union of the function's signature.
-function fft(x::Vector{T}) where T<:BigFloats
+
+function generic_fft(x::Vector{T}) where T<:AbstractFloats
     n = length(x)
-    ispow2(n) && return fft_pow2(x)
+    ispow2(n) && return generic_fft_pow2(x)
     ks = range(zero(real(T)),stop=n-one(real(T)),length=n)
     Wks = exp.((-im).*convert(T,π).*ks.^2 ./ n)
     xq, wq = x.*Wks, conj([exp(-im*convert(T,π)*n);reverse(Wks);Wks[2:end]])
     return Wks.*conv(xq,wq)[n+1:2n]
 end
 
 
-function fft!(x::Vector{T}) where T<:BigFloats
-    x[:] = fft(x)
+function generic_fft!(x::Vector{T}) where T<:AbstractFloats
+    x[:] = generic_fft(x)
     return x
 end
 
-# add rfft for BigFloat, by calling fft
+# add rfft for AbstractFloat, by calling fft
 #  this creates ToeplitzMatrices.rfft, so avoids changing rfft
 
-rfft(v::Vector{T}) where T<:BigFloats = fft(v)[1:div(length(v),2)+1]
-function irfft(v::Vector{T},n::Integer) where T<:BigFloats
+generic_rfft(v::Vector{T}) where T<:AbstractFloats = generic_fft(v)[1:div(length(v),2)+1]
+
+function generic_irfft(v::Vector{T},n::Integer) where T<:AbstractFloats
     @assert n==2length(v)-1
-    r = Vector{Complex{BigFloat}}(undef, n)
+    r = Vector{Complex{real(T)}}(undef, n)
     r[1:length(v)]=v
     r[length(v)+1:end]=reverse(conj(v[2:end]))
-    real(ifft(r))
+    real(generic_ifft(r))
 end
 
-ifft(x::Vector{T}) where {T<:BigFloats} = conj!(fft(conj(x)))/length(x)
-function ifft!(x::Vector{T}) where T<:BigFloats
-    x[:] = ifft(x)
+generic_ifft(x::Vector{T}) where {T<:AbstractFloats} = conj!(generic_fft(conj(x)))/length(x)
+function generic_ifft!(x::Vector{T}) where T<:AbstractFloats
+    x[:] = generic_ifft(x)
     return x
 end
 
-function conv(u::StridedVector{T}, v::StridedVector{T}) where T<:BigFloats
+function conv(u::StridedVector{T}, v::StridedVector{T}) where T<:AbstractFloats
     nu,nv = length(u),length(v)
     n = nu + nv - 1
     np2 = nextpow(2,n)
     append!(u,zeros(T,np2-nu)),append!(v,zeros(T,np2-nv))
-    y = ifft_pow2(fft_pow2(u).*fft_pow2(v))
+    y = generic_ifft_pow2(generic_fft_pow2(u).*generic_fft_pow2(v))
     #TODO This would not handle Dual/ComplexDual numbers correctly
     y = T<:Real ? real(y[1:n]) : y[1:n]
 end
@@ -60,7 +62,7 @@ end
 # c_radix2.c in the GNU Scientific Library and four1 in the Numerical Recipes in C.
 # However, the trigonometric recurrence is improved for greater efficiency.
 # The algorithm starts with bit-reversal, then divides and conquers in-place.
-function fft_pow2!(x::Vector{T}) where T<:BigFloat
+function generic_fft_pow2!(x::Vector{T}) where T<:AbstractFloat
     n,big2=length(x),2one(T)
     nn,j=n÷2,1
     for i=1:2:n-1
@@ -96,97 +98,129 @@ function fft_pow2!(x::Vector{T}) where T<:BigFloat
     return x
 end
 
-function fft_pow2(x::Vector{Complex{T}}) where T<:BigFloat
+function generic_fft_pow2(x::Vector{Complex{T}}) where T<:AbstractFloat
     y = interlace(real(x),imag(x))
-    fft_pow2!(y)
+    generic_fft_pow2!(y)
     return complex.(y[1:2:end],y[2:2:end])
 end
-fft_pow2(x::Vector{T}) where {T<:BigFloat} = fft_pow2(complex(x))
+generic_fft_pow2(x::Vector{T}) where {T<:AbstractFloat} = generic_fft_pow2(complex(x))
 
-function ifft_pow2(x::Vector{Complex{T}}) where T<:BigFloat
+function generic_ifft_pow2(x::Vector{Complex{T}}) where T<:AbstractFloat
     y = interlace(real(x),-imag(x))
-    fft_pow2!(y)
+    generic_fft_pow2!(y)
     return complex.(y[1:2:end],-y[2:2:end])/length(x)
 end
 
-
-function dct(a::AbstractVector{Complex{BigFloat}})
-	N = big(length(a))
-    c = fft([a; flipdim(a,1)])
+function generic_dct(a::AbstractVector{Complex{T}}) where {T <: AbstractFloat}
+	N = length(a)
+	twoN = convert(T,2) * N
+    c = generic_fft([a; flipdim(a,1)])
     d = c[1:N]
-    d .*= exp.((-im*big(pi)).*(0:N-1)./(2*N))
-    d[1] = d[1] / sqrt(big(2))
-    lmul!(inv(sqrt(2N)), d)
+    d .*= exp.((-im*convert(T, pi)).*(0:N-1)./twoN)
+    d[1] = d[1] / sqrt(convert(T, 2))
+    lmul!(inv(sqrt(twoN)), d)
 end
 
-dct(a::AbstractArray{BigFloat}) = real(dct(complex(a)))
+generic_dct(a::AbstractArray{T}) where {T <: AbstractFloat} = real(generic_dct(complex(a)))
 
-function idct(a::AbstractVector{Complex{BigFloat}})
-	N = big(length(a))
-    b = a * sqrt(2*N)
-    b[1] = b[1] * sqrt(big(2))
-    shift = exp.(-im * 2 * big(pi) * (N - big(1)/2) * (0:(2N-1)) / (2 * N))
+function generic_idct(a::AbstractVector{Complex{T}}) where {T <: AbstractFloat}
+	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)
     b = [b; 0; -flipdim(b[2:end],1)] .* shift
     c = ifft(b)
     flipdim(c[1:N],1)
 end
 
-idct(a::AbstractArray{BigFloat}) = real(idct(complex(a)))
+generic_idct(a::AbstractArray{T}) where {T <: AbstractFloat} = real(generic_idct(complex(a)))
 
-dct!(a::AbstractArray{T}) where {T<:BigFloats} = (b = dct(a); a[:] = b)
-idct!(a::AbstractArray{T}) where {T<:BigFloats} = (b = idct(a); a[:] = b)
+generic_dct!(a::AbstractArray{T}) where {T<:AbstractFloats} = (b = generic_dct(a); a[:] = b)
+generic_idct!(a::AbstractArray{T}) where {T<:AbstractFloats} = (b = generic_idct(a); a[:] = b)
+
+# These lines mimick the corresponding ones in FFTW/src/dct.jl, but with
+# AbstractFloat rather than fftwNumber.
+for f in (:dct, :dct!, :idct, :idct!)
+    pf = Symbol("plan_", f)
+    @eval begin
+        $f(x::AbstractArray{<:AbstractFloats}) = $pf(x) * x
+        $f(x::AbstractArray{<:AbstractFloats}, region) = $pf(x, region) * x
+    end
+end
 
 # dummy plans
 struct DummyFFTPlan{T,inplace} <: Plan{T} end
 struct DummyiFFTPlan{T,inplace} <: Plan{T} end
-struct DummyrFFTPlan{T,inplace} <: Plan{T} end
-struct DummyirFFTPlan{T,inplace} <: Plan{T} end
 struct DummyDCTPlan{T,inplace} <: Plan{T} end
 struct DummyiDCTPlan{T,inplace} <: Plan{T} end
+struct DummyrFFTPlan{T,inplace} <: Plan{T}
+	n	::	Integer
+end
+struct DummyirFFTPlan{T,inplace} <: Plan{T}
+	n	::	Integer
+end
 
 for (Plan,iPlan) in ((:DummyFFTPlan,:DummyiFFTPlan),
-                     (:DummyrFFTPlan,:DummyirFFTPlan),
+                     # (:DummyrFFTPlan,:DummyirFFTPlan),
                      (: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}()
     end
 end
 
+# Specific for rfft and irfft:
+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)
 
-for (Plan,ff,ff!) in ((:DummyFFTPlan,:fft,:fft!),
-                      (:DummyiFFTPlan,:ifft,:ifft!),
-                      (:DummyrFFTPlan,:rfft,:rfft!),
-                      (:DummyirFFTPlan,:irfft,:irfft!),
-                      (:DummyDCTPlan,:dct,:dct!),
-                      (:DummyiDCTPlan,:idct,:idct!))
+
+for (Plan,ff,ff!) in ((:DummyFFTPlan,:generic_fft,:generic_fft!),
+                      (:DummyiFFTPlan,:generic_ifft,:generic_ifft!),
+                      (:DummyrFFTPlan,:generic_rfft,:generic_rfft!),
+#                      (:DummyirFFTPlan,:generic_irfft,:generic_irfft!),
+                      (:DummyDCTPlan,:generic_dct,:generic_dct!),
+                      (:DummyiDCTPlan,:generic_idct,:generic_idct!))
     @eval begin
-        *(p::$Plan{T,true}, x::StridedArray{T,N}) where {T,N} = $ff!(x)
-        *(p::$Plan{T,false}, x::StridedArray{T,N}) where {T,N} = $ff(x)
+        *(p::$Plan{T,true}, x::StridedArray{T,N}) where {T<:AbstractFloats,N} = $ff!(x)
+        *(p::$Plan{T,false}, x::StridedArray{T,N}) where {T<:AbstractFloats,N} = $ff(x)
         function LAmul!(C::StridedVector, p::$Plan, x::StridedVector)
             C[:] = $ff(x)
             C
         end
     end
 end
 
+# Specific for irfft:
+*(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 LAmul!(C::StridedVector, p::DummyirFFTPlan, x::StridedVector)
+	C[:] = generic_irfft(x, p.n)
+	C
+end
 
-
-
-plan_fft!(x::Vector{T}) where {T<:BigFloats} = DummyFFTPlan{Complex{BigFloat},true}()
-plan_ifft!(x::Vector{T}) where {T<:BigFloats} = DummyiFFTPlan{Complex{BigFloat},true}()
-
-#plan_rfft!{T<:BigFloats}(x::Vector{T}) = DummyrFFTPlan{Complex{BigFloat},true}()
-#plan_irfft!{T<:BigFloats}(x::Vector{T},n::Integer) = DummyirFFTPlan{Complex{BigFloat},true}()
-plan_dct!(x::Vector{T}) where {T<:BigFloats} = DummyDCTPlan{T,true}()
-plan_idct!(x::Vector{T}) where {T<:BigFloats} = DummyiDCTPlan{T,true}()
-
-plan_fft(x::Vector{T}) where {T<:BigFloats} = DummyFFTPlan{Complex{BigFloat},false}()
-plan_ifft(x::Vector{T}) where {T<:BigFloats} = DummyiFFTPlan{Complex{BigFloat},false}()
-plan_rfft(x::Vector{T}) where {T<:BigFloats} = DummyrFFTPlan{Complex{BigFloat},false}()
-plan_irfft(x::Vector{T},n::Integer) where {T<:BigFloats} = DummyirFFTPlan{Complex{BigFloat},false}()
-plan_dct(x::Vector{T}) where {T<:BigFloats} = DummyDCTPlan{T,false}()
-plan_idct(x::Vector{T}) where {T<:BigFloats} = DummyiDCTPlan{T,false}()
+# We override these for AbstractFloat, so that conversion from reals to
+# complex numbers works for any AbstractFloat (instead of only BlasFloat's)
+AbstractFFTs.complexfloat(x::StridedArray{Complex{<:AbstractFloat}}) = x
+AbstractFFTs.realfloat(x::StridedArray{<:Real}) = x
+# We override this one in order to avoid throwing an error that the type is
+# unsupported (as defined in AbstractFFTs)
+AbstractFFTs._fftfloat(::Type{T}) where {T <: AbstractFloat} = T
+
+plan_fft!(x::StridedArray{T}, region) where {T <: AbstractFloats} = DummyFFTPlan{Complex{real(T)},true}()
+plan_ifft!(x::StridedArray{T}, region) where {T <: AbstractFloats} = DummyiFFTPlan{Complex{real(T)},true}()
+
+# plan_rfft!(x::StridedArray{T}) where {T <: AbstractFloat} = DummyrFFTPlan{Complex{real(T)},true}()
+# plan_irfft!(x::StridedArray{T},n::Integer) where {T <: AbstractFloat} = DummyirFFTPlan{Complex{real(T)},true}()
+plan_dct!(x::StridedArray{T}, region) where {T <: AbstractFloats} = DummyDCTPlan{T,true}()
+plan_idct!(x::StridedArray{T}, region) where {T <: AbstractFloats} = DummyiDCTPlan{T,true}()
+
+plan_fft(x::StridedArray{T}, region) where {T <: AbstractFloats} = DummyFFTPlan{Complex{real(T)},false}()
+plan_ifft(x::StridedArray{T}, region) where {T <: AbstractFloats} = DummyiFFTPlan{Complex{real(T)},false}()
+plan_rfft(x::StridedArray{T}, region) where {T <: AbstractFloats} = DummyrFFTPlan{Complex{real(T)},false}(length(x))
+plan_irfft(x::StridedArray{T}, n::Integer, region) where {T <: AbstractFloats} = DummyirFFTPlan{Complex{real(T)},false}(n)
+plan_dct(x::StridedArray{T}, region) where {T <: AbstractFloats} = DummyDCTPlan{T,false}()
+plan_idct(x::StridedArray{T}, region) where {T <: AbstractFloats} = DummyiDCTPlan{T,false}()
 
 
 function interlace(a::Vector{S},b::Vector{V}) where {S<:Number,V<:Number}

test/fftBigFloattests.jl

@@ -14,7 +14,7 @@ end
     @test norm(ifft(fft(c))-c) < 200norm(c)eps(BigFloat)
 
     p = plan_dct(c)
-    @test norm(dct(c) - p*c) == 0
+    @test norm(FastTransforms.generic_dct(c) - p*c) == 0
 
     pi = plan_idct!(c)
     @test norm(pi*dct(c) - c) < 1000norm(c)*eps(BigFloat)