137 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

Daniel M
on 28 Nov 2019

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

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

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.

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!