NUFFT

Examples illustrating the nufft methods in the Julia package MIRT.

The nufft functions in this package are wrappers around NFFT.jl.

The plots on this page are simply sanity checks about of the approximation error for that package and they help verify the correctness of the wrapper.

This page comes from a single Julia file: 2-nufft.jl.

You can access the source code for such Julia documentation using the 'Edit on GitHub' link in the top right. You can view the corresponding notebook in nbviewer here: 2-nufft.ipynb, or open it in binder here: 2-nufft.ipynb.

Setup

Packages needed here.

using Plots; default(markerstrokecolor = :auto, label="")
using MIRTjim: prompt
using MIRT: nufft_init, nufft_errors, dtft
using InteractiveUtils: versioninfo

The following line is helpful when running this file as a script; this way it will prompt user to hit a key after each figure is displayed.

isinteractive() && prompt(:prompt);

1D example

This code illustrates how to call the NUFFT.

Ω = range(-1,1,301) * π # digital frequencies (radians/sample)
N = 16 # # of samples
(nufft, nufft_adjoint, Anufft) = nufft_init(Ω, N)
x = randn(ComplexF32, N) # random 1D signal
y1 = nufft(x) # one way to compute NUFFT
y2 = Anufft * x # another way to compute NUFFT
@assert y1 == y2
yd = dtft(Ω, x) # exact (slow) nonuniform discrete Fourier transform
maximum(abs, yd - y1) # worst error for this 1D signal

8.668276793545973e-7

This plot shows that the NUFFT is very close to the exact nonuniform DFT.

xticks = ([-π,0,π], ["-π", "0", "π"]); xlims = (-π,π); xlabel="Ω"
p1 = plot(ylabel="Real part"; xlabel, xlims, xticks)
plot!(Ω, real(y1), label="real(NUFFT)", line=:dash)
plot!(Ω, real(yd), label="real(DTFT)", line=:dot)
p2 = plot(ylabel="Imag part"; xlabel, xlims, xticks)
plot!(Ω, imag(y1), label="imag(NUFFT)", line=:dash)
plot!(Ω, imag(yd), label="imag(NUFFT)", line=:dot)
p3 = plot(Ω, real(yd - y1), label="real(Error)"; xlabel, xlims, xticks)
p4 = plot(Ω, imag(yd - y1), label="imag(Error)"; xlabel, xlims, xticks)
plot(p1, p2, p3, p4)

prompt()

Plot worst-case errors vs $N$

Plot worst-case error over all frequencies $Ω$ between $0$ and $2π/N$ for various $N$.

Even for $N=512$ the error is below 4e-6.

Nlist = 2 .^ (4:9)
Nlist = [Nlist..., 191, 385] # test odd N values
errfunn = N -> maximum(nufft_errors( ; N)[2])
elist = errfunn.(Nlist)
pn = scatter(Nlist, elist,
     xtick=Nlist, color=:red, xlabel="N", ylabel="error", yaxis=:log10)

prompt()

Plot error vs NFFT m.

For the default m = 4 the worst-case error is below 4e-6.

mlist = 3:7
errfunm = nfft_m -> maximum(nufft_errors( ; nfft_m)[2])
worst_m = errfunm.(mlist)
pm = scatter(mlist, worst_m,
    xlabel="m", ylabel="error", color=:magenta, yaxis=:log10)

prompt()

Plot error vs NFFT sigma.

For the default σ = 4 the worst-case error is below 4e-6.

slist = [1.5; 2:6] # σ
errfuns = nfft_sigma -> maximum(nufft_errors( ; nfft_sigma)[2])
worst_s = errfuns.(slist)
ps = scatter(slist, worst_s,
    xlabel="σ", ylabel="error", color=:blue, yaxis=:log10)

prompt()

Reproducibility

This page was generated with the following version of Julia:

using InteractiveUtils: versioninfo
io = IOBuffer(); versioninfo(io); split(String(take!(io)), '\n')

12-element Vector{SubString{String}}:
 "Julia Version 1.10.5"
 "Commit 6f3fdf7b362 (2024-08-27 14:19 UTC)"
 "Build Info:"
 "  Official https://julialang.org/ release"
 "Platform Info:"
 "  OS: Linux (x86_64-linux-gnu)"
 "  CPU: 4 × AMD EPYC 7763 64-Core Processor"
 "  WORD_SIZE: 64"
 "  LIBM: libopenlibm"
 "  LLVM: libLLVM-15.0.7 (ORCJIT, znver3)"
 "Threads: 1 default, 0 interactive, 1 GC (on 4 virtual cores)"
 ""

And with the following package versions

import Pkg; Pkg.status()

Status `~/work/MIRT.jl/MIRT.jl/docs/Project.toml`
  [6e4b80f9] BenchmarkTools v1.5.0
  [e30172f5] Documenter v1.7.0
  [d3d80556] LineSearches v7.3.0
  [98b081ad] Literate v2.19.0
  [7035ae7a] MIRT v0.18.2 `~/work/MIRT.jl/MIRT.jl`
  [170b2178] MIRTjim v0.25.0
  [91a5bcdd] Plots v1.40.8
  [b77e0a4c] InteractiveUtils
  [9a3f8284] Random

This page was generated using Literate.jl.