Getting imaginary values for unequal spaced Gaussian function using NUFFT

Why does it give imaginary values for Gaussian function while using NUFFT function?

I am using the following script

clear;clc;
sigma=0.1;
a=1/(2*sigma*sigma);
xmax=50*sigma;
xrange=2*xmax;
n=2^8;
dx=xrange/n;
x=-xmax:dx:xmax-dx;
x1=x(1:10);
x2=x(21:120);
x3=x(135:178);
x4=x(197:250);
x=[x1 x2 x3 x4]; %comment this for equal spacing
f_x=exp(-a*x.*x);%original function
f=linspace(-1/2,1/2 -1/n, n);
f_k=nufft(f_x,[],f);
imag(f_k)

I tried using https://in.mathworks.com/matlabcentral/fileexchange/25135-nufft-nfft-usfft. Still getting imaginary values.

[Bjorn Gustavsson]

What's f_x?

[Athira B S]

f_x=exp(-a*x.*x);%original function

Answer 1

[Bjorn Gustavsson]

If you take a look at what you work with:

subplot(3,1,2)
plot(real(f_k))
hold on       
plot(imag(f_k))
plot(abs(f_k)) 
subplot(3,1,3)
plot(angle(f_k))
subplot(3,1,1)
plot(x,f_x,'.-')

You'll see that the cut around zero in x removes almost the entire peak of your Gaussian. Assymetrically. Then there's no wonder any Fourier-expansion will not be "schoolbook ideal".

HTH

[Athira B S]

Thanks @Bjorn Gustavsson. I understand your point. I have re-written the script in such a way that it does not remove the peak. But still getting imaginary values. Should I be expecting this?

clear;clc;
sigma=0.1;
a=1/(2*sigma*sigma);
xmax=50*sigma;
xrange=2*xmax;
n=2^8;
dx=xrange/n;
x1=-xmax:2*dx*1:0;
x2=dx:dx:xmax-dx; 
x=[x1 x2];
f_x=exp(-a*x.*x);
n = length(x);
f=linspace(-1/2,1/2 -1/n, n);
f_k=nufft(f_x,[],f);
imag(f_k)

[Bjorn Gustavsson]

Well....

...from the help to nufft we have this nugget of information:

If t is specified as [], the sample points in the transform are

0:(N-1).

Meaning that if you don't supply the information of what t - values your samples are at a uniform sampling is assumed. Therefore obviously your second example also will have significant imaginary components since it is interpreted as both assymetric and shifted from the mid-point. If you send in your x you'll be both more successful and happier with the results.

HTH

[Athira B S]

I used,

f_k=nufft(f_x,x,f);

Still getting imaginary components!

[Bjorn Gustavsson]

Why does that bother you? You have a nonuniform sampling and will not get as pretty behaviour as for the continuous Fourier-transform. Try to set the width to something wider and see that it behaves reasonably sensibly:

a=0.01/(2*sigma*sigma);
f_x=exp(-a*x.*x);
subplot(3,1,1)
plot(x,f_x,'.-')
f_k=nufft(f_x,x,f);
subplot(3,1,2)
hold off
plot(abs(f_k),'.-')
hold on 
plot(real(f_k),'.-')
plot(imag(f_k),'.-')
subplot(3,1,3)
plot(angle(f_k),'.-')