How to plot the amplitude ratio and phase differences for two time series signals ?
23 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Susan
am 17 Feb. 2022
Kommentiert: Mathieu NOE
am 1 Mär. 2022
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
Akzeptierte Antwort
Mathieu NOE
am 18 Feb. 2022
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 Kommentare
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Spectral Measurements finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!