When concatenating sine waves, how do I phase shift in order to prevent "pops" caused by sudden jumps in amplitude?

Question

I'm writing a synthesizer using numpy and python. I've run into a problem where concatenating sine waves causes 'pops' or 'clicks' when they are outputted to audio, say, via scipy.io.wavfile.write.

Here is some example code that illustrates the problem.

import matplotlib.pyplot as plt
import numpy as np

sr = 22050
ts = 1.0/sr
duration = 1
t = np.arange(0,duration,ts)

freq = 4.20
sine1 = 1*np.sin(2*np.pi*freq*t+0)

freq = 6.66
sine2 = 1*np.sin(2*np.pi*freq*t+0)

plot_tone = np.concatenate([sine1, sine2])
plot_t = np.arange(0,duration*2,ts)

plt.figure(figsize = (20, 14))
plt.plot(plot_t, plot_tone, 'green')
plt.ylabel('Amplitude')
plt.xlabel('Time (s)')
plt.annotate('tone1 end', (1, sine1[-1]), fontsize=14)
plt.annotate('tone2 start', (1, sine2[0]), fontsize=14)
plt.savefig('concat_sine_without_offset.png')
plt.show()

In the above output, you can see the jump in amplitude. If these frequencies were in the range of human hearing and I outputted as a wav file, the pop would be apparent, especially when concatenating repeatedly.

After some additional searching, I came across the Stackoverflow answer linked below. Based on my understanding of the answer, I attempted to calculate the phase offset, which is illustrated with the code below. https://stackoverflow.com/a/36466582/13959910

import matplotlib.pyplot as plt
import numpy as np

sr = 22050
ts = 1.0/sr
duration = 1
t = np.arange(0,duration,ts)

freq = 4.20
sine1 = 1*np.sin(2*np.pi*freq*t+0)

phase_increment = 2*np.pi*freq/sr
phase_offset = phase_increment * len(sine1)
print(f'{phase_offset=}')

freq = 6.66
sine2 = 1*np.sin(2*np.pi*freq*t+phase_offset)

plot_tone = np.concatenate([sine1, sine2])
plot_t = np.arange(0,duration*2,ts)

plt.figure(figsize = (20, 14))
plt.plot(plot_t, plot_tone, 'green')
plt.ylabel('Amplitude')
plt.xlabel('Time (s)')
plt.annotate('tone1 end', (1, sine1[-1]), fontsize=14)
plt.annotate('tone2 start', (1, sine2[0]), fontsize=14)
plt.savefig('concat_sine_with_offset.png')
plt.show()

I believe this gets me what I need, which is great. The 'end' and 'start' annotations line up, the concatenation of the two sine waves appears to be smooth.

However, when I try this with sines that have differing amplitudes, the problem returns.

e.g. in the second set of code, if you replace sine1 like so - sine1 = 0.75*np.sin(2*np.pi*freq*t+0) - there will be a sudden jump in amplitude.

To summarize: how could I phase shift a sine wave such that it begins near to the amplitude of the preceding sine wave? Or perhaps more fundamentally, are there recommended methods for preventing large jumps in amplitude when concatenating sine waves? Any help is greatly appreciated!

Answer 1

Phase continuous transitions can be accomplished by implementing the waveform in phase versus time instead of frequency versus time. Frequency is the time derivative of phase, so instead of stepping between two frequency values with a frequency vs time waveform, the phase vs time waveform would be ramps with the slope of the ramp setting the desired frequency. Note that a numerically controlled oscillator (NCO) does this directly and is simply a phase accumulator followed by a look-up table as one approach to compute the trigonometric $\sin(\theta)$ or $\cos(\theta)$ computation (or both for a complex $e^{j\theta}$ output!), converting instantaneous phase to a sinusoid. The rate of accumulation is set by a Frequency Control Word which sets the ramp rate of the phase versus time. In typical NCO implementations, the phase rolls over as a cyclical counter at $2\pi$ which is of no consequence to phase continuity in a sine wave and prevents any overflow issues. The complete details for a numerically controlled oscillator implementation are given here.

The NCO is ideal for this application as it is relatively simple in implementation and will naturally maintain phase continuity at any given change in frequency. Additionally we can implement arbitrary frequency modulation, and it can be easily adapted for a direct phase and amplitude modulation as well as explained in the referenced link. An example NCO implementation in Python is as follows, fcw is the frequency as an vector that update on every sample, thus we can change the frequency to any arbitrary value at any time; sr is the sample rate:

import numpy as np
def nco(fcw, sr):
    phase =  0
    phase_result = []
    for fcw_samp in fcw:
        ph_step = 2*np.pi* fcw_samp * 1/sr
        phase += ph_step
        phase_result.append(phase)
    return np.cos(phase_result)
    

