# Calculate phase shift i bode plot for 2 order system data

18 views (last 30 days)
Olga Rakvag on 2 Mar 2023
Edited: Olga Rakvag on 5 Mar 2023
Hi, I could not get right alignment for the data, but know that the experimental data should be almost the same as simulated. Can't figure out my mistake. I used formula 3 from https://www.analog.com/en/analog-dialogue/articles/phase-response-in-active-filters-2.html Is it a way to use Vout Vin instead? I know that magnitude can be computed from dB=20 log (V_out/V_in). Added formula for magnitude as function of voltage ratio, V_out/V_in. The system looks like
clc;clear all; close all;
%% variables
t = logspace(0, 0.1, 1000);
R = [1 100 560];
R_L = 1E+5;
c = 1E-6;
L = 0.1;
R_S = 65;
Tab=table();
for ii = 1:numel(R)
file=(['data_' num2str(R(ii)) 'ohm.txt']);
data.Properties.VariableNames = {['freq_' num2str(R(ii)) ' '], ['dB_' num2str(R(ii)) ' '], ['Vout_' num2str(R(ii)) ' '], ['Vin_' num2str(R(ii)) ' ']}; % Names of colum
f1=data.(['freq_' num2str(R(ii)) ' ']*2*pi);
f=2*pi*f1;
dB=data.(['dB_' num2str(R(ii)) ' ']);
Vout=data.(['Vout_' num2str(R(ii)) ' ']);
Vin=data.(['Vin_' num2str(R(ii)) ' ']);
A = [-1/(R_L*c) 1/c;-1/L -(R(ii)+R_S)/L];
B = [0; 1/L];
C = [1 0];
D = 0;
w_n=sqrt((R(ii)+R_S+R_L)/(c*L*R_L))
K=R_L/(R(ii)+R_L+R_S)
ksi=((c*R_L*(R(ii)+R_S)+L)/(R(ii)+R_S+R_L))*w_n/2
t_p=pi/(w_n*(sqrt(1-ksi^2)))*1000% *1000 for [ms]
OS=exp(-(ksi*pi/(sqrt(1-ksi^2))))*100
Tab(ii,:)=table(R(ii), K, w_n, ksi, t_p, OS);
num=K
den= [1 (1/w_n^2) (2*ksi/w_n)]
hp=tf(num, den)
model=ss(A,B,C, D);
theta=-atan(1/ksi*(2*((f)/(w_n))+sqrt(4-ksi^2)))-atan(1/ksi*(2*((f)/(w_n))-sqrt(4-ksi^2)));
thdeg=(theta*180)/pi; %i degrees
%teoretical vs experimental datas
figure (1)
disp(Tab)
h = bodeplot(model); hold on
setoptions(h,'MagVisible','off'); hold on
semilogx(f,thdeg,'--','Linewidth', 2) ;hold on, grid on
Legend{ii}=strvcat(['Simulated \$R =' num2str(R(ii)) '\Omega\$ '],[ 'Experimental \$R =' num2str(R(ii)) '\Omega\$']);
hold on
end
f = 300×1
5.0000 17.5890 31.9860 47.1530 62.8450 78.9380 95.3580 112.0550 128.9940 146.1450
Vout = 300×1
0.3060 0.3340 0.3330 0.3320 0.3330 0.3330 0.3330 0.3330 0.3330 0.3330
Vin = 300×1
0.3040 0.3340 0.3340 0.3340 0.3370 0.3400 0.3430 0.3470 0.3520 0.3570
w_n = 3.1633e+03
K = 0.9993
ksi = 0.1059
t_p = 0.9987
OS = 71.5638
num = 0.9993
den = 1×3
1.0000 0.0000 0.0001
hp = 0.9993 ----------------------------- s^2 + 9.993e-08 s + 6.696e-05 Continuous-time transfer function.
Var1 K w_n ksi t_p OS ____ _______ ______ ______ _______ ______ 1 0.99934 3163.3 0.1059 0.99875 71.564
f = 300×1
5.0000 17.5890 31.9860 47.1530 62.8450 78.9380 95.3580 112.0550 128.9940 146.1450
Vout = 300×1
0.3070 0.3330 0.3330 0.3350 0.3340 0.3320 0.3330 0.3320 0.3320 0.3330
Vin = 300×1
0.3060 0.3330 0.3330 0.3370 0.3370 0.3380 0.3410 0.3440 0.3470 0.3530
w_n = 3.1649e+03
K = 0.9984
ksi = 0.2623
t_p = 1.0286
OS = 42.5805
num = 0.9984
den = 1×3
1.0000 0.0000 0.0002
hp = 0.9984 ----------------------------- s^2 + 9.984e-08 s + 0.0001657 Continuous-time transfer function.
Var1 K w_n ksi t_p OS ____ _______ ______ _______ _______ ______ 1 0.99934 3163.3 0.1059 0.99875 71.564 100 0.99835 3164.9 0.26225 1.0286 42.58
f = 300×1
5.0000 17.5890 31.9860 47.1530 62.8450 78.9380 95.3580 112.0550 128.9940 146.1450
Vout = 300×1
0.3050 0.3330 0.3330 0.3340 0.3330 0.3330 0.3320 0.3310 0.3300 0.3290
Vin = 300×1
0.3030 0.3300 0.3290 0.3280 0.3250 0.3220 0.3160 0.3110 0.3050 0.2980
w_n = 3.1721e+03
K = 0.9938
ksi = 0.9867
t_p = 6.0959
OS = 5.1716e-07
num = 0.9938
den = 1×3
1.0000 0.0000 0.0006
hp = 0.9938 ----------------------------- s^2 + 9.938e-08 s + 0.0006221 Continuous-time transfer function.
Var1 K w_n ksi t_p OS ____ _______ ______ _______ _______ __________ 1 0.99934 3163.3 0.1059 0.99875 71.564 100 0.99835 3164.9 0.26225 1.0286 42.58 560 0.99379 3172.1 0.98671 6.0959 5.1716e-07
legend(Legend, 'Interpreter', 'Latex','Location', 'best')
hold off
%phase shift line
yline(-90,'r-', 'LabelHorizontalAlignment','left', 'Linewidth', 1.5)

