Huge difference between the result of fft function Matlab and analytical Fourier transform of the same function

Hey,

I am trying to find a way to obtain the numerical fourier transform of a function (it is not a signal and I only want to obtain the numerical fourier transform of a function). For a test code, I tried to see what is the result of fft matlab for a Gaussian function and compared it with the analytical fourier transform of this Gaussian. I have attached the plot of both results. Why the amplitude of fft result is that huge compared to analytical result? I can use fftshift to shift the result of fft to center but still the issue of amplitudes are there. I am wondering isn't it because it is a Gaussian function and so not periodic over time? Does fft only works if the function we have is spanned from -infinity to +infinity (basically a signal)? or we can use it if we want ot calculate the fourier transform of a function which spans from t-1 to t-2?

[Paul]

If the code is posted, I'm sure we can sort it out.

Also, by definiion the DFT (as computed by fft) can only be applied for signals that have have finite duration, so NOT from -inf to inf. If the signal of interest is of infinite duration, then fft can be used if we're willing to live with a windowed version of the signal, with the windowing fucntion being zero where the underlying signal is small.

Not sure what distinction is being made between signals and functions. Mabye you're using the term "function" to mean a finite duration signal? I've not heard of that nomenclature, but maybe it's used.

[Shaily_T]

Thanks for your response! I mean, my function starts at a certain point and ends at another certain point in time. For example, if we have a sin(2*pi*t) function, it is a periodic function that spans over the whole time range but a Gaussian function will start at a certain point and end at another certain point. My question is, can fft be used for both cases? By your second response, do you mean what I have done and shown in the code now? Multiplying fftshift by ts. Because now, I can get a similar result to the analytical form.

Here is my code for a Gaussian function.

I also have confusion which I appreciate your help for. I checked the result of fftshift by incorporating your guidance for multiplying it by ts, and its result is the same as using the direct definition of the fourier transform by allying trapz function. However, my question is, why is the result of fft before implementing fftshift strange? I have attached the plot. Indeed it shows two peaks in the frequency domain, which I’m pretty sure is not the case for the fourier transform of a Gaussian. I have seen some explanations, but still, I don’t understand why.

%% test for fourier of simple function

N = 5E4;

ts = 0.1;

t = (-N/2:N/2-1)*ts;

f = exp(-pi*(t.^2));

% plot(t,f)

deltaf = 1/(N*ts);

w = ((-N/2:N/2-1))*deltaf;

fourierf = exp(-pi*w.^2);

plot(w,fourierf)

fft_f = fft(f);

shift_fft_f = ts*fftshift(fft_f);

plot(w,abs(shift_fft_f))

Answer 1

[Paul]

Can't say for sure w/o seeing the code, but I supsect the fftshifted curve will be close to the blue curve if you mulitply the fftshifted result by the sampling period.

[Paul]

Hi Shaily_T,

Addressing the issues raised in this comment here.

"For example, if we have a sin(2*pi*t) function, it is a periodic function that spans over the whole time range but a Gaussian function will start at a certain point and end at another certain point."

A Gaussian (of the form exp(-a*t^2), a>0) is non-zero from -inf to inf, isn't it? It decays rapidly, but only approaches zero asymptotically. In other words, the Gaussian signal is of infinite duration (as is sin(t)).

"My question is, can fft be used for both cases? "

The FFT algorithm, which computes the Discrete Fourier Transform (DFT), is only applicable to discrete-time signals of finite duration, i.e., signals x[n] that are zero for n larger/smaller than an upper/lower bound. So no, fft can't be applied to sin(t) or exp(-a*t^2) (note that sin(t) is a different animal because it doesn't have convergent Continuous Time Fourier Transform (CTFT) integral, though it does have a Fourier transform as a sum of scaled Dirac deltas). However, because exp(-a*t^2) becomes small for some value of t s.t.. abs(t) > constant (c), we can multiply exp(-a*t^2) with a rectangular window of unit height and width from -c to c, take samples of that windowed signal, which is of finite duration, and then take the fft of those samples. As long as we choose c large enough, the FFT of those samples can be used to approximate samples of the CTFT of exp(-a*t^2).

"By your second response, do you mean what I have done and shown in the code now? Multiplying fftshift by ts."

