sadly I'm not too well aquinted with discretazation methods. At the moment I struggle to reproduce a bilinear transform with frequency prewarping.

I have the following transfer function in the s-domain:

$$\begin{align} \\ G(s)=\frac{-ms^2}{s^2+2abs+b^2} \\ \end{align}$$

There m,a,b are just some constants. Now I would like to transform this transfer function into the z-domain. In Matlab I am using the c2d-function. I compared all the available methods with a sampling frequency of 512Hz.

Comparison of discretization methods in matlab

One can clearly see, that with a pole being around 100Hz that is already too close to the Nyquist frequency, so that the Tustin method has to be wrapped in order to give a more accurate result.

Now I tried to reproduce the results from the matlab c2d-function by doing the bilinear transform myself.

I first tried it without prewarping using:

$$\begin{align} \\ s=A\cdot\frac{z-1}{z+1}, A=\frac{2}{T} \\ \end{align}$$

Comparing my result with Matlab does look really nice:

Tustin wo prewarping

Then I tried it with prewarping, using:

$$\begin{align} \\ s=A\cdot\frac{z-1}{z+1}, A=\frac{b}{\tan(b\cdot\frac{T}{2})}, b=2\pi\cdot100 \\ \end{align}$$

This results in:

Tustin w/ and wo prewarping

With the warping the result doesn't look the same as with the matlab c2d-function. So where is the difference??

In matlab I use:

G(z)=c2d(S,1/512,['Method','tustin', 'PrewarpFrequency', b]);

And I have derived:

$$\begin{align} \\ G(z)=\frac{-mA^2\cdot z^2+2mA^2\cdot z-mA^2}{(A^2+2abA+b^2)\cdot z^2+(2b^2-2A^2)\cdot z+(A^2-2abA+b^2)} \\ \end{align}$$

If I put in $A=2/T$ I achieve to get the same result as the Matlab c2d-function. But with $A=b/\tan(\frac{b\cdot T}{2})$ I don't get the same result as the c2d-function in Matlab.

Putting in some values for the constants I get he following:

Tustin wo prewarp using c2d yields: $\frac{-6.781 z^2+13.56 z - 6.781}{z^2-0.8456 z + 0.8669}$

Tustin with prewarp using c2d yields: $\frac{-7.997 z^2+15.99 z-7.997}{z^2-0.615 z+0.8217}$

My self derived Trustin-transform with prewarp yields: $\frac{-6.216 z^2+12.43 z-6.216}{z^2-0.6266 z+0.8599}$

So apparently I am doing something different than the matlab function. Maybe someone can enlighten me by pointing out where I am at odds. I also failed to recreate the discretized transfer function using the matched-method from c2d. Basically my problem is that I need to build nice discrete transfer functions while being unable to use matlab. So I would like to understand how the algorithm works. Thank you in advance for any help.

  • $\begingroup$ @DanBoschen There are many different discretization methods for c2d. It could be that matlab keeps adding more and more methods since the previous time I checked it did not have the least-squares method. $\endgroup$
    – fibonatic
    Feb 21, 2020 at 2:06
  • $\begingroup$ @fibonatic I see- so c2d does include Tustin! Thanks $\endgroup$ Feb 21, 2020 at 3:46

1 Answer 1


There might be some undocumented features for the c2d function. Namely, if I follow the documented way of specifying pre-warping (using c2dOptions) one would get:

opts = c2dOptions('Method', 'tustin', 'PrewarpFrequency', b);
G = c2d(S, T, opts);

When I used this I do get the same result as when I would manually do pre-warping using:

z = tf('z', h);
sp = b / tan(b * T / 2) * (z - 1) / (z + 1);
G = -m * sp^2 / (sp^2 + 2 * a * b * sp + b^2);

However, the notation that you use is not specified in the documentation for as far I could find. It does seem to give a better approximation of the continuous transfer function in this case. However, if I use the following transfer function

$$ P(s) = \frac{s}{s^2 + 0.1\,s + 1}, $$

with a sample time of $T = 1$ your method behaves in an unexpected way as shown in the figure below. In that figure D1 uses pre-warping using c2dOptions, D2 uses your method of pre-warping and D3 uses manual pre-warping (this completely overlaps with D1 in the figure).

enter image description here

It can be noted that D2 does follow the overal shape of $P(s)$ better, but it does not match it at the specified pre-warp frequency and the phase is much lower (when zoomed out to a bigger frequency spectrum it looks like a delay). On the other hand D1/D3 exactly matches $P(s)$ at the pre-warp frequency as expected, but the overall shape does deviate more the further away you get from the pre-warp frequency.

  • $\begingroup$ Thank you very much for your great answer! So, I understand while in my specific situation the documented c2d tustin method of matlab is worse in matching the continous function, there are also examples for the other way around. So I can also just stick with what I have. Do you further know, if there is a way to derive what matlab is doing with the undocumented tustin method? Or how I would implement the matched method? Because I failed with the additional poles and the DC gain. $\endgroup$ Feb 21, 2020 at 9:09
  • $\begingroup$ @MatthiasLa I do not know where to find more information about the undocumented behavior of c2d. Maybe you could try contact Matlab itself. As for the matched method I am not sure what you mean with additional poles. The gain can I think be obtained by setting their low frequency behavior (frequency of zero) equal to each other. So evaluating the continuous at $s=0$ and the discrete at $z=1$. $\endgroup$
    – fibonatic
    Feb 24, 2020 at 11:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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