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

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.

• I think I agree with you that a bilinear transform would not apply but without looking at the paper may be out of context: the bilinear transform of a analog integrator is the cascade of an accumulator which is z/(z-1) with a 2 point moving average resulting in (z+1)/(z-1). I don't see such a structure in your block diagram but would point out that the Tustin (bilinear) integrator structure has -90 degree shift over all frequencies which can therefore help to better match the analog response that it was mapped from. The use of Ts without the s subscript for the sampling period is confusing. – Dan Boschen Dec 12 '19 at 16:05
• But to "characterize" frequency warping I would simply compare G(s) to G(z). – Dan Boschen Dec 12 '19 at 16:06
• That's a good point, I will edit my block diagram. – Ben Dec 12 '19 at 16:11
• Also how does the Tustin create algebraic loops? It is just the Backward Euler cascaded with a 2 point moving average, isn't it? – Dan Boschen Dec 12 '19 at 16:12
• 2 Forward Euler integrators will create an algebraic loop, A tustin integrator is an average between a forward euler and a backward euler, therefore 2 tustin integrators will create an algebraic loop. I will try to check out the math to make sure that there really is an algebraic loop. – Ben Dec 12 '19 at 16:20

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.

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.

• If you try to compensate the 9th harmonic for example, the ratio between $f_{fundamental} and$f_{h9} will not be exactly 9 because of frequency warping, – Ben Dec 12 '19 at 16:21
• I like Tim's point--- just get into digital and not be concerned with the analog----- meaning start with and stay in the digital (z domain) and not try to copy some analog analysis. There was a day when there was so much more established in the analog domain than direct digital design that mapping made sense, but now copying the analog and mapping from it only leads to the burden of the limitations the analog design constrains you to. Mapping makes sense when you are trying to model an analog system in the digital domain but not to create a digital design. – Dan Boschen Dec 12 '19 at 16:36
• So to compute the 9th harmonic-- just compute it using G(z) and you will have the exact result. – Dan Boschen Dec 12 '19 at 16:41
• Yeah, I will try computing the 9th harmonic using G(z) – Ben Dec 12 '19 at 17:19

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.

• It's maybe worth mentioning that the backward Euler method and the impulse invariant method are generally not the same. Just in the special case of transforming an ideal integrator they result in the same discrete-time system (an accumulator). – Matt L. Dec 13 '19 at 14:42
• Thanks Matt - is misleading as written – Dan Boschen Dec 13 '19 at 14:50