Non-linear system - calculating phase shift between two sine waves - Fourier and Hilbert

Question

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:

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?

Answer 1

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.

Answer 2

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.