Search code examples
pythonscipysignalssignal-processing

Synchronizing two noisy sinusoids out of phase

I currently have two noisy sinusoids which are out of phase and I am looking for a solution to set them in phase. Firstly, I was trying to just adjust the sinusoids via finding the zero-crossing values with an upward trend, however, as there is some noise in the signal, this solution is not suitable. I made a small mock-up example:

time = np.arange(0, 60, 0.005)
amp = 2
F = 0.25
signal1 = amp*np.sin(2*np.pi*F*time + np.pi)
signal2 = amp*np.sin(2*np.pi*F*time + np.pi/3)

mu, sigma = 0.2, 0.1 # mean and standard deviation
noise1 = np.random.normal(mu, sigma, len(signal1))
noise2 = np.random.normal(mu, sigma, len(signal2))


noisy_signal1 = signal1 + noise1
noisy_signal2 = signal2 + noise2

plt.figure()
plt.plot(time, noisy_signal1)
plt.plot(time, noisy_signal2)
plt.show()

Of course that in reality, I don't know which is the phase shift, but I do know the amplitude and the frequency of each sinusoidal beforehand. What I want to achieve, is to go from the plots you can see when running the code, to this: enter image description here

EDIT1: Using the correlation function might be a suitable solution. However, you will have the problem that since the signal is noisy, the maximum correlation could be found in several periods ahead of one of the signals, and therefore, losing a lot of data points when doing the shift. I am looking specially for a solution in which I am identifying the phase and just shifting one of the sinuoids.

Solution

  • Thanks for asking this question. I'm looking into applications of cross-correlation and this was my first try at writing some code.

    Cross-correlation repeatedly convolves two signals and provides a (long) array as output, each expressing how well correlated the two signals are at any given time difference. For sinusoidal signals like yours a few cycles would probably be sufficient rather than correlating the entire array, and decimation might also be reasonable to improve execution time. Also I'm using the numpy correlate function rather than the scipy one, the latter using FFTs to work faster on large arrays. This should be a fascinating area of exploration.

    To shift the phase of the second signal, first the signals are cross-correlated, then I find the point with the maximum value. The second array is rotated by the index of that maximum value to synchronize the the two arrays so that they have maximum correlation. Rotation assumes full cycles so shifting and filling with zeroes would be more universal.

    It's just 3 lines of code (numpy is wonderful), inserted right before the plot. There may be edge cases I haven't thought about, but this works exactly the way you wanted it.

    corr = np.correlate(noisy_signal1, noisy_signal2, 'full')
corr_max_idx = np.argmax(corr)
noisy_signal2 = np.roll(noisy_signal2, corr_max_idx)

    To get a very similar result losing less than a half cycle, you just need to find the remainder of the argmax calculation when dividing by the frequency in samples. We know what the frequency is from your example program, but I calculate it anyway using FFT. While FFT isn't the recommended way of calculating a single frequency, I found the signal too noisy to reliably calculate from zero crossings, as there are several per cycle.

    Here's what I came up with. I changed a bit of the example code so the sample rate is available to me for later, although I ended up not needing it.

    import numpy as np
from scipy import fft
from matplotlib import pyplot as plt

sr = 200
time = np.arange(0, 60, 1/sr)

amp = 2
F = 0.25
omega = 2 * np.pi * F

signal1 = amp*np.sin(omega * time + np.pi)
signal2 = amp*np.sin(omega * time + np.pi/3)

mu, sigma = 0.2, 0.1 # mean and standard deviation
noise1 = np.random.normal(mu, sigma, len(signal1))
noise2 = np.random.normal(mu, sigma, len(signal2))


noisy_signal1 = signal1 + noise1
noisy_signal2 = signal2 + noise2

# Get frequency in samples
fft1 = fft.fft(noisy_signal1)
fft2 = fft.fft(noisy_signal2)
max_f = np.argmax(np.abs(fft1))
f = len(noisy_signal1) / max_f

# Perform correlation to determine best fit for phase difference (in samples)
corr = np.correlate(noisy_signal1, noisy_signal2, 'full')
delta_p = int(round(np.argmax(corr) % f))

# Chop phase difference (maximum 1/2 cycle)
fi = int(round(f,0))
if delta_p < fi/ 2:
    noisy_signal1 = noisy_signal1[delta_p:]
    noisy_signal2 = noisy_signal2[:-delta_p]
    time = time[delta_p:]
else:
    noisy_signal1 = noisy_signal1[:delta_p - fi]
    noisy_signal2 = noisy_signal2[fi - delta_p:]
    time = time[:delta_p - fi]

plt.figure()
plt.plot(time, noisy_signal1)
plt.plot(time, noisy_signal2)
plt.show()

    To get rid of correlation, which is computationally expensive, we can also use the FFT that was used to calculate the frequency, changing the two lines for the correlation to:

    # Do the same calculation with the FFT
p1 = np.angle(fft1[max_f])
p2 = np.angle(fft2[max_f])
delta_p = int((f * ((p2 - p1) / (2 * np.pi))) % f)

    Another possibility that requires neither correlation or FFT, if the frequency is known, is to use a single bin FFT then measure the difference between the phase angles to get the phase difference, e.g.

    mix = np.exp(-1j * omega * time)
p1 = np.angle(noisy_signal1 * mix)
p2 = np.angle(noisy_signal2 * mix)
delta_p = int((f * (np.mean(p2 - p1) / (2 * np.pi))) % f)