# Frequency prewarping of a bilinear transform (Tustin transform)

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.

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:

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:

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.

• @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. – fibonatic Feb 21 '20 at 2:06
• @fibonatic I see- so c2d does include Tustin! Thanks – Dan Boschen Feb 21 '20 at 3:46

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

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.

• 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. – Matthias La Feb 21 '20 at 9:09
• @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$. – fibonatic Feb 24 '20 at 11:20