Frequency warping when integrators are replaced with backward-euler and forward-euler integration

Question

Resonant controllers are used in the power industry. The transfer function is $G_{res}(s) = \frac{K_i s}{s^2 +2\omega_c s+ \omega_o^2}$

The "textbook" discrete implementation is depicted in the image The first integrator is a forward Euler and the second one a backward integrator. It is not possible to use Tustin integrators because this would create algebraic loops. Note that $\omega_o$ is not a constant, it will vary slowly and will be fed by a PLL or FLL.

Here's my question. How can one characterize the frequency warping?

In one paper, they evaluate frequency warping by applying the bilinear transform to $G_{res}(s)$, I feel like this is wrong.

Answer 1

The important thing here is that there is no conventional frequency warping with the forward or backward Euler methods. Frequency warping would mean that the discrete-time (DT) and continuous-time (CT) frequency responses are related by

$$H_d(f)=H_c(\phi(f))\tag{1}$$

with some warping function $\phi(f)$ that defines a mapping from the continuos-time frequency on the real line $(-\infty,\infty)$ to the discrete-time frequency in the interval $[-f_s/2,f_s/2]$, where $f_s$ is the sampling frequency.

Such a mapping is what happens when using the bilinear transform (Tustin's method). This means that the amplitude characteristic of the CT remains unchanged apart from a compression. This is why filters with optimum amplitude (e.g., Butterworth, Chebyshev, Cauer) remain optimal after transformation to the discrete domain.

However, for the Euler methods there is no frequency warping as in $(1)$ because the $j\omega$-axis of the complex $s$-plane is not mapped to the unit circle in the complex $z$-plane. We can ask which contour in the $s$-plane is actually mapped to the unit circle in the $z$-plane, i.e., which part of the CT transfer function becomes the frequency response of the DT filter. This is easily answered by looking at the mappings:

$$\textrm{forward Euler: }s\leftrightarrow \frac{1}{T}(z-1)\tag{1}$$

$$\textrm{backward Euler: }s\leftrightarrow \frac{1}{T}(1-z^{-1})\tag{2}$$

If in $(1)$ and $(2)$ we replace $z$ by $e^{j\omega}$ (i.e., the unit circle), we see that the corresponding regions in the $s$-plane are circles with radius $1/T$, centered at $s=-1/T$ for the forward Euler method, and centered at $s=1/T$ for the backward Euler method. I.e., only the DC value of the CT frequency response gets mapped to the DC value of the DT frequency response, the rest of the DT frequency response is not obtained from the CT frequency response but from values of the CT transfer function that don't lie on the $j\omega$ (frequency) axis.

This is why we cannot actually speak of frequency warping for the Euler methods. Yet, the frequency response of the the DT system can be easily predicted by simply evaluating the CT transfer function along the above mentioned circles in the complex plane.

This is shown in the figure below. The blue curve is the CT frequency response, and the other curves are the CT transfer function evaluated along the two circles in the $s$-plane mentioned above ("circle 1" corresponds to the forward Euler method, and "circle 2" corresponds to the backward Euler method). The resulting curves are exactly identical to the DT frequency responses implemented by the respective version of the Euler method.

Answer 2

To the extent this helps, here are some interesting magnitude and phase plots of the backward Euler (here using the method of impulse invariance mapping which for 1/s has the same result as backward Euler but as Matt L points out in the comments is not generally the case) vs Tustin (method of Bilinear Transform mapping) integrators. Notice the Euler on the left has dominant phase distortion while the Tustin on the right has more dominant amplitude distortion (no phase distortion from an analog integrator at all!).

The analog frequency response for a pure integrator as $1/s$ is -20 dB/decade in magnitude ($1/f$) and -90° for phase.

Answer 3

Here's my question, how can one characterize the frequency warping ?

I believe that you'd just have to calculate the poles in the z domain as a function of the parameters. Frequency warping with Tustin or forward or backwards Euler works because there's a consistent mapping from the value of your "fake s" to the value of the actual $s$. I'm almost certain that you just plain lose that here.

I'm known to march to my own drummer, but I gave up on the whole frequency warping thing a long time ago, unless I'm forced into it. I just move things into the $z$ domain as soon as possible, and do my calculations there.

In one paper, they evaluate frequency warping by applying the bilinear transform to Gres(s), I feel like this is wrong.

Probably, but it may be close enough. Consider that they're tweaking the value of $\omega_0$ anyway, so if the frequency calculations are off by a bit then the PLL will bring the real resonant frequency of the real filter into line quickly enough.