How to replicate bode plot

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

Answer 1

[Star Strider]

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)

AV = 

%first method: bode() function
[n,d]=numden(AV)

n = 

d = 

% n=fliplr(coeffs(collect(n,s),s))
% d=fliplr(coeffs(collect(d,s),s))
np = sym2poly(n)

np = 1×3

1.0e+36 * 0.2052 7.3283 0

dp = sym2poly(d)

dp = 1×4

1.0e+46 * 0.0000 0.0000 0.0033 3.2139

sys=tf(np, dp)

sys =
 
               2.052e35 s^2 + 7.328e36 s
  ---------------------------------------------------
  1.217e32 s^3 + 1.909e38 s^2 + 3.347e43 s + 3.214e46
 
Continuous-time transfer function.

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!

.

[Lucas]

What a wonderful morning when it works as you want !

Thanks a lot

[Star Strider]

As always, my pleasure!