Interpolating frequencies

This is more of a math problem than a programming problem, but it has now come up for me twice in two different programming projects.

I would like to create a function, w_sweep(t, f0, t0, f1, t1), that samples at time t a sine wave that sweeps through frequencies from f0 when t == t0 to f1 when t == t1. (You can assume frequencies are in Hz and times are in seconds. I don't care what happens when t is outside the range t0..t1.)

I started by writing a simpler function that uses a constant frequency:

float w_steady(float t, float freq) {
  return std::sin(2.0f * PI * freq * t);
}

Intuitively, it seems like w_sweep would just involve linear interpolation of the frequency.

// INCORRECT SOLUTION
float w_sweep(float t, float f0, float t0, float f1, float t1) {
  const float freq = std::lerp(f0, f1, (t - t0)/(t1 - t0));
  return std::sin(2.0f * PI * freq * t);
}

At first look, the result seems correct: The frequency is indeed f0 in the vicinity of t0 and then it starts to sweep toward f1. But by the time you reach t1, the actual waveform will have overshot f1 by quite a bit.

I think what's happening is that, instead of sweeping through the frequencies linearly with t, we're actually going at t*t, since t is a factor in the argument and a factor in freq. Another way to think of it is that each sample is drawn from a different wave function. They're all sine waves at different constant frequencies. But they aren't really correlated (except at time 0, when all of them will give the same value).

If I wanted to create sequential samples with a loop, I could solve the problem by accumulating the position along the sign wave.

float omega = 0;
for (float t = t0; t <= t1; t += dt) {
  samples.push_back(std::sin(omega));
  const float freq = std::lerp(f0, f1, (t - t0)/(t1 - t0));
  omega += 2.0f * PI * freq * dt;
}

That works, but I'd prefer a w_sweep that it can be sampled at arbitrary times.

Solution

You have to think about it in terms of phase. For a constant frequency the phase ist just the product of time and frequency. But in general it's the integral over the frequency function across time.

The easiest way to show how this works is with a diagram, but its also possible with formal math. Unfortunately SO markdown neither supports graphs nor latex markdown. But if work it through, the result is:

float w_sweep(float t, float f0, float t0, float f1, float t1) {
  float const freq = std::lerp(f0, f1, (t - t0)/(t1 - t0));
  return std::sin(PI * (f0 + freq) * (t-t0));
}

The graphical approch can be described as follows: The phase is area under the frequency function in time, covered since t0. For our example it is best to split it into a rectangular area at the base f0*(t-t0) and a triangular area on top 0.5*(freq-f0)*(t-t0). Multiply the sum by 2π and you get PI * (f0 + freq) * (t-t0).