Yes, the DFT samples (as computed by fft) of the sampled signal need to be muliplied by Ts to approximate samples of the CTFT of the continuous time signal for frequencies between -pi/Ts to pi/Ts (rad/sec). The reason for this is that the DFT samples are basically samples of the Discrete Time Fourier Transform (DTFT). The DTFT is basically a sum of the samples spaced by Ts of the CTFT integrand. If we want to approximate the CTFT using rectangular integration, we need to multiply the DTFT sum by Ts for the DTFT to approximate the CTFT. Because the DFT samples are samples of the DTFT, the DFT samples need to be muliplied by Ts as well.

"However, my question is, why is the result of fft before implementing fftshift strange? "

It's not really strange when plotted against the correct frequency vector and when we understand what the output of fft really means.

First, let's get the CTFT of the Gaussian signal

%% test for fourier of simple function
% Gaussian signal
syms t w real
y(t) = exp(-t^2/4/sym(pi))/2/sym(pi)

y(t) = 

% Its CTFT, w in rad/sec
Y(w) = fourier(y(t),t,w)

Y(w) = 

I'm using the default fourier parameters, so the transform isn't unitary, but that doesn't matter for this discussion.

Now do the fft stuff.

First, we need a sampled signal

N = 5E4;   % number of samples
Ts = 0.1;  % sampling period
nvals = (-N/2):(N/2-1);            % discrete time vector, N even
tvals = (-N/2:N/2-1)*Ts;           % sample of continuous time 
yvals = exp(-tvals.^2/4/pi)/2/pi;  % samples of the signal

This sampling assumes that y(t) is basically zero for t < tvals(1) and t > tvals(end). Let's check

vpa(y([tvals(1) tvals(end)]),10)

ans = 

Indeed, y(t) is quite small at the endpoints of the sampled time vector. So small that they are zero in double precision

yvals([1 end])

ans = 1×2

0 0

Now to the frequency domain

% Take the fft
Yfft = fft(yvals);
% The frequency vector corresponding to the elements of Yfft (rad/sec)
wvals = (0:N-1)/N*2*pi/Ts;
% Plot
figure;
subplot(211)
plot(wvals,abs(Yfft));ylabel('Magnitude')
subplot(212)
plot(wvals,angle(Yfft));ylabel('Phase (rad)')
xlabel('w (rad/sec)')

We see two issues.

