Skip to content

nufft precision #262

@PascalCarrivain

Description

@PascalCarrivain

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

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions