function y=blp(n,fbegin,fspace,dur,fs)
%BLP Band-limited pulse generator
%	y=BLP(n,fbegin,fspace,dur,fs) generates a band-limited pulse waveform consisting
%  of 'n' sinusoidal components beginning at frequency 'fbegin' Hz and spaced at
%  intervals of 'fspace' Hz. The duration of the waveform is 'dur' seconds, and
%  the sampling frequency is 'fs'.
%
%  Based on Moore p. 272.
 
%  Created 1998 Mar 4
%  Ed.Doering@Rose-Hulman.Edu
%  Department of Electrical and Computer Engineering
%  Rose-Hulman Institute of Technology
%  Terre Haute, IN 47803-3999
%
%  Copyright (c) 1998 by Ed Doering.
 

% Create the time vector (see NOTE below as to why the 1e-6 is there)
tt=(0:1/fs:dur)+1e-6;

% Create scaled time vectors  
theta=2*pi*fbegin*tt;
beta=2*pi*fspace*tt;

% Find B(n,theta)
% NOTE: Likely to get a 'divide by zero' error. Moore describes what to do about
% this (p. 272), but I have not implemented it here. So far I have only observed the
% problem for the very first element of beta (which is exactly zero), so I have 
% done the 'quick fix' by adding a small value to the time vector tt to eliminate the
% exact zero.
B=sin(n*beta/2)./(sin(beta/2));

% Finish generating the output waveform
y = sin(theta+((n-1)*beta)/2).*B;
y=y/max(y); %normalize to unit amplitude

%figure;plot(y)
return