Different results applying DFT (formula only) and FFT on Matlab?

Question

I have this function in time-domain and in frequency domain after fourier transform:

$$s_1(t) = (t-2)e^{-t}u(t-2) $$

$$S_1(f)=\frac{e^{-2(1+j2\pi f)}}{(1+j2\pi f)^2} $$

First I create a time vector and a frequency vector ranging from [0s,8s) and [-50hz,50hz) respectively and evaluate them on my functions, using these lines:

s1=(t1-2).*exp(-t1).*heaviside(t1-2);  %s1(t)
S1f=(exp(-2*(1+(2*pi*f1*1i)))./((1+(2*pi*f1*1i)).^2));  %S1(f)

What I get when plotting magnitude/phase is this:

After that, I wanted to compare those results using FFT command in Matlab, so I did this:

S1k=fft(s1);
figure % creates a figure
subplot(2,1,1)     %creates a grid of 2 plots in one figure, selecting the   stem as the first plot
stem(k1,abs(fftshift(S1k,2)),'red') %plots magnitude of S1f
title('Magnitude vs Frequency')
subplot(2,1,2)      %selects the phase plot as the second one in the grid
plot(k1,angle(fftshift(S1k,1)),'blue') %plots magnitude of S1f

title('Phase vs Frequency')

And the result, to my surprise, is this:

As you can tell, there are plenty of differences in both plots, though "shape" its similar at least in the magnitude plot, but with different values in the y-axis.

What can be the problem? I'm sure the fourier transform I did by hand is good, but yet results are different.

Why? Isn't FFT and DFT similar except for the speed of calculation?

Any hints will be appreciated.

Answer 1

The Fast Fourier Transform (FFT) is a class of algorithms that efficiently implement the Discrete Fourier Transform (DFT). However, what you computed by hand is the Continuous-Time Fourier Transform (CTFT), which is substantially different from the DFT. The DFT is applied to finite length discrete-time sequences, whereas the CTFT is applied to continuous-time functions of possibly infinite duration:

$$\begin{align}X[k]&=\sum_{n=0}^{N-1}x[n]e^{-j2\pi kn/N}\quad&\text{(DFT)}\\ X(j\omega)&=\int_{-\infty}^{\infty}x(t)e^{-j\omega t}dt\quad&\text{(CTFT)}\end{align}$$

However, you can use the DFT to approximately compute the CTFT, but there will always be at least one of the following two errors (and usually both):

1. Truncation error: if the signal is too long (or has infinite duration as in your example), it must be truncated to fit inside the DFT window.

2. Aliasing error: the continuous-time signal must be sampled, and sampling will usually introduce aliasing unless the signal is band-limited and the sampling theorem is satisfied.

Since the signal in your example is neither band-limited nor time-limited you will get both errors mentioned above.

More information on approximating the CTFT by the DFT can be found in this answer.

The following Octave script shows how the approximation works for the example in your question.

% choose parameters (sampling frequency and time interval)
Fs = 50;        % sampling frequency
T1 = 2;         % lower and upper
T2 = 20;        % integration limits

Ts = 1/Fs;
DT = T2 - T1;
N = round(DT/Ts);   % DFT length

% sampled time-domain function on [T1,T2]
t = (0:N-1)*DT/(N-1) + T1;
xd = (t-2).*exp(-t);

% DFT approximation of CTFT
Xd = Ts * fft(xd);
fd = (0:N-1)*Fs/N;

% phase correction for offset T1
Xd = Xd .* exp(-1i*T1*2*pi*fd);

% exact expression for CTFT
N2 = round(N/2);
fc = fd(1:N2);
Hc = exp(-2*(1+1i*2*pi*fc))./(1+1i*2*pi*fc).^2;

% plot results
subplot(2,1,1), plot(fc,20*log10(abs(Xd(1:N2))),fc,20*log10(abs(Hc)))
    title('magnitudes (dB)'), xlabel('f'), legend('DFT','CTFT'),
    axis([0,Fs/2,-110,-10]), grid on
subplot(2,2,3), plot(fc,20*log10(abs(Xd(1:N2))) - 20*log10(abs(Hc)))
    title('magnitude error (dB)'), xlabel('f'), grid on
subplot(2,2,4), plot(fc,unwrap(angle(Xd(1:N2))) - unwrap(angle(Hc)))
    title('phase error (rad)'), xlabel('f'), grid on

The plot below shows the magnitudes of the CTFT and its approximation by the DFT, as well as the magnitude and phase approximation errors.