How to plot the amplitude ratio and phase differences for two time series signals ?

Greetings,

I'd like to plot the amplitude ratio and phase differences of two signals. I used the following code to calculate and plot the amplitude ratio and phase differences by modifying the code provided here:

https://www.mathworks.com/matlabcentral/answers/91647-how-do-i-calculate-the-amplitude-ratio-and-phase-lag-for-two-sinusoidal-signals-in-matlab

clear all
% set the lag
lag = pi/2;
% amplitude for sinusodial 1
amp_x = 2;
% amplitude for sinusodial 2
amp_y = 4;
% Sampling frequency
Fs = 1024;
% Time vector of 0.5 second
t = 0:1/Fs:0.5*(Fs-1)/Fs;
% number of points
npts = length(t);
% Create a sine wave of 20 Hz.
x = amp_x * (sin(2*pi*t*20) + ...
0.25*randn(1,npts)) + 5;
% Create a sine wave of 20 Hz
% with a phase of pi/2.
y = amp_y * (sin(2*pi*t*20 + lag)...
+ 0.25*randn(1,npts));
% remove bias
x = x - mean(x);
y = y - mean(y);
% plot the signal
figure(1)
plot(t,x,t,y);
xlabel('Time (s)');
ylabel('Amplitude');
legend('Signal x(t):sin(2*pi*t*20)',...
'Signal y(t):sin(2*pi*t*20+lag)');
% take the FFT
X=fft(x);
Y=fft(y);
% Calculate the numberof unique points
NumUniquePts = ceil((npts+1)/2);
figure(2)
subplot(211);
f = (0:NumUniquePts-1)*Fs/npts;
plot(f,abs(X(1:NumUniquePts)));
title('X(f) : Magnitude response');
ylabel('|X(f)|')
subplot(212)
plot(f,abs(Y(1:NumUniquePts)));
title('Y(f) : Magnitude response')
xlabel('Frequency (Hz)');
ylabel('|Y(f)|')
figure(3)
subplot(211)
plot(f,angle(X(1:NumUniquePts)));
title('X(f) : Phase response');
ylabel('Phase (rad)');
subplot(212)
plot(f,angle(Y(1:NumUniquePts)));
title('Y(f) : Phase response');
xlabel('Frequency (Hz)');
ylabel('Phase (rad)');


% determine the phase difference
px = angle(X);
py = angle(Y);
phase_lag = py - px
% determine the amplitude scaling
amplitude_ratio = abs(Y)/abs(X)

Now, I'd like to plot the phase difference along with the phase of both X and Y and see how it looks different.

figure
polarplot(angle(X),abs(X),'--')
hold on
polarplot(angle(Y),abs(Y),'--')
hold on
polarplot(angle(phase_lag),rho,'--')

What would be the "rho" at the last line? is it the amplitude_ratio? Is there any better ways to see these phase differences/lags? I'd like to plot the amplitude ratio vs time too. How should I do that?

TIA

[Susan]

Any suggestions?

Answer 1

[Mathieu NOE]

hello Susan

I modified your code : take the amplitude ratio and phase lag at the X and Y fft frequency indexes corresponding to the max amplitude of the fft data

new code :

clear all
% set the lag
lag = pi/2;
% amplitude for sinusodial 1
amp_x = 2;
% amplitude for sinusodial 2
amp_y = 4;
% Sampling frequency
Fs = 1024;
% Time vector of 0.5 second
t = 0:1/Fs:0.5*(Fs-1)/Fs;
% number of points
npts = length(t);
% Create a sine wave of 20 Hz.
x = amp_x * (sin(2*pi*t*20) + ...
0.25*randn(1,npts)) + 5;
% Create a sine wave of 20 Hz
% with a phase of pi/2.
y = amp_y * (sin(2*pi*t*20 + lag)...
+ 0.25*randn(1,npts));
% remove bias
x = x - mean(x);
y = y - mean(y);
% plot the signal
figure(1)
plot(t,x,t,y);
xlabel('Time (s)');
ylabel('Amplitude');
legend('Signal x(t):sin(2*pi*t*20)',...
'Signal y(t):sin(2*pi*t*20+lag)');
% take the FFT
X=fft(x);
Y=fft(y);
% Calculate the numberof unique points
NumUniquePts = ceil((npts+1)/2);
figure(2)
subplot(211);
f = (0:NumUniquePts-1)*Fs/npts;
X=X(1:NumUniquePts);
Y=Y(1:NumUniquePts);
% find peak fft amplitude points
[valx,indx] = max(abs(X));
[valy,indy] = max(abs(Y));
% NB : the indx and indy indexes should be the same otherwise the signals have
% different base frequencies , which leads to non constant phase lag
if indx~=indy
    error('indx and indy indexes should be the same')
