function y = pulses(f, T, Fs)
   // Pulse train generator
   // This function creates a band-limited pulse train 
   // of given frequency and duration. It is equivalent 
   // to a series of Dirac deltas after being filtered
   // by an ideal filter with a cut off frequency Fs/2.
   //
   // Usage: 
   //        y = pulses(f, T, Fs)
   //
   // where
   //        f:  frequency in Hz
   //        T:  duration in s
   //        Fs: sampling rate
   //        y:  output filtered pulse train 
   // 
   // Output is equivalent to
   //
   //         N
   // y(k) = SUM cos(2*pi*n*k*f/Fs)
   //        n=1
   //
   // where N is the number of harmonics up to Fs/2, i.e., 
   // floor(Fs/2/f)
   //
   // Theory
   // This function uses the following equivalent short
   // formula, which is much more efficient:
   //
   //                           sin((N+1)*pi*k*f/Fs) 
   // y(k) = cos(N*pi*k*f/Fs) * -------------------- - 1
   //                              sin(pi*k*f/Fs)
   //
   // NOTE: In order to have the same energy as a pulse 
   //       train of frequency f it is necessary first to
   //       divide by the number of harmonics below Fs/2, 
   //       i.e., Fs/2/fo, and then by the number of
   //       samples within a period, i.e., Fs/fo
   //
   // Reference:
   //       Moore, F.R. (1990), "Elements of computer music". 
   //       Prentice Hall. Englewoog Cliffs. p. 271.
   //
   // -----------------------------------------
   // Author: Federico Miyara
   // Date:   2011-11-29
   //         2012-04-01
   //         2012-07-30
   //         2022-08-10

   // Time vector
   t = [0:1/Fs:T];

   // Number of harmonics to generate (harmonics up to 
   // Nyquist frequency
   M = floor(Fs/2/f);

   // Number of samples per period (may be non-integer)
   L = Fs/f;

   // Now we should apply the short formula for a sum of 
   // harmonics of equal amplitude. However, instead of 
   // applying directly the formula shown in the header
   // we will have to detect where sin(pi*f*t) = 0 to 
   // avoid 0/0's. But it is a little more complicated 
   // than that, since due to arithmetic errors sin(pi) 
   // is very small but not 0 so we have to capture
   // instants where sin(x) is too small and replace the 
   // quotient by its limit. A further difficulty is that
   // at alpha very close to pi sin(alpha) is so small
   // that arithmetic errors may cause the quotient 
   // sin(K*alpha)/sin(alpha) to be very different from K.

   u = sin(%pi*f*t);
   v = sin((M+1)*%pi*f*t);
   w = cos(M*%pi*f*t);

   // What follows is to prevent abnormal results due to 
   // errors in the computation of sin(x) when x is very
   // close to multiples of pi. If sin(x) is too small 
   // then sin((M+1)*x)/sin(x) is replaced by M+1 if M is 
   // even or by M+1 alternating with -(M+1) if M is odd. 
   // In the latter case the sign is according to a 
   // cosine of frequency f/2 
   toosmall = abs(u) < 0.000001;

   x = v./u;
   if 2*floor(M/2)==M
      x(find(toosmall)) = M+1;
   else
      x(find(toosmall)) = (M+1)*cos(%pi*f*t(find(toosmall)));
   end

   // Pulse train
   y = w .* x - 1;

   // Constant sqrt(2/K/L) is to normalize energy
   y = sqrt(2/M/L) * y;

endfunction

