How to make sure that a sinusoid wave with different subsequent frequencies matches up?

Question

A very basic formula to transform a sequence of notes into a signal would be:

$$ x(t) = A\cos( 2 \pi (2^{current\_note}/12) t) $$ where $current\_note$ is $notes[\lfloor t / period\_per\_note\rfloor]$ and $a[k]$ is the indexing operation into a string; in this case into a sequence of numerical values representing notes on a chromatic scale.

Listen to an example

However, when making a signal like this, whenever the frequency changes (when moving from one note to the next), the function is discontinuous, often starting at a completely different amplitude because in general, the periods of the two sinusoids with different frequencies end up being very different (not being a multiple of the same 'base frequency')-

This can also clearly be heard in the example, as 'clicks'.

Now the question: How can we make sure that the sinusoid with the 'new' frequency starts at the same value of x as the one with the old frequency?

I am guessing that this somehow involves keeping track of the current 'phase', so that we can time-shift the new frequency to make sure their amplitudes match up. But I have no idea how to approach this, neither from a mathematical nor from a programmatic perspective.

Answer 1

Since the OP asked for a real code example, there is one below. It's in Matlab, but optimized for implementation on a real time processor.

The idea is to implement the oscillator through a complex phasor multiplication, which only requires one complex multiply per sample and no transcendental functions (which tend to be very expensive).

The only potential risk of drift or numerical issues is the magnitude of the state phasor, which needs to stay unity. This can be compensated by occasionally dividing the state phasor by it's magnitude. Once every few minutes is plenty for 32-bit floating point. This requires two "expensive" functions, sqrt() and inverse. If that's a problem it can be made much quicker by using the Taylor expansion of $\frac{1}{\sqrt{x}}$ around $x = 1$.

%% Parameters
frequencies_in_Hz = [200 400 800 100];
durations_in_seconds = [0.333 0.41675 0.2 0.5011];
sample_rate_in_Hz = 44100;

%% setup the state variables. We do this all in single precision
duration_in_samples = round(durations_in_seconds*sample_rate_in_Hz);
output = zeros(sum(duration_in_samples),1,'single');
% initialize phasor and state
phasors =  single(exp(-1i*2*pi*frequencies_in_Hz/sample_rate_in_Hz));
state = single(1 + 0i); % real part 1, imaginary part 0

%% Now comes the real time part:
cnt = 1;
% loop over all frequencies
for freq = 1:length(frequencies_in_Hz)
  ph = phasors(freq);
  % innner loop over all samples per freqeuncies
  for i = 1:duration_in_samples(freq)
    % phasor multiply
    state = state *ph;
    % take real part for cosine, or imaginary for sine
    output(cnt) = real(state);
    cnt = cnt + 1;
  end
  % amplitude drift correction through a taylor expansion of 1/sqrt(x) aroud
  % x = 1 (the magintude of the phasor will be very close to one)
  g = single(.5)*(single(3)-real(state)*real(state)-imag(state)*imag(state) );
  state = state*g;

end
%% show the result
plot(output);

Answer 2

You need to make sure that the phase of continuous.

The easiest way would be something like this.

x(n) = cos(phase);
phase = phase + 2*pi*current_frequency*sample_time;
current_frequency= update_frequency(n);
n = n + 1;

Note, that this is very inefficient code and only for illustrative purposes.