nufft package
Core implementation of the non-uniform Fourier transform (NUFFT), based on [Fessler2003].
- nufft.compute_interp_coeffs(om, Nd, Jd, Kd, n_shift, varargin)
Compute the interpolation coefficients involved in the direct NUFFT operator (adapted from original code associated with [Fessler2003]).
- Parameters
om (double[:, d]) – Non-uniform coordinates of the target of the data points in the \(d\)-dimensional Fourier domain.
Nd (int[1, d]) – Input tensor size along each dimension.
Jd (int[1, d]) – Size of the interpolation kernel along each dimension.
Kd (int[1, d]) – Size of the Fourier domain along each dimension.
n_shift (int[1, d]) – Phase-shift along each dimension in Fourier space (expressed in number of samples in the Fourier domain).
varargin ({string} or {string, cell, cell}) – User-defined type of kernel, to be selected among the options
'minmax:kb'
(used by default),'minmax:tuned'
and'kaiser'
. Possible options are documented below.ktype (
'minmax:kb'
) – Minmax interpolator with excellent KB scaling ('minmax:kb'
is recommended, and used by default).ktype (
'minmax:tuned'
) – Minmax interpolator, somewhat numerically tuned.ktype (
'kaiser'
) – Kaiser-bessel (KB) interpolator (using default minmax best alpha, m).ktype, alpha, m (
'kaiser'
, double[:], double[:]) – Kaiser-bessel (KB) interpolator using the specified values for the parameters alpha and m).
- Returns
st – Structure containing the NUFFT scaling coefficients (
st.sn
) and Auxiliary parameters to build the NUFFT interpolation matrix (st.alpha
,st.beta
,st.uu
).- Return type
struct
- nufft.createG(om, J, N, K, nshift, varargin)
Create the gridding matrix
G
to perform the 2D NUFFT.- Parameters
om (double[M, 2]) – Normalized coordinates of the target of the data points in the \(d\)-dimensional Fourier domain (uv frequencies) (
M
data points).J (int[1, 2]) – Size of the interpolation kernel along each dimension.
N (int[1, 2]) – Image size along each dimension.
K (int[1, 2]) – Size of the Fourier domain along each dimension.
nshift (int[1, 2]) – Phase-shift along each dimension in Fourier space (expressed in number of samples in the Fourier domain).
ktype (string (varargin)) – Selected interpolator (among
'minmax:kb'
,'minmax:tuned'
and'kaiser'
, seenufft.compute_interp_coeffs_extended()
for associated kernel documentation). Recommended default is'minmax:kb'
.
- Returns
G (sparse complex[M, prod(K)]) – Sparse de-gridding interpolation matrix involed in the 2D NUFFT.
scale (double[:, :] [N]) – Scaling weights involed in the 2D NUFFT, of size
N
.
- nufft.createG_calib(D, om, N, K, J, S, nshift, varargin)
Create the gridding matrix
G
to perform the 2D NUFFT including DDE kernels.- Parameters
D (double[:, :]) – DDEs kernels for a single time instant
[S2, na]
.om (double[M, 2]) – Normalized coordinates of the target of the data points in the \(d\)-dimensional Fourier domain (uv frequencies).
N (int[1, 2]) – Image size along each dimension.
K (int[1, 2]) – Size of the Fourier domain along each dimension.
J (int) – Size of the square interpolation kernel along a single dimension.
S (int) – Size of the square DDE kernels (in the spatial Fourier domain) along a single dimension. In the following,
S2 = S^2
.nshift (int[1, 2]) – Phase-shift along each dimension in Fourier space (expressed in number of samples in the Fourier domain).
ktype (string (varargin)) – Selected interpolator (among
'minmax:kb'
,'minmax:tuned'
and'kaiser'
, seenufft.compute_interp_coeffs_extended()
for associated kernel documentation). Recommended default is'minmax:kb'
.
- Returns
G (sparse complex[M, prod(K)]) – Sparse de-gridding interpolation matrix involed in the 2D NUFFT.
scale (double[:, :] [N]) – Scaling weights involed in the 2D NUFFT, of size
N
.ll (int[:]) – Position of the nonzero values in each row of
G
.v (double[Q2, M]) – convolutions between the DDE kernels (
Q2 = (2 * S - 1)^2
).V (double[Q2, J2, M]) – values contained in
G
(J2 = J^2
).
- nufft.kaiser_bessel(x, J, alpha, kb_m, K_N)
Generalized Kaiser-Bessel function for
x
in support [-J/2,J/2] shape parameter “alpha” (default 2.34 J) order parameter “kb_m” (default 0) see (A1) in lewitt:90:mdi, JOSA-A, Oct. 1990- Parameters
x (double[M, 1]) – Arguments.
J (int[:]) – Number of neighbors.
alpha – Parameters of the interpolation kernel.
kb_m – Parameters of the interpolation kernel.
K_N – Parameters of the interpolation kernel.
- Returns
kb ([M, 1]) – KB function values, if
x
is numbers, orstring
forkernel(k,J)
, ifx
is astring
, or inline function, ifx
isinline
.alpha – Parameters of the interpolation kernel.
kb_m – Parameters of the interpolation kernel.
- nufft.kaiser_bessel_ft(u, J, alpha, kb_m, d)
Fourier transform of generalized Kaiser-Bessel function, in dimension
d
.- Parameters
u – frequency arguments or ‘string’ to return function string or ‘handle’ to return function handle
J (int) – Default 6.
alpha (double) – Shape parameter. Default
2.34 * J
.kb_m (int) – Order parameter. Default 0.
d (int) – Number of dimensions considered. Default 1.
- Returns
y – Transform values.
- Return type
complex[M, 1]
Note
Original code taken from [Fessler2003], available at https://github.com/JeffFessler/mirt. See (A3) in lewitt:90:mdi, JOSA-A, Oct. 1990.
- nufft.nufft_T(N, J, K, tol, alpha, beta, use_true_diric)
Precompute the matrix
T = [C' S S' C]\inv
used in NUFFT. This can be precomputed, being independent of frequency location.- Parameters
N (int[:]) – Signal length.
J (int[:]) – Number of neighbors.
K (int[:]) – FFT length.
tol (double) – Tolerance for smallest eigenvalue.
alpha (double[L+1, 1]) – Fourier coefficient vector for scaling.
beta (double[:]) – Scale
gamma = 2 * pi / K
by this for Fourier series.use_true_diric (bool) – Use true Diric function (default is to use sinc approximation).
- Returns
T – Precomputed matrix.
- Return type
double[J, J]
- nufft.nufft_alpha_kb_fit(N, J, K, varargin)
Return the alpha and beta corresponding to LS fit of L components to optimized Kaiser-Bessel scaling factors (
m=0
,alpha=2.34*J
). This is the best method known by J. Fessler currently for choosing alpha.- Parameters
N (int[:]) – Signal length.
J (int[:]) – Number of neighbors.
K (int[:]) – FFT length.
varargin (option, [1]) – Specify the
Nmid
midpoint:floor(N/2)
or default(N-1)/2
- Returns
alphas (complex[:]) – Parameters of the interpolation kernel.
beta (double[:]) – Parameters of the interpolation kernel.
Note
Original code taken from [Fessler2003], available at https://github.com/JeffFessler/mirt.
- nufft.nufft_coef
NUFFT interpolation coefficient vectors for a given kernel function.
- Parameters
om (double[M, 1]) – Digital frequency omega in radians.
J (int) – Number of neighbors used per frequency location.
K (_type_) – FFT size (should be >= N, the signal length).
kernel (function handle) – Kernel function handle.
- Returns
coef (double[J, M]) – Coef. vector for each frequency.
arg (double[J, M]) – Kernel argument.
Note
Original code taken from [Fessler2003], available at https://github.com/JeffFessler/mirt.
- nufft.nufft_diric(k, N, K, use_true_diric)
“Regular fourier” Dirichlet-function WITHOUT phase
nufft_diric(t) = sin(pi N t / K) / ( N * sin(pi t / K) ) \approx sinc(t / (K/N))
- Parameters
k (double[…]) – Sample locations (unitless real numbers)
N (int[…]) – Signal length.
K (int[…]) – FFT length.
use_true_diric (bool) – Use true Diric function (default is to use sinc approximation).
- Returns
f – Corresponding function values.
- Return type
double[…]
Warning
Caution: - matlab’s version is different:
sin(N * x / 2) / (N * sin(x / 2))
-nufft_diric()
isK
-periodic for oddN
but2K
-periodic for evenN
.Note
Original code taken from [Fessler2003], available at https://github.com/JeffFessler/mirt.
- nufft.nufft_offset(om, J, K)
Offset for NUFFT.
- Parameters
om (double[M, 1]) – omega (radians), typically in \([-\pi, \pi)\) (not essential!)
J (int[1, d]) – Number of neighbors used for NUFFT interpolation.
K (int[1, d]) – FFT size along each domension.
- Returns
k0 – Prepared for
mod(k0 + [1:J], K)
(+ 1
for matlab).- Return type
int[M,1]
Note
Original code taken from [Fessler2003], available at https://github.com/JeffFessler/mirt.
- nufft.nufft_r(om, N, J, K, alpha, beta, use_true_diric)
Make NUFFT
r
vector.- Parameters
om (double[M, d]) – Digital frequency omega in radians.
N (int[1, d]) – Signal length.
J (int[1, d]) – Number of neighbors used per frequency location.
K (int[1, d]) – FFT size (should be
> N
)alpha (complex[0:L]) – Fourier series coefficients of scaling factors.
beta (double[:]) – Scale
gamma=2pi./K
by this in Fourier series, typically isK./N
[Fessler2003] or 0.5 (Liu)use_true_diric (bool) – Optional flag for debugging purposes (defaults to false).
- Returns
rr (double[J, M]) –
r
vector for each frequencyarg (double[J, M]) – Dirac argument for
t=0
.
Note
Original code taken from [Fessler2003], available at https://github.com/JeffFessler/mirt.
- nufft.nufft_scale(Nd, Kd, alpha, beta, Nmid)
Compute scaling factors for NUFFT.
- Parameters
Nd (int[1, d]) – Input size (along each dimension).
Kd (int[1, d]) – Fourier size (along each dimension).
alpha (cell{d, 1}) – Parameters of the interpolation kernel.
beta (cell{d, 1}) – Parameters of the interpolation kernel.
Nmid (int[1, d]) – Midpoint:
floor(Nd/2)
or default(Nd-1)/2
.
- Returns
sn – Scaling factors.
- Return type
double[Nd]
Note
Original code taken from [Fessler2003], available at https://github.com/JeffFessler/mirt.
- nufft.nufft_sinc(x)
J. Fessler’s version of “sinc” function, because matlab’s
sinc()
is in a toolbox.- Parameters
x (double[:, …]) – Input array.
- Returns
y – Evaluation of
sinc(x)
.- Return type
double[:, …]
Note
Original code taken from [Fessler2003], available at https://github.com/JeffFessler/mirt.