Skip to content

Genericfft #64

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Apr 12, 2019
Merged

Genericfft #64

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
c00b38b
make fft more generic
daanhb Apr 8, 2019
0752928
Rename all fft routines to generic_fft
daanhb Apr 11, 2019
bfec169
Improve rfft and irfft
daanhb Apr 11, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
170 changes: 102 additions & 68 deletions src/fftBigFloat.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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}
Expand Down
2 changes: 1 addition & 1 deletion test/fftBigFloattests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down