end

plot(f,abs(X));
title('X(f) : Magnitude response');
ylabel('|X(f)|')
subplot(212)
plot(f,abs(Y));
title('Y(f) : Magnitude response')
xlabel('Frequency (Hz)');
ylabel('|Y(f)|')

% figure(3)
% subplot(211)
% plot(f,angle(X));
% title('X(f) : Phase response');
% ylabel('Phase (rad)');
% subplot(212)
% plot(f,angle(Y));
% title('Y(f) : Phase response');
% xlabel('Frequency (Hz)');
% ylabel('Phase (rad)');

% determine the phase difference
px = angle(X(indx));
py = angle(Y(indy));
phase_lag = py - px
% determine the amplitude scaling
amplitude_ratio = abs(Y(indy))/abs(X(indx))

[Susan]

Thank you so much for your response. I'd like to plot sth like this fig. Do you know how I can do that? This fig shows the comparison between phase of one signal (blue crosses) and the phase differences from this signal with another signal (red dots). I've attached both signals. I've also like to plot the amplitude ratio vs time .Any idea would be greatly appreciated.

[Mathieu NOE]

hello again

sorry I don't know how I can use the provided mat files as t, x, and y are missing

so using the same file as before , I did some mods to show you the polarplot and also the amplitude ratio vs time (smoothing the abs of x and y and then making the ratio of the two. Notice that if you do a sample by sample division you will get a very unstationnary amplitude ratio. Both traces are shown on figure 4 for info

hope it helps

NB : the polar plot need the radius info , which I don't know where to pick from as the signal has only one dominant peak frequency + random noise, if I were to use the abs(X) (fft amplitude) as the radius , then we would see one single blue spot with a large radius, and the rest very close to the origin point (0,0), which is not like the picture you mentionned. As a laternative I plotted all phase points on a radius = 1 value.

clear all
% set the lag
lag = pi/2;
% amplitude for sinusodial 1
amp_x = 2;
% amplitude for sinusodial 2
amp_y = 4;
% Sampling frequency
Fs = 1024;
% Time vector of 0.5 second
t = 0:1/Fs:0.5*(Fs-1)/Fs;
% number of points
npts = length(t);
% Create a sine wave of 20 Hz.
x = amp_x * (sin(2*pi*t*20) + ...
0.25*randn(1,npts)) + 5;
% Create a sine wave of 20 Hz
% with a phase of pi/2.
y = amp_y * (sin(2*pi*t*20 + lag)...
+ 0.25*randn(1,npts));

% remove bias
x = x - mean(x);
y = y - mean(y);
% plot the signal
figure(1)
plot(t,x,t,y);
xlabel('Time (s)');
ylabel('Amplitude');
legend('Signal x(t):sin(2*pi*t*20)',...
'Signal y(t):sin(2*pi*t*20+lag)');
% take the FFT
X=fft(x);
Y=fft(y);
% Calculate the numberof unique points
NumUniquePts = ceil((npts+1)/2);
figure(2)
subplot(211);
f = (0:NumUniquePts-1)*Fs/npts;
X=X(1:NumUniquePts);
Y=Y(1:NumUniquePts);
% find peak fft amplitude points
[valx,indx] = max(abs(X));
[valy,indy] = max(abs(Y));
% NB : the indx and indy indexes should be the same otherwise the signals have
% different base frequencies , which leads to non constant phase lag
if indx~=indy
    error('indx and indy indexes should be the same')
end

plot(f,abs(X));
title('X(f) : Magnitude response');
ylabel('|X(f)|')
subplot(212)
plot(f,abs(Y));
title('Y(f) : Magnitude response')
xlabel('Frequency (Hz)');
ylabel('|Y(f)|')

% determine the phase difference
px = angle(X(indx)); % in radians
py = angle(Y(indy)); % in radians
phase_lag = py - px   % in radians

figure(3)
% rho = smoothdata(abs(X),'gaussian',10); % not very meaningfull ? 
rho = ones(size(X));
polarplot(angle(X),rho,'b*',phase_lag,1,'rd') 
title('X(f) : Phase plot');


% determine the amplitude scaling
amplitude_ratio = abs(y)./(abs(x)+eps); % computed sample by sample
N = 150;
abs_x_s = smoothdata(abs(x),'gaussian',N);
abs_y_s = smoothdata(abs(y),'gaussian',N);
amplitude_ratio_smooth = abs_y_s./abs_x_s;
figure(4)
semilogy(t,amplitude_ratio,t,amplitude_ratio_smooth);
title('amplitude ratio');

[Susan]

Thank you so much for your help! I truly appreciate it.

[Mathieu NOE]

As always, my pleasure !