I am in the process of identifying a time shift between two sinusoidal signals in Python. A little research on this has revealed that correlation between the two input signals can solve the problem. However, when working with correlations, I do not understand the result of np.correlate
and scipy.signal.correlate()
as the result is different from the a priori result.
Let
t = np.arange(0, 10, 0.01)
sig1 = np.sin(t)
sig2 = np.cos(t)
I'd expect pi/2 = 1.571 difference between both signals. However, the following code determines |1.490|.
import scipy
import numpy as np
t = np.arange(0, 10, 0.01)
sig1 = np.sin(t)
sig2 = np.cos(t)
# Scipy
#x = scipy.signal.correlate(sig1, sig2, mode='full')
#tx = signal.correlation_lags(len(sig1), len(sig2))
# Numpy
x = np.correlate(sig1, sig2, 'full')
tx = np.linspace(-t[-1], t[-1], 2*len(t)-1)
ndx = np.argmax(np.abs(x))
print(tx[ndx])
>>> -1.4900000000000002
Do I miss something with the modes of the correlation?
However, the cosine similarity or least squares optimization gives the correct result. Also the simple comparison of the extrema of both signals works correctly.
Full code including the plot:
import scipy
import numpy as np
import matplotlib.pyplot as plt
fig, ax = plt.subplots(3, figsize=(8,6))
t = np.arange(0, 10, 0.01)
sig1 = np.sin(t)
sig2 = np.cos(t)
# Plot 1: Both input signals
ax[0].plot(t, sig1, label='sig1')
ax[0].plot(t, sig2, label='sig2')
ax[0].legend()
ax[0].set_xlabel('time in s')
## Apply correlation between the two given signals
# Scipy
#x = scipy.signal.correlate(sig1, sig2, mode='full')
#tx = signal.correlation_lags(len(sig1), len(sig2))
# Numpy
x = np.correlate(sig1, sig2, 'full')
tx = np.linspace(-t[-1], t[-1], 2*len(t)-1)
# Plot 2: Correlation between two input signals
ax[1].plot(tx, x, label='correlation', c='k')
ndx = np.argmax(np.abs(x))
ax[1].plot(tx[ndx], x[ndx], 'rx')
ax[1].legend()
# Plot 3: Both input signals, sig2 on corrected time values
ax[2].plot(t, sig1, label='sig1')
t_corr = t - tx[ndx] # Correct the time
ax[2].plot(t_corr, sig2, label='sig2')
ax[2].set_xlabel(r'Shift: %.5f' %(tx[ndx]))
ax[2].legend()
print(tx[ndx]) # Display the shift
Your code is correct, the problem is the correlation. This method assumes that you have an infinite signal. Then, the shorter is the signal, the larger will be the error.
I have executed your code with a longer t
vector and the result has been |1.569999999999709|, which is more similar to pi/2 = 1.571.
t = np.arange(0, 10000, 0.01)
If you are sure that the signals are sinusoidal, I would recommend you to use another method of the ones that you have cited.