I thought up another way to apply angle modulation to a synthesized sinusoid than
by calculating sines or cosines (or doing table lookups). From a bit of searching
it seems that this might be novel (then again it might just be that this already
exists with a name or acronym I don't know!). It might prove an interesting
variation on standard FM synthesis. I am working on an AudioStream implementation,
using floats. Watch this space!

There are two key realizations involved, one is that by using a complex-valued
signal you can use geometry in the complex plane to distort angles directly without
needing to explicitly track or represent the phase.

Such modulation necessarily can only modulate at the same frequency
as the waveform being modulated, as each complete turn (1 cycle) puts the complex
plane back where it started. All you can do is some stretching and squashing of the
waveform (but that's probably all you need).

Unit complex-valued tones run anti-clockwise around the unit-circle, mathematically
described by e^(jwt). A complex multiplication takes two such signals and adds
their phases (without ever making the phase explicit - only 6 floating point adds
and multiplies are needed to perform complex multiplication, no transcendental
operations or table lookups). This property of avoiding having to map between phase
and signal and back is part of the power of complex numbers.

The second realization is that once you are using a complex signal, shifting up and
down in frequency is just multiplication, and thus one oscillator can easily generate
harmonics on the fly cheaply.

This means you can modulate a signal of frequency 2f (with 2f modulation), shift
it down to f, and you have effectively modulated f with 2f, despite the limitation
mentioned above. You can then modulate the modulated signal again at frequency f.
You can even shift it down again to DC, giving a signal starved of the fundamental.

No aliasing happens in these frequency shifts, another powerful property of
complex signals.

The most convenient operations are frequency doubling, and shifts by the fundamental
which you have already generated. These are each a single complex multiplication or
a squaring. If you store some of these derived frequencies in variables you can then
use them to shift by larger amounts.

A subharmonic would need to be computed first, then doubled/shifted to derive the
notional fundamental.

The geometry operation that seems most straight-forward for angle distortion
to my mind is this (pseudocode):

z = z - o       ; shift unit circle by some offset (a complex value itself, note)
z = z / abs(z)  ; normalize back to the origin, warping angles.
What this does is take the notionally rotating signal (on the unit circle), and shift
the circle. Then the normalization step distorts the angles of the signal by forcing
the origin back to where is should be without moving the circle, then normalizing the
radii back to unit length. The larger the magnitude of the offset the more distortion, and
if its more than unit length the angles become limited to a fraction of a half-circle.

Click image for larger version. 

Name:	CSMod.png 
Views:	8 
Size:	74.0 KB 
ID:	21316

In full, pseudo-coded, using real components, where z = x+iy:
x -= ox           ; do the offsetting
y -= oy

mag2 = x*x + y*y  ; calculate abs(z)
mag = sqrt(mag2)

x /= mag          ; normalize
y /= mag
In fact the square-root operation doesn't need to be done to produce distortion of angles,
although it is needed to keep the magnitude of the resulting waveform tamed. Omitting it
strategically can save cycles, especially if sqrt isn't a hardware operation.

If you string a bunch of the angle-distorting and freqency-shifting steps together, the last
normalization should probably be a proper one using sqrt.

The final, real-valued, result is taken as the imaginary part at the end, i.e. the value y.

This is convenient if you confine the offsets to be real-valued (oy == 0), since the
output will then avoid a DC-component, as there is symmetry about the real axis - flip
vertically and reverse time and you get the same result.

Just as with FM synthesis stringing operations in series allows complex mixes of harmonics
to be generated, but the mathematics is different (and probably much more tricky).

Through experimentation with such strings of operators with a variety of frequency changing
steps its clear that a rich variety of partials ought to be possible:

Click image for larger version. 

Name:	CSMod_waveforms_spectra.png 
Views:	3 
Size:	185.4 KB 
ID:	21317
The top plot shows the variation in harmonic levels as the modulation depth (= offset sizes)
is varied - clearly rich structure, and the fundamental (black) drops in places as do most
harmonics (even = solid, odd = dashed). There was a string of 6 operators and various
frequency shifts in this experiment.

[Some example waveforms at different modulation depths are given in the bottom plot]

In processing a block of samples you can use complex multiplication to implement the
actual main oscillator from sample to sample, using sin/cos only at the start of the block.
Too many repeated complex multiplications might propagate rounding errors, so its sensible
to reset the value per block for long term stability.