The code below illustrates the analysis for one component of the signal. It can be adapted to include both components of the signal.
Note the two key differences between this code and the code in the question, both of which are noted below.
First, compute the exact CTFT
syms t w f
g(t) = rectangularPulse(-0.5,0.5,t/2);
G(w) = simplify(fourier(g(t),t,w))
which is the expected result stated in the question.
Now compute the samples of the signal g(t) at the speciifed sampling period over the specified time duration.
fs = 10000; Ts = 1/fs;
tval = -2:1/fs:2; % might be bettter to use linspace
gval = rectpuls(tval/2);
Note that rectpuls() results in gval = 1 at tval = -1, but gval = 0 at tval = 1 (which I believe is a desired result).
[tval(10001) gval(10001)]
[tval(30001) gval(30001)]
Consequently, gval has an even number of non-zero samples for parameters of the question
nnz(gval)
This fact will be important next.
Now compute the DTFT of gval using freqz() over the range -30 to 30 rad/sec. freqz() assumes that the first point in gval corresponds to the discrete-time sample n = 0, so we have to shift the output of freqz() usnig the DTFT time-shifting property to account for the fact that the first element of gval actually corresponds to n = -20000 (-20000/fs = -2)
wdtft = linspace(-1,1,1001)*30; % -30 rad/s to 30 rad/sec
dtft = freqz(gval,1,wdtft*Ts).*exp(1j*wdtft*Ts*20000);
dtft has a small, non-zero imaginary part because nnz(gval) is even, i.e., gval is not an odd sequence.
max(abs(imag(dtft)))
To compare to G(w), we'll just look at the real part.
dtft = real(dtft);
Now compute the DFT of gval, the elements of which are samples of the DTFT. Shift the DFT in the frequency domain for the same reason we shifted the DTFT.
dft = fft(gval).*exp(1j*(0:(numel(gval)-1))/numel(gval)*2*pi*20000);
wdft = (0:(numel(dft)-1))/numel(dft)*2*pi*fs;
Get rid of the small imaginary part.
max(abs(imag(dft)))
dft = real(dft);
Now compare the CTFT and DTFT over the frequency range -30 to 30 rad/sec, and the DFT over 0-30 rad/sec. Could use fftshift() if we wanted to get the image of the DFT onto negative frequencies. Note that the code in the question did not use fftsthift() even though f was defined over from -fn to fn. Also, the DTFT and DFT have to be scaled by the sample period to approximate the CTFT; this scaling also needs to be added ot the code in the question.
figure
hold on
fplot(G(w),[-30 30])
plot(wdtft,dtft*Ts)
plot(wdft(wdft<30),Ts*dft(wdft<30),'ko'),grid
xlabel('rad/sec')
legend('CTFT','T*DTFT','T*DFT')
Note that starting with the third point of the DFT every other point is zero, which would result in a triangular look if we plotted the abs() of the three curves and connected the samples of the DFT.
