Skip to content

Commit c7934db

Browse files
committed
Define the nufft() command, which just calls nufft1()
1 parent 9a3ebf4 commit c7934db

File tree

2 files changed

+6
-21
lines changed

2 files changed

+6
-21
lines changed

src/FastTransforms.jl

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
__precompile__()
22
module FastTransforms
33

4-
using Base, ToeplitzMatrices, Compat
4+
using Base, ToeplitzMatrices, Compat, LambertW
55

66
import Base: *
77
import Compat: view
@@ -11,8 +11,8 @@ export leg2cheb, cheb2leg, leg2chebu, ultra2ultra, jac2jac
1111
export gaunt
1212
export paduatransform, ipaduatransform, paduatransform!, ipaduatransform!, paduapoints
1313
export plan_paduatransform!, plan_ipaduatransform!
14-
export nufft, nufft1, nufft2, nufft3
15-
export nufft_plan nufft1_plan, nufft2_plan, nufft3_plan
14+
export nufft, nufft1, nufft2
15+
export nufft_plan, nufft1_plan, nufft2_plan
1616

1717
# Other module methods and constants:
1818
#export ChebyshevJacobiPlan, jac2cheb, cheb2jac

src/nufft.jl

Lines changed: 3 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,3 @@
1-
21
function nufft1_plan{T<:AbstractFloat}( x::AbstractVector{T}, ϵ::T )
32
(s_vec, t_idx, gam) = FindAlgorithmicParameters( x )
43
K = FindK(gam, ϵ)
@@ -15,19 +14,19 @@ In = speye(eltype(c), N, N)
1514
p( c ) = (v.*(N*conj(ifft(In[:,t_idx]*conj(Diagonal(c)*u),1))))*ones(K)
1615
end
1716

17+
18+
nufft_plan{T<:AbstractFloat}( x::AbstractVectorP{T}, ϵ::T ) = nufft1_plan( x, ϵ )
19+
nufft{T<:AbstractFloat}( c::AbstractVector, x::AbstractVector{T}, ϵ::T ) = nufft_plan(x, ϵ)(c)
1820
nufft1{T<:AbstractFloat}( c::AbstractVector, x::AbstractVector{T}, ϵ::T ) = nufft1_plan(x, ϵ)(c)
1921
nufft2{T<:AbstractFloat}( c::AbstractVector, ω::AbstractVector{T}, ϵ::T ) = nufft2_plan(ω, ϵ)(c)
20-
nuftt3{T<:AbstractFloat}( c::AbstractVector, x::AbstractVector{T}, ω::AbstractVector{T}, ϵ::T ) = nufft3_plan(x, ω, ϵ)(c)
2122

2223
FindK{T<:AbstractFloat}::T, ϵ::T) = Int( ceil(5.0*γ.*exp(lambertw(log(10.0/ϵ)./γ/7.0))) )
2324

2425
function FindAlgorithmicParameters{T<:AbstractFloat}( x::AbstractVector{T} )
25-
2626
N = size(x, 1)
2727
s_vec = round(N*x)
2828
t_idx = mod(round(Int64, N*x), N) + 1
2929
γ = norm(N*x - s_vec, Inf)
30-
3130
return (s_vec, t_idx, γ)
3231
end
3332

@@ -39,14 +38,11 @@ N = size(x, 1)
3938
(s_vec, t_idx, γ) = FindAlgorithmicParameters( x )
4039
er = N*x - s_vec
4140
scl = exp( -1im*pi*er )
42-
4341
# Chebyshev polynomials of degree 0,...,K-1 evaluated at er/gam:
4442
TT = Diagonal(scl)*ChebyshevP(K-1, er/γ)
45-
4643
# Chebyshev expansion coefficients:
4744
cfs = Bessel_coeffs(K, γ)
4845
u = zeros(eltype(cfs), N, K)
49-
5046
# Construct them now:
5147
for r = 0:K-1
5248
for j = 1:N
@@ -56,19 +52,15 @@ for r = 0:K-1
5652
end
5753
end
5854
end
59-
6055
# rowspace vectors:
6156
v = ChebyshevP(K-1, 2.0*collect(0:N-1)/N - ones(N) )
62-
6357
return (u, v)
6458
end
6559

66-
6760
function Bessel_coeffs(K::Int64, γ::Float64)
6861
# Calculate the Chebyshev coefficients of exp(-2*pi*1im*x*y) on [-gam,gam]x[0,1]
6962

7063
cfs = complex(zeros( K, K ))
71-
7264
arg = -γ*pi/2
7365
for p = 0:K-1
7466
for q = mod(p,2):2:K-1
@@ -77,37 +69,30 @@ for p = 0:K-1
7769
end
7870
cfs[1,:] = cfs[1,:]/2.0
7971
cfs[:,1] = cfs[:,1]/2.0
80-
8172
return cfs
82-
8373
end
8474

85-
8675
function ChebyshevP{T<:AbstractFloat}(n::Int64, x::AbstractVector{T})
8776
# Evaluate Chebyshev polynomials of degree 0,...,n at x:
8877

8978
N = size(x,1)
9079
Tcheb = zeros(eltype(x), N, n+1)
91-
9280
# T_0(x) = 1.0
9381
for j = 1:N
9482
Tcheb[j, 1] = 1.0
9583
end
96-
9784
# T_1(x) = x
9885
for k = 2:min(n+1,2)
9986
for j = 1:N
10087
Tcheb[j, 2] = x[j]
10188
end
10289
end
103-
10490
# 3-term recurrence relation:
10591
twoX = 2.0*x
10692
for k = 2:n
10793
for j = 1:N
10894
Tcheb[j, k+1] = twoX[j]*Tcheb[j, k] - Tcheb[j, k-1]
10995
end
11096
end
111-
11297
return Tcheb
11398
end

0 commit comments

Comments
 (0)