function y = blitsaw(n, p, m)
% BLITSAW  A bandlimited sawtooth generator using closed-form, explicit
% calculation.
%
%   Y = BLITSAW(N,P,M) generates a bandlimited sawtooth signal of period
%   P (samples), length N, and with M harmonics (odd value less or
%   equal to P).
%  
%   The code below implements the BLIT algorithm of Stilson and
%   Smith, followed by a summation and filtering operation to produce
%   a sawtooth waveform.  After experimenting with various approaches
%   to calculate the average value of the BLIT over one period, I
%   found that an estimate of C2_ = 1.0 / period (in samples) worked
%   most consistently.  A "leaky integrator" is then applied to the
%   difference of the BLIT output and C2_.
%
%   By Gary P. Scavone, McGill University, 2005.

y = zeros(1, n);
c2_ = 1 / p;
rate_ = pi * c2_;
m_ = 2 * m + 1;
a_ = m_ / p;
phase_ = 0.0;
state_ = 0.0;

epsilon = 0.000001;
for c = 1:n
  den = sin( phase_ );
  if abs( den ) < epsilon
    y(c) = a_;
  else
    y(c) = sin( m_ * phase_ ) / ( p * den );
  end

  y(c) = y(c) + state_ - c2_;
  state_ = y(c) * 0.995;

  phase_ = phase_ + rate_;
  if ( phase_ >= pi ) phase_ = phase_ - pi; end
end
