How to concatenate sine waves without phase jumps

Question

I need to make a python script which generates sine waves of a given frequency and plays them using pyaudio (blocking mode), I also need to be able to change this frequency while it runs, modulate it and plot it using pyqtgraph. For now I have a thread generating chunks of data and my approach to 'connect' those sines was to get the fft and then calculate the angle (numpy.angle), store it in a variable and use it as the phase offset to the next chunk, but I'm not getting the results I expected, maybe I'm missing something or mixing them up.

import matplotlib.pyplot as plt
import numpy as np
import pyaudio

#-----------------------
CHUNK = 1024
RATE =  44100
CHANNELS = 2
FORMAT = pyaudio.paFloat32
#-----------------------

samples = int(CHUNK)
t = np.arange(samples) / RATE
con = 0


def generate_sine(a: float = 0.5, freq: float = 440.0):

    global con

    sine = a * np.sin(2.0 * np.pi * freq * t + con)

    # get the angle of the wave

    phase = np.angle(np.fft.fft(sine))

    # update ref var to generate subsequent sines
    # begining where the last ended

    con = phase[-1]

    return sine


def play_sine(data):

    pa = pyaudio.PyAudio()

    stream = pa.open(format=FORMAT,
                         channels=CHANNELS,
                         rate=RATE,
                         input=False,
                         output=True,
                         frames_per_buffer=CHUNK)

    stream.write(np.array(data).astype(np.float32).tostring())

    stream.close()

if __name__ == '__main__':

    f = 80

    chunks = generate_sine(freq=f)

    for i in range(0,4):

        chunks = np.concatenate((chunks, generate_sine(freq=f)))

    #for i in range(0,10):

    #play_sine(chunks)

    plt.plot(chunks)

    plt.show()

Demo plot

You can see in the imagem linked that there are discontinuities around x=1024, x=2048 and so on.

Answer 1

You are generating your signal with

    a * sin(2πf * t + con)

where t ranges over [0 .. CHUNK/RATE).

When you start the next chunk, t is reset to zero. To generate a continuous waveform, you need to modify con to generate the same resulting phase value as the previous sample.

Using an FFT isn't going to work because the signal you are generating isn't an exact multiple of the sample window, plus you are actually interested in the phase at the end of the sample window, not the phase at the beginning.

Instead, you simply need the generating function's phase value at t=t_end, modulo 2π.

Ie, you could simply use:

con = 2.0 * np.pi * f * CHUNK/RATE + con

but that value will grow, possibly causing eventual numerical issues if you concatenate a lot of chunks together, where the frequencies are high. Since the sine function is periodic, you just need to normalize the end phase into the 0 to 2π range:

con = math.fmod(2.0 * np.pi * f * CHUNK/RATE + con, 2.0 * np.pi)

---

If you modify your generation function to be:

    a * sin(2π * (f * t + con))

then con represents the fraction of a full cycle carried forward from the previous chuck, and you can avoid the modulo division by 2π, which might improve accuracy slightly.

con = math.modf(f * CHUNK/RATE + con)[0]

---

Attempting a clearer explanation:

Note: This technique only works because you exactly know the generation equation for the previous chunk, and are generating the following chunk. If either of those conditions change, you'll need a different technique to match the chunks.

The previous chunk is generated using the sin() of the following sequence:

2πf₁*0/RATE+con₁, 2πf₁*1/RATE+con₁, ..., 2πf₁*1022/RATE+con₁, 2πf₁*1023/RATE+con₁

It should be clear that, for a smooth transition to the next chunk, that chunk should start with the sin() of 2πf₁*1024/RATE+con₁

The next chunk starts with the sin() of 2πf₂*0/RATE+con₂.

Therefore, we will have a smooth transition if:

2πf₂*0/RATE + con₂ = 2πf₁*1024/RATE + con₁

or

     0      + con₂ = 2πf₁*1024/RATE + con₁

or

              con₂ = 2πf₁*1024/RATE + con₁

which can be written in your generate_sine function as:

con = 2.0 * np.pi * f * CHUNK/RATE + con

which is the equation that came "out of nowhere" in my above answer. From that point, since the sin() function is 2π periodic, I'm just performing modulo 2π reductions, to keep the argument to sin() from growing without bound, causing numerical inaccuracies.

Hope that makes things clearer.

Answer 2

The following solution is not 100% satisfactory but seems to work well if you are able to discard a bit of each snippet's signal at the joints.

It uses Hilbert transform to reconstruct the phase. Unfortunately near the edges of the signal there tend to be artifacts (see for example very end of the red signal in the plot), which I could only get rid of by cutting them off.

>>> import numpy as np
>>> from scipy.signal import hilbert
>>> 
# create two sinewaves and join them discarding some signal near the joint
# this is the black curve in the plot
>>> A = np.sin(np.linspace(0, 30, 10000))
>>> B = np.cos(np.linspace(0, 15, 10000))
>>> Cjump = np.concatenate([A[:-1000], B[1000:]])
>>> 
# join the same sinewaves using Hilbert
>>> AP = np.unwrap(np.angle(hilbert(A)))
>>> BP = np.unwrap(np.angle(hilbert(B)))
>>> CP = np.concatenate([AP[:-1000], BP[1000:] - BP[1000] + AP[-1000]])
>>> C = np.cos(CP)
# this is the red curve in the plot

Why Hilbert transform? Because it can to some extent handle modulated signals: example (blue orig, green from reconstructed phase):

Code to generate signals:

# original (blue):
>>> A = np.sin(np.convolve(np.random.uniform(0.2, 1, 10).repeat(1000), np.ones(10)/1000, 'same').cumsum())
>>> AP = np.unwrap(np.angle(hilbert(A)))
# recon (green):
>>> AR = np.cos(AP)