Phase lag between two signals
93 views (last 30 days)
Show older comments
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
3 Comments
Daniel M
on 28 Nov 2019
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.
Accepted Answer
Daniel M
on 28 Nov 2019
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
1 Comment
Daniel M
on 28 Nov 2019
Edited: Daniel M
on 28 Nov 2019
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
More Answers (1)
Linas Svilainis
on 2 Dec 2019
I am not sure whether it is the phase you are looking for. Would be better if I knew the application.
Nevertherless:
- 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.
- 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.
0 Comments
See Also
Categories
Find more on Switches and Breakers in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!