Matlab - Bode plot of discrete and continuous Function

I am rather new to Matlab and I just cant make sense of what I see in the bode plot of the continuous and discrete version of the same function. The bode plot of the continuous function looks as expected. However the bode plot of the discrete version has a phase offset of +90 degrees and the gain stays the same at lower frequencies. It's basically a lag compensator with an integrator. This is not the final result I am going for, but the easiest example I could think of to make the problem as clear as possible. Please note that I am using a rather high sampling rate of ~12 MHz, which is what I need in the final hardware design. The high sampling rate seems to be part of the problem. With lower sampling rates the frequency at which the bode plots differ gets shifted to the left.

Can anyone explain this behaviour and maybe how I can achieve a result that is closer to the continuous bode plot using the high sampling rate?

Code of the example:

    FS = 48000*256;
    T_sample = 1/FS;
    
    accu = tf([0 1], [1 0]);
    lag = tf([10 1], [100 1]);
    
    lag_d = c2d(lag, T_sample, "zoh");
    accu_d = c2d(accu, T_sample, "zoh");
    
    bode(lag*accu, lag_d*accu_d);

[Benjamin Pommer]

Had you been finding an answer to that problem? I am struggeling with the same thing?

[Luca Ferro]

it doesn't make any sense to bode plot a discrete transfer function.

Bode is defined in the Laplace domain (s) while discrete transfer functions are in the z domain.

It's like putting a cat in a pool and expecting it to behave like a fish.

Answer 1

[Paul]

One problem is that the discretization of the product is not the product of the discretization. But with T_sample so small, that's probably not the issue. It appears that bode is having a problem with numerical calculations in the tf form. Instead, use zpk. zpk is preferred over tf in most (all?) cases.

FS = 48000*256;
T_sample = 1/FS;
    
accu = tf([0 1], [1 0]);
lag = tf([10 1], [100 1]);

htf = lag*accu;
hzpk = zpk(lag)*zpk(accu);
hzpkd = c2d(hzpk,T_sample,'zoh');
    
%lag_d = c2d(lag, T_sample, "zoh");
%accu_d = c2d(accu, T_sample, "zoh");
    
bode(htf, hzpk, hzpkd);

[Pascal Gahinet]

Paul is correct, this is a classic illustration of the pitfalls of the TF representation. Here lag_d and accu_d both have a pole close to z=1. When forming lag_d*accu_d, the denominator is close to (z-1)^2 and rounding errors on repeated eigenvalues move the two poles to a complex pair with natural frequency 5e-3:

>> damp(lag_d*accu_d)

Pole Magnitude Damping Frequency Time Constant

(rad/seconds) (seconds)

1.00e+00 1.00e+00 1.00e+00 5.00e-03 2.00e+02

1.00e+00 1.00e+00 1.00e+00 5.00e-03 2.00e+02

This is why the Bode plot is not what you'd expect.

The standard recommendation is to use the state-space representation (@ss) for all operations on LTI models. Here for example:

lag_d = c2d(ss(lag), T_sample, "zoh");

accu_d = c2d(ss(accu), T_sample, "zoh");

bode(lag*accu, lag_d*accu_d);

For more information see

https://www.mathworks.com/help/control/ug/using-the-right-model-representation.html

[Paul]

Hi Pascal,

The reason for the observed behavior is a bit more nuanced than I orginally appreciated

FS = 48000*256;
T_sample = 1/FS;
accu = tf([0 1], [1 0]);
lag = tf([10 1], [100 1]);
lag_d = c2d(lag, T_sample, "zoh");
accu_d = c2d(accu, T_sample, "zoh");

Keep lag_d and accu_d in tf form and compute the magnitude of their product at low frequency using freqresp, which yields the expected result.

w = 1e-4;
h1 = lag_d*accu_d;
db(abs(freqresp(h1,w)))

ans = 80.0000

If we compute the tf product and convert the result to zpk, we get the incorrect result

h2 = zpk(lag_d*accu_d);
db(abs(freqresp(h2,w)))

ans = 52.0377

which is the same result as from bode keeping everything in tf form on input

db(bode(lag_d*accu_d,w))

ans = 52.0377

Internal to bode the input, for this example, is converted to zpk and then the freqresp is computed from that (I checked).

Therefore, in this particular example, the problem isn't with tf form per se (calling freqresp on the tf form with a full frequency vector yields the correct result). Rather, the problem is with the little bit of inaccuracy in the conversion of the product from tf to zpk, which is probably due to the (nearly?) repeated root. I'm sure that internal conversion in bode to zpk is done for good reason in general, but I believe it's detrimental in this particular instance.

If we convert to zpk and then take the product, we get the correct result

h3 = zpk(lag_d)*zpk(accu_d);
db(abs(freqresp(h3,w)))

ans = 79.9996

The difference between h2 and h3 is

format long e
[h2.k h3.k]

ans = 1×2

1.0e+00 * 8.138020833333334e-09 8.138020833333334e-09

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[h2.z{:},h3.z{:}]

ans = 1×2

9.999999918619791e-01 9.999999918619791e-01

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[h2.p{:},h3.p{:}]

ans = 2×2

9.999999995930989e-01 9.999999991861979e-01 9.999999995930989e-01 1.000000000000000e+00

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

So very small differences in computed pole locations cause ~30 dB difference in gain at low frequency.

Of course, using ss from the outset is a good thing to do. In retrospect, there was time in the past, possibly the distant past, that I thought bode always converted the input to ss before computing the response.

And I believe starting from zpk in cases like this, where only products are involved, should be fine as well. Would be very interested to know if you disagree with this assertion.