Calculate phase shift i bode plot for 2 order system data
11 views (last 30 days)
Show older comments
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)
%load files
file=(['data_' num2str(R(ii)) 'ohm.txt']);
data = readtable(file, 'HeaderLines',1);
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);
%calculates phase radians
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
legend(Legend, 'Interpreter', 'Latex','Location', 'best')
hold off
%phase shift line
yline(-90,'r-', 'LabelHorizontalAlignment','left', 'Linewidth', 1.5)
0 Comments
Accepted Answer
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);
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.
9 Comments
William Rose
on 5 Mar 2023
@Olga Rakvag, you are welcome. I don;t understand why Zumbahlen gave that strange formula for phase. The derivation of the formulas for K, , and ζ, for your filter circuit, is below.
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!