% har_sqw.m - harmonics of a square wave
% HJSIII - 98.07.07

clear
clf

% bipolar square wave +/-5 V at 50 Hz, 50% duty cycle
A0 = 5;
f0 = 50;
tau = 1 / f0;

% sampling rate for display
fs = 5000;
dt = 1 / fs;

% loop to fill buffers
for is=1:121,
  t = (is-11) * dt;
  time(is) = t;
  if mod( (t+tau), tau ) < 0.5*tau,
    sqw(is) = A0;
  else
    sqw(is) = -A0;
  end    
  v1(is) = (4*A0/pi/1) * sin( 1 * 2*pi*f0*t );
  v3(is) = (4*A0/pi/3) * sin( 3 * 2*pi*f0*t ) + v1(is);
  v5(is) = (4*A0/pi/5) * sin( 5 * 2*pi*f0*t ) + v3(is);
  v7(is) = (4*A0/pi/7) * sin( 7 * 2*pi*f0*t ) + v5(is);
  v9(is) = (4*A0/pi/9) * sin( 9 * 2*pi*f0*t ) + v7(is);
  v11(is) = (4*A0/pi/11) * sin( 11 * 2*pi*f0*t ) + v9(is);
  v13(is) = (4*A0/pi/13) * sin( 13 * 2*pi*f0*t ) + v11(is);
  v15(is) = (4*A0/pi/15) * sin( 15 * 2*pi*f0*t ) + v13(is);
end
  
% plot
subplot(4,2,1)
  plot( time, sqw, time, v1 )  
  title( '1 harmonic' )
  axis( [ time(1) time(121) -1.5*A0 1.5*A0 ] )
  axis off
subplot(4,2,3)
  plot( time, sqw, time, v3 )  
  title( '1+3 harmonics' )
  axis( [ time(1) time(121) -1.5*A0 1.5*A0 ] )
  axis off
subplot(4,2,5)
  plot( time, sqw, time, v5 )  
  title( '1+3+5 harmonics' )
  axis( [ time(1) time(121) -1.5*A0 1.5*A0 ] )
  axis off
subplot(4,2,7)
  plot( time, sqw, time, v7 )  
  title( '1+3+5+7 harmonics' )
  axis( [ time(1) time(121) -1.5*A0 1.5*A0 ] )
  axis off



subplot(4,2,2)
  plot( time, sqw, time, v9 )  
  title( '1+3+5+7+9 harmonics' )
  axis( [ time(1) time(121) -1.5*A0 1.5*A0 ] )
  axis off
subplot(4,2,4)
  plot( time, sqw, time, v11 )  
  title( '1+3+5+7+9+11 harmonics' )
  axis( [ time(1) time(121) -1.5*A0 1.5*A0 ] )
  axis off
subplot(4,2,6)
  plot( time, sqw, time, v13 )  
  title( '1+3+5+7+9+11+13 harmonics' )
  axis( [ time(1) time(121) -1.5*A0 1.5*A0 ] )
  axis off
subplot(4,2,8)
  plot( time, sqw, time, v15 )  
  title( '1+3+5+7+9+11+13+15 harmonics' )
  axis( [ time(1) time(121) -1.5*A0 1.5*A0 ] )
  axis off