Phase lag between two signals

I'm trying to obtain the phase lag between two signals in Matlab.

The signals were obtained experimentally over time.

I'm using cpsd function but I get wrong values.

Im attaching an excel file with the signal.

Any help would be greatly appreciated

[Daniel M]

Use xcorr()

[U A]

This gives me the phase lag in a numeric way.

I need the phase lag in degrees

[Daniel M]

They look to be 180 degrees out of phase with each other visually, would you agree? I'm guessing you need a programmatic way to assess that.

Answer 1

[Daniel M]

Here's what I came up with using xcorr.

It calculates that Signal 2 lags behind Signal 1 by 179.6 degrees.

clearvars
close all
clc

[data,headers] = xlsread('SignalsPhaselag.xlsx'); % read data

% time and sampling frequency
t = data(:,1);
fs = 1/mean(diff(t));
x1 = detrend(data(:,2));
x2 = detrend(data(:,3));

% plot raw data
figure
plot(t,x1,t,x2)
xlabel('Time [s]')
legend('x1','x2')

[r,lags] = xcorr(x1,x2); % get the cross-correlation
figure
stem(lags,r)
xlabel('Lags [samples]')
ylabel('Cross-correlation')

% Do some frequency analysis using function below
[~,X1,F1] = plotFFT(x1,fs,[],[],1); title('x1')
[~,X2,F2] = plotFFT(x2,fs,[],[],1); title('x2')

% Find the peaks of the frequency spectrum
[pks1,locs1] = findpeaks(X1,F1,'SortStr','descend');
[pks2,locs2] = findpeaks(X2,F2,'SortStr','descend');
hold on
plot(locs2(1),pks2(1),'rv') % plot the peak

% Check if they have the same fundamental 
assert(isequal(locs1(1),locs2(1))); % i.e. peak at the same frequency

fundamentalFreq = locs1(1); 

[~,maxlags] = max(r); % position of maximum lag

delay = lags(maxlags)/fs; % convert to time (not in number of samples)

x3 = circshift(x2,lags(maxlags)); % shift the sample

% View the raw data and the shifted sample
% Notice how x3 lies on top of x1 nicely
figure
plot(t,x1,t,x2,t,x3)
xlabel('Time [s]')
legend('x1','x2','x3')

% Get the delay in terms of degrees
delayInDegrees = rad2deg(delay*fundamentalFreq*2*pi);

fprintf('\nSignal 2 lags behind Signal 1 by ~%.1f degrees\n',abs(delayInDegrees))



%%%%%%%%%%%%%%%%%%%% Helper function %%%%%%%%%%%%%%%%%

function [H,YY,ff] = plotFFT(x,fs,begtime,endtime,demean,plotType)
% [H] = plotFFT(x,fs,begtime,endtime,demean,plotType)
% This function takes a signal, x, and the sample rate, fs, and
% plots the FFT between begtime and endtime (defaults to
% beginning and end of signal). Times are in seconds, not
% indices.
%
% x must be a signal from one channel
% demean: 1 or 0 to detrend data. 0 is default
% plotType: 'plot' (default), 'semilogy', 'semilogx', 'loglog'
%
% Returns a handle to a figure, as well as YY, the single sided FFT, and ff
% the single sided frequency vector.

T = 1/fs;               % sampling period

if min(size(x)) ~= 1
    error('x must be a vector')
end
if nargin < 3 || isempty(begtime) || begtime < 0
    begtime = 0;
end
if nargin < 4 || isempty(endtime) || endtime > (length(x)-1)*T
    endtime = (length(x)-1)*T;
end
if begtime >= endtime
    error('begtime must be before endtime')
end
if nargin < 5 || isempty(demean)
    demean = 0;
elseif ~(demean==1 || demean==0)
    error('demean must be true or false')
end
if nargin < 6 || isempty(plotType)
    plotType = 'plot';
elseif ~any(strcmp(plotType,{'plot','semilogy','semilogx','loglog'}))
    error('Incorrect plotType')
end

% Get indices from beg/end times
begind = floor(begtime*fs+1);
endind = floor(endtime*fs+1);

X = x(begind:endind);   % new signal
if demean
    X = detrend(X);
end

Y = fft(X);             % fourier transform
L = length(X);          % length of signal
t = (0:L-1)*T;          % time vector
f = (0:L-1)/L * fs;     % frequency vector

% if bitget(L,1)
%     % L is odd
%     f((L+1)/2+1:end) = f((L+1)/2+1:end)-fs;
% else % L is even
%     f(L/2+1:end) = f(L/2+1:end)-fs;
% end

% shouldn't matter if odd or even, this should handle it right
f(ceil(L/2)+1:end) = f(ceil(L/2)+1:end)-fs;

% get positive frequencies only
ff = f(1:ceil(L/2));
YY = abs(Y(1:ceil(L/2)));

H = figure;
subplot(2,1,1)
plot(t,X)
xlabel('Time [s]')
ylabel('Signal')
axis tight

subplot(2,1,2)
switch plotType
    case 'plot'
        plot(ff,YY)
    case 'semilogy'
        semilogy(ff,YY)
    case 'semilogx'
        semilogx(ff,YY)
    case 'loglog'
        loglog(ff,YY)
end

xlabel('Frequency [Hz]')
ylabel('|P(f)|')
axis tight
grid on

% H = tightfig(H);

end

[Daniel M]

Using the cpsd function, the frequency of the peak is the same as doing an FFT (because fundamentally they do the same thing). So you can use either one.

But you still need to get the delay in the time-domain (using xcorr), and use the fundamental frequency and sampling rate to translate that into a shift in degrees.

[pxy,f] = cpsd(x1,x2,length(x1),0,length(x1),fs);

figure
plot(f,abs(pxy))
xlabel('Frequency [Hz]')
ylabel('cpsd')

[~,maxf] = max(abs(pxy));

if isequal(fundamentalFreq,f(maxf))
    fprintf('The fundamental frequencies are equal.')
end

Answer 2

[Linas Svilainis]

I am not sure whether it is the phase you are looking for. Would be better if I knew the application.

Nevertherless:

1. The signals you have are not pure sinusoid. Hovewer, if you know the exact frequency (most probably it is you who is injecting the probing sinusoid), you can use sine wawe correlation: SWCtruncated You get complex amplitude U1 and U2 for signal 1 and signal 2. Taking angle(U1) and angle(U2) gives you their phase in rad. Subtract and scale to geds by dividing by pi and multiplying by 180.
2. If you get time delay estimate and still know the frequency phase lag is needed for, then you can use subsample delay estimators and then convert to phase by scaling to period: GetTOFfftPhase or GetTOFcos(MySignal,RefSignal) Those estimate delay with subsample resolution and your SNR is good so resulting accuracy should be high. In order to get time, divide by signal sampling frequency. In order to convert time to phase in deg, divide by period and multiply by 360.