Calculate phase shift i bode plot for 2 order system data

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

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)

Answer 1

[William Rose]

@Olga Rakvag,

[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.

[William Rose]

@Olga Rakvag,

You asked if there is a way to use Vin, Vout data from the files for the phase analysis.

No, there is not. Vin and Vout in each of the three text files are the input and output voltage amplitudes, at each frequency. This data does not contain information about the phase difference.

[William Rose]

@Olga Rakvag, Matlab Answers says you edited your question, but I don;t know what you changed. It is helpful, when you edit a previous post, to specify what you have edited, at the top.

[William Rose]

@Olga Rakvag, the data in Vin, Vout columns in the three text files are consistent with the data in the dB column of the respective file. The dB value is 20*log10(Vout/Vin), up to the precision of the Vin. Vout data.

If adding the factor of 2*pi does not solve the discrepancy between the phase response from bodeplot(model) and the phase response obtained from the formula for theta, then it is likely that there is an inconsitency between the formulas for A,B,C,D (used by bodeplot(model)) and the formulas for ksi and w_n (used by the formula for theta). The document you shared from Analog Devices does not include any of those formulas. I would review the sources for the formulas to make sure they are correct and consistent.

By the way, your code includes the line below:

den= [1 (1/w_n^2) (2*ksi/w_n)]

which looks wrong. I think it should be

den= [1  (2*ksi/w_n) (1/w_n^2)]

den is not used to compute the plotted phase angle, as far as I can tell. Therefore this mistake, if it is a mistake, is not the cause of the phase discrepancy.

[William Rose]

@Olga Rakvag,

[edit: Change " is the natural frequency" to " is the natural frequency"]

The PDF you cited (Zumbahlen2) has an error. Zumbahlen2 says, regarding his equation for phase, "In Equation 3, α, the damping ratio of the filter, is the inverse of Q (that is, Q = 1/α)." This is wrong. The damping ratio is 1/(2Q). (See here and here, for example.) Zumbahlen2 does not provide a formula for alpha, so it is not clear exactly what he meant by alpha. He certainly did not mean the decay rate alpha=1/tau, used by some authors. The decay rate alpha has units of 1/time, and Zumbahlen's alpha, in equation 3, must be dimensionless.

Another thing: Zumbahlen2 says the phase can be "approximated by" equation 3:

Why use an approximation, when the exact equation is simpler? The exact equation is

where ζ is the (true) damping ratio, and is the natural frequency. Maybe the discrepancy between phase curves which you observed is due to the fact that equation 3 is just an approximation. Maybe the phase curve pairs will match, if you use the exact equation for phase. To use the exact equation above, you will have to determine the equation for ζ for your filter, in terms of the filter component values.

[Olga Rakvag]

See your point, it is actually what I wanted to implement at the beginning! Good idea! The right equation! But it somehow did not work out with exact equation * i replaced theta in for loop with your propoced theta=-atan((ksi*(f/w_n).^2)/(1-(f/w_n).^2)), got 3 straight lines to my own dissapointment. I used ksi from for loop ksi=((c*R_L*(R(ii)+R_S)+L)/(R(ii)+R_S+R_L))*w_n/2

[William Rose]

@Olga Rakvag,

The formula you used is not the one I proposed, which is

You have (f/w_n).^2 in the numerator, but it should be f/w_n. Also, atan(y/x) only returns values between -pi/2 and +pi/2. Therefore, to get angles such as -pi, use atan2(y,x):

theta=-atan2(2*zeta*f/w_n,1-f.^2/w_n^2);

where zeta = damping ratio:

zeta=0.5*(R(ii)+R_S)*R_L*c*w_n/(R(ii)+R_S+R_L); % damping ratio

I computed K=gain, w_n=natural frequency, and zeta=damping ratio, for the circuit in the diagram you supplied. My formulas for K and w_n are the same as the forulas in your code. However, my zeta is slightly different from your ksi. Your ksi has an additional "L" term in the numerator. Using my formula, I get plots that match well:

Good luck.

[Olga Rakvag]

Awesome! I couldn't believe my eyes when I saw that your formula worked! Thank's a lot! :-)

[William Rose]

@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.

[Olga Rakvag]

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