-
Notifications
You must be signed in to change notification settings - Fork 27
Open
Description
Dear all,
I am trying to implement a Python version of the algorithm from the paper A nonuniform fast Fourier transform based on low rank approximation.
I think nufft.jl from FastTransforms.jl uses the same algorithm.
I first did a comparison between nufft and its naive implementation for a given x
and uniform samples samples
.
As expected I get the same result for k=1
.
julia> x = complex([-0.1, 0.2, 0.3, 0.01, -1.0, 4.0])
6-element Vector{ComplexF64}:
-0.1 + 0.0im
0.2 + 0.0im
0.3 + 0.0im
0.01 + 0.0im
-1.0 + 0.0im
4.0 + 0.0im
julia> samples = collect(0:6-1) / 6
6-element Vector{Float64}:
0.0
0.16666666666666666
0.3333333333333333
0.5
0.6666666666666666
0.8333333333333334
julia> sum(x .* exp.((2.0 * im * pi * k) .* samples)) / sqrt(6)
0.9553009996854398 - 0.8838834764831844im
julia> f = getindex(nufft2(x, samples, eps()) / sqrt(6), k + 1)
0.9553009996854397 + 0.8838834764831844im
However, for nonuniform samples I get an error:
julia> samples = [0.05, 0.051, 0.3, 0.47, 0.65, 0.9]
6-element Vector{Float64}:
0.05
0.051
0.3
0.47
0.65
0.9
julia> sum(x .* exp.((2.0 * im * pi * k) .* samples)) / sqrt(6)
1.557891331279937 - 0.49922137697035407im
julia> f = getindex(nufft2(x, samples, eps()) / sqrt(6), k + 1)
-0.03047518500012294 - 1.3431665848287448im
Do I miss something or the error is expected ?
Do you have a way to estimate the error for each entry of the output ?
Thank you.
Best,
Pascal
Metadata
Metadata
Assignees
Labels
No labels