a saw-shaped IFFT in matlab

I have a problem with the inverse Fourier transform of a function that is supposed to be smooth, but in practice, it becomes saw-shaped. The function in the Fourier domain is:

f(m,n)=i*m*(heaviside(m)+heaviside(-m))/(m^2+n^2)

Applying ifft2 to this function must result in a smooth function as:

F(x,y)=x/(x^2+y^2)

where m and n are the variables in fourier domain corresponding to x and y in space domain, respectively and i is sqrt(-1).

However, doing this leads to a saw-shaped function in the x-direction.

Here is the code I have written:

clear all

N=pow2(8);

Lx=0.01;

XX=linspace(-Lx,Lx,N);

[X,Y]=meshgrid(XX,XX);

dx=XX(2)-XX(1);

% the grid in fourier domain is chosen in such a way to avoid the

% singularity at the origin

mm=2*pi/N/dx*linspace(-N/2,N/2,N);

[m,n]=meshgrid(mm);

fmn=1i*m.*(heaviside(m)+heaviside(-m))./(m.^2+n.^2);

% the numerical ifft2 of the spectrum

fxy=1/dx^2*fftshift(ifft2(ifftshift(fmn)));

% the analytical ifft2 of the spectrum

Fxy=X./(X.^2+Y.^2);

subplot(2,1,1)

mesh(X,Y,abs(fxy))

title('the output by ifft2')

subplot(2,1,2)

mesh(X,Y,abs(Fxy))

title('the expected (analytical) inverse transform')

I would appreciate any suggestion in this regard.

Mohammad

[Walter Roberson]

Formally speaking your function is not defined at m=0, but if you assign any non-infinite value to heaviside(0) it becomes -i*m/(m^2+n^2) . It is not obvious to me that the 2D inverse fourier of that would be of the form you indicate.

[[e-mail]]

Thank you Walter for your remark. I don't think the problem is the heaviside function, because it's just a simple example that I'm working on to figure out the problem. However, the practical function (not included here! but typically similar to this example) does not have heaviside function in it, but its ifft2 output has the same behavior. If you run the code you will see the ifft2 output is not smooth as it is expected to be (the expected output is also plotted when running the code).

Answer 1

[Bjorn Gustavsson]

I think that the problem you have is that your frequency coordinates are defined without 0 when you calculate fmn - this means that your fmn doesn't have any DC-component, then you send that into the ifft2 function after properly doing ifftshift. However the ifft2 function will treat first row and column as the DC/0-frequency component of an FFT. So when you send ifft2 something you intend to be symmetric, ifft2 doesn't agree with you so you'll get something other than you expect. To give ifft2 something it will transform right you'd have to give it the 0-frequency components, and to do that you'll have to decide how to handle the divergence to infinity...

HTH