clearvars
freq = 50;
% https://www.mathworks.com/help/matlab/ref/fft.html
Fs = 1000;            % Sampling frequency                    
T = 1/Fs;             % Sampling period       
L = 1500;             % Length of signal
t = (0:L-1)*T;        % Time vector
f = Fs*(0:(L/2))/L;

%

% Beautiful square: theta redefined
theta = linspace(-2*pi, 2*pi, 1e5);
x1 = square(theta);
x2 = square(theta + deg2rad(15))/2 + square(theta - deg2rad(15))/2;
figure()
hold on
plot(theta/pi, x1)
plot(theta/pi, x2)
xlabel('\theta (\times\pi rad)')
ylabel('MMF (Ni/2)')
hold off

%

% FFT calculation
theta = 2*pi*freq*t;
% theta = linspace(-2*pi, 2*pi, 1e5)

figure()
hold on

x1 = square(theta);
Y = fft(x1);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
plot(f, P1) 

x2 = square(theta + deg2rad(15))/2 + square(theta - deg2rad(15))/2;
Y = fft(x2);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
plot(f, P1) 
xlabel('f (Hz)')
ylabel('|P1(f)|')

%

% Original plot
figure()
hold on
plot(theta/pi, x1)
plot(theta/pi, x2)
xlim([0, 4])
hold off
xlabel('\theta (\times\pi rad)')
ylabel('MMF (Ni/2)')
hold off