How to manually replicate the bode() gain plot from a transfer function

28 views (last 30 days)
Hello,
I have been trying to replicate a bode gain plot from a given transfer function (see code below).
I am not seeing where I went wrong, but my manual plots looks like this:
But the actual bode plots should look like this:
I thought that I was setting up my amplitude equation exactly the same as how the bode() function operates, but apparently not. Can you see where I made the mistake?
Here is my code:
syms C R_1 R_2 R_3 L s V_o i omega
%Equations
R_1 = 10;
R_2 = 10;
R_3 = 10;
L = 0.001;
C = 2*10^-6;
Z_1 = R_1+L*s;
Z_2 = 1/(C*s)+R_2;
Z_3 = R_3;
Z_23 = 1/(1/(Z_2)+1/(Z_3));
Z_tot = Z_1+Z_23;
V_i = Z_1*i+V_o;
V_o = (V_i/(Z_1))/((1/((1/(C*s))+R_2))+1/R_3+1/(Z_1));
Transfer_func(s) = vpa(simplify(V_o/V_i),5);
amp = (sqrt((5.903*10^24)^2+(1.1806*10^20*omega)^2))/(sqrt((1.1806*10^25+(-2.3613*10^16*omega^2))^2+(9.4447*10^20*omega)^2))
amp_dB = 20*log(amp)
fplot(amp_dB,[0.001 1*10^6])
figure
subplot(2,1,1)
fplot(20*log(abs(Transfer_func)), [0.001 1*10^6])
legend('|H( j\omega )|')
title('Gain (dB)')
grid
subplot(2,1,2)
fplot(180/pi()*angle(Transfer_func), [0.001 1*10^6])
grid
title('Phase (degrees)')
xlabel('Frequency (rad/s)')
legend('\phi( j\omega )')
H = tf([118059162071741125000 5902958103587056517120000],[23611832414348225 944473296573929026712 11805916207174113034240000]);
figure
bode(H,{0.001,1*10^6})

Accepted Answer

Star Strider
Star Strider on 4 Jun 2020
But the actual bode plots should look like this:’
It does.
Harkening back to yesterday:
syms C R_1 R_2 R_3 L s V_o i omega
%Equations
R_1 = 10;
R_2 = 10;
R_3 = 10;
L = 0.001;
C = 2*10^-6;
Z_1 = R_1+L*s;
Z_2 = 1/(C*s)+R_2;
Z_3 = R_3;
Z_23 = 1/(1/(Z_2)+1/(Z_3));
Z_tot = Z_1+Z_23;
V_i = Z_1*i+V_o;
V_o = (V_i/(Z_1))/((1/((1/(C*s))+R_2))+1/R_3+1/(Z_1));
Transfer_func(s) = vpa(simplify(V_o/V_i), 5);
Transfer_func(omega) = subs(Transfer_func, {s},{1j*omega});
figure
subplot(2,1,1)
fplot(20*log10(abs(Transfer_func)), [0.001 1E6])
set(gca, 'XScale','log')
legend('|H( j\omega )|')
title('Gain (dB)')
grid
subplot(2,1,2)
fplot(180/pi*angle(Transfer_func), [0.001 1E6])
set(gca, 'XScale','log')
grid
title('Phase (degrees)')
xlabel('Frequency (rad/s)')
legend('\phi( j\omega )')
produces:
Note the slight changes between your code and mine.
  7 Comments
Star Strider
Star Strider on 9 Jun 2020
I suspect Shane Palmer is encountering this for the first time and is experimenting in order to understand how it all works. (I did something similar when I first encountered them, so I understand.)
Shane Palmer
Shane Palmer on 9 Jun 2020
Star Strider is correct, I am experimenting and learning how to work with matlab and transfer functions simultaneously. I originally started the code symbolically, and then went back in afterwards to add the input values. I don't have a solid answer for you on whether the H is 100% correct, but the manual plot of the bode charts match the matlab bode() command for the charts, so I would say yes as far as the process goes.

Sign in to comment.

More Answers (0)

Categories

Find more on Get Started with Control 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!