% m


% T. Scott Dattalo, 06MAR00

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 = '010110100000';
%bitstring =  '111111000000';
%bitstring =  '11110000';

[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)) - 48;
end

phase = [0:(P-1)] .* bits * T/P;
phi = [0:(P-1)] .* bits * T/P;

harmonics = 20*P;
g = zeros(size(1:harmonics));
ph  = zeros(size(1:harmonics));

Tb = T/P;
%d = tau/T;
tr = Tb/1000;
tau = Tb - tr;
tau = Tb;
sqwave = zeros(size(1:N))+ sum(bits)*(tau+tr)/T;


for n=1:harmonics
  w = n*2*pi/T;
  Rn = 2*sin(n*pi*tau)/(n*pi);
  %Rn = ( -cos(w*(tau/2+tr))/tr + cos(w*tau/2)/tr)/(w*w*T/2);

  for m=1:P
    if bits(m) ~= 0
      sqwave = sqwave + Rn*cos(w*(t-tau/2-phase(m)));
      g(n) = g(n) + Rn*exp(j*w*phase(m));
    end
  end
end
ti = sprintf('Square Wave, Bit Stream: %s',bitstring);
subplot(211)
title(ti)
plot(t,sqwave)
subplot(212)
title('Harmonics')
h = harmonics/P
%plot([1:h]/P, abs(g(1:h)),'*')
plot( abs(g(1:h)),'*')