phase difference between two sines with wt and 2wt frequencies

Hi all,

I have two data vectors (x,y) that very fairly follow these equations:

how can I calculate ϕ from the data? At the moment I have been fitting two dummy sine functions with different ϕ until it matches the data shape (x,y). I also tried to fit sine functions but I don't know how to relate the phase angle I get for x and y. I think I'm missing something really obvious here.

Regards,

Answer 1

[Star Strider]

I usually let the Symbolic Math Toolbox do the ‘heavy lifting’ in situations like these (particularly because it doesn’t make the stupid algrbra errors I’m prone to making):

syms  Ax Ay phi t w x y  
Y = Ay * sin(w*t);
X = Ax * cos(2*w*t + phi);

phi_sln = solve(X == Y, phi)

producing:

phi_sln =
  acos((Ay*sin(t*w))/Ax) - 2*t*w
- acos((Ay*sin(t*w))/Ax) - 2*t*w

This assumes that you have some idea of what and and ω are, however that should be apparent from the magnitudes of the data.

A numerical parameter estimation approach:

t = linspace(0,5).';                                                            % Column Vector

Ax = 2;                                                                         % Define Parameter
Ay = 3;                                                                         % Define Parameter
w = 5;                                                                          % Define Parameter
phi = pi/4;                                                                     % Define Parameter

y = Ay*sin(w*t) + randn(size(t))*0.05;                                          % Create Data
x = Ax*cos(2*w*t + phi) + randn(size(t))*0.05;                                  % Create Data

objfcn = @(b,t) [b(1).*sin(b(2).*t), b(3).*cos(2*b(2).*t + b(4))];              % Objective Function

ypks = islocalmax(y);                                                           % Estimate ‘w’
estw = 2*pi/(mean(diff(t(ypks))));                                              % Estimate ‘w’

B = fminsearch(@(b) norm([y, x] - objfcn(b,t)), [max(y); estw; max(x); rand]);  % Estimate Parameters

yxest = objfcn(B, t);                                                           % Fitted Data

figure
plot(t, y,    t, x)
hold on
plot(t, yxest, '--') 
hold off
grid
legend('y', 'x', 'y_{est}', 'x_{est}')

Note that ‘x’ and ‘y’ have to be essentially noise-free for this to work, so some filtering may be necessary to remove high-frequency and low-frequency noise if any exist in your signals.

EDIT — (26 Dec 2020 at 17:36)

A much more robust estimator, especially in the presence of noise, for the radian frequency ‘w’ is:

yzx = find(diff(sign(y)));                                                      % Estimate ‘w
estw = pi/(mean(diff(t(yzx))));                                                 % Estimate ‘w’

I am leaving my original code unchanged, however this version is likely a better option.

.

[Daniel Mella]

Hi Star Strider,

Thank you for your fast reply and sorry for the late response. It took me a while to understand your code since I am not very familiar with handles in Matlab.

The first part of your code was really helpful since I didn't know you could do that. Unfortunately, I am not sure that doing solve(X==Y,phi) is the right approach since X is not equal to Y. In theory, we can solve by using only solve(X,phi) but phi is dependent on the relative phase between X and Y.

The second part doesn't give me accurate results. The values of B using my data are [0.0004 , 18.11 , 0 , 0.282], which do not match with the amplitudes or frequency of the signal.

Since I am pretty sure I am missing something important. I uploaded two data examples I am dealing with. The time step is the same for both signals. The goal is to find the phase angle between X and Y.

Thank you for your help. You have gave me new insights into Matlab I was not aware about.

PD. I edited my original question to change a cosine for a sine.

[Star Strider]

My pleasure!

The first part was just a symbolic exploration. The objective function in the second part only uses those relationships to determine the frequency ‘w’, since according to them, the second function frequency is twice the first functon frequency.

With respect to the ‘data.mat’ file, does ‘x1’ match ‘y1’ and ‘x2’ match ‘y2’? If not, what should I do with these data? I need to know what you want to do.

[Daniel Mella]

Thanks!

The data I attached is just to give a few examples I am trying to analyse. (x1,y1) and (x2,y2) are the displacement of a structure under different conditions. I am trying to explain the difference between the trajectory traced under these two conditions. The first one (x1,y1) have a trajectory that resembles an eight, while (x2,y2) has a more complex trajectory. Given the harmonic nature of the displacement in x and y under both conditions, I thought that the difference could be explained in terms of phase difference between x and y.

In other words, if (x,y) follow this relationship

Then, the difference between (which should be close to zero) and should explain the difference in the observed trajectory. Now, I am not sure anymore.

[Star Strider]

My pleasure!

The problem was using fminsearch on your data. It worked fine with my example data, however it failed with the complexity of your actual data.

The (slightly revised) code:

D = load('data.mat');

t = D.ts;
x = D.x1;
y = D.y1;

objfcn = @(b,t) [b(1).*sin(b(2).*t), b(3).*sin(2*b(2).*t + b(4))];              % Objective Function

yzx = find(diff(sign(y)));                                                      % Estimate ‘w
estw = pi/mean(diff(t(yzx)));                                                   % Estimate ‘w’

Binit = [max(y); estw; max(x); rand];

% B = fminsearch(@(b) norm([y, x] - objfcn(b,t)), Binit);  % Estimate Parameters
B = lsqcurvefit(objfcn, Binit, t, [y x]);                                       % Estimate Parameters

yxest = objfcn(B, t);                                                           % Fitted Data

