Merge pull request #59 from daanhb/master

Fix idct for complex-valued bigfloat vectors

Author: MikaelSlevinsky

src/fftBigFloat.jl

@@ -110,7 +110,7 @@ function ifft_pow2(x::Vector{Complex{T}}) where T<:BigFloat
 end
 
 
-function dct(a::AbstractArray{Complex{BigFloat}})
+function dct(a::AbstractVector{Complex{BigFloat}})
 	N = big(length(a))
     c = fft([a; flipdim(a,1)])
     d = c[1:N]
@@ -121,14 +121,14 @@ end
 
 dct(a::AbstractArray{BigFloat}) = real(dct(complex(a)))
 
-function idct(a::AbstractArray{Complex{BigFloat}})
-    N = big(length(a))
+function idct(a::AbstractVector{Complex{BigFloat}})
+	N = big(length(a))
     b = a * sqrt(2*N)
     b[1] = b[1] * sqrt(big(2))
-    b .*= exp.((im*big(pi)).*(0:N-1)./(2*N))
-    b = [b; 0; conj(flipdim(b[2:end],1))]
+    shift = exp.(-im * 2 * big(pi) * (N - big(1)/2) * (0:(2N-1)) / (2 * N))
+    b = [b; 0; -flipdim(b[2:end],1)] .* shift
     c = ifft(b)
-    c[1:N]
+    flipdim(c[1:N],1)
 end
 
 idct(a::AbstractArray{BigFloat}) = real(idct(complex(a)))

test/fftBigFloattests.jl

@@ -17,10 +17,16 @@ end
     @test norm(dct(c) - p*c) == 0
 
     pi = plan_idct!(c)
-    @test norm(pi*dct(c) - c) < 500norm(c)*eps(BigFloat)
+    @test norm(pi*dct(c) - c) < 1000norm(c)*eps(BigFloat)
 
     @test norm(dct(c)-dct(map(Float64,c)),Inf) < 10eps()
 
     cc = cis.(c)
     @test norm(dct(cc)-dct(map(Complex{Float64},cc)),Inf) < 10eps()
+
+    c = big.(rand(100)) + im*big.(rand(100))
+    @test norm(dct(c)-dct(map(ComplexF64,c)),Inf) < 10eps()
+    @test norm(idct(c)-idct(map(ComplexF64,c)),Inf) < 10eps()
+    @test norm(idct(dct(c))-c,Inf) < 1000eps(BigFloat)
+    @test norm(dct(idct(c))-c,Inf) < 1000eps(BigFloat)
 end