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

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

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

db(bode(lag_d*accu_d,w))

So 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 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 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)))

The difference between h2 and h3 is

format long e
[h2.k h3.k]
[h2.z{:},h3.z{:}]
[h2.p{:},h3.p{:}]

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

Of course, using ss from the outset is a good thing to do.

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.