3

Edit: From the suggested use of the xcorr function, which gave a result of zero, this data is more likely to represent a non-linear system than I had first thought. In this case, the expressions for a and b are simplified and do not necessarily compare entirely with the real data. Moreover, the value for the phase shift d can no longer be assumed to be a constant.

I have two arrays of raw data that correspond to two sinusoidal signals at the same (known) frequency over time, where one is phase shifted by an unknown amount delta, d. These signals can be described as:

a = A*sin(wt)
b = B*sin(wt+d)

Where w = 2*pi*f and A, B are the amplitudes of a and b, respectively.

The goal here is to evaluate the phase shift over time for the given signals. This could either be as an instantaneous value or as an averaged value over a number of cycles (which is small compared to the total number of cycles, i.e. at f = 150Hz a 10 second test is equivalent to 1500 cycles).

I am aware that there are many methods that can be used to evaluate such a problem and that many other users have already asked/answered questions on this subject. Below the methods with links to the original posts and the code that I have implemented. The problems I am having at the moment are:

  1. How to interpret the phase shift after it has been evaluated (for my scenario the phase shift is usually quoted as a single, positive number at a point in time that stays fairly constant over many cycles)
  2. Differences in value between different methods
  3. Are there any additional methods that I've missed?

For my data with a and b both sized [1 x 399] and covering 20 cycles. This is a small portion of the expected data and I'm using it to test out the different methods. It would be reasonable to expect a phase shift of at least 10 degrees although it is not known exactly:

Fourier Transform

Source: How do I calculate the phase shift between two sinusoidal signals?

fft_a = fft(a);
fft_a = fft_a(2:end);       %Drop the first point
angle_a = angle(fft_a);     %Angle the result

fft_b = fft(b);
fft_b = fft_b(2:end);       %Drop the first point
angle_b = angle(fft_b);     %Angle the result

ps1 = rad2deg(abs(angle_a - angle_b));     %Phase shift calculation

Using this method - what can I do to eliminate the spikes at the beginning/end? Do I need to set the the transform to perform over a certain number of points?

And how can I then plot the phase shift against time as surely this result is now in the frequency domain?

Hilbert Transform

Source: Identifying phase shift between signals

ha = hilbert(a);                 %Hilbert transform
hb = hilbert(b);

ps2 = rad2deg(angle(hb./ha));    %Phase shift calculation

Here the resultant phase shift over time is both negative and positive and seems to oscillate like a waveform. It's fairly consistent in shape but again, I do not know how to interpret the results. After taking the abs value of the phase shift, the result now varies between 0 and 30 degrees.

Graph: Hilbert Phase Shift Time Plot

For both these methods, how can I use the results to quote the phase shift as being a single value over a set period of cycles? Taking the mean doesn't seem like an exact method.

For a quick calculation of the phase shift I used:

ps3 = rad2deg(acos(dot(a,b)/(norm(a)*norm(b))))    %Norm product method

ps3 =

    11.8289

I'm fairly new to this type of analysis and a bit of a novice at Matlab so any corrections/guidance would be greatly appreciated.

Edit: With the knowledge that this is a non-linear system, what methods can be used to best evaluate the phase shift?

Community
  • 1
  • 1
Natalie
  • 31
  • 5
  • Have you looked at `fminsearch`? Just subtract one signal from the other using a parameterizable phase offset, optimizing for the smallest sum of squared differences – Rody Oldenhuis Dec 12 '16 at 13:00
  • ...fairly sensitive to noise though. Plus, if the signals have different amplitudes, you need to normalize those first – Rody Oldenhuis Dec 12 '16 at 13:06
  • Is this an alternative method to using the `xcorr` function? Also it would seem like I could evaluate the phase offset for each pair of points rather than across the entire arrays....thank you @rodyoldenhuis for your suggestion I shall definitely try it out. – Natalie Dec 12 '16 at 14:45

2 Answers2

3

What about cross-correlation?

with cross-correlation you should find an absolute maxima in the shift, since when you shift one of the two signal of the right amount, the two will become the same.

a = A*sin(wt)
b = B*sin(wt+d)
[c, time_c] = xcorr(a,b)
[m, i_m] = max(c)
d_ext = time_c(i_m) * T * w % T is the sampling time

it should work for a constant d, i'm not sure for a variable one.

bracco23
  • 2,181
  • 10
  • 28
  • I just tried using the `xcorr` function and it gave me an answer of zero.... does that mean that a phase shift does exist but it's being averaged to zero by using this function? I'm not too familiar with this function or how it operates to be certain of why I got a zero answer. – Natalie Dec 12 '16 at 14:35
  • The function only implement the cross-correlation, you can find more on the [wikipedia page](https://en.wikipedia.org/wiki/Cross-correlation). What exactly is zero? have you tried the cose above? It worked fine when i tried it with a constant shift – bracco23 Dec 12 '16 at 14:44
  • I think this has more to do with the type of data that I'm using rather than which method is 'better'...and yes, you're right in that it should work fine for a constant value for `d`. – Natalie Dec 12 '16 at 15:03
  • There is an interesting paragraph on non-linear systems and zero-valued results on the Wikipedia page you referenced: "...This problem arises because some quadratic moments can equal zero and this can incorrectly suggest that there is little "correlation" (in the sense of statistical dependence) between two signals, when in fact the two signals are strongly related by nonlinear dynamics. " Therefore it's sensible to conclude that there is a non-linear phase shift here. Any thoughts on how to approach this? – Natalie Dec 12 '16 at 15:31
  • It should be noted that `xcorr` is part of the [Signal Processing Toolbox](https://nl.mathworks.com/help/signal/index.html) – Rody Oldenhuis Dec 12 '16 at 16:14
0

What I mean is this: 

clc

% Noiseless example data
x = 0:0.1:40*pi;

p1 = deg2rad(+5.74);
p2 = deg2rad(-4.28);

y = sin(x + p1);
z = sin(x + p2);

input = rad2deg(p1 - p2) % display injected phase difference


% Retrieve values using no assumptions 

options = optimset('TolX'   , 0,...
                   'TolFun' , 0, ...
                   'display', 'off');

p1 = fminsearch(@(q) sum( (y - sin(x + q)).^2 ), 0, options);
p2 = fminsearch(@(q) sum( (z - sin(x + q)).^2 ), 0, options);

computed = rad2deg(p1 - p2) % display computed phase difference

For this example, the output is: 

input =
    1.002000000000000e+001
computed =
    1.002000000000000e+001

But, as noted, this method is quite simplistic and the accuracy will decrease much faster than one of the other methods when noise comes into play.

Rody Oldenhuis
  • 37,726
  • 7
  • 50
  • 96