add inufft1 and inufft2

Author: MikaelSlevinsky

README.md

@@ -79,8 +79,11 @@ then increment/decrement operators are used with linear complexity (and linear c
 
 The NUFFTs are implemented thanks to [Alex Townsend](https://github.com/ajt60gaibb):
  - `nufft1` assumes uniform samples and noninteger frequencies;
- - `nufft2` assumes nonuniform samples and integer frequencies; and,
- - `nufft3 ( = nufft)` assumes nonuniform samples and noninteger frequencies.
+ - `nufft2` assumes nonuniform samples and integer frequencies;
+ - `nufft3 ( = nufft)` assumes nonuniform samples and noninteger frequencies;
+ - `inufft1` inverts an `nufft1`; and,
+ - `inufft2` inverts an `nufft2`.
+Here is an example:
 ```julia
 julia> n = 10^4;
 

docs/src/index.md

@@ -64,6 +64,14 @@ nufft2
 nufft3
 ```
 
+```@docs
+inufft1
+```
+
+```@docs
+inufft2
+```
+
 ```@docs
 plan_nufft1
 ```
@@ -76,6 +84,14 @@ plan_nufft2
 plan_nufft3
 ```
 
+```@docs
+plan_inufft1
+```
+
+```@docs
+plan_inufft2
+```
+
 ```@docs
 paduatransform
 ```

src/FastTransforms.jl

@@ -19,12 +19,12 @@ export plan_normleg12cheb2, plan_cheb22normleg1
 
 export gaunt
 
+export nufft, nufft1, nufft2, nufft3, inufft1, inufft2
+export plan_nufft, plan_nufft1, plan_nufft2, plan_nufft3, plan_inufft1, plan_inufft2
+
 export paduatransform, ipaduatransform, paduatransform!, ipaduatransform!, paduapoints
 export plan_paduatransform!, plan_ipaduatransform!
 
-export nufft, nufft1, nufft2, nufft3
-export plan_nufft, plan_nufft1, plan_nufft2, plan_nufft3
-
 export SlowSphericalHarmonicPlan, FastSphericalHarmonicPlan, ThinSphericalHarmonicPlan
 export sph2fourier, fourier2sph, plan_sph2fourier
 export sphones, sphzeros, sphrand, sphrandn, sphevaluate
@@ -55,6 +55,7 @@ include("ChebyshevUltrasphericalPlan.jl")
 include("ultra2cheb.jl")
 include("cheb2ultra.jl")
 include("nufft.jl")
+include("inufft.jl")
 
 include("cjt.jl")
 

src/inufft.jl

@@ -0,0 +1,104 @@
+doc"""
+Pre-computes an inverse nonuniform fast Fourier transform of type `N`.
+
+For best performance, choose the right number of threads by `FFTW.set_num_threads(4)`, for example.
+"""
+immutable iNUFFTPlan{N,T,S,PT} <: Base.DFT.Plan{T}
+    pt::PT
+    TP::Toeplitz{T}
+    ϵ::S
+end
+
+doc"""
+Pre-computes an inverse nonuniform fast Fourier transform of type I.
+"""
+function plan_inufft1{T<:AbstractFloat}(ω::AbstractVector{T}, ϵ::T)
+    N = length(ω)
+    p = plan_nufft1(ω, ϵ)
+    pt = plan_nufft2(ω/N, ϵ)
+    c = p*ones(Complex{T}, N)
+    r = conj(c)
+    avg = (r[1]+c[1])/2
+    r[1] = avg
+    c[1] = avg
+    TP = Toeplitz(c, r)
+
+    iNUFFTPlan{1, eltype(TP), typeof(ϵ), typeof(pt)}(pt, TP, ϵ)
+end
+
+doc"""
+Pre-computes an inverse nonuniform fast Fourier transform of type II.
+"""
+function plan_inufft2{T<:AbstractFloat}(x::AbstractVector{T}, ϵ::T)
+    N = length(x)
+    pt = plan_nufft1(N*x, ϵ)
+    r = pt*ones(Complex{T}, N)
+    c = conj(r)
+    avg = (r[1]+c[1])/2
+    r[1] = avg
+    c[1] = avg
+    TP = Toeplitz(c, r)
+
+    iNUFFTPlan{2, eltype(TP), typeof(ϵ), typeof(pt)}(pt, TP, ϵ)
+end
+
+function (*){N,T,V}(p::iNUFFTPlan{N,T}, x::AbstractVector{V})
+    A_mul_B!(zeros(promote_type(T,V), length(x)), p, x)
+end
+
+function Base.A_mul_B!{T}(c::AbstractVector{T}, P::iNUFFTPlan{1,T}, f::AbstractVector{T})
+    pt, TP, ϵ = P.pt, P.TP, P.ϵ
+    cg(TP, c, f, 50, 100ϵ)
+    conj!(A_mul_B!(c, pt, conj!(c)))
+end
+
+function Base.A_mul_B!{T}(c::AbstractVector{T}, P::iNUFFTPlan{2,T}, f::AbstractVector{T})
+    pt, TP, ϵ = P.pt, P.TP, P.ϵ
+    cg(TP, c, conj!(pt*conj!(f)), 50, 100ϵ)
+    conj!(f)
+    c
+end
+
+doc"""
+Computes an inverse nonuniform fast Fourier transform of type I.
+"""
+inufft1{T<:AbstractFloat}(c::AbstractVector, ω::AbstractVector{T}, ϵ::T) = plan_inufft1(ω, ϵ)*c
+
+doc"""
+Computes an inverse nonuniform fast Fourier transform of type II.
+"""
+inufft2{T<:AbstractFloat}(c::AbstractVector, x::AbstractVector{T}, ϵ::T) = plan_inufft2(x, ϵ)*c
+
+function cg{T}(A::ToeplitzMatrices.AbstractToeplitz{T}, x::AbstractVector{T}, b::AbstractVector{T}, max_it::Integer, tol::Real)
+	n = length(b)
+	n1, n2 = size(A)
+	n == n1 == n2 || throw(DimensionMismatch(""))
+    nrmb = norm(b)
+    if nrmb == 0 nrmb = one(typeof(nrmb)) end
+	copy!(x, b)
+    r = zero(x)
+    p = zero(x)
+    Ap = zero(x)
+    # r = b - A*x
+    copy!(r, b)
+    A_mul_B!(-one(T), A, x, one(T), r)
+	copy!(p, r)
+	nrm2 = r⋅r
+    for k = 1:max_it
+        # Ap = A*p
+        A_mul_B!(one(T), A, p, zero(T), Ap)
+		α = nrm2/(p⋅Ap)
+        @inbounds @simd for l = 1:n
+            x[l] += α*p[l]
+            r[l] -= α*Ap[l]
+        end
+		nrm2new = r⋅r
+        cst = nrm2new/nrm2
+        @inbounds @simd for l = 1:n
+            p[l] = muladd(cst, p[l], r[l])
+        end
+		nrm2 = nrm2new
+        if sqrt(abs(nrm2)) ≤ tol*nrmb break end
+	end
+    return x
+end

src/nufft.jl

@@ -1,5 +1,5 @@
 doc"""
-Pre-compute a nonuniform fast Fourier transform of type `N`.
+Pre-computes a nonuniform fast Fourier transform of type `N`.
 
 For best performance, choose the right number of threads by `FFTW.set_num_threads(4)`, for example.
 """
@@ -14,7 +14,7 @@ immutable NUFFTPlan{N,T,FFT} <: Base.DFT.Plan{T}
 end
 
 doc"""
-Computes a nonuniform fast Fourier transform of type I:
+Pre-computes a nonuniform fast Fourier transform of type I:
 
 ```math
 f_j = \sum_{k=1}^N c_k e^{-2\pi{\rm i} (j-1)/N \omega_k},\quad{\rm for}\quad 1 \le j \le N.
@@ -37,7 +37,7 @@ function plan_nufft1{T<:AbstractFloat}(ω::AbstractVector{T}, ϵ::T)
 end
 
 doc"""
-Computes a nonuniform fast Fourier transform of type II:
+Pre-computes a nonuniform fast Fourier transform of type II:
 
 ```math
 f_j = \sum_{k=1}^N c_k e^{-2\pi{\rm i} x_j (k-1)},\quad{\rm for}\quad 1 \le j \le N.
@@ -59,7 +59,6 @@ function plan_nufft2{T<:AbstractFloat}(x::AbstractVector{T}, ϵ::T)
 end
 
 doc"""
-Computes a nonuniform fast Fourier transform of type III:
 
 ```math
 f_j = \sum_{k=1}^N c_k e^{-2\pi{\rm i} x_j \omega_k},\quad{\rm for}\quad 1 \le j \le N.
@@ -93,21 +92,19 @@ function (*){N,T,V}(p::NUFFTPlan{N,T}, x::AbstractVector{V})
     A_mul_B!(zeros(promote_type(T,V), length(x)), p, x)
 end
 
-function Base.A_mul_B!{T}(y::AbstractVector{T}, P::NUFFTPlan{1,T}, c::AbstractVector{T})
+function Base.A_mul_B!{T}(f::AbstractVector{T}, P::NUFFTPlan{1,T}, c::AbstractVector{T})
     U, V, p, t, temp, temp2, Ones = P.U, P.V, P.p, P.t, P.temp, P.temp2, P.Ones
 
-    # (V.*(N*conj(ifft(In[:,t]*conj(Diagonal(c)*U),1))))*ones(K)
-
     broadcast!(*, temp, c, U)
     conj!(temp)
     fill!(temp2, zero(T))
     recombine_rows!(temp, t, temp2)
     p*temp2
     conj!(temp2)
     broadcast!(*, temp, V, temp2)
-    A_mul_B!(y, temp, Ones)
+    A_mul_B!(f, temp, Ones)
 
-    y
+    f
 end
 
 function Base.A_mul_B!{T}(Y::Matrix{T}, P::NUFFTPlan{1,T}, C::Matrix{T})
@@ -117,7 +114,7 @@ function Base.A_mul_B!{T}(Y::Matrix{T}, P::NUFFTPlan{1,T}, C::Matrix{T})
     Y
 end
 
-function A_mul_B_col_J!{T}(Y::Matrix{T}, P::NUFFTPlan{1,T}, C::Matrix{T}, J::Int)
+function A_mul_B_col_J!{T}(F::Matrix{T}, P::NUFFTPlan{1,T}, C::Matrix{T}, J::Int)
     U, V, p, t, temp, temp2, Ones = P.U, P.V, P.p, P.t, P.temp, P.temp2, P.Ones
 
     broadcast_col_J!(*, temp, C, U, J)
@@ -128,9 +125,9 @@ function A_mul_B_col_J!{T}(Y::Matrix{T}, P::NUFFTPlan{1,T}, C::Matrix{T}, J::Int
     conj!(temp2)
     broadcast!(*, temp, V, temp2)
     COLSHIFT = size(C, 1)*(J-1)
-    A_mul_B!(Y, temp, Ones, 1+COLSHIFT, 1)
+    A_mul_B!(F, temp, Ones, 1+COLSHIFT, 1)
 
-    Y
+    F
 end
 
 function broadcast_col_J!(f, temp::Matrix, C::Matrix, U::Matrix, J::Int)
@@ -144,30 +141,28 @@ function broadcast_col_J!(f, temp::Matrix, C::Matrix, U::Matrix, J::Int)
     temp
 end
 
-function Base.A_mul_B!{T}(y::AbstractVector{T}, P::NUFFTPlan{2,T}, c::AbstractVector{T})
+function Base.A_mul_B!{T}(f::AbstractVector{T}, P::NUFFTPlan{2,T}, c::AbstractVector{T})
     U, V, p, t, temp, temp2, Ones = P.U, P.V, P.p, P.t, P.temp, P.temp2, P.Ones
 
-    # (U.*(fft(Diagonal(c)*V,1)[t+1,:]))*ones(K)
-
     broadcast!(*, temp, c, V)
     p*temp
     reindex_temp!(temp, t, temp2)
     broadcast!(*, temp, U, temp2)
-    A_mul_B!(y, temp, Ones)
+    A_mul_B!(f, temp, Ones)
 
-    y
+    f
 end
 
-function Base.A_mul_B!{T}(y::AbstractVector{T}, P::NUFFTPlan{3,T}, c::AbstractVector{T})
+function Base.A_mul_B!{T}(f::AbstractVector{T}, P::NUFFTPlan{3,T}, c::AbstractVector{T})
     U, V, p, t, temp, temp2, Ones = P.U, P.V, P.p, P.t, P.temp, P.temp2, P.Ones
 
     broadcast!(*, temp2, c, V)
     A_mul_B!(temp, p, temp2)
     reindex_temp!(temp, t, temp2)
     broadcast!(*, temp, U, temp2)
-    A_mul_B!(y, temp, Ones)
+    A_mul_B!(f, temp, Ones)
 
-    y
+    f
 end
 
 function reindex_temp!{T}(temp::Matrix{T}, t::Vector{Int}, temp2::Matrix{T})
@@ -189,17 +184,17 @@ function recombine_rows!{T}(temp::Matrix{T}, t::Vector{Int}, temp2::Matrix{T})
 end
 
 doc"""
-Pre-compute a nonuniform fast Fourier transform of type I.
+Computes a nonuniform fast Fourier transform of type I.
 """
 nufft1{T<:AbstractFloat}(c::AbstractVector, ω::AbstractVector{T}, ϵ::T) = plan_nufft1(ω, ϵ)*c
 
 doc"""
-Pre-compute a nonuniform fast Fourier transform of type II.
+Computes a nonuniform fast Fourier transform of type II.
 """
 nufft2{T<:AbstractFloat}(c::AbstractVector, x::AbstractVector{T}, ϵ::T) = plan_nufft2(x, ϵ)*c
 
 doc"""
-Pre-compute a nonuniform fast Fourier transform of type III.
+Computes a nonuniform fast Fourier transform of type III.
 """
 nufft3{T<:AbstractFloat}(c::AbstractVector, x::AbstractVector{T}, ω::AbstractVector{T}, ϵ::T) = plan_nufft3(x, ω, ϵ)*c
 

test/nuffttests.jl

@@ -45,19 +45,26 @@ using FastTransforms, Base.Test
 
     for n in N, ϵ in (1e-4,1e-8,1e-12,eps(Float64))
         c = complex(rand(n))
+        err_bnd = 500*ϵ*n*norm(c)
 
-        ω = collect(0:n-1) + rand(n)
+        ω = collect(0:n-1) + 0.25*rand(n)
         exact = nudft1(c, ω)
         fast = nufft1(c, ω, ϵ)
-        @test norm(exact - fast, Inf) < 500*ϵ*n*norm(c)
+        @test norm(exact - fast, Inf) < err_bnd
 
-        x = (collect(0:n-1) + 3rand(n))/n
+        d = inufft1(fast, ω, ϵ)
+        @test norm(c - d, Inf) < err_bnd
+
+        x = (collect(0:n-1) + 0.25*rand(n))/n
         exact = nudft2(c, x)
         fast = nufft2(c, x, ϵ)
-        @test norm(exact - fast, Inf) < 500*ϵ*n*norm(c)
+        @test norm(exact - fast, Inf) < err_bnd
+
+        d = inufft2(fast, x, ϵ)
+        @test norm(c - d, Inf) < err_bnd
 
         exact = nudft3(c, x, ω)
         fast = nufft3(c, x, ω, ϵ)
-        @test norm(exact - fast, Inf) < 500*ϵ*n*norm(c)
+        @test norm(exact - fast, Inf) < err_bnd
     end
 end