function [y coeff] = square_wave(x,M,duty)
% SQUARE_WAVE
%
%      Finite Fourier Series evaluation of square wave
%
%   y = square_wave(x,M,width)
%
%      where x are the values to evaluate (-0.5<x<0.5)
%      M is the number of terms to use
%      duty is the duty cycle (0-1)
%
if (nargin<3)
   duty = 0.5;
end

p = (1:M)';
coeff = sinc(duty*p);
c = coeff*ones(size(x));
d = cos(2*pi*p*x);
y = duty+2.0*duty*sum(c.*d);
end