Add examples folder showing full transform from points to coefficients

[dlfivefifty]

I'm always confused which points to use for which transform, especially with the SphericalHarmonic/TriangleHarmonics case. I think an examples folder would help here.

If I figure out TriangleHarmonics, I'll add the example.

• [x] Chebyshev T and U
• [x] Jacobi (commit https://github.com/JuliaApproximation/FastTransforms.jl/commit/65ffcdb1fb27d2deffa5167bde79c518250e6a28)
• [x] Spherical harmonics
• [x] Triangle harmonics
• [ ] nfft
• [x] Padua

[MikaelSlevinsky]

An example that uses the nufft for a semi-Lagrangian time stepper would be neat too

[ajt60gaibb]

@danfortunato and I could happily write an example using the nufft for a semi-Lagrangian time stepper on the surface of the surface.

[pnavaro]

In order to compute advection, i used nufft to interpolate a function after a displacement of alpha = 0.5.

I made this code with FINUFFT.jl, do you know how to do it with FastTransforms ?

using Plots, FFTW, FINUFFT
"""
    interpolate_with_nufft_1d( x, f, alpha)

Interpolate the 1d function f after a displacement alpha

"""
function interpolate_with_finufft_1d!( f, lx, alpha)
    nj = length(f)
    f̂  = fft(f) ./ nj
    x̂  = range(0, stop=2π, length=nj+1)[1:end-1] |> collect
    x̂  .-= alpha/lx * 2π
    opts = finufft_default_opts()
    opts.modeord = 1
    tol = 1e-12
    nufft1d2!(x̂, f, 1, tol, f̂, opts)
end

nj = 64
f  = zeros(ComplexF64, nj)
x  = range(-2π, stop=2π, length=nj+1)[1:end-1] |> collect
dx = 4π / nj
f .= exp.(-x.^2 )
lx = 4π
alpha = 0.5
plot(x .+ alpha, real(f);  marker=:o)
interpolate_with_finufft_1d!( f, lx, alpha)
plot!(x, real(f);  marker=:o)

[MikaelSlevinsky]

Looking at the documentation, there are a couple differences in the mathematical definitions. In particular, for the type-II transform, this package is explicit about the "negative 2pi" in the frequency, and the summation starts at frequency 0 rather than at a negative integer. See https://finufft.readthedocs.io/en/latest/math.html vs. http://mikaelslevinsky.github.io/FastTransforms.jl/latest/

Otherwise, the API is slightly different (but probably subject to change).

nufft2(f, x, eps()) is sufficient for one-time use, but you can use plan_nufft2 and a mul! to separate pre-computation from execution.

[pnavaro]

I am very interested to use FastTransforms plans. Unfortunately i don't know how to reorder the coefficients computed with FFTW. Here i should get the same function. I'm lost on what computes fftw and FastTransforms.

using FastTransforms
n  = 64
f  = zeros(ComplexF64, n)
x  = range(-π, stop=π, length=nj+1)[1:end-1] |> collect
f .= exp.(-x.^2 )
plot(x, real(f);  marker=:o)
f̂  = fft(f) 
x̂ = (collect(0:n-1) .+ rand(n))./n
f .= inufft2(f̂, x̂, eps())
plot!(x̂ .* 2π .- π, real(f);  marker=:o)

[MikaelSlevinsky]

Hi Pierre,

Whenever I get lost in the NUFFT's, I find it's easiest to understand in my notebook first. Perhaps you may find the following comments and code helpful.

Although separating pre-computation into a plan and execution into a mul! is in principle computationally optimal, I haven't benchmarked this library versus the Flatiron's NUFFT, and it's possible that although the Julia wrappers to FINUFFT don't explicitly store pre-computation at the Julia level, some planning may be cached at the C level (provided that, e.g. parameters don't change between executions).

In this tutorial, we consider how to use the FFT in Julia to return Fourier coefficients to a function that is sampled at n equispaced points on [-1,1). That is:

    x_j = -1+2j/n,    for    j = 0,…,n-1.

For simplicity, we only consider n to be odd. Our goal is to recover the Fourier coefficients in the trigonometric interpolating polynomial:

    p_n(x) = Σ_{k=-n÷2}^{n÷2} c_k e^{i k π x}.

Note that in Julia, array indexing starts at 1, so c_{-n÷2} will actually be c[1].

If we sampled p_n at the equispaced points, we would find:

    p_n(x_j) = Σ_{k=-n÷2}^{n÷2} c_k e^{i k π (-1+2j/n)},

             = Σ_{k=-n÷2}^{n÷2} e^{2 π i j k / n} e^{-i k π} c_k,

             = Σ_{k=-n÷2}^{n÷2} e^{2 π i j k / n} (-1)^k c_k,

             = Σ_{k=0}^{n÷2} e^{2 π i j k / n} (-1)^k c_k

               + Σ_{k=n÷2}^{n-1} e^{2 π i j (k - n) / n} (-1)^{k-n} c_{k-n}.

In the last line, we changed the summation index by k^{new} = k^{old} + n. Then:

    p_n(x_j) = Σ_{k=0}^{n÷2} e^{2 π i j k / n} (-1)^k c_k

               + e^{-2 π i j} Σ_{k=n÷2}^{n-1} e^{2 π i j k / n} (-1)^{k-n} c_{k-n},

             = Σ_{k=0}^{n÷2} e^{2 π i j k / n} (-1)^k c_k

               + Σ_{k=n÷2}^{n-1} e^{2 π i j k / n} (-1)^{k-n} c_{k-n},

             = Σ_{k=0}^{n-1} e^{2 π i j k / n} d_k.

Now, we recognize that the evaluation of the trigonometric polynomial at equispaced points may be expressed as the conjugate transpose of the DFT.

Furthermore, the coefficients d_k are related to c_k by two simple steps. Firstly, there is the alternation in sign, and secondly, the first and second halves of the c_k coefficients are swapped. In summary, all the steps are invoked in reverse order, and FFTW has a convenient function fftshift for swapping the coefficients.

When n is even, we have to make a choice to keep either the highest positive Fourier mode or highest negative Fourier mode. Below, we implement the keeping of the highest negative Fourier mode.

using Base.MathConstants, FFTW

function vals2coeffs(f::Vector)
    n = length(f)
    c = fftshift(fft(f)/n)
    s = iseven(n÷2) ? 1 : -1
    @inbounds for k = 1:n
        c[k] = s*c[k]
        s *= -1
    end
    c
end

function coeffs2vals(c::Vector)
    n = length(c)
    s0 = iseven(n÷2) ? 1 : -1
    s = s0
    @inbounds for k = 1:n
        c[k] = s*c[k]
        s *= -1
    end
    f = n*ifft(ifftshift(c))
    s = s0
    @inbounds for k = 1:n
        c[k] = s*c[k]
        s *= -1
    end
    f
end

Consider the function f = x -> (√2 + e*im) * exp(-2*π*im*x) + (√3+π*im) * exp(-3*π*im*x) + (√5 + catalan*im) * exp(2*π*im*x) + (√6+eulergamma*im) * exp(3*π*im*x). If everything is correct, we should be able to exactly recover the Fourier coefficients of f with as few as n = 7 equispaced points.

f = x -> (√2 + e*im) * exp(-2*π*im*x) + (√3+π*im) * exp(-3*π*im*x) + (√5 + catalan*im) * exp(2*π*im*x) + (√6+eulergamma*im) * exp(3*π*im*x)

n = 7

x = -1 .+ 2*(0:n-1)/n

fvals = map(f, x)

c = vals2coeffs(fvals)

ftest = coeffs2vals(c)

norm(fvals-ftest)/norm(fvals)

Next, we show how to use the type-II NUFFT to evaluate the trigonometric interpolating polynomial at a different set of points. Again, using the trigonometric interpolating polynomial:

    p_n(x) = Σ_{k=-n÷2}^{n÷2} c_k e^{i k π x},

we reverse the summation index, resulting in:

    p_n(x) = Σ_{k=-n÷2}^{n÷2} c_{-k} e^{-i k π x},

and apart from the common factor of e^{π i x (n÷2)}, we may express evaluation at a set of points x_j ∈ [-1,1) by the type-II NUFFT as follows:

    p_n(x) = e^{π i x (n÷2)} Σ_{k=-n÷2}^{n÷2} c_{-k} e^{-i (k+n÷2) π x}

           = e^{π i x (n÷2)} Σ_{k=0}^{n-1} c_{n÷2-k} e^{-2 π i k (x/2)}.

           = e^{π i x (n÷2)} Σ_{k=0}^{n-1} [e^{π i k} c_{n÷2-k}] e^{-2 π i k [(x+1)/2]}.

If the points x ∈ [-1,1), then [(x+1)/2] ∈ [0,1), ready for the NUFFT.

using FastTransforms

function my_nufft_interp(cfs::Vector, x::Vector)
    n = length(cfs)
    c = cfs[end:-1:1]
    s = 1
    for k = 1:n
        c[k] *= s
        s *= -1
    end
    vals = nufft2(c, (x .+ 1) ./ 2, eps())
    for i = 1:n
        vals[i] *= exp(π*im*x[i]*(n÷2))
    end
    vals
end

y = 2 .* rand(n) .- 1 # -1 .+ 2*(0:n-1)/n .+ 1 ./ n # shifted points

fshiftedvals = map(f, y)

ffastshiftedvals = my_nufft_interp(c, y)

norm(fshiftedvals-ffastshiftedvals)/norm(fshiftedvals)

The inner workings of my_nufft_interp could be optimized with plans.

[pnavaro]

Thank you very much Mikael.