Nonuniform fast Fourier transforms


#1

FastTransforms.jl v0.2.2 includes nonuniform fast Fourier transforms thanks to Alex Townsend:

  • nufft1 assumes uniform samples and noninteger frequencies;
  • nufft2 assumes nonuniform samples and integer frequencies;
  • nufft3 ( = nufft) assumes nonuniform samples and noninteger frequencies;
  • inufft1 inverts an nufft1; and,
  • inufft2 inverts an nufft2.

Basic usage is like so:

julia> Pkg.add("FastTransforms")

julia> using FastTransforms

julia> n = 10^4;

julia> c = complex(rand(n));

julia> ω = collect(0:n-1) + rand(n);

julia> nufft1(c, ω, eps());

julia> p1 = plan_nufft1(ω, eps());

julia> @time p1*c;
  0.002383 seconds (6 allocations: 156.484 KiB)

The new algorithms are described in the following pre-print:

D. Ruiz—Antolín and A. Townsend. A nonuniform fast Fourier transform based on low rank approximation, arXiv:1701.04492, 2017.


#2

@MikaelSlevinsky: Are you aware of https://github.com/tknopp/NFFT.jl? It seems that nufft just allows for 1D NFFTs whereas NFFT.jl allows for d-dimensional transforms.

I am pretty interested in a performance comparison of both approaches.


#3

And from an interface point of view: Isn’t nufft2 the adjoint of nufft2 (modulo conjugation). The NFFT is quite often used as an operator when solving linear systems of equations and there you usually want that pair. In my (internal) software on MRI reconstruction we overload A_mul_B! and Ac_mul_B! for this purpose.


#4

Hi Tobias, there is some support for 2D nufft’s of types I-I and II-II, but the multi-dimensional support is not fully fleshed out. You are right, nufft1 and nufft2 are adjoints, modulo conjugation.

I’m interested in a performance comparison, too, and the first step is actually getting a release of this approach.

There are four levels to the 1D interface:

  • the functions nufft1/2/3(args...) plan & execute;
  • the functions plan_nufft1/2/3(args...) construct plans;
  • if p is an NUFFTPlan, then p*c applies the plan; and,
  • A_mul_B!(out, p, in) applies the plan to the in array in-place on the out array.

For a fair comparison, methods should be compared at the same tolerance ϵ.

Note that a hierarchical approach to the Chebyshev–Legendre transform accelerated the execution phase (of any previous methods in FastTransforms.jl) by a factor of ~20. With some effort, Dutt & Rokhlin’s hierarchical nonuniform FFT method could be implemented based on the HierarchicalMatrices.jl framework.


#5

I just saw that the paper actually makes a comparison with the NFFT.jl, which is great.

I wonder if it would be good to share an interface. NFFT.jl could have different implementations that could be switchable.


#6

Yes, in Figure 2.2 we did a quick comparison to NFFT.jl for the 1D NUFFT of type-I:


#7

Maybe a meta-package a la DifferentialEquations.jl is in order? Transforms.jl is available it looks like…there could even be an org JuliaTransforms.


#8

No, in this case I do not really see a need for that. It should be possible to have a single NFFT package that allows to switch different algorithms in its backend.
FastTransforms.jl is more high level and can then integrate the low level packages (where NFFT is just one example)


#9

I think that’s what I meant: a package that wraps the different algorithms for the different transforms, with a consistent naming scheme and consistent way to specify algorithms.


#10

Talking about a common interface maybe this package could be interesting: AbstractOperators.jl.

This offers a common interface that allows to easily combine linear and nonlinear transformations and use them in inverse problems (through the package StructuredOptimization).

In case you find this interesting feel free to contribute!