function y = blit(n, p, m)
% BLIT  A bandlimited impulse generator using closed-form, explicit
% calculation.
%
%   Y = BLIT(N,P,M) generates a bandlimited impulse train of period
%   P (samples), length N, and with M harmonics (odd value less or
%   equal to P).  The output is normalized by scaling by P / M;
%
%   By Gary P. Scavone, McGill University, 2005.

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

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

  phase_ = phase_ + rate_;
  if ( phase_ >= pi ), phase_ = phase_ - pi; end
  % The phase range can be 0 to pi because the result from pi to 2*pi is
  % the same ... the negative signs cancel in sin / sin.
  %if ( phase_ >= 2*pi ), phase_ = phase_ - 2*pi; end
end
