%% specrect.m plot the spectrum of a square wave
close all 
Ts=1/100;          % time interval between adjacent samples
t=0:Ts:20;         % create a time vector
x=[t <= 2];        % rectangular pulse $1[0 \leq t \leq 2]$
plotspect(x,t)     % call plotspect to draw spectrum
xlim([-5,5])       % look only from f = -5 to f = 5 Hz

%% SymbFourier
syms t f;                              
syms a positive
s = exp(-a*t)*heaviside(t);
S = fourierf(s)                 % Find S(f)


%%===============================================

% plotspect(x,t) plots the spectrum of the signal x
% whose values are sampled at time (in seconds) specified in t

function plotspect(x,t)
N=length(x);                               % length of the signal x
Ts = t(2)-t(1);                            % find the sampling interval
ssf=((-N/2):(N/2-1))/(Ts*N);               % frequency vector
fx=Ts*fft(x(1:N));                         % do DFT/FFT
fxs=fftshift(fx);                          % shift it for plotting
subplot(2,1,1); 
set(plot(t,x),'LineWidth',1.5);            % plot the waveform
xlabel('Seconds');                         % label the axes
subplot(2,1,2);          
set(plot(ssf,abs(fxs)),'LineWidth',1.5);   % plot magnitude spectrum
xlabel('Frequency [Hz]'); ylabel('Magnitude')   % label the axes
end

%% fourierf.m
function G = fourierf(g)
syms f
G = simplify(subs(fourier(g),'w',2*pi*f));
end