Well, at least partially...

It finally dawned on me yesterday how magic sine waves can be sythesized.
What's frustrating, is that the answer is so simple that it's obvious; so
obvious that it's been overlooked in the discussions we had off and on for
the last couple of years (at least by me). Lancaster is probably a
chuckling lurker :).

There are probably a few ways to go about it, but this is the one I tried.
Imagine how you'd go about synthesizing a sine wave if you had an A/D
converter. The most common way is to produce `N' samples for one sine wave
cycle and repeatedly write these to an A/D converter. The sine wave
frequency is determined by the rate at which the samples are written to
the A/D converter.

Now imagine replacing the A/D converter with a PWM generator that has the
same resolution. So for example, if you had a 12-bit A/D converter the PWM
generator would be capable of producing 4095 (2^12 - 1) different duty
cycles. Theoretically, the low-pass filtered output of the A/D converter
and the PWM are the same. Practically however, it's difficult to design a
filter that would have the proper cut-off characteristics. So to ease the
filter requirements, the PWM synthesized sine wave is usually made to be a
much higher frequency than the frequency at which samples are written to
the A/D converter.

Now when we, or at least I, begin to think about `magic sine waves' the
idea that the pwm frequency can be reduced is forgotten. Furthermore, we
constrain ourselves to the low-resolution 8 or 10 bit pwm generators that
are available in the various pic's. What the magic sine wave algorithm
does, AFAICT, is find a suitable time discretization for a high
resolution PWM synthesized sine wave.

So, there are two parts to the synthesis. Part 1 invovles finding the high
resolution PWM representation of the sine wave. Part 2 involves finding an
optimum way to make time discrete. The first part is easy and the attached
octave (or matlab) program will do this. The second part can be determined
using brute force search techniques. I don't have an elegant solution
though, but hopefully some others have some insight.


Scott

PS, this octave program is meant to be tool and otherwise not too useful
(in other words, I don't bother to explain hardly anything...)

% `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 = 32;

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 a pwm sine wave with edges based on floating point
% numbers:

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

% Now calculate the harmonic content
harmonics = 20*pulses;
mag = zeros(1:harmonics);
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

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

  m=2;
  for m=1:(pulses)
    tau = (1 + sine_samples(m))/2*tp;
    d = tau/T;
    sqwave = sqwave + 2*sin(n*pi*d)/(n*pi)*cos(w*(u-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

sqwave = sqwave + 0.5 + 2;
dc = sum(sine_samples);
mag(2) = 0;
ti = sprintf('PWM Synthesis');
title(ti)
subplot(111)
%plot(t,sine_samples)
%subplot(212)
% plot the harmonics:
plot([0:harmonics-1]/pulses,abs(mag))
% uncomment to plot the square wave train re-synthesized from its
% harmonics
%plot(ts,p,u,sqwave)