% EE3032 W5D3 Winter 2019-'20, Dr. Durant

close all % close any open figures

%% I. Show that the phasors (Fourier series coefficients) given constitute a square wave.
dt = 0.01;
t = -3:dt:3;

T0 = 1; % s, fundamental period
omega0 = 2*pi/T0; % radians/s, fundamental frequency
N = [1   3   5    7 ];
X = [1 -1/3 1/5 -1/7] * 4/pi; % will derive this in a few days; in general, phasors X can be complex
% Try adding more components to N and X to continue the pattern and better approximate x(t).

figure, subplot(2,1,1) % start a new figure and plot in the top half
x = zeros(size(t));
for idx = 1:length(N)
    part = abs(X(idx)) * cos(N(idx)*omega0*t + angle(X(idx)));
    plot(t,part), hold on % when plotting, don't clear what was there before
    x = x + part; % add parts of x(t) to which we will apply superposition
end
plot(t,x), title('x(t) and components using phasors X')

%% II. Apply transfer function theory and superposition to determine the
% system output given a square wave input.
tau = 0.1; % s, time constant in h(t) = (1/tau)*exp(-t/tau)u(t)
a = 1/tau; % in class, we forgot to take the reciprocal, so we were using large time constants when we intended short ones and vice versa
H = 1 ./ (a + 1j * omega0*N);
disp('Magnitudes and angles of H:'), disp(abs(H)), disp(angle(H))

subplot(2,1,2) % switch to the bottom half of the current figure
y = zeros(size(t));
for idx = 1:length(N)
    part = abs(X(idx)) * (abs(H(idx))) ...
        * cos(N(idx)*omega0*t + angle(X(idx)) + angle(H(idx)) );
    plot(t,part), hold on
    y = y + part;
end
plot(t,y), xlabel('Time (s)'), title('y(t) and components, found using phasors Y = HX')