AN INTRODUCTION TO FM

		Bill Schottstaedt

In frequency modulation we modulate the frequency Ð "modulation" hereis just a latinate word for "change".  Vibrato and glissando arefrequency modulation.  John Chowning tells me that he stumbled on FM when he sped up vibrato to the point that  it was creating audible sidebands (perceived as a timbral change) rather than faster warbling (perceived as a frequency change).  We can express this (the vibrato, not the neat story) as:

    cos(wct + f(t))
 
where the c subscript stands for "carrier" and f(t) means "some arbitrary function added to the carrier".  Since cos takes an angle as its argument, f(t) "modulates" (that is, changes) the angle passed to the cosine.  Hence the generic name "angle modulation".  Phase and angle seem like synonyms, hence the name "phase modulation".  Since we can always translate from instantaneous radial frequency 

    wct + f(t) 
 
evaluated at a particular time to the corresponding frequency (differentiate and divide by 2 pi), this same fomula can viewed as "frequency modulation".  Textbooks usually mutter "strictly speaking we can only take the cosine of an angle", and they put an explicit integral in when they mean frequency modulation:

    cos(wct + B ò0t f(t) dt)   
			(please forgive the odd layout of these formulas -- Edit is 
			inadequate in this regard).
 
The terms "angle modulation", "phase modulation", and "frequency modulation" (not to mention "exponential modulation" and "angular modulation") are used almost interchangeably in the radio literature, but the textbooks on the subject normally distinguish between PM and FM as follows: if f(t) is directly proportional to the modulating signal then it's PM; if f(t) is directly proportional to the derivative of the modulating signal, then it's FM.  Of course, you can't tell which is in use either from the waveform or the mathematical expression of the waveform Ð you have to know what the modulating signal was.  That is a roundabout way of saying that in computer music applications there is no essential difference between frequency and phase modulation.  However, the fact that phase and frequency modulation are just two ways of interpreting the same thing should not mislead us into thinking "phase doesn't matter" Ð the initial phases of the parts of the FM formula do make a difference, as we'll see below. 

Just to beat this dead horse very thoroughly, here are two clm instruments, one performing textbook phase modulation, the other performing frequency modulation.  I have tried to make the innards explicit at each step, and match the indices so that the instruments produce identical results given identical parameters.  Also, to lay a different controversy to rest, it should be obvious from these two functions that there is no difference in run-time computational expense between PM and FM.

(definstrument pm (beg end freq amp mc-ratio index)
  (let ((carrier-phase 0.0)
        (carrier-phase-incr (* freq (/ two-pi sampling-rate)))
        (modulator-phase 0.0)
        (modulator-phase-incr (* freq mc-ratio (/ two-pi sampling-rate)))
        (modulation 0.0)
        (val 0.0))
   (Run
     (loop for i from beg to end do
       (setf modulation (* index (sin modulator-phase)))
       (setf val (* amp (sin (+ carrier-phase modulation))))
       ;; no integration in phase modulation
       (incf carrier-phase carrier-phase-incr)
       (incf modulator-phase modulator-phase-incr)
       (outa i val)))))

(definstrument fm (beg end freq amp mc-ratio index)
  (let ((carrier-phase 0.0)
        (carrier-phase-incr (* freq (/ two-pi sampling-rate)))
        (modulator-phase (* .5 pi))
        (modulator-phase-incr (* freq mc-ratio (/ two-pi sampling-rate)))
        (fm-index (* freq mc-ratio index (/ two-pi sampling-rate)))
         ;; ^ fix up fm index to take integration into account
        (val 0.0)
        (modulation 0.0))
  (Run
     (loop for i from beg to end do
       (setf modulation (* fm-index (sin modulator-phase)))
       (incf carrier-phase modulation)
       ;; here is the fm integration ^
       (setf val (* amp (sin carrier-phase)))
       (incf carrier-phase carrier-phase-incr)
       (incf modulator-phase modulator-phase-incr)
       (outa i val)))))

(with-sound () (pm 0 10000 1000 .25 0.1 4))
(with-sound () (fm 0 10000 1000 .25 0.1 4))