And demonstrating the OP's example:

sr = 22050
ts = 1.0/sr
duration = 2
t = np.arange(0,duration,ts)
freq1 =4.20
freq2 = 6.66
t_change = 1
fcw = np.zeros_like(t)
fcw[t<=t_change] = freq1
fcw[t>t_change] = freq2
result = nco(fcw, sr)

Results in the following plot:

The amplitude itself can be changed completely independently with a similar “Amplitude Control Word” which scales its output. This can be tapered as desired to avoid abrupt transitions. The taper can be in the vector containing the amplitudes for each sample (similar to the vector of frequency control words above), or the resulting signal can be passed through an FIR low pass filter with the cutoff set above the highest frequency desired.

Below shows an example of this with the python NCO function modified to support amplitude control:

def nco(fcw, acw, sr):
    # fcw: array-like container of frequency control values, one for each sample
    # acw: array-like container of amplitude control values, one for each sample
    # sr: sample rate
    phase =  0
    phase_result = []
    for fcw_samp in fcw:
        ph_step = 2*np.pi* fcw_samp * 1/sr
        phase += ph_step
        phase_result.append(phase)
    return acw * np.cos(phase_result)

And demonstrating with the OP's example where a Gaussian filtered transition was use to change the amplitude:

import scipy.signal as sig
sr = 22050
ts = 1.0/sr
duration = 2
t = np.arange(0,duration,ts)
freq1 =4.20
freq2 = 6.66
amp1 = 0.75
amp2 = 1
t_change = 1
taper_dur = 0.2
b=10
fcw = np.zeros_like(t)
acw = np.zeros_like(t)
fcw[t<=t_change] = freq1
fcw[t>t_change] = freq2
acw[t<=t_change] = amp1
acw[t>t_change] = amp2
nsamps = int(taper_dur * sr) * 2
coeff = sig.gaussian(nsamps, 500)
acw2 = sig.filtfilt(coeff, np.sum(coeff), acw) 
result = nco(fcw, acw2, sr)

Resulting in the following plot:

---

Mathematical Description of Changing Frequency with Continuous Phase

Below provides a complete mathematical relationship for precomputing the phase waveforms needed for frequency change with continuous phase.

Specific to the OP's example, the first phase ramp is given by $2\pi 4.20$ radians/sec starting at $t=0$ and transitions to $2\pi 6.66$ radians per second at $t=t_1$.

This makes it then very easy to keep the phase continuous at the transitions by simply accumulating the phase as given by the desired slope (which is the target frequency), similar to the operation of the NCO. Note in this case, the offset is given by:

$$\theta_{offset} = \theta(t_1) - 2\pi f_2 t_1 = 2\pi f_1 t_1 - 2\pi f_2 t_1$$

From the resulting phase ramp we can plot the sine or cosine of that phase versus time as desired and it will be phase continuous.

This is demonstrated in Python below with the following results:

Python:

import matplotlib.pyplot as plt
import numpy as np
sr = 22050 
ts = 1.0/sr 
duration = 2 
t = np.arange(0,duration,ts)
t_change = 1
freq1 =4.20 
phase1 = 2*np.pi*freq1* t[t<=t_change]
freq2 = 6.66 
phase2 =  2*np.pi*freq2*t[t>t_change] + 2*np.pi*(freq1-freq2)*t_change
total_phase = np.concatenate((phase1, phase2))
out = np.cos(total_phase)

Considering more generalized cases with multiple steps, the phase offset to add for any given segment is:

$$\theta_n = \sum_{m=1}^{n-1}2\pi f_m t_m - \sum_{m=1}^{n-1}2\pi f_{m+1} t_m \tag{1} \label{1}$$

With the resulting phase ramp proceeding starting at $t_n$ as:

$$\theta(t) = 2\pi f_n t + \theta_n$$

This is demonstrated with a 3 segment ramp with reference to the figure below, where we note that each next phase ramp will start at the completed phase from the previous frequency. For example, the second phase offset given as $\theta_2$ below is similarly determined as above by solving for the intersection of linear line equations to be:

$$\theta_2 = \theta(t_2)-2\pi f_3 t_2$$ $$= 2\pi f_2 t_2 + \theta_1 - 2\pi f_3 t_2$$ $$ = 2\pi f_2 t_2 + 2\pi (f_1-f_2)t_1 - 2\pi f_3 t_2$$

Which is consistent with the generalized form given in $\ref{1}$ as:

$$ \theta_2 = 2\pi f_1 t_1 + 2\pi f_2 t_2 -2\pi f_2 t_1 - 2\pi f_3 t_2$$

With this, the third segment plotted below is given as:

$$\theta(t) = 2\pi f_3 t + \theta_2$$

Summary:

