Phase lag between two signals
216 views (last 30 days)
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 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.
[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
[r,lags] = xcorr(x1,x2); % get the 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');
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
% 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
% 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')
if nargin < 3 || isempty(begtime) || begtime < 0
begtime = 0;
if nargin < 4 || isempty(endtime) || endtime > (length(x)-1)*T
endtime = (length(x)-1)*T;
if begtime >= endtime
error('begtime must be before endtime')
if nargin < 5 || isempty(demean)
demean = 0;
elseif ~(demean==1 || demean==0)
error('demean must be true or false')
if nargin < 6 || isempty(plotType)
plotType = 'plot';
% Get indices from beg/end times
begind = floor(begtime*fs+1);
endind = floor(endtime*fs+1);
X = x(begind:endind); % new signal
X = detrend(X);
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;
% 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;
% H = tightfig(H);
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.
- 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.