What is the best way to get the phase between noisy signals

I have a few sinusoidal signals at approximately the same frequency. I would like to know the phase difference between the three signals. What is the best way to do this? It would also be helpful to attach an uncertainty to these phase measurements. I've attached a picture of the signals

• If it is "approximately the same frequency", then the average phase difference for three signals will be 120 degrees. – jippie Sep 3 '14 at 16:55
• Modified the question to hopefully make it a bit more clear – cpc333 Sep 3 '14 at 17:01
• Make a square wave from the sine waves using a Schmitt-trigger input. Then use a timer to measure delay between rising edges of the different signals. – jippie Sep 3 '14 at 17:01
• I'd like to measure it from old data. So no hardware solutions allowed. I suppose I could try to replicate that in software maybe? – cpc333 Sep 3 '14 at 17:03
• How much time can I have to process it? (Oh and what's the frequency?) – George Herold Sep 3 '14 at 18:59

The obvious answer would be to take the FFT of each signal. The phase angle of the highest (magnitude) bin gives you the phase of the fundamental sinewave buried in each signal. Compare and contrast as needed.

Use an appropriate window function to minimize edge effects.

I recommend transfer-function estimation techniques such as the Welch Periodogram, used in conjunction with coherence estimates, which can be used to put confidence intervals around your estimate. I recommend Bendat and Piersol, Random Data: Analysis and Measurement Procedures for a very thorough treatment.

The Welch periodogram involves taking the spectrum of overlapping segments of data to yield the periodogram.

Coherence pretty much indicates the noise and nonlinearities inherent in the relationship between input and output. A coherence of 1.0 reflects a perfectly linear noise-free system. Both noise and nonlinearity reduce coherence.

You can also plot one signal against the other, and calculate phase from the resulting Lissajou figure

Another option would be to examine the peaks in the cross correlation between the two signals.

That appears to be a matplotlib graph correct? This would imply you are using python to post-process. Equally you want to post-process already collected data.

What I would suggest is apply a moving average filter to each signal. While such a filter will "phase shift" an individual waveform with respect to the original, relative to each other (post-filter) there will be no difference as the same "shift" will have been applied

A simple

1/4x + 1/4(x-1) + 1/4(x-2) + 1/4(x-3) should do or a simple numpy.mean with a fix slice window and varying index.

Either way once you have lightly filtered the data to minimise zero crossing jitter you can utilise a couple of numpy functions:

zero_crossings = numpy.where(numpy.diff(numpy.sign(DATA)))[0]

This will return an array of all the zero crossings. You may need to increase the filtering if there are multiple very close together.

If you have matlab a similar scheme can be achieved using the filter command & there is a crossing.m file publically available from mathworks exchange.

NOTE this won't give a distinction for positive going or negative going crossings.

--EDIT--

Take this piece of python code:

#!/usr/bin/env python3

import numpy as np
from pylab import *
import matplotlib
matplotlib.pyplot.switch_backend('QT4Agg')

def filt(data, l=4):
tmp = np.zeros(len(data))
for i in range(l,len(data)):
tmp[i] = np.mean(data[i-l:i])
return tmp

t = np.arange(0,10*np.pi,0.1e-3)
s = np.sin(t) + np.random.normal(0,0.01,len(t))

crossings = np.where(np.diff(np.sign(s)))[0]
print(len(crossings))
s2 = filt(s,l=4)
crossings = np.where(np.diff(np.sign(s2)))[0]
print(len(crossings))
s3 = filt(s,l=10)
crossings = np.where(np.diff(np.sign(s3)))[0]
print(len(crossings))
s4 = filt(s,l=150)
crossings = np.where(np.diff(np.sign(s4)))[0]
print(len(crossings))

filt(s)
plot(t,s)
grid(True)
hold(True)
plot(t,s2)
plot(t,s3)
plot(t,s4,linewidth='2')
show()


This generates a 300,000 datapoint sin wave and couples onto it a normally distributed noise source with a sd of 0.01

As you can see a lot of zero crossing noise, an annoying amount.

The rest of the script does a really crude moving average (a better implementation would be a movingAverageFilter & this is easy in a matlab script via the "filter" command. (http://www.mathworks.co.uk/help/matlab/data_analysis/filtering-data.html)

Each iteration uses a different window size & equally prints out the total number of zero crossings detected.

[jrb@FluidMotion:~/tmp] 27s 1 \$ ./test.py
1126
272
120
12


The number of detected zero crossings has fallen from 1126 to 12

As mentioned there is an absolute phaseshift (would be significantly reduced if a better filter was used, but this was a 5min example), but for comparing relative differences it is valid as all datasets are shifted by the same amount.

Once you know the indexes of the zero crossings and you know the timestep of one sample, you are able to calculate the time difference between zero crossings. Couple that with the number of samples between three zero crossings of the same data & you know the number of samples per 360degrees

This took 150samples out of 300,000 (to produce a reasonable zero-crossing detection but the source had alot of random zero crossings. The OP data appears to only have a couple of zero-crossing events.

• What makes you think that "light" filtering is going to significantly reduce the zero-crossing jitter? Have you ever actually used this method? – Dave Tweed Sep 3 '14 at 20:36
• yes I have used this method. I actually have been using it daily for the last couple of weeks with both matlab and python implementations as I have been trying to track down the source of some orthoganality errors – JonRB Sep 3 '14 at 20:45
• as to what makes me think "light" filtering will work... It really comes downto how many multiple crossings there are to then drive how "light" the filtering needs to be. For instance some data I took last week, 10million datapoints, 20 cycles worth I needed 256samples moving average filter. Another set at 30,000 with only a couple of additional crosses needed 5 – JonRB Sep 3 '14 at 20:51
• My point is that even if you use a low-pass filter to eliminate the "multiple crossings", you're ignoring the fact that the low-frequency components of the noise -- including any sort of DC offset -- will still shift the timing of the zero-crossings significantly, making any sort of "phase" estimate moot. – Dave Tweed Sep 4 '14 at 3:54
• You would obviously have to remove the dc component ( data = data - np.mean(data) – JonRB Sep 4 '14 at 6:38

You could aggressively band-pass filter the signals, since it sounds like you know the approximate frequency you're interested in. Low pass to get rid of the noise, and high pass because different peaks of a single signal are at wildly differing amplitudes. Then compare the times of the zero crossings.