How to replicate bode plot
4 views (last 30 days)
Show older comments
Hi everyone,
I am trying to plot the magnitude and phase of a transfer function. For that I used the bode() function, but the result seems a bit weird to me. Therefore I thought I could plot the bode myself.
I know that a very similar question has already been asked: How to manually replicate the bode() gain plot from a transfer function. Unfortunately, even with the answers I do not manage to find the same result between: the bode() function and my "handmade script".
The bode() function result:
And my script result:
I think the problem comes from the first method (bode() function) with probably wrong arguments. As the second method (handmade) is directly inspired from the accepted answer of the question How to manually replicate the bode() gain plot from a transfer function.
Here is the script:
clc
clear
%variables
syms s ESR C LP LS k RS RP RB RDSON CD AV omega f
%Numerical values
LP=0.000083;
LS=0.000083;
k=0.999;
RS=0.07;
RP=0.07;
RB=0.01;
RDSON=0.11;
CD=0.000022;
C=100;
ESR=0.00028;
%Equations
ZL=ESR+(1/(s*C));
ZM=s*k*LP;
ZP=RB+RP+s*(1-k)*LP;
ZS=RDSON+(1/(s*CD))+RS+s*(1-k)*LS;
AV=(ZL*ZM)/((ZS+ZL)*(ZP+ZM)+ZP*ZM);
%first method: bode() function
[n,d]=numden(AV);
n=fliplr(coeffs(collect(n,s),s));
d=fliplr(coeffs(collect(d,s),s));
sys=tf(double(n),double(d));
bode(sys,{0.01 1E7})
grid
%second method: handmade
Transfer_func(f) = subs(AV,{s},{1j*2*pi*f});
figure
subplot(2,1,1)
fplot(20*log10(abs(Transfer_func)), [0.01 1E7])
set(gca, 'XScale','log')
legend('|H( j\omega )|')
title('Gain (dB)')
grid
subplot(2,1,2)
fplot(180/pi*angle(Transfer_func), [0.01 1E7])
set(gca, 'XScale','log')
grid
title('Phase (degrees)')
xlabel('Frequency (Hz)')
legend('\phi( j\omega )')
PS: I am quite a newbie to matlab, sorry if I made some huge mistakes
0 Comments
Accepted Answer
Star Strider
on 9 Nov 2022
The problem is that you wer not sending the same vectors to your tf call that you were using in the symbolic solution. Note that here I used the sym2poly function on ‘n’ and ‘d’ and then sent those results to tf.
With that change, the plots match, although the frequency scale of one plot is in rad/sec and the other is in Hz. That can be fixed with a bodeoptions call, and then the units match —
% clc
% clear
%variables
syms s ESR C LP LS k RS RP RB RDSON CD AV omega f
%Numerical values
LP=0.000083;
LS=0.000083;
k=0.999;
RS=0.07;
RP=0.07;
RB=0.01;
RDSON=0.11;
CD=0.000022;
C=100;
ESR=0.00028;
%Equations
ZL=ESR+(1/(s*C));
ZM=s*k*LP;
ZP=RB+RP+s*(1-k)*LP;
ZS=RDSON+(1/(s*CD))+RS+s*(1-k)*LS;
AV=(ZL*ZM)/((ZS+ZL)*(ZP+ZM)+ZP*ZM);
AV = vpa(simplify(AV, 500), 5)
%first method: bode() function
[n,d]=numden(AV)
% n=fliplr(coeffs(collect(n,s),s))
% d=fliplr(coeffs(collect(d,s),s))
np = sym2poly(n)
dp = sym2poly(d)
sys=tf(np, dp)
opts = bodeoptions;
opts.FreqUnits = 'Hz';
bode(sys,{0.01 1E7}, opts)
grid
%second method: handmade
Transfer_func(f) = subs(AV,{s},{1j*2*pi*f});
figure
subplot(2,1,1)
fplot(20*log10(abs(Transfer_func)), [0.01 1E7])
set(gca, 'XScale','log')
legend('|H( j\omega )|')
title('Gain (dB)')
grid
subplot(2,1,2)
fplot(180/pi*angle(Transfer_func), [0.01 1E7])
set(gca, 'XScale','log')
grid
title('Phase (degrees)')
xlabel('Frequency (Hz)')
legend('\phi( j\omega )')
Thank you for referencing my previous post!
.
2 Comments
More Answers (0)
See Also
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!