Here's some Python code I wrote once to test out the max correlation method. The idea here is to perform correlations between the two signals, and then whichever phase angle produces the maximum such correlation is the winner. Naturally, more noise would make this process more difficult: I would recommend a median filter for smoothing data before trying this, although there are loads of options.
import numpy as np
def float_equal(x, y) -> bool:
"""
This function measures whether two floating-point numbers are
equal, to within a certain tolerance
:param x:
:param y:
:return: Whether the two numbers are equal, to within machine
epsilon.
"""
tol = 0.0000001
return np.abs(x - y) < tol
# Define the array length, and the begin and end of the x location
# array.
array_length = 10000
begin_x_location = 0
end_x_location = 2 * np.pi
# x array of values from 0 to 2Pi.
signal1_x = np.arange(begin_x_location,
end_x_location,
end_x_location / array_length)
# Need a basic cosine.
signal1_y = np.cos(signal1_x)
# A shifted cosine.
signal2_y = np.cos(signal1_x - np.pi / 4)
# Add some noise to each y signal above, so that the correlations aren't
# perfect. Want the magnitude to be significantly less than the signal.
full_signal1 = signal1_y + np.random.normal(0, 0.01, array_length)
full_signal2 = signal2_y + np.random.normal(0, 0.01, array_length)
# Debugging how to do the correlations.
print(str(len(full_signal1[:array_length-1250])))
print(str(len(full_signal2[1250:])))
# Note that 1250 corresponds to pi/4, the shift. We expect this to be
# the maximum correlation.
corr = np.corrcoef(full_signal1[:array_length-1250],
full_signal2[1250:])
# The correlation at the known phase shift.
print(str(corr[0][1]))
print("Now we find max correlation coefficient:")
# Create a vector of correlation coefficients. Each coefficient is
# generated by shifting one array relative to the other. We need the
# two arrays to be the same length, so we simply trim off the unwanted
# numbers to get two equal-sized arrays.
corrs = np.array([np.corrcoef(full_signal1[:array_length-i],
full_signal2[i:])[0][1]
for i in range(array_length-100)])
# This is the max correlation coefficient.
print('Maximum correlation: ' + str(np.max(corrs)))
# This is where that maximum occurred. Use the float_equal function.
index_location = np.where(float_equal(corrs, np.max(corrs)))[0][0]
print('Maximum occurred at index: ' + str(index_location))
# Now convert index to x location and print out where that happened.
# Note that the answer is in radians.
x_location = (end_x_location - begin_x_location) * index_location \
/ array_length + begin_x_location
print('Maximum occurred at x value: ' + str(x_location))
```