I'm trying to understand the derivation of the bilinear transform for a set of continuous time state-space matrices. I've found plenty of websites which list steps to perform the conversion (here 1 or here 2 or here 3) - but haven't found any derivation or insight into how they came about. When I attempt to do it myself, I keep getting terms in z which I can't seem to get rid of.

To make the problem clear, if I define the bilinear transform as:

$$ s=\frac{\alpha\left(z-1\right)}{z+1} $$

Where alpha is usually given as 2/T.

The transfer function of a continuous-time state-space system can be given as:

$$ H(s) = \mathbf{C}(s\mathbf{I}-\mathbf{A})^{-1}\mathbf{B} + \mathbf{D} $$

How do I derive the expressions for the matrices of the discrete time version bilinear transformed model:

$$ H(z) = \mathbf{C_d}(z\mathbf{I}-\mathbf{A_d})^{-1}\mathbf{B_d} + \mathbf{D_d} $$

When I make the variable change, I keep getting extra terms in z which I don't know what to do with. I'm kinda hoping this is just me failing to have a strong grasp of matrix operations. Any insight (including if I am going about this completely the wrong way) would be really appreciated!

--- Edit with answer.

Thanks to Klaz for posting an answer to this first. I managed to work through the problem today and thought I would post my own alternate working here.

Given the typical state-space equations:

$$ s\mathbf{Q}(s) = \mathbf{A}\mathbf{Q}(s) + \mathbf{B}\mathbf{X}(s) $$ $$ \mathbf{Y}(s) = \mathbf{C}\mathbf{Q}(s) + \mathbf{D}\mathbf{X}(s) $$

Focus on the state update equation first and make the bilinear substitution:

$$ \frac{\alpha\left(z-1\right)}{z+1} \mathbf{Q}(z) = \mathbf{A}\mathbf{Q}(z) + \mathbf{B}\mathbf{X}(z) $$

We manipulate this expression to make the LHS $z\mathbf{Q}(z)$:

$$ \alpha\left(z-1\right) \mathbf{Q}(z) = \mathbf{A}\left(z+1\right)\mathbf{Q}(z) + \mathbf{B}\left(z+1\right)\mathbf{X}(z) $$

$$ z\left(\alpha\mathbf{I}-A\right)\mathbf{Q}(z) = \left(\alpha\mathbf{I}+A\right)\mathbf{Q}(z) + \mathbf{B}\left(z+1\right)\mathbf{X}(z) $$

$$ z\mathbf{Q}(z) = \left(\alpha\mathbf{I}-A\right)^{-1}\left(\alpha\mathbf{I}+A\right)\mathbf{Q}(z) + \left(\alpha\mathbf{I}-A\right)^{-1}\mathbf{B}\left(z+1\right)\mathbf{X}(z) $$

This instantly gives us $\mathbf{A_d}$ as:

$$ \mathbf{A_d} = \left(\alpha\mathbf{I}-A\right)^{-1}\left(\alpha\mathbf{I}+A\right) $$

We know that the linear $(z+1)$ term of $\mathbf{X}(z)$ can be applied wherever $\mathbf{Q}(z)$ is used. We apply this in the output expression to give the new output discrete equation as:

$$ z\mathbf{Q}(z) = \left(\alpha\mathbf{I}-A\right)^{-1}\left(\alpha\mathbf{I}+A\right)\mathbf{Q}(z) + \left(\alpha\mathbf{I}-A\right)^{-1}\mathbf{B}\mathbf{X}(z) $$ $$ \mathbf{Y}(z) = \mathbf{C}\left(z+1\right)\mathbf{Q}(z) + \mathbf{D}\mathbf{X}(z)$$ $$ \mathbf{Y}(z) = \mathbf{C}z\mathbf{Q}(z) + \mathbf{C}\mathbf{Q}(z) + \mathbf{D}\mathbf{X}(z)$$

We can now pull out $\mathbf{B_d}$ as:

$$ \mathbf{B_d} = \left(\alpha\mathbf{I}-A\right)^{-1}\mathbf{B} $$

Because we know an expression for $z\mathbf{Q}(z)$, we can formulate $\mathbf{Y}(z)$ as:

$$ \mathbf{Y}(z) = \mathbf{C}\left(\mathbf{A_d}\mathbf{Q}(z) + \mathbf{B_d} \mathbf{X}(z) \right) + \mathbf{C}\mathbf{Q}(z) + \mathbf{D}\mathbf{X}(z)$$

$$ \mathbf{Y}(z) = \mathbf{C}\left(\mathbf{A_d}+\mathbf{I}\right)\mathbf{Q}(z) + \left(\mathbf{C}\mathbf{B_d} + \mathbf{D}\right)\mathbf{X}(z) $$


$$ \mathbf{C_d} = \mathbf{C}\left(\mathbf{A_d}+\mathbf{I}\right) $$ $$ \mathbf{D_d} = \mathbf{C}\mathbf{B_d} + \mathbf{D} $$

  • $\begingroup$ I noticed that your solution is slight different from the result in 1, I also find different solution here. Obviously, moving the linear term $(z+1)$ from $B$ matrix to $C$ matrix would not change the total transfer function. But would inner state matrix $Q(z)$ lose its structure? It is no longer discretization of original $Q(s)$ but something else? $\endgroup$
    – HaoranLin
    Sep 14, 2020 at 3:35
  • $\begingroup$ I don't understand the part: We know that the linear $(z+1)$ term of $\mathbf{X}(z)$ can be applied wherever $\mathbf{Q}(z)$ is used. Why? $\endgroup$ Aug 10, 2021 at 15:35
  • $\begingroup$ First, if you're doing control systems then your knowledge of the parameters is inexact, so you may as well use something that makes the math easy. The forward Euler method works. I.e. $x_n = (\mathbf I + T_s \mathbf A)x_{n-1} + T_s \mathbf B$. Second, if you feel that exactitude is important, and since it's 2021 and we're surrounded by a sea of computers, use the exact method. Note that it's easy to get over the problem with $\mathbf A$ being singular. $\endgroup$
    – TimWescott
    Nov 15, 2021 at 16:07
  • $\begingroup$ @JuanGonzalezBurgos yes, this was very poor wording on my part. If 𝐘(𝑧)=𝐀(𝑧) and $\mathbf{A}(z)=(z+1)\mathbf{X}(z)$, because the expressions are linear, we can move the (𝑧+1) from the second expression into the first e.g. 𝐘(𝑧)=(𝑧+1)𝐀(𝑧) and $\mathbf{A}(z)=\mathbf{X}(z)$. This is equivalent to a change of variable - but I was lazy and didn't change the variable. $\endgroup$ Jul 27, 2022 at 9:35
  • $\begingroup$ @TimWescott thanks for the additional comments. The particular use-case surrounding my question didn't arise from control systems, but from attempting to build algorithms that would efficiently realise parallel-form (as opposed to cascade) discrete implementations of high-order filters. State-space representation is IMO widely under-appreciated for it's uses outside of control systems. Even regular digital IIR filters can have their performance dramatically improved by designing the filter in state space. $\endgroup$ Jul 27, 2022 at 9:43

2 Answers 2


I've had the same question last week, but I've managed to find how to derive it (getting rid of those $z$ terms is indeed tricky). I will give here detailed demonstration of how to arrive to the result given in 1 (with, in your notation, $\alpha = 2 \lambda$).

So we define our new discrete-time function transfer as

