Phase lag between two signals

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 个评论
显示 -1更早的评论隐藏 -1更早的评论

Daniel M 2019-11-28

编辑：Daniel M 2019-11-28

在 MATLAB Online 中打开

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 2019-12-2

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/493692-phase-lag-between-two-signals#answer_404244

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Phase lag between two signals

3 个评论
显示 1更早的评论隐藏 1更早的评论

采纳的回答

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

更多回答（1 个）

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

Phase lag between two signals

3 个评论 显示 1更早的评论隐藏 1更早的评论

采纳的回答

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

更多回答（1 个）

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

3 个评论
显示 1更早的评论隐藏 1更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论