First, the amplitude of the DFT has a peak at low frequency starting at w = 0 and a peak at what appears to be high frequency ending at w = 62.28. But recall that the DFT is actually evaluated at angles ccw around the unit circle starting at theta = 0. Halfway around the unit circle, theta = pi, corresponds to the Nyquist frequency of pi/Ts = 31.14 rad/sec (approximately), and frequencies "higher" than that are actually coming back around the unit circle towards dc. (Just redo the plots above, but with wvals*Ts as the x-axis values, and you'll see the plots from 0 to 2*pi (actually, one frequency sample less than 2*pi). This wrapping is exactly addressed by fftshift which reorders the DFT samples so they instead correspond to angles around the unit circle roughly from -pi to pi, instead of 0 to 2*pi, or frequencies from -pi/Ts to pi/Ts.

In short, those two peaks are perhaps more properly viewed as the right and left sides of a single peak centered at w = 0, as we'll see below.

But we also see that the DFT has a nonzero imaginary part, as exhibited by the phase. This happens because the DFT (at least as implemented in Matlab by fft) assumes the first sample of the input corresponds to n = 0, whereas in our signal the first sample corresponds to n = -N/2, In other words, we're taking the DFT of a time-shifted version of our signal, and time shifting corresponds to a phase shift in the frequency domain. So we have to adjust the output of fft to account for this effect

Yfft = Yfft.*exp(-1j*nvals(1)*wvals*Ts);
figure;
subplot(211)
plot(wvals,abs(Yfft));ylabel('Magnitude')
subplot(212)
plot(wvals,angle(Yfft));ylabel('Phase (rad)')

Now see the phase is zero as expected for low and "high" frequencies where abs(Yfft) isn't very small. Where abs(Yfft) is very small, the phase is just floating point trash. So take the real part

Yfft = real(Yfft);

BTW, it may have been better to apply ifftshift to yvals before applying fft and then we wouldn't have to do the phase adjustment, but I think the approach above is more instructive.

Now do fftshift

Yfft = fftshift(Yfft);

The corresponding frequencies (rad/sec) are now

wvals = nvals/N*2*pi/Ts;

Compare the (shifted) DFT, scaled by Ts, to the CTFT

figure
fplot(Y(w),[wvals(1) wvals(end)])
hold on
plot(wvals,Ts*Yfft) % scale by Ts
xlim([-1 1]);
xlabel('rad/sec')
legend('CTFT','Ts*DFT')

The curves are basically on top of each other.

[Shaily_T]

@Paul

Hi Paul,

Thank you so much for your time and your comprehensive and instructive response!

I think I got the main points and message! I still do have some questions from your instrutions to make sure I undrestood the whole idea and if you could comment on them I really appreciate it. I undrestand you have already been so generous in your response and it is really appreciated.

1. So based on your explanations, I can only use fft for a sin(t) if it is a finite signal and has a start and end point. Is that correct?
2. I think I do have a confusion regarding how we transfer time to frequency when calculating the fourier transform. My undrestanding is for example for fft because in the fourier transform definition there is exp(-i2*pi*t*k) so for going from time to frequency if ts is time intervals then 1/(N*ts) would be the frequency interval (because indeed it is frequency and not angular frequency). But I see you have a 2*pi factor as well. Could you please comment on it?
3. Indeed, all of your explanations to instruct me about the sample reordering and the fact that the result of the fft is properly viewed as the right and left sides of a single peak centered at w = 0 and the phase adjustment point that you mentioned we should do, all of these stuff are done by fftshift(fft). Is that correct? I mean indeed our correct answer for fourier transform of a function would be fftshift(fft) and if we use fftshift of fft and not fft itself, indeed we are incorporating all of the points you mentioned above to instruct me. Is that correct?

Also, I have a question which is in a similar spirit to my Q2 here. And it is few days that I am thinking

how to deal with that. I have posted it here https://www.mathworks.com/matlabcentral/answers/1802830-how-to-obtain-the-time-domain-from-the-frequency-domain-while-computing-the-numerical-inverse-fourie?s_tid=prof_contriblnk

I totally undrestand it is another question but I thought if it is possible for you to have a look at it as well. I really appreciate your time.

Again thank you!

[Paul]

Shaily_T,

These are good questions and I'll do my best to answer them here briefly. However, this forum isn't a substitute for a well written text book, a well-writen set of course notes, or a well-instructed class on these topics.

1. The DFT, which is what fft computes, is only applicable to discrete-time, finite duration signals (discussed above). sin(t) is not of finite duration, so it can't be sampled into a disctete-time, finite duration signal from which a DFT can be computed. We can immediately see a problem when we realize that the DFT samples of any finite duration signal are finite in magnitude (and duration), but the DTFT of x[n] = sin(n*Ts) is given in terms of Dirac deltas (link), as is the CTFT of sin(t) (entry 305). The output of fft will never be a Dirac delta. What one can do, and examples abound on this forum, is take samples of sin(t) over a finite duration, and then compute the DFT of those samples. In this case, the DFT yields frequency domain samples of the DTFT of samples of sin(t)*window(t), which in turn is an approximation ot the CTFT of sin(t)*window(t). But sin(t)*window(t) is not the same as sin(t). Alternatively, with an appropriate selection of Ts and appropriate selection of the finite number of samples of sin(n*Ts), the DFT of the samples can be related to Fourier series coeffients of sin(t), but, again, that's not the Fourier transform. Of course, all of these concepts are related to each other, but (IMO) it's important to understand the differences between them.
2. The 2*pi factor comes directly from the definiton of the DFT (link). Each sample of the DFT, X[k], is a weighted sum of powers of exp(-1j*2*pi/N*k), where k goes from 0 to N-1. We see that the unit circle is being sampled at angles (rad) of 2*pi/N*k. So my preference is to keep the 2*pi scaling because it reminds me of what's happening mathematically, but its just a preference. I've seen treatments that explicilty bring Ts into the fft sum, but it's really just multiplying by 1, i.e., exp(-1j*2*pi/Ts/N*k*Ts), which can be way of viewing the frequency domain samples as being spaced by 2*pi/Ts/N.
3. I'm afraid that's not totally correct. When using fft we have to keep in mind two things. First, fft as implemented in Matlab assumes that the first point in the input sequence is a sample at n = 0 (or t = 0 if viewed as a sample of a continuous-time signal). But, in our example problem, the first point in our sequence is sampled at a negative time, corresponding to n = -N/2. We could define the fft as a sum from -N/2 to N/2 -1 (N even), which would then fit nicely with the problem at hand (but wouldn't for others), but that's not how Matlab does it, so we have to compensate for what Matlab assumes relative to what we actually have. For this particular problem, with N even and n defined as -N/2 : (N/2-1) we could take advantage of ifftshift before taking the fft as a different option. Second, because the Matlab implementation of fft takes the sum from 0 to N-1 (all positive numbers) the corresponding frequency samples are also all positive and monotonically increaing from 0 to ~2*pi (or 2*pi/Ts is scaling by Ts). Then we use fftshift to just reorder the ouput from fft as if fft computed the DFT sum from -N/2 to N/2 - 1 (N even), and then we compute the corresponding frequency vector in whatever units we prefer. In short, we can always use fftshift if we choose (we don't have to), but whether or not we need to phase-adjust the output of fft is problem-dependent, and only matters if we care about the phase in the first place.

[Shaily_T]

@Paul

Paul,

Thank you so much for your explanations and clarifications! I appreciate that.

so, for Q2, having 2*pi or not having it is a preference and we should just make sure that we are consistent in keeping it or not keeping itthrough our calculations. Is that correct?

And for Q3, all you mentioned only affects the phase of the Fourier transform and the magnitude would not be affected. Is that correct?

[Paul]

Q2: Correct. The output of an N-point FFT are frequency domains samples with indepedent variable n = 0:(N-1) (the samples are zero outside this range thought not shown explicitly). The sample values n = 0:(N-1) can also be mapped to frequencies by multiplying n by 2*pi/N (rad/sample), or 1/N (cycles/sample), or 2*pi/N/Ts (rad/sec), or 1/N/Ts (Hz) depending on whatever is preferred. After applying fftshift, the independent variable is n = -N/2:(N/2 -1) for N even and n = ( (N-1)/2 ) : ( (N-1)/2 ) for N odd. However, keep in mind that if doing phase adjustment that the "frequency," w, used in exp(1j*n0*w) must be in rad (aka rad/sample)

Q3:Correct. Translation in the time domain doesn't affect the magnitude of the Fourier transform. However, keep in mind that magnitude and phase don't have clear meaning for signals with FTs that include Dirac deltas. But even in this case, the time translation still only results in a multiplication by a complex exponential of the FT of the unshifted signal.

syms t w
f(t) = cos(t);
F(w) = fourier(f(t))

F(w) = 

% show that the FT of the shifted signal is an exp*FT of the unshifted
% signal
simplify(fourier(f(t-0.5)) - exp(-1j*sym(0.5)*w)*F(w))

ans = 

0

Also, for the N-point DFT, instead of using the explicit phase adjustment for an N-point signal that doesn't start at n = 0 (as done above) we can use a circular shift on the input to fft. For example, suppose we have a signal with N = 8 values

rng(100)
N = 8;
x = rand(1,N);

that correspond to discrete-time values n of

n = -17:10;

If we want to use the DFT to get N samples of the DTFT we can do either

w = (0:(N-1))/N*2*pi;
X1 = fft(x).*exp(-1j*w*n(1));

or

X2 = fft(circshift(x,n(1)));
% show equal results
max(abs(X2-X1))

ans = 6.7642e-15

I suspect the latter is more numerically robust.

[Shaily_T]

@Paul

Thank you so much Paul! I really appreciate your time to clarify all of this things to me! Thank you!

Answer 2

[David Goodmanson]

Hi ST,

Your frequency grid runs from -5 to 5, which for an fft is -fs/2 to fs/2 (fs being the sampling frequency), so fs = 10. The sampling interval delta_t = 1/fs = 1/10. It appears that you are trying to approximate a continuous Riemann integral

Int g(t) e^(-2pi i f t) dt   

by using the fft. That approximation is a sum over indices (done by the fft) times the width of the intervals.

Sum  (stuff)  * delta_t

So if you multiply by 1/10, that takes the fft result down by exactly the right amount.

[Shaily_T]

@David Goodmanson

Hi DG,

Thanks for your nice explanation! so, fft function indeed do a summation over the indices. Is that correct?

Also, If possible could you please comment on my latter questions as well? I mean this part:"I also have confusion which I appreciate your help for. I checked the result of fftshift by incorporating your guidance for multiplying it by ts, and its result is the same as using the direct definition of the fourier transform by allying trapz function. However, my question is, why is the result of fft before implementing fftshift strange? I have attached the plot. Indeed it shows two peaks in the frequency domain, which I’m pretty sure is not the case for the fourier transform of a Gaussian. I have seen some explanations, but still, I don’t understand why."

Thank you so much!