or (to use Chowning's example):

(with-sound () (fm 0 10000 100 .25 1.0 4.0))
(with-sound () (pm 0 10000 100 .25 1.0 4.0))

Given our formula for FM, let's assume, for starters, that f(t) is a sinusoid:

    cos(wct + B sin wmt)

where the "m" stands for "modulator" and the B is usually called the "modulation index".  In clm syntax this is:

    (oscil carrier (* B (oscil modulator)))

where oscil is:

    (defun oscil (oscillator fm-input &optional (pm-input 0))
      (prog1 
        (sin (+ (osc-phase oscillator) pm-input))
        (incf (osc-phase oscillator) 
          (+ (osc-freq oscillator) fm-input))))

Since it is generally believed that the ear performs some sort of projection of the time domain waveform into the frequency domain (a Fourier Transform), and that timbre is at least partly a matter of the mix of frequencies present (the spectrum), our main interest in the FM formula is in the spectrum it produces.  To determine that spectrum, we have to endure some tedious mathematics.  By the trigonometric identity:

    cos(a + b) = cos a cos b - sin a sin b

we can substitute wct for "a" and B sin wmt for "b" and get:

    cos(wct + B sin wmt) =
        cos wct cos(B sin wmt) - sin wct sin(B sin wmt)

If we can get a Fourier transform of the two inner portions:

    cos(B sin wmt)  and   sin(B sin wmt)

we can use:

    cos A cos B = 1/2 (cos (A - B) + cos (A + B))

and

    sin A sin B = 1/2 (cos (A - B) - cos (A + B))

to get the final results.  "A" here is wct in the earlier formulas,  and "B" is either  
cos(B sin wmt) or sin(B sin wmt).

You can slog through the Fourier series expansion by hand, or look it up in Abramowitz and Stegun, "Handbook of Mathematical Functions" (formulas 9.1.42 and 9.1.43), or take my word for it that:

    cos(B sin wmt) = 
        J0(B) + 2J2(B)cos 2wmt + ... 
         + 2J2n(B)cos 2nwmt + ...

    sin(B sin wmt) = 
         2J1(B)sin wmt + 2J3(B)sin 3wmt + ... 
          + 2J2n-1(B)sin(2n-1)wmt ...

Here the J's refer to the Bessel functions which we will return to later. First, let's finish this expansion Ð we take these two sums and wct and plug them into our first expansion of the FM formula, and out pops:

    cos(wct + B sin wmt) =
       J0(B) cos wct  
          - J1(B)[cos(wc - wm)t - cos (wc + wm)t]
          + J2(B)[cos(wc - 2wm)t + cos(wc + 2wm)t]
          - J3(B)[cos(wc - 3wm)t - cos(wc + 3wm)t] + ...

or in a slightly more compact form:

    ¥
    S [Jn(B)cos(wc + nwm)t]
   n=-¥

Here we are using the fact that

    J-n(x) = (-1)n Jn(x)

We can change our point of view on the first part of the expansion given above, and ask for the amplitude of a given sideband before it is translated up to the carrier (by multiplication):

    Jn(B) = 2/p ò0p/2sin(B sin w)sin nw dw   (n odd)

    Jn(B) = 2/p ò0p/2cos(B sin w)cos nw dw   (n even)

So we end up with a spectrum made up of a "carrier" at wc and   symmetrically placed sidebands separated by wm.  The amplitudes  follow the Bessel functions.  I put carrier in quotes because in computer music we listen to the result of the modulation (this was Chowning's idea Ð see "The Synthesis of Complex Audio Spectra by Means of Frequency Modulation").  The Bessel functions are nearly 0 until the index (B) equals the order (n).  Thereupon they have a bump and then tail off in a sort of damped sinusoidal way (see Chowning's article for some graphs).  So as the index sweeps upward, energy is swept gradually outward into higher order side bands Ð this is the originally exciting, now extremely annoying "FM sweep".  It is worrisome that the simple case has so quickly become jejune.  The only important thing to get from these functions is that the higher the index, the more dispersed the spectral energy Ð normally a "brighter" sound.

There is a rule of thumb, Mr Carson's rule, about the overall bandwidth of the resultant spectrum (it follows from our description of the Bessel functions): Roughly speaking, there are fm-index+1 significant sidebands on each side of the carrier, so our total bandwidth is on the order of 

    2 * mc-ratio * carrier * (fm-index + 1)

This is a good approximation Ð 99% of the signal power is within its limits.

Here's a simple FM instrument:

(definstrument fm (beg end freq amp mc-ratio index 
                   &optional (index-env '(0 1 100 1)))
  (let* ((cr (make-oscil :frequency freq))
         (md (make-oscil :frequency (* freq mc-ratio)))
         (fm-index (* index mc-ratio freq (/ two-pi sampling-rate)))
         (ampf (make-env :envelope index-env :scaler amp 
                         :start beg :end end))
         (indf (make-env :envelope index-env :scaler fm-index 
                         :start beg :end end)))
    (Run
      (loop for i from beg to end do
        (outa i (* (env ampf) 
                   (oscil cr (* (env indf)
                                (oscil md)))))))))

I put an envelope on the fm-index so that we can try out the "dynamic spectra" ("dynamic" means "changing" here).  For now, don't worry too much about the actual side band amplitudes Ð these will not always match Chowning's description Ð we'll get around to an explanation eventually.

(with-sound () (fm 0 10000 100 .5 1.0 4.0))

is Chowning's first example.  Sure enough, it's a complex spectrum (that is, it has lots of components Ð try an index of 0 to hear a sine wave, if you're suspicious).  Since our modulating frequency to carrier frequency ratio (mc-ratio above) is 1.0, we get sidebands at harmonics of the carrier. If we use an mc-ratio of .25 and a carrier of 400:

(with-sound () (fm 0 10000 400 .5 0.25 4.0))

we end up with the same perceived pitch because the sidebands are still at multiples of 100 Hz.

(with-sound () (fm 0 10000 400 .5 1.1414 4.0))

has inharmonic sidebands.

Most real sounds seem to change over the course of a note, and it was at one time thought that most of this change was spectral.  To get a changing spectrum, we need only put an envelope on the fm-index:

(with-sound () (fm 0 10000 400 .5 1.0 5.0 
                    '(0 0 20 1 40 .6 90 .5 100 0)))

making a "brass-like" sound.  Similarly, Chowning suggests that

(with-sound () (fm 0 10000 900 .5 1/3 2.0 
                    '(0 0 6 .5 10 1 90 1 100 0)))

is a woodwind-like tone,

(with-sound () (fm 0 10000 500 .5 .2 1.5 
                    '(0 0 6 .5 10 1 90 1 100 0)))

is bassoon-like, and finally

(with-sound () (fm 0 10000 900 .5 2/3 2 
                    '(0 0 25 1 75 1 100 0)))

is clarinet-like.  Many simple extensions can be made to the fm instrument given above and something good can often be coaxed from it.  But, to clear up one source of confusion Ð if you looked at the spectrum of our first example, and compared it to the spectrum Chowning works out, you may wonder what's gone awry.  We have to return to our initial set of formulas.  If we consider that:

    sin(a + b) = sin a cos b + cos a sin b

so

    sin(wct + B sin wmt) =
       sin wct cos(B sin wmt) + cos wct sin(B sin wmt)

and using our previous formulas for the expansion of the cos(sin) and sin(sin) terms, with the elementary identity:

    sin a cos b = 1/2 (sin(a - b) + sin(a + b))

we see that we still have a spectrum symmetric around the carrier, and the amplitude and frequencies are just as they were before, but the initial phases of the side bands have changed. Our result is now 

    sin(wct + B sin wmt) =
       J0(B)sin wct
          + J1(B)[sin(wc + wm)t - sin(wc - wm)t]
          + J2(B)[sin(wc + 2wm)t + sin(wc - 2wm)t]
          + J3(B)[sin(wc + 3wm)t - sin(wc - 3wm)t]
          + ...

This is Chowning's version of the expansion.  Here we are also using:

    sin(-x) = -sin x  

(i.e sin(wmt - wct) = -sin(wct - wmt))

In general:

    cos(wct + B sin(wmt + qm) + y) =

       ¥
         S Jk(B)cos((wc + kwm)t + kqm + y)
       k=-¥

Our first reaction is, "so what if one's a sine and another's a cosine, or whatever Ð they'll sound the same", but we are being hasty.  What if (for example), the modulator has the same frequency as the carrier, and its index (B) is high enough that some significant energy appears at

    wc - 2wm = - wc?

Where does energy at a negative frequency go?  We once again fall back on elementary trigonometry:

    sin(-x) = -sin x

but

    cos(-x) = cos x

so the negative frequency component adds to the positive frequency component if it's a cosine, but subtracts if it's a sine.  The upshot of all this is that the results will be different Ð we get a different pattern of cancellations depending on the initial phases of the carrier and modulator.  Take the clm instrument:

(definstrument fm (beg end freq amp mc-ratio index 
                   carrier-phase mod-phase)
  (let ((cr (make-oscil :frequency freq 
                        :initial-phase carrier-phase))
        (md (make-oscil :frequency (* freq mc-ratio) 
                        :initial-phase mod-phase))
        (fm-index (* index mc-ratio freq (/ two-pi sampling-rate))))
    (Run
      (loop for i from beg to end do
        (outa i (* amp (oscil cr (* fm-index (oscil md)))))))))

and compare:

(with-sound () (fm 0 10000 100 .25 1.0 4 0 (* .5 pi)))

with

(with-sound () (fm 0 10000 100 .25 1.0 4.0 (* .5 pi) (* .5 pi)))

There is a noticeable difference!  The first case is that given by Chowning in his fig 1.4a, the second is Bate's fig 3 ("The Effects of Modulator Phase on Timbres in FM Synthesis").  (At least I think it is Ð Bate says he is changing the modulator's initial phase, but I can't reproduce his graph unless I change the carrier initial phase Ð the important point is that these phases matter).  Also, remember that we are thinking in terms of the pre-integration phase here (hence the (* .5 pi)) Ð

    sin(x + p/2) = cos x

and:

    ò cos x dx = sin x

so by starting with a cosine in the modulator, after the integration in oscil (and ignoring the constant for the moment) we have

    sin(wct + B sin wmt)

As Bate points out, by varying the relative phase (by very slowly changing the modulator frequency independently of the carrier), we get a changing spectrum because of these cancellations.  Here is a clm instrument that shows this effect:

(definstrument fm (beg end freq amp mc-ratio index 
                   carrier-phase mod-phase skew-func skew-amount)
  (let ((cr (make-oscil :frequency freq :initial-phase carrier-phase))
        (md (make-oscil :frequency (* freq mc-ratio) 
                        :initial-phase mod-phase))
        (skewf (make-env :envelope skew-func 
                         :start beg :end end
                         :scaler (* skew-amount mc-ratio freq 
                                    (/ two-pi sampling-rate))))
        (fm-index (* index mc-ratio freq (/ two-pi sampling-rate))))
    (Run
      (loop for i from beg to end do
        (outa i (* amp (oscil cr (* fm-index 
                                    (oscil md (env skewf))))))))))

(with-sound () (fm 0 40000 100 .25 1.0 4.0 0 0 
                   '(0 0 50 1 100 0) .02))

One subtle pitfall here, pointed out by Holm (1992), involves the constant offset caused by the integration.  If we integrate f(x) we get some F(x) + F(a) where a is the starting point of the integration.  In the case of sin(B*sin), a is 0, and F(x) is -cos(x), so the constant term is B (i.e. -B*(-cos(0)).  That is, we get a phase offset in the outer sin calculation equal (in this case) to the modulation index!   So the true carrier initial phase is a combination of the initial phase we specify and the modulation index times the initial value of the integrated wave.

The next question is obviously, "if we can get cancellations, can we twiddle the phases and get single sideband FM?".  We could then get away from the symmetric spectra FM always tries to produce, and edge toward phase quadrature waveshaping.  And as far as I can tell, the answer is, "yes, but..."  We can get good results with low indices if we take advantage of the fact that:

    cos A cos B = 1/2 (cos(A - B) + cos(A + B))

and

    sin A sin B = 1/2 (cos(A - B) - cos(A + B))

If we have a spectrum B made up entirely of sines (or entirely cosines), we can then multiply it by sin A (or cos A) then add the two resulting spectra, and the (A + B) parts cancel.  Here's an example:

(definstrument fm (beg end freq amp mc-ratio index cr0p cr1p md0p md1p)
  (let ((cr0 (make-oscil :frequency (* freq mc-ratio) 
                         :initial-phase cr0p))
        (cr1 (make-oscil :frequency (* freq mc-ratio) 
                         :initial-phase cr1p))
        (md0 (make-oscil :frequency (* freq mc-ratio) 
                         :initial-phase md0p))
        (md1 (make-oscil :frequency (* freq mc-ratio) 
                         :initial-phase md1p))
        (am0 (make-oscil :frequency freq :initial-phase 0))
        (am1 (make-oscil :frequency freq :initial-phase (* .5 pi)))
        (fm-index (* index mc-ratio freq (/ two-pi sampling-rate))))
    (Run
      (loop for i from beg to end do
        (outa i (* amp (+ (* (oscil am0) 
                             (oscil cr0 (* fm-index (oscil md0))))
                          (* (oscil am1) 
                             (oscil cr1 (* fm-index (oscil md1)))))))))))

(with-sound () (fm 0 10000 1000 .25 .1 1.0 
                   0 (* .5 pi) (* .5 pi) (* .5 pi)))

(cr0p = 0 for sin in sin(sin), cr1p=(* .5 pi) for cos in cos(sin), md0p and md1p=(* .5 pi) so that both will be sin upon integration) Ð we use cos*cos(sin)+sin*sin(sin) to cancel the upper sideband).

The cancellation is nearly perfect in:

(with-sound () (fm 0 10000 1000 .25 .1 .1 
                   0 (* .5 pi) (* .5 pi) (* .5 pi)))

Palamin et al in "A Method of Generating and Controlling Musical Asymmetrical Spectra" started at the same point (the expansion of simple FM), added a factor rn to each sideband, and then worked backwards to get the exponential modulation that would produce that spectrum:

    e(B/2(r-1/r)cos wmt) sin(wct + B/2(r+1/r)sin wmt)

         = S rnJn(B)sin(wct + nwmt)

This idea was mentioned by Moorer in "Signal Processing Aspects of Computer Music: A Survey", but he didn't go into much detail.  It does give very pretty asymmetric spectra, but I believe waveshaping does as well for less computation, is more versatile, and is easier to set up and control.

So far we've been using just a sinusoid for the modulator Ð what if we make it a more complicated signal?  Here again elementary trigonometry can be used to expand

    sin(wct + B1sin wm1t + B2sin wm2t)

The modulating signal is now made up of two sinusoids (don't despair; this is a terminating sequence).  Since sine is not linear, this is not the same thing as

    sin(wct + B1sin wm1t) + sin B2sin wm2t

In the second case we just get the superposition of the two simple FM spectra, but in the first case we get a more complex mixture involving all the sums and differences of wm1and wm2.  If the indices are small (the Bi's are much less than 1), then there isn't much difference between the two versions.  Once again, the expansion is tedious, but it really involves nothing more than putting the sum of the sines in place of the single sine in the previous case, then turning the algebraic crank.  See Le Brun "A Derivation of the Spectrum of FM with a Complex Modulating Wave".  The result can be expressed:

        S S Ji(B1) Jk(B2) sin(wct + iwm1t + kwm2t)
        i  k

You can chew up any amount of free time calculating the resulting side band amplitudes Ð see Schottstaedt "The Simulation of Natural Instrument Tones Using Frequency Modulation with a Complex Modulating Wave".  These extra modulating components flatten and spread out the spectrum somewhat.  In general:

              k
    cos(wct + S Bisin(wit + qi) + y) =
              i=1

                  k                   k
         S...S  (P Jki(Bi)) cos(wct + S ki(wit + qi) + y)
        kk   k1   i=1                  i=1

A clm instrument to produce something like this might be:

(definstrument fm (beg end freq amp mc-ratios indexes 
                   carrier-phase mod-phases)
  (let* ((cr (make-oscil :frequency freq :initial-phase carrier-phase))
         (n (length mc-ratios))
         (md-arr (make-array n :element-type 'osc))
         (fm-index-arr (make-array n :element-type 'short-float))
         (sum 0.0))
    (loop for i from 0 below n do
      (setf (aref md-arr i) 
               (make-oscil :frequency (* freq (aref mc-ratios i)) 
                           :initial-phase (aref mod-phases i)))
      (setf (aref fm-ind-arr i) (* (aref indexes i)
                                   (aref mc-ratios i)
                                   freq (/ two-pi sampling-rate))))
    (Run
      (loop for i from beg to end do
        (setf sum 0.0)
        (dotimes (k n) 
           (incf sum (* (aref fm-ind-arr k) 
                        (oscil (aref md-arr k)))))
        (outa i (* amp (oscil cr sum)))))))

But this is silly Ð I don't think all these modulators are buying us anything.  I say three is enough.  My favorite computer instrument, the FM violin, uses three sinusoidal components in the modulating wave Ð see v.ins in the clm area Ð and for more complex spectra these violins are then ganged together (see fmviolin.clm for many examples).  By using a few sines in the modulator, you get away from the "simple FM sweep" that has become anathema, and the broader, flatter spectrum is somewhat closer to that of a real violin.

The lumpy nature of the simple FM spectrum can be seen in Chowning's figure 1.8, but I find it easier to read the graph of Jn(10) as a function of n in Abramowitz and Stegun fig 9.3 Ð think of the n axis as frequency, and remember that we only  "see" n when it is an integer Ð otherwise the negative part of the graph might give us pause.  The main lump is near n=B for all B, a most unfortunate fact if we are trying to mimic the spectra of normal musical instruments.  Long time FM'ers have tied themselves in elaborate knots trying to palliate this sore thumb of a bump; parallel modulators are merely a simple and obvious bandage.

The fm-violin in v.ins is complicated (a working instrument rapidly becomes encrusted with hard-to-explain features).  A pared down version is:

(definstrument violin (beg end frequency amplitude fm-index)
  (let* ((frq-scl (in-Hz frequency))
         (maxdev (* frq-scl fm-index))
         (index1 (* maxdev (/ 5.0 (log frequency))))
         (index2 (* maxdev 3.0 (/ (- 8.5 (log frequency)) 
                                  (+ 3.0 (/ frequency 1000)))))
         (index3 (* maxdev (/ 4.0 (sqrt frequency))))
         (carrier (make-oscil :frequency frequency))
         (fmosc1 (make-oscil :frequency frequency))
         (fmosc2 (make-oscil :frequency (* 3 frequency)))
         (fmosc3 (make-oscil :frequency (* 4 frequency)))
         (ampf  (make-env :envelope '(0 0 25 1 75 1 100 0) 
                          :scaler amplitude :start beg :end end))
         (indf1 (make-env :envelope '(0 1 25 .4 75 .6 100 0) 
                          :scaler index1 :start beg :end end))
         (indf2 (make-env :envelope '(0 1 25 .4 75 .6 100 0) 
                          :scaler index2 :start beg :end end))
         (indf3 (make-env :envelope '(0 1 25 .4 75 .6 100 0) 
                          :scaler index3 :start beg :end end))
         (pervib (make-triangle-wave :frequency 5 
                                     :amplitude (* .0025 frq-scl)))
         (ranvib (make-randi :frequency 16 :amplitude (* .005 frq-scl)))
         (vib 0.0))
    (Run
     (loop for i from beg to end do
       (setf vib (+ (triangle-wave pervib) (randi ranvib)))
       (outa i (* (env ampf)
                  (oscil carrier
                         (+ vib 
                            (* (env indf1) (oscil fmosc1 vib))
                            (* (env indf2) (oscil fmosc2 (* 3.0 vib)))
                            (* (env indf3) (oscil fmosc3 (* 4.0 vib)))
                          )))))))))

(with-sound () (violin 0 10000 440 .1 2.5))

There is an extension of Carson's rule to the case of complex modulation. The quick-and-dirty form says "take the highest wm and use its index"  The fancier version redefines the modulation index to be
                                        K
    B = sqrt(1/K S Bi2/2) / wK
                 i=1

where Bi is the index of the i-th component and wK is the highest modulating frequency and K is the number of components present.  When the modulating signal gets very complex, you might use:

    B = peak-deviation * sqrt(power-of-modulating-signal)/ 
           highest-modulating-frequency

but our rule of thumb is rapidly getting out of hand, as it were.

We could, of course, use FM (or anything) to produce the modulating signal.  When FM is used, it is sometimes called "cascade FM":

    sin(wct + B1sin(wm1t + B2sin wm2t))

In clm nomenclature:

    (* A (oscil carrier 
                (* B (oscil modulator 
                            (* C (oscil cascade))))))

and we can get an extremely complex spectrum in no time.  Unfortunately, FM can produce energy at 0Hz (when wc = nwm for example), and that DC component becomes a constant offset in the lower of the two FM pairs. So your carrier no longer has any obvious relation to wc!  That is, we can expand our cascade formula into:

    sin(wct + S [Jn(B)cos(wc + nwmn)t])

but now wherever wc= -nwmn we get Jn(B) cos 0 = Jn(B).  This could be a disaster, because in most cases where we care about the perceived fundamental, we are trying to create harmonic spectra, and that is pretty hard if we can't predict what our modulator-carrier ratios will be.  Jan Mattox's fm-drum (ins.lisp) can ignore such niceties.  If you are using low indices and the top pair's mc-ratios are below 1.0 (in vibrato, for example), you have a good chance of getting usable results.  If you want cascade FM to work in other situations, make sure the top oscil has an initial phase of p/2 Ð the middle FM spectrum will then have only sines (not cosines), so the DC amplitude will (sometimes) be 0  Even the tiniest error in the top oscil's phase becomes a beat, inharmonicity, or a timbral change in the output, so cascade FM is almost impossible to keep under control.

The irascible reader may be grumbling about angels and pins, so here's an example of cascade FM to show how strong this effect is:

(definstrument cascade (beg end freq amp modrat modind casrat casind caspha)
  (let ((cr (make-oscil :frequency freq))
        (md (make-oscil :frequency (* freq modrat)))
        (ca (make-oscil :frequency (* freq casrat) :initial-phase caspha))
        (fm-ind0 (* modind modrat freq (/ two-pi sampling-rate)))
        (fm-ind1 (* casind (/ casrat modrat) freq (/ two-pi sampling-rate))))
    (Run
     (loop for i from beg to end do
       (outa i (* amp 
                 (oscil cr 
                        (* fm-ind0 
                           (oscil md (* fm-ind1 
                                        (oscil ca)))))))))))

Compare

(with-sound () (cascade 0 10000 400 .25 1.0 1.0 1.0 1.0 0))

with

(with-sound () (cascade 0 10000 400 .25 1.0 1.0 1.0 1.0 (* .5 pi)))

A similar trick comes up in "feedback FM" used in some synthesizers. Here the output of the carrier is low-pass filtered, then fed back into its fm-input location.  Cellon.ins, written by Stanislaw Krupowicz, is an example Ð the feedback FM occurs in the lines (edited here for clarity):

    (setf fm (one-zero low (* B (oscil fmosc fm))))
    (outa i (* amp (oscil carrier fm)))

or

    sin(y <= wct + B sin y) 

(ignoring the low pass filter for the moment).  This is expanded by Tomisawa as:

     ¥
     S (2/(nB))Jn(nB) sin nwct
     n=1

As Tomisawa points out this is very close to the formulas we've been looking at already, except that the argument to the Bessel function depends on the order, we have only multiples of the frequency wc in the expansion, and the elements of the sequence are multiplied by 2/nB.  The result is a much broader, flatter spectrum than one normally gets from FM.  If you just push the index up in normal FM, the energy is pushed outward in a lumpy sort of fashion, not evenly spread across the spectrum.  In effect we've turned the axis of the Bessel functions so that the higher order functions start at nearly the same time as the lower order functions.  The new function Jn(nB) decreases gradually.  For example if the index (B) is 1:

    J1(1) = .440 
    J2(2) = .353 
    J3(3) = .309
    J200(200) = .076

 and so on. Since the other part of the equation goes down as 1/n, we get essentially a sawtooth wave out of this equation (its harmonics go down as 1/n).  Tomisawa suggests that B should be between 0 and 1.5. Since we are dividing by B in the equation, we might worry that as B heads toward 0, all hell breaks loose, but luckily:

    lim 2/B (J1(B)) = 1
    B=>0 

and for all the other components

    lim 2/nB Jn(nB) = 0
    B=>0 

so, just as in normal FM, if the index is 0, we get a pure sine wave. This also seems natural from 

    sin(y <= wct + B sin y) 

with B=0 (but computer arithmetic doesn't care what "seems natural").

Apparently, if B>1.0, this formula becomes numerically unstable, Ð the result is a burst of noise.  Tomisawa says the problem "is considered to be caused by an error in the digital computation".  We may be rubbing against the behavior of Jn(nB) for  large n and B>1.0 Ð see Watson, "Theory of Bessel Functions", chapter 8.  Jn(nB) might not be a decreasing function of n, so there may be hundreds of noticeable side bands, most of them aliasing.  One of the standard ways to make noise (deliberately) with FM is to increase the index until massive aliasing is taking place.  However, this doesn't explain the odd blippy nature of the noise burst.  Whatever the problem really is, the low pass filter is cutting off the sum before it becomes unmanageable.

Speaking of noise, let's put that into our FM equations and see if anything useful pops out.  The fm-noise instrument in noise.ins performs one form of modulation with noise.  The actual modulation portion is:

    (oscil carrier (* fm-index (randh modulator)))

where randh is producing white noise.  You can produce a sequence of different noises, moving toward gaussian noise, simply by adding randh's together (this is known as the "central-limit theorem")  The probability density of a sum of n independent random variables (randh's in clm) tends to approach a gaussian density as n increases.  We start with "white" noise which has an equal probability of any value. Add another randh, thereby convolving the two flat "white" densities, and you get a /\ shape.  The more you add the closer you get to the bell shaped curve of gaussian noise.  When frequency modulation gets noise, the power spectral density of the output has the same form as the density function of the modulating wave, but now centered around the carrier.  For either kind of noise (white or gaussian), the bandwidth of the result is about 4 times the peak deviation (the random number frequency times its amplitude).  

The real value of noise, however, is not in making hiss.  One of the great improvements to the world of computer music came about when Heinrich Taube fed noise (as a sort of cascade FM) into the parallel modulators of an fm-flute, but not into the carrier.  The modulating signal becomes a sum of two or three narrow band noises (narrow because normally the amplitude of the noise is low), and these modulate the carrier.  It is my belief that you get the normal spectrum with each component smeared out by a copy of the noise band. 

In clm, the basic idea is:

    (oscil carrier 
           (* fm-index (oscil fm 
                              (* noise-index (randh noise)))))

In the realm of "anything" as the modulating signal, consider 

    (sin (+ sound-file (* index (sin (* 2 pi sound-file))))

where "sound file" is any recorded sound.  I call this "contrast-enhancement" in the clm package.  It makes a sound "crisper" Ð "Wait for Me!" uses it whenever a sound needs to cut through a huge mix of sounds.

We can use more than one sinusoidal component in our carrier, or multiple banks of carriers and modulators, and depend upon vibrato and "spectral fusion" to make the result sound like one voice.  In this cross between additive synthesis (the multiple carriers) and FM (the "formant" centered on each carrier), we get around many of the limitations of the Bessel functions.  There are numerous examples in fmviolin.clm.  Chowning's version is reson.ins and jcvoi.ins. One of the raspier versions of the fm-violin used a sawtooth wave as the carrier.

One of the more elaborate multi-carrier FM instruments is the voice instrument written by Le Brun and used in "Colony" and other pieces. In a simplified form it is:

(definstrument vox (beg end freq amp 
                    &optional (indexes '(.005 .01 .02)) 
                              (formant-amps '(.86 .13 .01)))
  (let* ((car-os (make-oscil :frequency 0))
         (evens (make-array 3 :element-type 'osc))
         (odds (make-array 3 :element-type 'osc))
         (amps (make-array 3 :element-type 'short-float 
                             :initial-contents formant-amps))
         (ampf (make-env :envelope '(0 0 25 1 75 1 100 0) 
                         :scaler amp :start beg :end end))
         (frmfs (make-array 3 :element-type 'seg-func))
         (indices (make-array 3 :element-type 'short-float 
                                :initial-contents indexes))
         (frq 0.0) (car 0.0) (frm 0.0) (frm-int 0) (frm0 0.0) 
         (even-amp 0.0) (odd-amp 0.0) (even-freq 0.0) 
         (odd-freq 0.0) (sum 0.0)
         (per-vib (make-triangle-wave :frequency 6 
                                      :amplitude (* freq .03)))
         (ran-vib (make-randi :frequency 20 
                              :amplitude (* freq .5 .02))))
    (dotimes (i 3)
      (setf (aref evens i) (make-oscil :frequency 0))
      (setf (aref odds i) (make-oscil :frequency 0)))
    (setf (aref frmfs 0) (make-env :envelope '(0 520 100 490) 
                                   :scaler 1.0 :start beg :end end))
    (setf (aref frmfs 1) (make-env :envelope '(0 1190 100 1350) 
                                   :scaler 1.0 :start beg :end end))
    (setf (aref frmfs 2) (make-env :envelope '(0 2390 100 1690) 
                                   :scaler 1.0 :start beg :end end))
    (Run
     (loop for i from beg to end do
       (setf frq (+ freq (triangle-wave per-vib) (randi ran-vib)))
       (setf car (oscil car-os (in-Hz frq)))
       (setf sum 0.0)
       (dotimes (k 3)
         (setf frm (env (aref frmfs k)))
         (setf frm0 (/ frm frq))
         (setf frm-int (floor frm0))
         (if (evenp frm-int)
             (progn
               (setf even-freq (in-Hz (* frm-int frq)))
               (setf odd-freq (in-Hz (* (+ frm-int 1) frq)))
               (setf odd-amp (- frm0 frm-int))
               (setf even-amp (- 1.0 odd-amp)))
           (progn
             (setf odd-freq (in-Hz (* frm-int frq)))
             (setf even-freq (in-Hz (* (+ frm-int 1) frq)))
             (setf even-amp (- frm0 frm-int))
             (setf odd-amp (- 1.0 even-amp))))
         (incf sum (+ (* (aref amps k) 
                      (+ (* even-amp 
                            (oscil (aref evens k) 
                              (+ even-freq (* (aref indices k) car))))
                         (* odd-amp 
                             (oscil (aref odds k) 
                               (+ odd-freq (* (aref indices k) car))))
                       )))))
      (outa i (* (env ampf) sum))))))

(with-sound () (vox 0 20000 110 .5 
                    '(0.02 0.01 0.02) '(.9 .09 .01)))

References

Abramowitz and Stegun, Handbook of Mathematical Functions, 1965.

Bate, "The Effect of Modulator Phase on Timbres in FM Synthesis", CMJ vol 14, no3 1990, 
	p38-45.

Chowning,  "The Synthesis of Complex Audio Spectra by Means of Frequency Modulation", 
	JAES 21:526-534, 1973

De Poli, "A Tutorial on Digital Sound Synthesis Techniques", CMJ vol7 no 4, 1983.

Gagliardi, "Introduction to Communications Engineering", Wiley Interscience, 1978.

Holm, "Understanding FM Implementations:  A Call for Common Standards", CMJ vol 16 no. 1.

Le Brun, "A Derivation of the Spectrum of FM with a Complex Modulating Wave", CMJ vol1 
	no 4 1977 p51-52

Moorer, "Signal Processing Aspects of Computer Music: A Survey" Proc IEEE  vol 65 1977.

Palamin, Palamin, Ronveaux "A Method of Generating and Controlling Asymmetrical Spectra"  	JAES vol 36, no 9, Sept 88, p671-685.

Schottstaedt,  "The Simulation of Natural Instrument Tones Using Frequency Modulation with 	a Complex Modulating Wave", CMJ vol 1 no 4 1977 p46-50

Taub and Schilling, "Principles of Communications Systems", McGraw-Hill, 1986.

Tomisawa, "Tone Production Method for an Electronic Musical Instrument" US Patent 		4,249,447, 1981.

Watson, "A Treatise on the Theory of Bessel Functions", Cambridge, 1922.

1991 (revised 1992 to include Holm's paper)