$$ \begin{array}{rcl} H_d(z) &=& H(\frac{\alpha\left(z-1\right)}{z+1}) \\ &=& \mathbf{C} \left[ \frac{\alpha\left(z-1\right)}{z+1}\mathbf{I}-\mathbf{A} \right]^{-1}\mathbf{B} + \mathbf{D} \\ &=& \frac{z+1}{\alpha} \mathbf{C} \left[ \left(z-1\right)\mathbf{I}- \frac{z+1}{\alpha}\mathbf{A} \right]^{-1}\mathbf{B} + \mathbf{D} \\ &=& \frac{z+1}{\alpha} \mathbf{C} \left[ z \left(\mathbf{I} - \frac{1}{\alpha} \mathbf{A} \right) - \left(\mathbf{I} + \frac{1}{\alpha} \mathbf{A}\right) \right]^{-1}\mathbf{B} + \mathbf{D} \\ \end{array} $$

For sake of notation, let $\mathbf{P} = \mathbf{I} - \frac{1}{\alpha} \mathbf{A}$ and $\mathbf{Q} = \mathbf{I} + \frac{1}{\alpha} \mathbf{A}$ (important remark useful later on: $\mathbf{P} + \mathbf{Q} = 2 \mathbf{I}$). Then

$$ \begin{array}{rcl} H_d(z) &=& \frac{z+1}{\alpha} \mathbf{C} \left[ z \mathbf{P} - \mathbf{Q} \right]^{-1}\mathbf{B} + \mathbf{D} \\ &=& \frac{z+1}{\alpha} \mathbf{C} \left[ z \mathbf{I} - \mathbf{P}^{-1} \mathbf{Q} \right]^{-1} \mathbf{P}^{-1} \mathbf{B} + \mathbf{D} \\ &=& \frac{z+1}{\sqrt{2 \alpha}} \mathbf{C} \left[ z \mathbf{I} - \mathbf{P}^{-1} \mathbf{Q} \right]^{-1} \left( \sqrt{\frac{2}{\alpha}} \mathbf{P}^{-1} \mathbf{B} \right) + \mathbf{D} \\ &=& \frac{z+1}{\sqrt{2 \alpha}} \mathbf{C} \left[ z \mathbf{I} - \mathbf{A_d} \right]^{-1} \mathbf{B_d} + \mathbf{D} \end{array} $$

where we defined

$ \mathbf{A_d} = \mathbf{P}^{-1} \mathbf{Q} $ and $ \mathbf{B_d} = \sqrt{\frac{2}{\alpha}} \mathbf{P}^{-1} \mathbf{B} $. Furthermore, we have

$$ \begin{array}{rcl} \frac{z+1}{\sqrt{2 \alpha}} \mathbf{C} &=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left( z + 1 \right) \\ &=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left( z \mathbf{I} + \mathbf{I} + \mathbf{A_d} - \mathbf{A_d} \right) \\ &=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left[ \left( z \mathbf{I} - \mathbf{A_d} \right) + \left( \mathbf{I} + \mathbf{A_d} \right) \right] \\ &=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left[ \left( z \mathbf{I} - \mathbf{A_d} \right) + \mathbf{P}^{-1} \left( \mathbf{P} + \mathbf{Q} \right) \right] \\ &=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left[ \left( z \mathbf{I} - \mathbf{A_d} \right) + 2 \mathbf{P}^{-1} \right] \end{array} $$

which gives us

$$ \begin{array}{rcl} H_d(z) &=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \left[ \left( z \mathbf{I} - \mathbf{A_d} \right) + 2 \mathbf{P}^{-1} \right] \left[ z \mathbf{I} - \mathbf{A_d} \right]^{-1} \mathbf{B_d} + \mathbf{D} \\ &=& \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \mathbf{B_d} + \sqrt{\frac{2}{\alpha}} \mathbf{C} \mathbf{P}^{-1} \left[ z \mathbf{I} - \mathbf{A_d} \right]^{-1} \mathbf{B_d} + \mathbf{D} \\ &=& \mathbf{C_d}(z\mathbf{I}-\mathbf{A_d})^{-1}\mathbf{B_d} + \mathbf{D_d} \end{array} $$ with $ \mathbf{C_d} = \sqrt{\frac{2}{\alpha}} \mathbf{C} \mathbf{P}^{-1} $ and $ \mathbf{D_d} = \mathbf{D} + \frac{1}{\sqrt{2 \alpha}} \mathbf{C} \mathbf{B_d} $. Replacing $\mathbf{P}$ and $\alpha = 2 \lambda$ gives the results in 1.

