Get lag with cross-correlation?

Question

Let's say have have two signals:

import numpy
dt = 0.001

t_steps = np.arange(0, 1, dt)
a_sig = np.sin(2*np.pi*t_steps*4+5)
b_sig = np.sin(2*np.pi*t_steps*4)

I want to shift the first signal to match the second signal. I know this can be completed using cross-correlation, as evidenced by Matlab, but how do I accomplish this with SciPy.

Inbeing

Answer 1

See some examples first. Assume we are in unit tests class already.

# Autocorrelation.
y1 = [1, 1, 0, 0, 1, -1, -1]
corr, lag = cross_corr(y1, y1)
self.assertEqual(lag, 0)

y1 = [1, 1, 0 ,1, -1, -1]
y2 = [1, 0, 1, 0, 0, 2]
corr, lag = cross_corr(y1, y2)
self.assertEqual(lag, -2)

here is my code.

import numpy as np    

def cross_corr(y1, y2):
  """Calculates the cross correlation and lags without normalization.

  The definition of the discrete cross-correlation is in:
  https://www.mathworks.com/help/matlab/ref/xcorr.html

  Args:
    y1, y2: Should have the same length.

  Returns:
    max_corr: Maximum correlation without normalization.
    lag: The lag in terms of the index.
  """
  if len(y1) != len(y2):
    raise ValueError('The lengths of the inputs should be the same.')

  y1_auto_corr = np.dot(y1, y1) / len(y1)
  y2_auto_corr = np.dot(y2, y2) / len(y1)
  corr = np.correlate(y1, y2, mode='same')
  # The unbiased sample size is N - lag.
  unbiased_sample_size = np.correlate(
      np.ones(len(y1)), np.ones(len(y1)), mode='same')
  corr = corr / unbiased_sample_size / np.sqrt(y1_auto_corr * y2_auto_corr)
  shift = len(y1) // 2

  max_corr = np.max(corr)
  argmax_corr = np.argmax(corr)
  return max_corr, argmax_corr - shift

Creamy

Is it possible for this to work with y1 and y2 with different length like in MATLAB. – Jarodjarosite

@AlexRamses Sure. It is possible. But the question is where do you want them to be aligned (position of lag=0)? Left end, right end, middle? It's a little bit tricky, so I didn't do that. – Creamy

Answer 2

Here is an example code to get the lag of cross-relation using SciPy.

from scipy.signal import correlate
from scipy.signal import correlation_lags

x = np.asarray([1,2,3,4])
y = np.asarray([.5,1,2,3])
lags = correlation_lags(x.size, y.size, mode="full")
lag = lags[np.argmax(correlation)]
print(lag)

Please see the following links to know more about the correlate and correlation_lags from Scipy.

[1] https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.correlation_lags.html

[2] https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.correlate.html

Colon

Answer 3

For completeness: Scipy.stats and Scipy.signal offers the functionality that you are looking for, as do several of the Scipy-derived packages such as astropy etc.

Exhibitionist

Which functions exactly should be used? – Inbeing

@Inbeing Without wishing to avoid the issue, I suggest looking at: en.wikipedia.org/wiki/Correlation_coefficient and the disambiguation page on Wikipedia to see which you could to use. Scipy has scipy.signal.correlation_lags, but lthe general scipy.signal page which sets out the various options. Astropy has similar functionality for example at specutils.readthedocs.io/en/stable/analysis.html but this is more directed at astronomy/astrophysics and uses specific astropy/specutils interfaces. However it is very powerful. I hope that this helps. – Exhibitionist