% `magic'sine wave synthesis
%   The purpose of this program is to find an optimum pulse stream
% that when low pass filtered will produce a sine wave. The optimization
% criterion is to minimize the harmonics for a given number of pulses.


% by  T. Scott Dattalo 04MAR00
clear;
clg;
pulses = 16;

T = 1;              % Period
f = 1/T;            % fundamental frequency

tp = T/pulses;      % Time window for each pulse

t = [0:(pulses-1)]/pulses;
sine_samples = sin(2*pi*f * t);

ts = zeros(1,4*(pulses-1));
p = ts;

% Create the pulse width modulated sine wave.
% Floating point numbers are used to represent the location of
% the edges of the pulses.

for i=1:pulses

  ts(4*i) = t(i);
  ts(4*i+1) = t(i);
  ts(4*i+2) = t(i) + (1 + sine_samples(i))/2 * tp + 1000*eps;
  ts(4*i+3) = ts(4*i+2);
  p(4*i+0) = 0;
  p(4*i+1) = 1;
  p(4*i+2) = 1;
  p(4*i+3) = 0;

end

% quantize time
time_steps = 10000;
quant = time_steps/T;
tq = floor(ts * quant) / quant;



% Now calculate the harmonic content
harmonics = 20*pulses;
mag = zeros(1:harmonics);
g = mag;
ph  = zeros(1:harmonics);
phase = [0:(pulses-1)] * T/(pulses);
N = 1000;                 %Number of samples
sqwave = zeros(1:N);
u = [1:N]/N * T;          %Time vector

tenth = harmonics/10;
i=1;
for n=1:harmonics
  w = n*2*pi/T;

  for m=1:(pulses)
    tau = ts(4*m+2) - tq(4*m+1);
    phase = ts(4*m+1);
    d = tau/T;
    Rn = 2*sin(n*pi*d)/(n*pi);
    %sqwave = sqwave + 2*sin(n*pi*d)/(n*pi)*cos(w*(u-tau/2-phase));
    %mag(n) = mag(n) + 2*sin(n*pi*d)/(n*pi)*cos(w*phase);
    sqwave = sqwave + Rn*cos(w*(u-tau/2-phase));
    g(n) = g(n) + Rn*exp(j*w*phase);

    %ph(n) = ph(n) + cos(w*phase);
  end

  if (n>i*tenth)
    printf("%d\n",i);
    i = i + 1;
  end

end

sqwave = sqwave + 0.5 ;
dc = sum(sine_samples);
mag(2) = 0;
ti = sprintf('PWM Synthesis');
title(ti)
subplot(111)
%plot(t,sine_samples)
%subplot(212)
plot([0:harmonics-1]/pulses,abs(g),'*')

% plot the pwm wave three different ways:
%  1 - using floating point numbers to represent edge locations
%  2 - quantized version of (1)
%  3 - Harmonic reconstruction of (2)
%plot(ts,p,u,sqwave,tq,p)