PS: the same idea can be used to prove the discrete-time state-space representation found using Generalized Bilinear Transform ($s \leftarrow \alpha \frac{z - 1}{\beta z + \left(1-\beta\right)}$) or First-order Holder methods.

  • $\begingroup$ Thanks Klaz for taking the time to post a solution! I managed to work through this one today after almost giving up. I think I went about it in a slightly different way and have added my own working to the original question. $\endgroup$ Nov 10, 2017 at 12:42

I wanted to add to this some specifics on the MATLAB implementation.

Continuous time system: $$\dot{x} = A x + B u$$

Laplace transform: $$sx = Ax + Bu$$

Tustin discretization: $$\frac{\alpha(z-1)}{z+1}x = Ax + Bu$$ where $\alpha = 2f_\text{s}$, $f_\text{s}$ being the sampling frequency.

When $z \ne -1$, multiplying through by $(z+1)$, and rearranging, $$ z\left[x - \frac{1}{\alpha}(Ax+Bu)\right] = x + \frac{1}{\alpha}(Ax + Bu)$$

Motivated by this, do the following coordinate transformation (this is what is done in MATLAB, see here): $$ \begin{aligned} x_\text{d} &= x - \frac{1}{\alpha}(Ax+Bu) \\ &= \left(I - \frac{1}{\alpha}A\right)x - \frac{1}{\alpha}Bu \end{aligned} $$

Then the inverse coordinate transformation, $$x = \left(I - \frac{1}{\alpha}A\right)^{-1} \left(x_\text{d}+\frac{1}{\alpha}Bu\right)$$

Substituting back, $$z x_\text{d} = \left(I + \frac{1}{\alpha}A\right) \left(I - \frac{1}{\alpha}A\right)^{-1}x_\text{d} + \left[\left(I + \frac{1}{\alpha}A\right) \left(I - \frac{1}{\alpha}A\right)^{-1} + I\right] \frac{1}{\alpha}Bu $$

Define $$ \begin{aligned} A_\text{d} &= \left(I + \frac{1}{\alpha}A\right) \left(I - \frac{1}{\alpha}A\right)^{-1} \\ B_\text{d} &= \frac{1}{\alpha}(A_\text{d} + I)B \end{aligned} $$

Then, $$x_\text{d}^{k+1} = A_\text{d}x_\text{d}^k + B_\text{d}u^k$$

To check if this is indeed what $\texttt{c2d}$ in MATLAB is doing, try the following:

% Construct a second order system to test
f = 5; % frequency
zet = 0.05; % damping ratio
w = 2*pi*f; % circular frequency

% state space representation
A = [0 1; -w^2 -2*zet*w]; 
B = [0; 1];
C = [1 0];
D = [];
sys = ss(A, B, C, D);

% Discrete time system using Tustin
Ts = 0.01; % sampling time (100Hz)
sysd = c2d(sys, 0.01, 'Tustin');

% Construct discrete time matrices using the formulas derived
alph = 2/Ts;
Ad = (eye(2)+A/alph)/(eye(2)-A/alph);
Bd = (1/alph)*(Ad + eye(2))*B;

% Compare
fprintf('Compare Ad:\n')
disp([Ad sysd.A])
fprintf('Compare Bd:\n')
disp([Bd sysd.B])

This gives $$A_\text{d} = \begin{bmatrix}0.9526 & 0.0096 \\-9.4865 & 0.9224 \end{bmatrix}$$ and $$B_\text{d} = \begin{bmatrix}0.0000 \\ 0.0096 \end{bmatrix}$$

  • $\begingroup$ This is a good answer, but I had already upvoted it. It's actually a really good answer. $\endgroup$ Mar 20 at 15:56

Your Answer

By clicking β€œPost Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.