Bilinear Transform (Tustin's Method) applied to the Derivative I hope that I have not misunderstood something terribly wrong, but the continuous derivative $D=d/dt$ can be considered a transfer function in Laplace space $D(s) = s$, right?

So when I try to discretize it using the bilinear transform (Tustin's method) I trivially get

$D(z) = \frac{2}{T} \frac{1-z^{-1}}{1+z^{-1}}$

When I apply this to a series containing one discrete impulse, the response oscillates at the Nyquist frequency. Even worse, the spectrum around $\omega=0$ is quadratic and not $\sim i\omega$ like it would be expected from the derivative. (Edit: the latter was just due to roundoff error, because the low-freq amplitude got swamped by the Nyquist-peak)

Although I know of course how well the bilinear transform works for discretizing all kinds of filters, it is hard for me to understand, why it is considered superior if it seems to fail so miserably for something as simple as the derivative, which can otherwise be easily represented by first order finite differences:

$D(z)=\frac{1-z^{-1}}{T}$

or even second order (symmetric) finite differences

$D(z)=\frac{1-z^{-2}}{2Tz^{-1}}$

I am sure it all has a very simple explanation, but I can't see it.

PS: What is all the more confusing: when I apply the bilinear D(z) to a step function, the result is (correctly) a single peak. Consequently the inverse of the bilinear D(z) applied to an impulse yields the step function, like it has to be. What is going on there?
It's a matter of perspective if the bilinear transform "fails so miserably" when trying to approximate a derivative. First of all, it's not true that the approximation is quadratic for small $\omega$, it's linear as it should be:

$$D(e^{j\omega})=\frac{2}{T}\frac{1-e^{-j\omega}}{1+e^{-j\omega}}=\frac{2j}{T}\tan\left(\frac{\omega}{2}\right)\approx \frac{j\omega}{T}\quad\textrm{for}\quad |\omega|\ll 1$$

So for small frequencies the discrete-time system behaves as a (discrete-time) differentiator. For larger frequencies we expect the effect of frequency warping, which is exactly what happens here: the infinite continuous-time frequency range is mapped to the interval $\omega\in[0,\pi]$.

In sum, the bilinear transform behaves as expected, but it is also clear that there are better ways to approximate a derivative in discrete-time than applying the bilinear transform to the continuous-time transfer function $H(s)=s$.