% sqseries.m
%
% The purpose of this program is to draw a the frequency
% components of a pulse train


% T. Scott Dattalo, 14FEB97

clear; clg;

T = 1;
N = 1000;                 %Number of samples
t = [1:N]/N * T;          %Time vector

%            0   4   8   c   0   4   8   c   0   4   8   c 
#bitstring = '00000001001101111110110010000000';
bitstring = '10010001110000111100010101101010';
[m,P] = size(bitstring);  % Number of pulses in the period of T

bits = zeros(1,P);
for m=1:P
  bits(m) = toascii(bitstring(m)) - toascii('0');
end


phase = [0:(P-1)] .* bits * T/P;
  
harmonics = 100;
mag = zeros(harmonics,1);
ph  = zeros(harmonics,1);

tau = T/P;
d = tau/T;
sqwave = zeros(size(t));

for n=1:harmonics
  w = n*2*pi/T;
  for m=1:P
    if bits(m) ~= 0
      sqwave = sqwave + 2*sin(n*pi*d)/(n*pi)*cos(w*(t-tau/2-phase(m)));
      mag(n) = mag(n) + 2*sin(n*pi*d)/(n*pi)*cos(w*phase(m));
      ph(n) = ph(n) + cos(w*phase(m));
    end
  end
end


ti = sprintf('SQUARE WAVE');
title(ti)
subplot(211)
plot(t,sqwave)
subplot(212)
plot([1:(harmonics)], mag,'*')