c2d Tustin not matching analytical method

I have a second order lowpass filter. Discrete transter function using c2d with tustin method does not match the analytical solution (see the code below with bode comparison). Analytical solution was obtained by subing s = 2/Ts * (z-1)/(z+1) in continuous transfer function which is apparantly consistent with Matlab c2d tsting option with the difference that the Matlab code is first converting TF to PZK and then fomulating c2d (see reference folders below) which is not an exact match to my fomulation but I did not expect this big of a difference. Thoughts appreciated.

/Applications/MATLAB_R2024b.app/toolbox/shared/controllib/engine/+ltipack/@tfdata/c2d.m   % ltipack.tfdata method
/Applications/MATLAB_R2024b.app/toolbox/shared/controllib/engine/+ltipack/@zpkdata/c2d.m  % ltipack.zpkdata method

%% Filter Specs
Ts = 0.01;

zeta   = 1;
wn_rps = 30 *2*pi;

%% C2D (tustin)
s = tf('s');

tf_cont = wn_rps^2 / (s^2 + 2*zeta*wn_rps*s + wn_rps^2);

tf_disc = c2d(tf_cont, Ts, 'Method','tustin');

%% Analytical (tustin)
% s = 2/Ts * (z-1)/(z+1)
K = 2/Ts;

num = wn_rps^2 * [1 2 1];
den = [   K^2 + 2*zeta*wn_rps*K +   wn_rps^2, ...
       -2*K^2                   + 2*wn_rps^2, ...
          K^2 - 2*zeta*wn_rps*K +   wn_rps^2 ];

tf_disc_analytical = tf(num, den, Ts);

%% Plot
figure, bode(tf_disc), hold on, bode(tf_disc_analytical) 

Answer 1

[Paul]

Hi Alborz,

%% Filter Specs
Ts = 0.01;
zeta   = 1;
wn_rps = 30 *2*pi;
%% C2D (tustin)
s = tf('s');
tf_cont = wn_rps^2 / (s^2 + 2*zeta*wn_rps*s + wn_rps^2);

This call signature is not shown on the doc page c2d. There should only be 2 or 3 input arguments.

tf_disc = c2d(tf_cont, Ts, 'Method','tustin')

tf_disc =
 
     0.3597 z + 0.3597
  ------------------------
  z^2 - 0.3037 z + 0.02305
 
Sample time: 0.01 seconds
Discrete-time transfer function.

A little-known "feature" of c2d is that it only looks at the first character of the 'method' argument, which is the third argument, and matches that to the first character of one of the allowable methods. In this case, the 'M' corresponds to 'matched' (and it appears that the fourth argument, in this case 'tustin', is ignored.

tf_disc = c2d(tf_cont, Ts, 'matched')

tf_disc =
 
     0.3597 z + 0.3597
  ------------------------
  z^2 - 0.3037 z + 0.02305
 
Sample time: 0.01 seconds
Discrete-time transfer function.

The correct conversion for Tustin would be

tf_disc = c2d(tf_cont, Ts, 'tustin')

tf_disc =
 
  0.2354 z^2 + 0.4708 z + 0.2354
  ------------------------------
   z^2 - 0.05923 z + 0.0008769
 
Sample time: 0.01 seconds
Discrete-time transfer function.

%% Analytical (tustin)
% s = 2/Ts * (z-1)/(z+1)
K = 2/Ts;
num = wn_rps^2 * [1 2 1];
den = [   K^2 + 2*zeta*wn_rps*K +   wn_rps^2, ...
       -2*K^2                   + 2*wn_rps^2, ...
          K^2 - 2*zeta*wn_rps*K +   wn_rps^2 ];
tf_disc_analytical = tf(num, den, Ts)

tf_disc_analytical =
 
  3.553e04 z^2 + 7.106e04 z + 3.553e04
  ------------------------------------
     1.509e05 z^2 - 8939 z + 132.4
 
Sample time: 0.01 seconds
Discrete-time transfer function.

minreal normalizes the leading coefficient of the denominator to unity to match the computed result above.

tf_disc_analytical = minreal(tf_disc_analytical)

tf_disc_analytical =
 
  0.2354 z^2 + 0.4708 z + 0.2354
  ------------------------------
   z^2 - 0.05923 z + 0.0008769
 
Sample time: 0.01 seconds
Discrete-time transfer function.

[Al]

Thanks Paul. Great catch. It seems there is an oppotunity to make c2d code more robust. I should have gotten an error in this case.