フィルターのクリア

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

32 ビュー (過去 30 日間)
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

採用された回答

Mathieu NOE
Mathieu NOE 2022 年 2 月 18 日
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 件のコメント
Susan
Susan 2022 年 2 月 18 日
Thank you so much for your help! I truly appreciate it.
Mathieu NOE
Mathieu NOE 2022 年 3 月 1 日
As always, my pleasure !

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeSpectral Measurements についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by