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

Question

Susan 2022-2-17

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1652400-how-to-plot-the-amplitude-ratio-and-phase-differences-for-two-time-series-signals

评论： Mathieu NOE 2022-3-1

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

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

Susan 2022-2-18

Any suggestions?

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Mathieu NOE 2022-2-18

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1652400-how-to-plot-the-amplitude-ratio-and-phase-differences-for-two-time-series-signals#answer_899260

在 MATLAB Online 中打开

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))

4 个评论
显示 2更早的评论隐藏 2更早的评论

Mathieu NOE 2022-2-18

在 MATLAB Online 中打开

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');