figure
plot(t, y,    t, x)
hold on
plot(t, yxest, ':', 'LineWidth',2) 
% plot(t(yzx(1:2:end)), y(yzx(1:2:end)), 'r+')
hold off
grid
legend('y', 'x', 'y_{est}', 'x_{est}')
title(sprintf('$y(t) = %7.4f\\cdot sin(%5.2f\\cdot t)$\n$x(t) = %7.4f\\cdot sin(%5.2f\\cdot t%+6.4f)$',[B(1:3); B(2)*2; B(4)]), 'Interpreter','latex')
xlim([-50 -40])                                                                 % Zoom In To See Detail

With the appropriate data substitutions in runs for each variable pair:

Phi_1 = -0.2290826914261480
Phi_2 = +0.3623175895494797

with the difference:

Phi_Dif = -0.591400280975628

or more accurately:

Phi_Dif = ±0.591400280975628

depending on how they’re compared.

Plots for —

and for —

demonstrate that the fits are appropriate with respect to frequency and phase, although slightly less so with respect to amplitude. Note — these are ‘zoomed’ to show detail that would not be apparent otherwise.

[Daniel Mella]

Hello! I have been playing around with your code but I think it gives different phase angles depending on its initial random value. I attached a screenshot where I runned it several times.

I did some small changes, like using "findpeaks" to calculate the amplitudes (Ax,Ay) and changing the objective function from [y,x] to [x,y]. With these small changes and selecting an initial value of the phase angle equal to zero, I got these results

Ax = -0.002
w = 5.4285
phi = -1.1565
Ay = -0.0021

Which is also different from previous results. Later, I decided to "fix" the the amplitudes (Ax,Ay) and frequency (w), and bound the phase angle from -pi to pi. These results are

Ax = 0.0024
w = 5.4328
phi = 1.59
Ay = 0.0033

Unfortunately, when I plot them I get quite different trajectories:

I noticed that the trajectory is inverted regarding the y-axis (the estimated trajectory is "looking" right while the original data is skewed towards left). So, I changed Ax to -Ax from the initial value and I got better results.

Here is the code I am using

%load data
load('data.mat')

%adjust ts. Start time at zero
ts = ts - ts(1);

%Amplitude of x and y. Mean peak values
[a,b] = findpeaks(x1);
Ax = mean(a);
[a,b] = findpeaks(y1);
Ay = mean(a);

%y main frequency (fx = 2*fy)
[p,f] = pwelch(y1,length(y1),0,[],1/(ts(end)-ts(end-1)));
fy_max = f(find(p==max(p)));
w = 2*pi*fy_max;

%---------1-----------%
%Objective function [x,y]
objfcn = @(b,t) [b(1).*sin(2*b(2).*t+b(3)),b(4).*sin(b(2).*t)];

%initial values
Bini = [-Ax;w;0;Ay];

%fit
B_1 = lsqcurvefit(objfcn, Bini, ts, [x1,y1])
yxest = objfcn(B_1, ts);

%plot estimate and original values
figure(1)
plot(yxest(:,1),yxest(:,2))
hold on
plot(x1,y1)

%---------2-----------%
%Fit by fixing Ax w and Ay. Bound phase angle from -pi to pi
%initial values
Bini_2 = [Ax;w;-pi;Ay];

%fit
B_2 = lsqcurvefit(objfcn, Bini, ts, [x1,y1],Bini_2-0.000001,[Ax;w;pi;Ay]+0.000001)
yxest = objfcn(B_2, ts);

figure(2)
plot(yxest(:,1),yxest(:,2))
hold on
plot(x1,y1)

[Star Strider]

Nonlinear paramerer estimation routines usding gradient-descent are extremely sensitive to the initial parameter estimates provided to them, and will converge on the minimum closest to the initial estimates, whether it is a local minimum or the global minimum. This is inherent in all gradient-descent algorithms, of which lsqcurvefit is one. The initial parameter estimates I use for this are calculated from the data themselves, and so will give a reasonably accurate (and reproducable) estimate for the parameters when lsqcurvefit converges.

I also suggest that you examine the results of the fit itself, such as the norm of the residuals, to see how well a particular set of parameters fits the data. Lower residual norm values indicate a better fit. That may help you assess the various parameter estimates.

It might be worthwhile of you to use the ga (genetic algorithm) function to search the entire parameter space for more appropriate parameters, if this is an important consideration. It would be relatively straightforward to modify my ‘objfcn’ to work with it. (I have posted several Answers that describe my approach to writing ga code for these problems.) Note that ga will likely require several runs to converge on a successful set of parameters, however for well-posed problems, it will find them.

[Daniel Mella]

Thank you for your reply. I have learn a lot !

At the end I decided to fit a sine function to x and y

where each parameter (amplitudes, frequency and phase angles) are given by the fit function using 'sin1'. Then, I removed the phase angle in the y function

which implies a change in as

I tested different until The trajectory matched. Unfortunatelly, this approach is manual and I haven't decided on a way to see how accurate the match is. For the (x1,y1) data, the phase angle is -0.235 (or is it -0.23? I cant say for sure) .

[Star Strider]

I would not put phase variables in both terms, since the actual phase difference is then not easily determined.

I provided a reliable way to determine the phase, based on the actual data and reasonably reliable nonlinear parameter estimation techniques. The parameter estimates are what the functions determined, and while the phase ϕ may be shifted by integer multiples of π, they remain reliable.

The results I got for ϕ I posted in my previous Comment. Numerical computations can vary slightly between runs. It is straightforward to use nlparci with the results of lsqcurvefit (it will be necessary to request additional outputs from lsqcurvefit, such as the residuals and the Jacobian) to see what the 95% confidence limits are for ϕ. That will give you a reasonable idea of its accuracy, and whether it is statistically different from 0 (if it is not, the confidence limits will have oppposite signs, and it is not necessary to include it in the model).

If my Answer helped you solve your problem, please Accept it!

.