function buildsq(speed)
% buildsq(speed)   Build up a square wave by adding terms of Fourier series
% 2004-09-28 Dan Ellis dpwe@ee.columbia.edu

if nargin < 1
  speed = 0;
end

% Time extent (normalized by pi)
tmax = 1.5;
tstep = 0.01;
nn = -tmax:tstep:tmax;

% Vertical scale limit for displays
ymax = 1.5;

% The function itself - starts at zero, accumulates terms
yy = zeros(1, length(nn));

yold = yy;

% Loop around adding odd harmonics, up to 49th

for k = 1:2:49
  
  % Sign and scale of this harmonic
  ak = ( (-1)^((k-1)/2) ) * 1/k;
  
  % The new term
  xx = ak*cos(pi*k*nn);
  
  % Add it in
  yy = yy + xx;
  
  % Display them both and scale
  subplot(211)
  plot(nn,xx);
  % Manually set axis range to avoid autoscaling
  axis([-tmax tmax -ymax ymax]);
  grid on
  title([num2str(k),'th harmonic']);
  
  % Make sure double buffering is on for this figure
  set(gcf, 'DoubleBuffer', 'On');
  
  subplot(212)
  if speed > 0
    % Fade in new term
    for w = 0:(1/speed):1
      plot(nn,yold,'c',nn,yold+w*xx,'b');
      axis([-tmax tmax -ymax ymax]);
      grid on
      drawnow;
    end
  else
    % Plot this and previous versions
    plot(nn,yold,'c',nn,yy,'b');
    grid on
  end
  
  title('Composite square wave approximation');

  pause
  
  yold = yy;
  
end