# Do you recognize this discrete realization of a low pass filter?

I have a simple 1st order lag filter implementation that performs well, shown in the block diagram below. This is legacy 'code' we've been using for many years, the original developers have long since moved on. I've been asked to doc the requirements (development of the difference equation) starting from the s domain and assuming $$\frac{Y(s)}{U(s)} = \frac{1}{\tau s+1}$$.

$$k1$$ is $$0.5+\frac{\tau}{T_s}$$, $$k2$$ is $$0.5$$. $$\tau$$ is the time constant, in sec, and $$T_s$$ is the time between samples, in sec.

My assumption was I would arrive at a bilinear transform of $$\frac{1}{s+1}$$ but it's just not happening as the topology is wrong. Any hints?

Thanks, Chris

• You didn't show us how you arrived at that topology. – Matt L. Aug 12 '20 at 8:03

This is a discrete time filter and you have a block diagram, so forget the s-domain, and just work directly with difference equations and the Z transform.

In your block diagram, let's call the signal going into the delay block $$w[n]$$ and the signal coming out of the division block $$v[n]$$.

Writing some initial difference equations:

\begin{align*} w[n] &= w[n-1] + v[n] \quad \Rightarrow \quad v[n] = w[n] - w[n-1]\\ \\ v[n] &= \dfrac{1}{k_1} u[n] - \dfrac{1}{k_1} w[n-1] \\ \\ y[n] &= w[n] - k_2 v[n]\\ \end{align*}

Performing some simple algebra, we get:

\begin{align*} y[n] &= \left(1 - k_2\right) w[n] +k_2 w[n-1]\\ \\ w[n] - \left(1 - \dfrac{1}{k_1}\right)w[n-1] &= \dfrac{1}{k_1}u[n]\\ \end{align*}

So take the Z transforms and manipulate:

\begin{align*} Y(z) &= \left(1 - k_2\right) W(z) +k_2 z^{-1}W(z)\\ \\ \dfrac{Y(z)}{W(z)} &= \left(1 - k_2\right) +k_2 z^{-1} \\ \\ \dfrac{Y(z)}{W(z)} &= \left(1 - k_2\right) \left(1 - \dfrac{k_2}{k_2-1}z^{-1}\right)\\ \\ W(z) - \left(1-\dfrac{1}{k_1}\right)z^{-1} W(z) &= \dfrac{1}{k_1}U(z)\\ \\ \dfrac{W(z)}{U(z)} &= \dfrac{1}{k_1} \cdot \dfrac{1}{1-\dfrac{k_1 -1}{k_1}z^{-1} }\\ \end{align*}

The final transfer function is then:

\begin{align*} \dfrac{Y(z)}{U(z)} &= \dfrac{Y(z)}{W(z)}\cdot\dfrac{W(z)}{U(z)}\\ \\ &= \dfrac{1-k_2}{k_1} \cdot \dfrac{1 - \dfrac{k_2}{k_2-1}z^{-1}}{1-\dfrac{k_1 -1}{k_1}z^{-1} }\\ \\ &= \dfrac{0.5}{0.5 + \frac{\tau}{T_s}}\cdot \dfrac{1 - (-1)z^{-1}}{1-\dfrac{\left(0.5 + \frac{\tau}{T_s}\right) -1}{0.5 + \frac{\tau}{T_s}}z^{-1} }\\ \end{align*}

You have a zero at the Nyquist frequency, $$-1$$ on the unit circle, so this is definitely looking like a lowpass filter.

You have a pole on the real axis just inside the unit circle, near $$1$$, at

$$z_p = \dfrac{\left(0.5 + \frac{\tau}{T_s}\right) -1}{0.5 + \frac{\tau}{T_s}} = 1 - \dfrac{1}{0.5 + \frac{\tau}{T_s}}$$

which should provide a strong response near DC. This is what you would expect from a lowpass filter.

You can do some algebra and take an inverse Z-transform to get the final difference equation:

\begin{align*} \left( 1-\dfrac{\left(0.5 + \frac{\tau}{T_s}\right) -1}{0.5 + \frac{\tau}{T_s}}z^{-1}\right)Y(z) = \dfrac{0.5}{0.5 + \frac{\tau}{T_s}}\cdot \left(1 - (-1)z^{-1} \right) U(z) \\ \\ y[n] - \dfrac{\left(0.5 + \frac{\tau}{T_s}\right) -1}{0.5 + \frac{\tau}{T_s}} y[n-1] = \dfrac{0.5}{0.5 + \frac{\tau}{T_s}} u[n] + \dfrac{0.5}{0.5 + \frac{\tau}{T_s}} u[n-1]\\ \\ y[n] = \left(1-\dfrac{1}{0.5 + \frac{\tau}{T_s}}\right) y[n-1] + \dfrac{0.5}{0.5 + \frac{\tau}{T_s}} u[n] + \dfrac{0.5}{0.5 + \frac{\tau}{T_s}} u[n-1] \\ \end{align*}

