fixed bug in IRDFT

Author: nantonel

src/linearoperators/IRDFT.jl

@@ -42,8 +42,8 @@ function IRDFT(x::AbstractArray{Complex{T},N}, d::Int, dims::Int=1) where {T<:Nu
 	dim_out = ()
 	idx = ()
 	for i = 1:N
-		dim_out = i == dims ? (dim_out...,          d) : (dim_out...,dim_in[i]   )
-		idx     = i == dims ? (idx...    ,2:dim_in[i]) : (idx...    ,Colon()     )
+		dim_out = i == dims ? (dim_out...,               d) : (dim_out...,dim_in[i]   )
+        idx     = i == dims ? (idx...    , 2:ceil(Int,d/2)) : (idx...    ,Colon()     )
 	end
 	At = plan_rfft(zeros(dim_out),dims)
 	IRDFT{T,N,dims,typeof(A),typeof(At),typeof(idx)}(dim_in,dim_out,A,At,idx)

src/utilities/block.jl

@@ -36,7 +36,7 @@ blockeltype(x::AbstractArray) = eltype(x)
 blocklength(x::Tuple) = sum(blocklength.(x))
 blocklength(x::AbstractArray) = length(x)
 
-blockvecnorm(x::Tuple) = sqrt(blockvecdot(x, x))
+blockvecnorm(x::Tuple) = sqrt(real(blockvecdot(x, x)))
 blockvecnorm(x::AbstractArray{R}) where {R <: Number} = vecnorm(x)
 
 blockmaxabs(x::Tuple) = maximum(blockmaxabs.(x))
@@ -55,8 +55,7 @@ blockset!(y::Tuple, x) = blockset!.(y, x)
 blockset!(y::AbstractArray, x) = (y .= x)
 
 blockvecdot(x::T1, y::T2) where {T1 <: Tuple, T2 <: Tuple} = sum(blockvecdot.(x,y))
-blockvecdot(x::AbstractArray{R1}, y::AbstractArray{R2}) where {R1 <: Number, R2 <: Number} = real(vecdot(x, y))
-# inner product must be always real see section 4.2 of TFOCS manual
+blockvecdot(x::AbstractArray{R1}, y::AbstractArray{R2}) where {R1 <: Number, R2 <: Number} = vecdot(x, y)
 
 blockzeros(t::Tuple, s::Tuple) = blockzeros.(t, s)
 blockzeros(t::Type, n::NTuple{N, Integer} where {N}) = zeros(t, n)

test/test_block.jl

@@ -48,8 +48,8 @@ blockset!(yb,xb)
 
 z = blockvecdot(x,x2)
 zb = blockvecdot(xb,x2b)
-@test z == real(vecdot(x,x2))
-@test zb == real(vecdot(xb[1],x2b[1])+vecdot(xb[2],x2b[2])+vecdot(xb[3],x2b[3]))
+@test z == vecdot(x,x2)
+@test zb == vecdot(xb[1],x2b[1])+vecdot(xb[2],x2b[2])+vecdot(xb[3],x2b[3])
 
 y = blockzeros(x)
 yb = blockzeros(xb)

test/test_linear_operators.jl

@@ -207,19 +207,43 @@ op = RDFT(n,n)
 @test is_full_column_rank(op) == false
 
 ####### IRDFT ############
-n = 10
-op = IRDFT(Complex{Float64},(n,),19)
-x1 = rfft(randn(19))
-y1 = test_op(op, x1,irfft(randn(n),19), verb)
-y2 = irfft(x1,19)
+n = 12
+op = IRDFT(Complex{Float64},(div(n,2)+1,),n)
+x1 = rfft(randn(n))
+y1 = test_op(op, x1,irfft(randn(div(n,2)+1),n), verb)
+y2 = irfft(x1,n)
 
 @test all(vecnorm.(y1 .- y2) .<= 1e-12)
 
-n,m,l = 4,10,5
-op = IRDFT(Complex{Float64},(n,m,l),19,2)
-x1 = rfft(randn(n,19,l),2)
-y1 = test_op(op, x1, irfft(x1,19,2), verb)
-y2 = irfft(x1,19,2)
+n = 11
+op = IRDFT(Complex{Float64},(div(n,2)+1,),n)
+x1 = rfft(randn(n))
+y1 = test_op(op, x1,irfft(randn(div(n,2)+1),n), verb)
+y2 = irfft(x1,n)
+
+@test all(vecnorm.(y1 .- y2) .<= 1e-12)
+
+n,m,l = 4,19,5
+op = IRDFT(Complex{Float64},(n,div(m,2)+1,l),m,2)
+x1 = rfft(randn(n,m,l),2)
+y1 = test_op(op, x1, irfft(x1,m,2), verb)
+y2 = irfft(x1,m,2)
+
+@test all(vecnorm.(y1 .- y2) .<= 1e-12)
+
+n,m,l = 4,18,5
+op = IRDFT(Complex{Float64},(n,div(m,2)+1,l),m,2)
+x1 = rfft(randn(n,m,l),2)
+y1 = test_op(op, x1, irfft(x1,m,2), verb)
+y2 = irfft(x1,m,2)
+
+@test all(vecnorm.(y1 .- y2) .<= 1e-12)
+
+n,m,l = 5,18,5
+op = IRDFT(Complex{Float64},(div(n,2)+1,m,l),n,1)
+x1 = rfft(randn(n,m,l),1)
+y1 = test_op(op, x1, irfft(x1,n,1), verb)
+y2 = irfft(x1,n,1)
 
 @test all(vecnorm.(y1 .- y2) .<= 1e-12)