Controversy in calculating the Fourier transform of a function numerically and analytically

Hey,

I am trying to calculate the inverse fourier transform of a complex function numerically. And I am trying to start with calculating the fourier transform of a simple gaussian to make sure my code is correct for simple case. The issue is for a Gaussian function (F = exp(-t^2 / 2*(c^2))), the analytical result for the fourier transform is F.F = 2^(1/2)*pi^(1/2)*exp(-2*c^2*f.^2*pi^2)*abs(c). And it is a real function. But when I calculate the fourier transform numerically using trapz function, the result has both real and imaginary parts (although the imaginary part is too small) but the abs value of the result give the same result as the analytical fourier transform. But for a Gaussian centered at zero there is not an imaginary part in the fourier transform. I don't know why I get this imaginary part when I use the definition of fourier transform directly to calculate fourier transform numerically. Can anyone help why there is a controvercy? I need to know to make sure for my general function it works correctly.

I have copied my code and the imaginary part of the fourier transform plot here.

t=-2:.001:2;
c= 0.5;
x = exp(-t.^2/(2*c^2));

% analytical result of fourier transform
% f=-5:.01:5;
% fourierfunc = 2^(1/2)*pi^(1/2)*exp(-2*c^2*f.^2*pi^2)*abs(c);
% plot(f,fourierfunc)

% numerical result of fourier transform
k=0;
for f=-5:.001:5
    k=k+1;
    fourierx(k)=trapz(t,x.*exp(-1i*2*pi*f*t));
end
f=-5:.001:5;
plot(f,real(fourierx))

[Walter Roberson]

http://matlab.wikia.com/wiki/FAQ#Why_is_0.3_-_0.2_-_0.1_.28or_similar.29_not_equal_to_zero.3F

[Shaily_T]

Thank you!

Answer 1

[Image Analyst]

That's basically quantization noise. Since the computations are done digitally/numerically, not analytically/symbolically, it's hard to make sure you always get exactly zero. If you'd taken a course on linear algebra, numerical analysis, or something similar, they would have covered that.

However here and virtually everywhere in MATLAB if you have values that are around 10^-16 or smaller than your other values, you essentially consider it zero.

[Shaily_T]

Thanks for your response! So, it means if I compute the fourier transform numerically I'll always get this imaginary part even if the fourier transform of that function is a real function. Is that correct? And my code for the numerical calculation of the fourier transform is fine. Right? I mean we can use trapz or integral to calculate the integration of complex functions or numbers as well. Is that correct?

[Walter Roberson]

So, it means if I compute the fourier transform numerically I'll always get this imaginary part even if the fourier transform of that function is a real function. Is that correct?

"Always" is possibly an exaggeration. "Typically" is more likely than "Always" for this situation. That is, it might be the case that you could carefully construct a signal that your calculation left with no residual imaginary portion.

I mean we can use trapz or integral to calculate the integration of complex functions or numbers as well

trapz() of sampling of a function approaches the integral of the function as the number of samples approaches infinity.

[Shaily_T]

Thank you for your explanation!

For my first question, what if I have a function which is determined and I cannot change it? Because you mentioned: " you could carefully construct a signal that your calculation left with no residual imaginary portion".

My main question is, is it harmfull for the result of fourier transform then?

[David Goodmanson]

Hello ST,

Whether it is created by trapz or any other way, if you have an array y whose ratio of [imaginary part] to [real part] is down around 10^-15 or 10^-16, then the function is real for all practical purposes. After you have checked that to be true (and only after you have checked that to be true), you are perfectly justified in using the command y = real(y) and proceeding from there.

Leaving the tiny imaginary part is not harmful, but it can be annoying. If you plot y, you will get a warning that only the real part will be plotted. If you input y into another function, sometimes that other function only takes real arguments, not complex arguments. Also, complex y takes twice as much memory as real y, half of which is used up saving tiny imaginary parts that do no good.

[Walter Roberson]

Gaussian is your test function. but if you define

 f = exp(-(t-t1)^2 / 2*(c^2))

in order to shift the peak to t1, then the fourier transform is complex valued.

You should probably be assuming that your actual function is not symmetric around 0 and that it probably will have some phase in some component.

There are some real-valued functions that have a purely real fourier transform (such as Heaviside) but most are going to give you a complex-valued transform

[David Goodmanson]

Hi Walter,

I don't believe that the heaviside function has a real fourier transform. As continuous function its transform is

N/(2*pi*i*f)

where N is a real normalization factor depending on which transform convention is used. Going to the fft does not alter that conclusion.

However, any real function, symmetric about the origin, has a real (and symmetric) fourier transform.

[Shaily_T]

Thank you so much @David Goodmanson @Walter Roberson for your nice explanations! I appreciate that!