33# ---------------------------------------------------------
44
55#=
6- Examples illustrating the `nufft` option in the Julia package
7- [`MIRTjim`](https://github.com/JeffFessler/MIRTjim.jl).
6+ Examples illustrating the `nufft` methods
7+ in the Julia package
8+ [`MIRT`](https://github.com/JeffFessler/MIRT.jl).
89
910The `nufft` functions in this package
1011are wrappers around
@@ -38,7 +39,7 @@ This page was generated from a single Julia file:
3839
3940using Plots; default (markerstrokecolor = :auto , label= " " )
4041using MIRTjim: prompt
41- using MIRT: nufft_errors
42+ using MIRT: nufft_init, nufft_errors, dtft
4243using InteractiveUtils: versioninfo
4344
4445# The following line is helpful when running this file as a script;
@@ -47,6 +48,36 @@ using InteractiveUtils: versioninfo
4748isinteractive () && prompt (:prompt );
4849
4950
51+ #=
52+ ## 1D example
53+
54+ This code illustrates how to call the NUFFT.
55+ =#
56+ Ω = range (- 1 ,1 ,301 ) * π # digital frequencies (radians/sample)
57+ N = 16 # # of samples
58+ (nufft, nufft_adjoint, Anufft) = nufft_init (Ω, N)
59+ x = randn (ComplexF32, N) # random 1D signal
60+ y1 = nufft (x) # one way to compute NUFFT
61+ y2 = Anufft * x # another way to compute NUFFT
62+ @assert y1 == y2
63+ yd = dtft (Ω, x) # exact (slow) nonuniform discrete Fourier transform
64+ maximum (abs, yd - y1) # worst error for this 1D signal
65+
66+ # This plot shows that the NUFFT is very close to the exact nonuniform DFT.
67+ xticks = ([- π,0 ,π], [" -π" , " 0" , " π" ]); xlims = (- π,π); xlabel= " Ω"
68+ p1 = plot (ylabel= " Real part" ; xlabel, xlims, xticks)
69+ plot! (Ω, real (y1), label= " real(NUFFT)" , line= :dash )
70+ plot! (Ω, real (yd), label= " real(DTFT)" , line= :dot )
71+ p2 = plot (ylabel= " Imag part" ; xlabel, xlims, xticks)
72+ plot! (Ω, imag (y1), label= " imag(NUFFT)" , line= :dash )
73+ plot! (Ω, imag (yd), label= " imag(NUFFT)" , line= :dot )
74+ p3 = plot (Ω, real (yd - y1), label= " real(Error)" ; xlabel, xlims, xticks)
75+ p4 = plot (Ω, imag (yd - y1), label= " imag(Error)" ; xlabel, xlims, xticks)
76+ plot (p1, p2, p3, p4)
77+
78+ #
79+ prompt ()
80+
5081#=
5182## Plot worst-case errors vs ``N``
5283
0 commit comments