close all

pkg load signal  % is used for square


% square generates square wave with period 2pi and values between -1
% and 1

fs = 44100;  % sampling rate
fs = 44100*2;
L = 150;    % Frequency
t = linspace(0,2*pi*L,fs);

tend = 1000; % For plotting, given in samples

y = 0.5*square(t);

N1 = 4;  % first sum goes from 1:2:2*4+1 = 1:2:9
N2 = 16;
N3 = 32;
N4 = 64;
N5 = 128;

SN = 0.0;

for n = 0:N1,
  SN = SN + 4.0/( (2*n+1)*pi ) * sin( (2*n+1) * t );
end

S1 = 0.5*SN;   % remove problems with clipping when playing the sound

for n = N1+1:N2,
  SN = SN + 4.0/( (2*n+1)*pi ) * sin( (2*n+1) * t );
end

S2 = 0.5*SN;

for n = N2+1:N3,
  SN = SN + 4.0/( (2*n+1)*pi ) * sin( (2*n+1) * t );
end

S3 = 0.5*SN;

for n = N3+1:N4,
  SN = SN + 4.0/( (2*n+1)*pi ) * sin( (2*n+1) * t );
end

S4 = 0.5*SN;

for n = N4+1:N5,
  SN = SN + 4.0/( (2*n+1)*pi ) * sin( (2*n+1) * t );
end

plot(t(1:tend),y(1:tend))
hold
%keyboard
plot(t(1:tend),S1(1:tend),'r')
%keyboard
plot(t(1:tend),S2(1:tend),'g')
%keyboard
plot(t(1:tend),S3(1:tend),'k')
%keyboard
% if you want to play the signals just use
%  sound(S1,fs) where fs is the sampling frequency
% if you change the sampling frequency the sound gets higher or lower
%also listen to the difference
%
% sound(y - S4,fs) but there is an artefact created by the sampling
% and the numerical accuracy, because the max position cannot be found
% exactly by the sampling, so the sound changes its amplitude. That
% should get better with higher sampling frequency