So the filter is a scaled average of the two most recent input values with decaying feedback of the output value.

This could have been implemented with a much more straightforward block topology than the wonky one you were given.

In fact this difference equation is a 2 sample Moving Average filter, followed by an Exponentially Weighted Moving Average (EWMA) filter:

\begin{align*} x[n] &= 0.5u[n] + 0.5u[n-1] \quad \text{(Moving Average)}\\ \\ y[n] &= \left(1-\dfrac{1}{0.5 + \frac{\tau}{T_s}}\right) y[n-1] + \dfrac{1}{0.5 + \frac{\tau}{T_s}} x[n] \quad \text{(EWMA)} \\ \end{align*}

From this answer: https://dsp.stackexchange.com/a/40465/28112 the -3 dB cutoff frequency of your EWMA filter is:

$$f_{3dB} = \dfrac{1}{2\pi T_s} \cos^{-1}\left[\dfrac{\left(\dfrac{1}{0.5 + \frac{\tau}{T_s}}\right)^2}{2\left(1-\dfrac{1}{0.5 + \frac{\tau}{T_s}}\right)}\right]$$

Two other observations for you:

$$k_2$$ sets the balance of the weight given to the inputs $$u[n]$$ vs. $$u[n-1]$$

$$k_1$$ sets the $$\alpha$$ of the EWMA filter, which sets the balance of the weight given to the previous output vs. the new input.

Update

So now that we know where the digital filter's pole is, let's work through from the continuous time filter's pole, to see if it matches, using the bilinear transform relationship.

\begin{align*} s_p = -\dfrac{1}{\tau} &= \dfrac{2}{T_s} \cdot \dfrac{z_p -1}{z_p+1}\\ \\ z_p+1 &= -2\dfrac{\tau}{T_s}\left(z_p-1\right) \\ \\ z_p\left(1 + 2\dfrac{\tau}{T_s}\right) &= -1 + 2\dfrac{\tau}{T_s} \\ \\ z_p\left(1 + 2\dfrac{\tau}{T_s}\right) &= -2 + \left(1 + 2\dfrac{\tau}{T_s}\right) \\ \\ z_p\left(0.5 + \dfrac{\tau}{T_s}\right) &= \left(0.5 + \dfrac{\tau}{T_s}\right) -1 \\ \\ z_p & = \dfrac{\left(0.5 + \dfrac{\tau}{T_s}\right) -1}{0.5 + \dfrac{\tau}{T_s}} = \dfrac{k_1 -1}{k_1} \end{align*}

So that checks out. Someone did apparently start from a single pole filter of the form

$$H_a(s) = g\cdot\dfrac{1}{1+\tau s}$$

$$H_a(s)$$ has a zero as $$s \rightarrow \infty$$ (along any direction from the origin of the s-plane). Looking again at the bilinear transform relationship:

\begin{align*} s_0 &= \dfrac{2}{T_s} \cdot \dfrac{z_0 -1}{z_0+1}\\ \\ \lim_{z_0 \to -1} s_0 &\rightarrow -\infty\\ \\ z_0 &= -1\\ \end{align*}

Which also agrees with the zero we found analyzing the digital filter.

Update 2

And just to give some insight as to why $$k_2 = 0.5$$

\begin{align*} s_0 &= \dfrac{2}{T_s} \cdot \dfrac{z_0 -1}{z_0 +1}\\ \\ &= \dfrac{2}{T_s} \cdot \dfrac{\frac{k_2}{k_2-1} -1}{\frac{k_2}{k_2-1} +1} \\ &= \dfrac{2}{T_s} \cdot \dfrac{k_2 - k_2 + 1}{k_2 + k_2-1}\\ \\ &= \dfrac{2}{T_s} \cdot \dfrac{ 1}{2k_2 - 1}\\ \\ &= \dfrac{1}{T_s} \cdot \dfrac{ 1}{k_2 - 0.5}\\ \end{align*}

• That's awesome Andy, thanks so much, I owe you a beer or coffee or both! I'm a bit dismayed by the answer though, in that, if I read correctly, there really isn't a meaningful S domain equivalent of this filter. Our development philosophy is to design, simulate and analyze (for stability, performance) in the S domain and implement in the Z within limitations (such as using a bilinear transform and limits on tau) that ensure those analysis translate. Was really hoping for a 'this is a 1st order B – c hecker Aug 13 '20 at 14:46
• You're welcome. It's not that this couldn't have come from a continuous time/s-domain model first, but working out that lost rationale is going to be near impossible. You could make up anything that comes to the right answer, but you'll never know for sure. From the z-domain analysis, the final difference equation breaks out so nicely into a MA and EWMA filter, I'd have a hard time believing it wasn't designed with those filters in mind. – Andy Walls Aug 13 '20 at 15:11