How to use fft on the continuous time signal

Hi I have the system defined by the continuous impulse response h_f:

T_L = 12.5*10^-9; % 12.5 ns
T_S = 0.2*10^-9; % 0.2 ns
T = 1*10^-10; % unit time
Fs = 1/T; % sampling frequency
tau = 1*10^-9; % tau
sigma = 0.1*10^-9; % here we make sigma 0.1ns
t_0 = 0.1*10^-9; % there we make t_0 = 0.1ns
t = -1000*T:T:1000*T; % this is the time array

h_f = exp(-t./tau).*(t>=0); % this is the h function

And I want to compute both the analytical and the numerical solution of it. So here is my code for analytical answer,

% $h_f(t) = exp(-t/\tau)u(t)$ has the fourier transform of
% $\frac{1}{1/\tau+jw}$. So that we can plot the magnitude and phase of
% both.

% theoretical
w = Fs/2001*(-1000:1000);
tau = 1*10^-9;
a = 1/tau;
H_f = 1 ./ (a + 1i*w); % Fourier transform of h_f
magnitude = abs(H_f); % Magnitude
phase = angle(H_f); % Phase
figure;
subplot(2,1,1);
plot(w,magnitude);
ylabel('|H_f| from analytical solution');
subplot(2,1,2);
plot(w,phase);
ylabel('\angle H_f from analytical solution');
xlabel('w (rad/s)');
sgtitle('Magnitude and Phase of the Analytical Solution')

And here is my code for numerical solution using the fft()

% numerical
H_f_numerical = fftshift(fft(ifftshift(h_f)));% the shifted time and then shifted frequency
magnitude_numerical = abs(H_f_numerical)/2001; % Magnitude of numerical Fourier transform
phase_numerical = unwrap(angle(H_f_numerical)); % Phase of numerical Fourier transform
% Normalize the numerical Fourier transform magnitude
figure;
subplot(2,1,1);
plot(w, magnitude_numerical);
ylabel('|H_f| from numerical solution');
subplot(2,1,2);
plot(w, phase_numerical);
ylabel('\angle H_f from numerical solution');
xlabel('w (rad/s)');
sgtitle('Magnitude and Phase of the Numerical Solution')

But when you plot them you can see that they look very different in the magnitude and the phase ( the magnitude at -5 and 5 *10^9 are not exactly 0.2, and the peak height is not 1) (The phase is even worse that it makes no sense) , so I was wondering where I do the code wrong? How to use fft for this continuous time fourier transform?

[Walter Roberson]

The symbolic calculation agrees with your theoretical version.

sympref('heavisideatorigin', 1);
tau = sym(10)^-9;
syms t
h_f = exp(-t/tau) * heaviside(t)

h_f = 

HF = fourier(h_f)

HF = 

absHF = abs(HF)

absHF = 

angHF = angle(HF)

angHF = 

figure
fplot(absHF, [-5 5]*1e9)

figure
fplot(angHF, [-5 5]*1e9)

However... when you are working with control systems, you usually take the laplace transform, not the fourier transform.

Answer 1

[Matt J]

Here is the proper way to set up the time and frequency sampling axes:

tau = 1*10^-9; % tau

%Set up frequency and time sampling parameters and axes

N=1e5;
T=40*tau/(N-1);
dF=1/N/T;

ax = -ceil((N-1)/2):floor((N-1)/2); 
t=ax*T;
w=2*pi*ax/(N*T);

h_f = exp(-t./tau).*(t>=0); % this is the h function

% theoretical
a = 1/tau;
H_f = 1 ./ (a + 1i*w); % Fourier transform of h_f
magnitude = abs(H_f); % Magnitude
phase = angle(H_f); % Phase

Note also that in the numerical calculation, the FFT must be weighted by T to approximate a continuous Fourier transform:

% numerical
H_f_numerical = fftshift(fft(ifftshift(h_f)))*T;% the shifted time and then shifted frequency
magnitude_numerical = abs(H_f_numerical); % Magnitude of numerical Fourier transform
phase_numerical = angle(H_f_numerical); % Phase of numerical Fourier transform

Now we visualize:

skip=1:round(N/1e4):N;
tiledlayout(2,1);
nexttile
plot(w, magnitude_numerical,w(skip),magnitude(skip),'x');   xlim([-1,1]*1e11)
ylabel('|H_f|');

nexttile
plot(w, phase_numerical, w(skip), phase(skip),'x');    xlim([-1,1]*1e11)
ylabel('\angle H_f');
xlabel('w (rad/s)');

sgtitle('Magnitude and Phase')
leg = legend('Numerical','Analytical'); leg.Layout.Tile='East';

Warning: Hardware-accelerated graphics is unavailable. Displaying fewer markers to preserve interactivity.

Note that the numerical version of the phase response will be sensitive to undersampling and time truncation, because the amplitude spectrum tails are nearly zero. The phase is therefore not stable there numerically.

[Zeyuan]

This is great, thank you!