Discretization method for a simple first order system

Question

I know there are multiple ways of discretizing a continuous system. In many occasion, I needed a discrete model for a simple first order system (RC circuit, inertial load, etc.) and most of the times I just went by intuition and converted it to something that ended up being an exponential averaging filter. I realized that I knew no formal method of passing from the continuous representation of the system to its discrete version that I just made by intuition.

Here's the approach I usually do. Consider a simple RC circuit with $$ i=\frac{V_R}{R} = \frac{V-V_c}{R} =C\frac{dV_c}{dt} $$ Rearrange and we get $$ dV_c = \frac{V-V_c}{RC}dt$$

Then I take that as is and make model like that

with $\tau=RC$, this would have the following discrete transfer function

$$ H(z) = \frac{T_s/\tau}{1 -(1-T_s/\tau)z^{-1}} $$

Which is, from what I know, an exponential averager with $\alpha = T_s/\tau$

Now. My problem is that I don't know what mathematical method can transform the continuous version of the RC circuit to what I modeled.

$$ \frac{1}{\tau s+1} \Rightarrow?\Rightarrow \frac{T_s/\tau}{1 -(1-T_s/\tau)z^{-1}}$$

I used Matlab to look at these methods : zoh, impulse invariance, foh, tustin, least-square. None gives exactly this result.

What's the name of the discretization method I used by intuition (if any)?

Answer 1

What you're doing is slightly different from the common mappings from the $s$-domain to the $z$-domain, but it's a peculiar mix of the well-known backward and forward Euler methods.

Using $x$ and $y$ for the input and output, respectively, your continuous-time system is

$$y'(t)=\frac{1}{\tau}\big[x(t)-y(t)\big]\tag{1}$$

The way you translate $(1)$ into discrete time is

$$\frac{y[n]-y[n-1]}{T_s}=\frac{x[n]-y[n-1]}{\tau}\tag{2}$$

The unconventional bit here is the delayed output $y[n-1]$ on the right-hand side of $(2)$. Backward Euler simply is

$$\frac{y[n]-y[n-1]}{T_s}=\frac{x[n]-y[n]}{\tau}\tag{3}$$

Forward Euler would result in

$$\frac{y[n+1]-y[n]}{T_s}=\frac{x[n]-y[n]}{\tau}\tag{4}$$

which is equivalent to

$$\frac{y[n]-y[n-1]}{T_s}=\frac{x[n-1]-y[n-1]}{\tau}\tag{6}$$

Now the difference with $(2)$ is the delayed input signal $x[n-1]$ on the right-hand side of $(6)$.

Because the input and output are not aligned in time on the right-hand side of $(2)$, unlike forward and backward Euler, your discretization method cannot be described by a mapping $s\Rightarrow g(z)$.

Answer 2

The somewhat annoying answer this question is: there is no single correct answer. You can't do this correctly and you have to pick whatever discretization artifacts is least objectionable for your specific application.

The underlying reason here is that a first order lowpass is unlimited in time (impulse response) and it frequency (transfer function), so it can't meet the requirements of the sampling theorem.

Let's look at a simple example: a 1kHz low pass filter sampled at $fs = 48 kHz$ with a sampling period of $T = 1/f_s$

The continuous transfer function and impulse response would be $$H(\omega) = \frac{1}{1+j\omega/\omega_c} \;\; h(t) = u(t)e^{-\omega_c t}$$

We can sample this both in the time one (impulse invariance), i.e. $$h_1[n] = e^{-\omega_c nT} $$

We can also just sample it in the frequency domain as $$H_2[k] = \frac{1}{1+j2\pi k /(N\omega_c)} $$

where N is the FFT length of your frequency vector.

The third method would be the bilinear transform: map the s-plane poles and zeros into the z-plane. You get a third transfer function

$$H_3[k] = g\frac{1-ze^{-j2\pi k/N}}{1-pe^{-j2\pi k/N}}$$ where $z$ is the mapped zero and $p$ the mapped pole and $g$ some gain coefficient.

All three methods give different results. Just for grins we will add a 4th one where we use the bilinear transform but move the zero from $z=-1$ a little inwards to, say $z = -.7$.

Let's see what we get

The purple line "sample in freq" represents the behavior of the continuous system. The blue line "sample in time" shows clear aliasing in the frequency domain. The spectrum must be real at the Nyquist frequency and hence the phase of any transfer function that's not zero at Nyquist must be a whole multiple of $\pi$. This results in significant phase error even at fairly low frequencies.

The bilinear version gets around this constraint by having a magnitude of 0 at Nyquist so it does a good job matching the phase, but the amplitude at high frequencies is way off. The "fudged" version is a compromise. It matches the amplitude very well up to 20 kHz and has less phase error than the "impulse invariance" version.

We see similar differences when we look at the time domain.

All impulses have the general exponential decay shape but there are significant differences as well. In particular, the "sample in freq" curve has a lot of ringing and is non-causal. So what's the closest in the frequency domain is the worst in the time domain (and vice versa).