%% Fourier Series Representation of Square Wave
%% Acknowledgements
% While all of the code is my own unless otherwise noted, parts of this notebook were dervied from mathematics presented in:
% Tolga Tasdizen's image processing course, Fall 2019, University of Utah.
% Brad Osgood's course notes from The Fourier Transform and its Applications, EE261, Stanford. (Permission for use in derivative works is granted from Stanford Engineering Everywhere under the Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License.)
%% Set up parameters:

num_harmonics = 1000;
n_arr = 1:num_harmonics; % stores the range of basis frequencies 
T_0 = 2*pi; % period 
T_1 = 1.5; % width of the "nonzero part" of a period of the square wave
omega_0 = 2*pi/T_0; % angular frequency of square waves (radians / second)

%
%% Define the square wave:
% We will define our square wave as a periodic function f(t), with period T_0, and where T_1 controls how wide the nonzero part of the function is.
% f(t) = \left\{
% \begin{array}{ll}
%       1 & | t | \lt T_1\\
%       0 & T_1 \lt | t | \lt \frac{T_0}{2}\\
% \end{array}
% \right.

plt_num_gridpts = 2500;
domain_min = -3*pi;
domain_max = 3*pi;
t_arr = linspace(domain_min,domain_max,plt_num_gridpts);
sq_wave = zeros(size(t_arr));
for i=1:size(t_arr,2)
    t = mod(t_arr(i),T_0); % modulus operator, needed for periodicity
    if T_1 < abs(t) && abs(t) < (T_0 - T_1)
        sq_wave(i) = 0;
    else 
        sq_wave(i) = 1;
    end
end

%

clf;
set(gcf,'position',[0,0,1920,1080])
plot(t_arr,sq_wave,'LineWidth',3)
ylim([-0.1 1.1])
%title("Square Wave") 
xlabel('t') 
ylabel('f(t)') 

%
%% Compute the Fourier cofficients for the square wave:
% Our square wave has an angular frequency of \omega_0 = \frac{2\pi}{T_0} radians per second. Because our square wavef(t) is an even function, we can approximate it as follows:
% f(t) \approx \frac{a_0}{2} + \sum_{n=1}^{N}a_n\cos(2\pi n t)
% Where the Fourier cofficients for n > 0 are given by:
% a_n = \frac{4\sin(\omega_0 n T_1)}{\omega_0 n T_0}
% Finally, a_0, which is sometimes called the DC component, is given by:
% a_0 = \frac{2 T_1}{T_0}

c_n = 2*sin(omega_0*n_arr*T_1)./(omega_0*n_arr*T_0); % complex Fourier cofficients for n > 0
c_0 = 2*T_1/T_0; % n=0, zero frequency, a.k.a. "DC component"

a_n = 2*real(c_n); % "Cosine series" Fourier cofficients for n > 0 
a_0 = 2*c_0; % for n = 0

%
%% Sample a Fourier basis:
% Below, we sample the cosine function:
% At spatial locations t given in the array t_arr
% At the frequencies F = n cycles per second (alternatively, at angular frequencies \omega = 2\pi F = 2 \pi n radians per second), where the values of n are given in the array n_arr
% The aforementioned samples are stored in the 2D array cos_basis. cos_basis(i,j) stores the value of the cos(2\pi nt), where the value of n is given by n_arr(i) and the value of t is given by t(j)

nt = n_arr' * t_arr; % outer product, equivalent to nt(i,j) = n(i) * t_arr(j) for all valid i,j
cos_basis = cos(2*pi*nt/T_0); % cos_basis[i,j] = cos(2*pi*n[i]*t_arr[j]) for all valid i,j

%
%% Visualize the harmonic basis (or at least part of it)
% What this means plotting cos(\frac{2\pi n t}{T_0})for n = 1,2,3,4,5.
% I don't plot higher frequencies here because it would make the figure too messy.

clf;
hold on;
axis([-0 T_0/2 -1.1 1.1]);
plot(t_arr,cos_basis(1,:),'LineWidth',3);
plot(t_arr,cos_basis(2,:),'LineWidth',3);
plot(t_arr,cos_basis(3,:),'LineWidth',3);
plot(t_arr,cos_basis(4,:),'LineWidth',3);
plot(t_arr,cos_basis(5,:),'LineWidth',3);
legend({'cos(ωₙt)','cos(2ωₙt)','cos(4ωₙt)','cos(6ωₙt)','cos(8ωₙt)'},'Location','east');
title("First five harmonics")
xlabel('t')  
hold off;

%
%% Approximate f(x) using a truncated Fourier basis:

scaled_basis = a_n' .* cos_basis; % componentwise multiply using the " .* " operator
clf;
set(gcf,'position',[0,0,1600,720])

%
%% Plot of approximated function, sinusoids, and spectrum, for Fourier series with N terms
% Move the slider below to change N:

N = 10; % <- slider, goes up to N=10
f_approx = a_0/2 + sum(scaled_basis(1:N,:),1);
subplot(1,2,1); plot(t_arr,f_approx,'LineWidth',3); xlabel('t'); 
subplot(1,2,2); plot(t_arr,scaled_basis(1:N,:),'LineWidth',3); 

%

clf;
set(gcf,'position',[0,0,1280,720])

%

N = 1000; % <- slider, goes up to N=1000
f_approx = a_0/2 + sum(scaled_basis(1:N,:),1);
plot(t_arr,f_approx,'LineWidth',1); xlabel('t'); ylabel('Approx. of f(t)');

%%
%
%
%% Licensing
% This work is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License.