In summary I have given two approaches to creating continuous phase frequency transitions: as an NCO which automatically maintains phase continuity through frequency transitions, and a computationally based approach for use when the formula for the phase waveform is desired.

Answer 2

IMO the best way to implement this is a rotating phasor. Recall that

$$ e^{jx} = \cos(x) + j \cdot \sin(x) $$

and

$$ e^{j\omega(n+1) } = e^{j\omega n } e^{j\omega}$$

That means we can calculate the a sine wave as the imaginary part of a complex phasor and that the next value of the phasor is the current value multiplied with $e^{j\omega}$.

Below is the Matlab code that implements your example. This approach has a lot of advantages.

1. It doesn't require any transcendent functions like sine or cosine.
2. It's very fast. You only need one complex multiply per sample. No branching or function calls required
3. You can change the current frequency at any sample simply by changing the multiplier from $e^{j\omega_1}$ to $e^{j\omega_2}$
4. It has extremely low memory footprint. No need for look up tables, coefficient tables etc.
5. It doesn't need to be constraint to $[0,2\pi]$ like a phase accumulator (which is otherwise also a good option)
6. It's a "quadrature" oscillator: i.e. you get both sine and cosine for the price of one (if you need that).
7. It's great for real time and frame based implementation. It's very fast and state keeping between frames is trivial.

MATLAB code

%% variable frequency oscillator using phasor rotation

% parameters
fs = 22050; % sample rate in Hz
freq1  = 4.2; % first frequency in Hz
n1 = 22050; % length of first frequency in samples
freq2 = 6.66;
n2 = 22050;

% state variables
y = zeros(n1+n2,1); % output vector
W1 = exp(1i*2*pi*freq1/fs); % first frequency rotation
W2 = exp(1i*2*pi*freq2/fs); % 2nd  frequency rotation
phasor = 1; % initial state of the phasor

%% create the sine waves
for i = 1:n1
  y(i) = imag(phasor);
  phasor = phasor*W1;
end
for i = n1+(1:n2)
  y(i) = imag(phasor);
  phasor = phasor*W2;

end

%% plot the result

plot((0:(n1+n2-1))/fs,y);
grid on
xlabel('Time in seconds');

Answer 3

Use phase from the start, and do not just multiply by time t - this is vital

For each sine wave, use a variable to hold the current phase. Then for each time step, multiply 2*np.pi*freq by the time step interval and add that incremental change to the previous phase value. You'll also need to check for the sum exceeding 2*np.pi and subtract 2*np.pi when it does. With this solution, you're guaranteed a stepless frequency change.

There's an even more important benefit from this too, which is essential to know for real-world DSP. As time t becomes larger, the value will start to run out of floating-point precision for the time step interval. Your sine wave will become progressively less accurate as the program runs, until at some point the value loses precision entirely and the sine wave will simply stop moving. The lower the frequency and the shorter the sample rate, the sooner this'll happen; but maths says it will always happen eventually. This is a very common bug for novice DSP implementations, and it's often not obvious because everything may seem to work fine until you've had your code running for a few hours and then it mysteriously just freezes.

Incremental phase changes do have a similar issue of course, where the increment may be too small to represent with your floating-point precision. Using double-precision floating-point is usually the way to go here. The advantage you have with this though is that if there is a precision problem, you're guaranteed to see it within one sine-wave period, so it won't come as some kind of surprise.

Ramp amplitudes if they can change

Your code has a fixed amplitude of 1 for both sine waves. (By the way, if the amplitude is 1 then you don't need to multiply by it.) If you want changes to amplitude to be stepless, one solution is to ramp the amplitude from the previous value to the new value. You can choose a suitable ramp period for your application.

And/or change settings at zero crossing

The other simple strategy for amplitude changes is to hold off changing until the next zero-crossing of the sine wave. Multiplying anything by zero results in zero, so you're guaranteed that this will be stepless. This works well for fast sine waves. For slow sine waves where the sine wave period may be longer than a reasonable ramp transition, it might not be the best strategy though.

Changing settings at the zero crossing can also work for frequencies too, but there's generally less reason to need that. Still though, if you want the frequency and magnitude to both change in step with each other, you can change them both when the phase crosses 2*np.pi, and you're good. By the way, don't expect the phase to ever exactly hit 2*np.pi, so a zero-crossing transition is never going to be perfectly stepless, but it's still likely to be better than most other strategies.

Answer 4

One solution I see is to concatenate the two, different amplitude sines at their zero crossings, making sure they have the same phase (e.g. if one goes up, the other has to go up, too). This way, the discontinuity should be minimal, but there is no other way to do it without discontinuities, since the amplitudes differ. The other two answers already explain at large how to match the phase, I'm just saying where. It might pay off to calculate the derivatives and see the point were the slopes have the minimal difference; or not.