sine wave glissando from one pitch to another in Numpy

Question

I have been working on a program where I need to slowly and smoothly change the pitch of a sine wave from one pitch to another. I am able to get an array of the frequency the pitch should be at any given moment (for instance, [440, 526.5, 634.2 794.8, 880], though much, much longer) but it seems I am unable to actually apply that frequency to a wave. My best attempt is:

numpy.sin(2*math.pi*x*freq/self.sample_rate)

where "freq" is the array of frequencies and x is an enumeration array ([0,1, 2, 3, 4...]). This method sort of works, however it makes the frequency go above the expected frequency, and then back down. I have been working on this problem for a very long time and have been unable to make any progress on finding a more appropriate method. Any advice? Was I clear enough in expressing my dilemma?

Thank you.

Answer 1

The issue is that as you ramp through the frequencies, each frequency effectively has a different phase for the given time. When you scroll through these phases quickly and continuously, they drive the sine wave at higher frequency (or lower is also possible).

Imagine, for example, that you changed the frequency instantaneously -- to do this you'd have to supply the phase correction p_1 = p_0 + 2*pi*t*(f_0-f_1) to make the phases match up at time t. As you do this is little steps, you also have to make a similar phase correction, with each phase correction adding to the previous.

Here's the resulting figure, with the code below. The top figure is the frequency the middle is without the phase correction, and the bottom has the continuously corrected phase.

from pylab import *

sample_rate = .001
f0, f1 = 10, 20
t_change = 2

times = arange(0, 4, sample_rate)

ramp = 1./(1+exp(-6.*(times-t_change)))
freq = f0*(1-ramp)+f1*ramp
phase_correction = add.accumulate(times*concatenate((zeros(1), 2*pi*(freq[:-1]-freq[1:]))))

figure()
subplot(311)
plot(times, freq)
subplot(312)
plot(times, sin(2*pi*freq*times))
subplot(313)
plot(times, sin(2*pi*freq*times+phase_correction))

show()