1
$\begingroup$

I want to compute the phase shift between two 1-D signals of same frequency, but before I'm trying to compute the time shift between. The cross-correlation function seems to be ideal for that but I'm confused on how to interpret scipy cross-correlation. Let's take two sinus with a frequency f0 = 200 Hz, a sample frequency fs = 10000 Hz, playing during 0.1s and with a phase difference of pi.

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# We build the two sinus
t = 0.1
f = 200
fs = 10000
x = np.arange(0, t, step = 1/fs)
sin1 = np.sin(2 * np.pi * f * x)
sin2 = np.sin(2 * np.pi * f * x + np.pi)

# We compute the scipy cross-correlation    
scorr = signal.correlate(sin1, sin2)
plt.plot(scorr)
plt.title('Scipy Correlation')
plt.show()

With this last graphic output I get:

enter image description here

Now a few questions:

  1. We have a x-axis spanning on 2 * fs, as the function is hermitian, I guess that we have the hermitian symmetry? In this case we could just shift it around zero (like fftshift) and only consider the positive axis, right ?

  2. If 1) is ok, does my x time vector could fit the x-axis of my cross-correlation ?

  3. Finally, and not necessarily related to previous questions, how to read x and y axis ? There is not much practical documentation on cross-correlation product, the only thing I know is that we have to look where the function takes its maximum in order to get the time lag between the two signals. For me the y-axis is just the result of the product of the two signals as in the formula (cross-correlation) (but I don't get why the product of two sinus with amplitude 1 could ouput 500 ...), and the x-axis gives the indice corresponding to the time difference ( and in this case, the indice where the function takes its max corresponds to the time shift I am searching for, hence the utility of plotting the cross-correlation with a fitted x time vector). Is it right or do I misunderstand something ?

$\endgroup$
0

1 Answer 1

1
$\begingroup$
  1. No, the output is len(x)*2-1 long, an odd number
  2. I don't understand the question
  3. The x axis is the delay in samples, and the y axis is the cross-correlation. The number of x samples is odd, and the middle sample represents 0 delay.

If you cross-correlate the sin with itself, you will see a peak at sample 999, which is the middle sample, which represents 0 delay. With the signal you've shown, the peak is at sample 1024, which represents a time delay of 1024-999 = +25 samples. There is also a peak at sample 974, which represents a time delay of 974-999 = −25 samples. Sine waves are repetitive and will show multiple peaks as they line up with themselves at different lags.

If you do it again with a white noise signal (randn()), you will see only a single peak and it will be clearer.

$\endgroup$
4
  • $\begingroup$ 1. Effectively, but something is not clear as Scipy documentation says that the function outputs a N-dimensional array from two N-dimensional input arrays. 2. I just wonder if it's possible to set a x-axis with time values and not samples (as we set a frequency x-axis for a FFT). $\endgroup$
    – terzan5
    Aug 18, 2021 at 13:42
  • 1
    $\begingroup$ @terzan5 1. Yes, N-dimensional, not N-length. If you are feeding it 1-dimensional arrays, you will get a 1-dimensional array back. 2. If delay is 25 samples and fs is 10 kHz, then delay is 25 samples / (10000 samples/sec) = 2.5 ms $\endgroup$
    – endolith
    Aug 18, 2021 at 14:10
  • $\begingroup$ I was computing it and that's ok. The delay is effectively 2.5 ms which is exactly 1/(2 * f0) = T0/2 which significates that the signals have opposite phases. Thanks endolith ! $\endgroup$
    – terzan5
    Aug 18, 2021 at 14:13
  • 1
    $\begingroup$ @terzan5 You can also do interpolation to get a better estimate of the peak. gist.github.com/endolith/376572 $\endgroup$
    – endolith
    Aug 18, 2021 at 14:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.