I have a system that's sampling two 50Hz sine waves at 1kHz. What's the best way to estimate the phase shift between the two?
The measurements are somewhat noisy and the sampling period is not very even. I know the sample times to a high resolution but can't control them to be evenly spaced. The phase shift needs to be estimated on chunks of data covering approximately one second.
The measurement noise makes the zero-crossing method unreliable; if the noise causes a measured signal to cross zero twice at one "real" zero-crossing then the algorithm breaks down somewhat.
At present I'm doing something like this (python code):
t, v, a = [NumPY arrays of samples] p_inst = v * a s = numpy.sign(p_inst) s_avg = numpy.trapz(s, x=t) / (t[-1] - t) shift_fraction = s_avg / 2 + 0.5 shift_angle = (1 - shift_fraction) * numpy.pi
s is the sign of the product of the two sine waves. Finding the average of this using
numpy.trapz gives a value representing how often the two signals have the same sign, with 1 meaning "all of the time", -1 meaning "none of the time" and zero meaning "half of the time". This is converted to a fraction of the time the signals have the same sign (
shift_fraction) and this is converted to a phase-shift angle.
The biggest problem with this is that it doesn't tell you whether the phase difference is leading or lagging.
Does anyone have a better method?
Edit For those suggesting finding the correlation peak, am I right in thinking that, with 20 samples per cycle, that method will necessarily limit the resolution of the estimate to 18 degrees? What would you think of instead calculating
numpy.correlate(v, -a), finding the first zero crossing and then interpolating to estimate the actual zero point? I think that, for sinusoidal signals, the correlation should also be sinusoidal, and so the small angle approximation should make the interpolation reasonably good.
Can anyone compare the correlation and demodulation methods?
Further edit And, for those giving the demodulation method, am I right in thinking that if I don't care about the actual phase difference, just the power factor, then I can do this:
s = sin(t * 2 * pi * f) c = sin(t * 2 * pi * f) a_s_d = s * a v_s_d = s * v a_c_d = c * a v_c_d = c * v s_a = numpy.trapz(a_s_d, t) c_a = numpy.trapz(a_c_d, t) s_v = numpy.trapz(v_s_d, t) c_v = numpy.trapz(v_c_d, t) # cos ( a - b ) = cos a * cos b + sin a * sin b power_factor = c_a * c_v + s_a * s_v
Or can I not rely on
0 < s_a, c_a, s_v, c_v < 1?