Examples of FINUFFT?

[erny123]

I’m trying to learn some basic Non-Uniform FFT.

I want to see if I can get a frequency spectrum of a laser pulse with a NUFFT type 1 using the FINUFFT.jl library.

Here’s the example:

using FFTW
using FINUFFT

c = 2.998e8; #m/s speed of light
wl = 525.0e-9; #m wavelength 
pl = wl/c; #period
wlfreq = c/wl; #2pi*Hz frequency

fwhm = 500.0e-15; #temporal fwhm
tau0 = fwhm/sqrt(2*log(2)); #gaussian fwhm
tlbw = 0.44/fwhm #transform limited bandwidth
	
maxf = c/wl +tlbw
minf = c/wl-tlbw
df = (maxf - minf)/1.0e3
	
	
#simulation
dt = pl/5.0;
t0 = -2.0*fwhm
tend = 2.0*fwhm
ttot = ceil(Int64,(tend-t0)/dt)

E = [exp(-t^2/tau0^2)*cos(2.0*pi*wlfreq*t) for t in t0:dt:tend];
t = [i for i in t0:dt:tend];
lent = length(E)
Nt = length(t);

The pulse looks like:

Taking the fft:

freq = FFTW.fftshift(FFTW.fftfreq(lent,1/dt))
Ef = FFTW.fftshift(FFTW.fft(E));
If = abs2.(Ef); #intensity spectrum

Looks like:

Which is exactly what I expect.

Now I want to us a NUFFT type 1 according to:
https://finufft.readthedocs.io/en/latest/math.html

I decide to take nonuniform points in the form of a tangent curve so as to get most of the FWHM of the pulse:

#nonuniform indices to extract sampling points
#using tangent curve
maxpi = pi*22.0/100.0;
dpi = 2.0*maxpi/(lent/2);
l = [i for i in -maxpi:dpi:maxpi];
tt = @. tan(1.95*l);
ttt = @. -tt[1]+tt;
ttt = @. ttt*lent/8.81+1;
t4 = @. ceil(Int64,ttt);

#extract nonuniform points
Enu = convert.(ComplexF64,E[t4]); #field at nonuniform points
tnu = t[t4]; #time at noniniform points
lentnu = length(Enu)

This looks like:

wmax = freq[end]; #max freq get from previous FFTW, around 1.5e15
Nw = length(freq); #number of total points|
dw = 2*wmax/Nw; # spacing of target k grid|
w = dw .* [i for i in ceil(Int64,-Nw/2):ceil(Int64,Nw/2-1)]; #    % a particular uniform M-grid of this spacing|
#normalize
Enunorm = Enu ./ maximum(E)

fhat = FINUFFT.nufft1d1(tnu, Enunorm, 1, 1.0e-12, Nw);   #% type 1, requesting M modes

And the intensity spectrum looks like:

Not sure what to make of it. The amplitudes change but by an extremely small bit?
The main pulse frequency is off by more than 2.

[rafael.guerra]

Fyi, the code posted doesn’t run: what are pl and wlfreq?

[erny123]

Oops sorry, edited my original post

[erny123]

Ok, looks like I still need more sampling points to get a better NUFFT spectrum.

For example increasing the sampling frequency to 1/10 the period of the wavelength rather than 1/5 in the original signal:

#simulation
dt = pl/10.0;
t0 = -2.0*fwhm
tend = 2.0*fwhm
ttot = ceil(Int64,(tend-t0)/dt)

E = [exp(-t^2/tau0^2)*cos(2.0*pi*wlfreq*t) for t in t0:dt:tend];
t = [i for i in t0:dt:tend];
lent = length(E)
Nt = length(t);

And then changing the nonuniform tangent sampling to match the number of points from indices 1 to length(t):

#nonuniform points
maxpi = pi*22.0/100.0;
dpi = 2.0*maxpi/(lent/2);
l = [i for i in -maxpi:dpi:maxpi];
tt = @. tan(1.95*l);
ttt = @. -tt[1]+tt;
ttt = @. ttt*lent/8.82+1;
t4 = @. ceil(Int64,ttt);

and I had to normalize the timed samples and the pulse amplitudes;

wmax = freq[end]; #max freq get from previous FFTW, around 1.5e15
Nw = length(freq); #number of total points|
dw = 2*wmax/Nw; # spacing of target k grid|
w = dw .* [i for i in ceil(Int64,-Nw/2):ceil(Int64,Nw/2-1)]; #    % a particular uniform M-grid of this spacing
#normalize
Enunorm = Enu ./ maximum(E)
tnunorm = tnu ./ tend

fhat = FINUFFT.nufft1d1(tnunorm, Enunorm, 1, 1.0e-12, Nw);   #% type 1, requesting M modes

With the rest of the code unchanged, my spectrum is now much cleaner:

The only issue I have now is that the frequency isn’t on the center frequency of approximately 5.7e14

So I’m not sure whether the frequencies I chose are exactly the way to do it

[erny123]

Ok figured it out. The nufft1 takes time samples normalized from [-pi , pi] and the amplitudes need to be complex normalized to 1. So the time sampling needs to be normalized as follows:

#normalize
Enunorm = Enu ./ maximum(E)
tnunorm = tnu ./ tend .* (pi)

With this:

Which is extremely close to the actual value.