FM synthesis using Numpy

[Bas]

You got your math wrong. What you are calculating is:
sin(2*pi*(1000+15*sin(2*pi*6*t))*t) = sin(2*pi*1000*t +
2*pi*15*sin(2*pi*6*t)*t)
The 'instantaneous frequency' can be calculated by differentiating the
argument of the sine and dividing by 2pi:
x = sin(phi(t)) -> f_inst = d (phi(t)) / dt / (2*pi)
So in your case:
f_inst = 1000 + 15*sin(2*pi*6*t) + 2*pi*t*6*15*cos(2*pi*6*t)
the last term explains the effect you hear.

What you want is:
f_inst = f0 + df*cos(2*pi*fm*t)
Integrating this and multiplying by 2pi gives the phase:
phi(t) = 2*pi*f0*t + sin(2*pi*fm*t)*df/fm

So you can achieve the frequency modulation by using phase modulation
(these two are related). You can do this with your own code by

phi = oscillator(t, freq=6, amp=15/6)
tone = oscillator(t, freq=1000, amp=0.1, phase=phi)

cheers,
Bas