Comments (7)
@danfortunato and I could happily write an example using the nufft for a semi-Lagrangian time stepper on the surface of the surface.
from fasttransforms.jl.
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.
from fasttransforms.jl.
An example that uses the nufft for a semi-Lagrangian time stepper would be neat too
from fasttransforms.jl.
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)
from fasttransforms.jl.
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.
from fasttransforms.jl.
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)
from fasttransforms.jl.
Thank you very much Mikael.
from fasttransforms.jl.
Related Issues (20)
- Why are second kind Chebyshev points the first kind extrema?
- jac2cheb: Crash on M1 (with Rosetta)
- Remove DSP dependency HOT 4
- Should multiplication of `FTPlan` take `StridedArray`s as arguments instead of `Array`s? HOT 1
- Move Chebyshev functionality to Chebyshevtransforms.jl
- What algorithm does cheb2leg() use? HOT 2
- Potentially suboptimal performance in r2r transforms? HOT 1
- sph methods slow? HOT 4
- Chebyshev transforms do not work with the FFTW v1.6. HOT 2
- Add a slow path in ultra2ultra transforms?
- Reduction in accuracy in `cheb2leg ∘ leg2cheb` on v0.15 HOT 4
- Regression in `ultra2ultra` with identical orders HOT 1
- Write test for Normalization in `cheb2leg` and `leg2cheb` HOT 4
- `cheb2ultra` seems unreasonably slow HOT 2
- Julia rewrite of libFastTransforms HOT 2
- 1d transforms with multi-dimensional regions are slow HOT 3
- Use Toeplitz-dot-Hankel for very large dimensional transforms HOT 4
- Add multidimensional interface for lib transforms HOT 3
- Avoid reexporting dependencies?
- Move out forwardrecurrence and clenshaw
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from fasttransforms.jl.