William Rose on 2 Mar 2023
Edited: William Rose on 2 Mar 2023
[edit: In case it's not obvious: Add the 2*pi factor in both places where f/w_n appears, in your formula for theta.]
I think you are asking why the solid curves in your plot are significantly different from the corresponding dashed lines in the same plot. The legend in the plot is mis-sized, so it is not possible to determine exactly what curve is what, from the plot as it appears in your post.
I think you need to add a factor of 2*pi in your equation for theta (dashed line plot), as explained below:
The dashed lines in the plot are generated by
semilogx(f,thdeg,'--','Linewidth', 2);
Unrecognized function or variable 'f'.
where the phase angle estimate thetadeg=(180/pi)*theta, where theta is given by
theta=-atan(1/ksi*(2*((f)/(w_n))+sqrt(4-ksi^2)))-atan(1/ksi*(2*((f)/(w_n))-sqrt(4-ksi^2)));
You say you were using equation 3, which is
Notice how equation 3 has , but you have f/w_n. Replace f with 2*pi*f.
Your ksi corresponds to α in equation 3:
ksi=((c*R_L*(R(ii)+R_S)+L)/(R(ii)+R_S+R_L))*w_n/2
I do not have the necessary information to check whether this formula for ksi is correct.
Your w_n corresponds to in equation 3.
w_n=sqrt((R(ii)+R_S+R_L)/(c*L*R_L))
I do not have the necessary information to check whether this formula for w_n is correct.
The solid curves in the plot are generated by .
h = bodeplot(model);
where model is a state space model defined by
model=ss(A,B,C, D);
where A,B,C,D are given by
A = [-1/(R_L*c) 1/c;-1/L -(R(ii)+R_S)/L];
B = [0; 1/L];
C = [1 0]; D = 0;
I do not know if the equations you have provided for A,B,C,D are a correct interpretation of the second order low pass active filter which is described by equation 3 in the paper you cited.
Olga Rakvag on 5 Mar 2023
Edited: Olga Rakvag on 5 Mar 2023
Zumbahlen for sure has written the right formula. Because his fomula for 1 order system worked for me. You see, it is not only because of his formel, but I calculated ksi wrong, the one you caled for zeta! But thank you a lot, for helping me to see where I was mistaken and find the solution. ;-)

### Categories

Find more on Get Started with Phased Array System Toolbox in Help Center and File Exchange

### Community Treasure Hunt

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

Start Hunting!

Translated by