% Square Wave Fourier Series Representation

clc             % clear the command window
clear           % clear all variables
close all       % close any existing windows
format compact	% prevent extra blank lines in the output

% Define and plot the square wave function
% file: square_wave_func.m
% function [f] = square_wave_func(t, T)
% % Function to evaluate a square wave of period T
%   if (t <= T/2)
%         f = 1;
%     else
%         f = -1;
%     end
% end
T = 1;
t = 0 : T/1000 : T;
for i = 1 : 1001
    f(i) = square_wave_func(t(i), T);
end
plot (t, f);

% Plot Fourier Series harmonics
% file: square_wave_amp.m
% function [B] = square_wave_amp(i)
% % Function to evaluate a square wave Fourier Series term amplitude
% B = 2 /(i*pi) * (1 - cos(i*pi));
% end
% file: square_wave_harmonic.m
% function [H] = square_wave_harmonic(i, t, w0)
% % Function to evaluate a Fourier Series term (harmonic) for the square wave
% B = square_wave_amp(i);
% H = B * sin(i*w0*t);
% end
w0 = 2*pi / T;
H1 = square_wave_harmonic(1, t, w0);
H3 = square_wave_harmonic(3, t, w0);
H5 = square_wave_harmonic(5, t, w0);
figure      % open new figure window (to prevent previous from being lost)
plot (t, H1,'Color',[1 0 0]);
hold on;
plot (t, H3,'Color',[0 0 1]);
hold on;
plot (t, H5,'Color',[0 1 0]);

% Plot Fourier Series with different numbers of terms
% file: square_wave_series.m
% function [F] = square_wave_series(t, n, w0)
% % Function to evaluate a square wave Fourier Series to a given number of
% % terms
% F = 0;
% for i = 1 : n
%     F = F + 4/((2*i-1)*pi)*sin((2*i-1)*w0*t);
% end
% end
F1 = square_wave_series(t, 1, w0);
F2 = square_wave_series(t, 2, w0);
F3 = square_wave_series(t, 3, w0);
figure      % open new figure window (to prevent previous from being lost)
plot (t, F1,'Color',[1 0 0]);
hold on;
plot (t, F2,'Color',[0 0 1]);
hold on;
plot (t, F3,'Color',[0 1 0]);
F20 = square_wave_series(t, 20, w0);
figure      % open new figure window (to prevent previous from being lost)
plot (t, F20);
F100 = square_wave_series(t, 100, w0);
figure      % open new figure window (to prevent previous from being lost)
plot (t, F100);

% Plot square wave spectrum
i = 1 : 13;
B = square_wave_amp(i);
figure
bar(i, B);

disp 'Hit Enter to close all windows and quit'
pause
close all