# Convert filtering source C-code into difference equation

I'm struggling to convert three simple lines of code into a difference equation to calculate the frequency response.

The C-code is as simple (and legacy), as

float temp = 0.0f;
const a = 6;
void filter(float *inOut)
{
temp += *inOut;       // add *input* sample to accum
*inOut = temp / 2^a;  // scale accum and assign *output*
temp -= *inOut;       // subtract *output* from accum
}


inOut is the data to be filtered. Let's assume, at each function call it's a one to calculate the step response.

There's no documentation, but it seems to be a version of the recursive moving average filter, i.e. $$y = (1-\alpha)y_\text{old}+\alpha x$$, but I can't prove it.

I've re-implemented the code in Matlab and plotted the step response, looks like a moving average filter.

I'd like to write down the difference equation at the $$k-$$th timestep, this is my attempt:

$$\text{temp}_k = \text{temp}_{k-1} + \text{in}_k \qquad(1)\\ \text{out}_k = \text{temp}_k / 2^a \qquad(2)\\ \text{temp}_{k} = \text{temp}_k-\text{out}_k \qquad(3)$$

But what now? How can I eliminate $$\text{temp}$$? I'm expecting to get a difference equation with variables out (or y) and in (or x) and some coefficients. Then, I'd like to use z-transform to get the complex transfer function of that filter for further analysis.

I think, re-ordering the last line to become the first one and shift the indices properly should do the trick, but I don't find a solution yet.

Any hints are greatly appreciated, thanks!

## EDIT 1:

Looks like eqn (3), $$\text{temp}_{k+1} = \text{temp}_k-\text{out}_k$$ is wrong. This is not timestep $$k+1$$ but $$k$$, does that make sense? I've edited it.

## EDIT 2:

My calculation steps so far:

Shift eqn. (3) one timestep. Is this allowed? $$\text{temp}_{k} = \text{temp}_k-\text{out}_k \\ \Rightarrow \text{temp}_{k-1} = \text{temp}_{k-1}-\text{out}_{k-1}$$

Substitute into (1)

$$\text{temp}_k = (\text{temp}_{k-1}-\text{out}_{k-1}) + \text{in}_k \qquad(4)$$

Transform, timeshift and substitute (2) into (4)

$$\text{out}_k = \text{temp}_k / 2^a \Leftrightarrow \text{temp}_k = \text{out}_k 2^a \\ \Rightarrow \text{out}_{k} 2^a = (\text{out}_{k-1} 2^a-\text{out}_{k-1}) + \text{in}_k$$

Rearrange and solve for $$\text{out}_k$$

$$\text{out}_{k} 2^a = (\text{out}_{k-1} 2^a-\text{out}_{k-1}) + \text{in}_k \\ \Leftrightarrow \text{out}_{k} 2^a = \text{out}_{k-1}( 2^a-1) + \text{in}_k \\ \Leftrightarrow \text{out}_{k} = \text{out}_{k-1}(1-2^{-a}) + 2^{-a}\text{in}_k$$

This is the expected difference equation. z-transform and transform to transfer function form (with X = input and Y = output)

$$Y = z^{-1}Y(1-2^{-a})+2^{-a}X \\ \Leftrightarrow \frac{Y}{X} = \frac{2^{-a}}{1-z^{-1}(1-2^{-a})}$$

Big question: Is this correct? May I timeshift my difference equations and substitute them as shown?

Thanks!

• Your variable in (line 5) is undefined. Where is the input coming from ? Commented Nov 25, 2022 at 12:42
• @Hilmar, sorry, fixed. inOut variable was meant. It's always one at call time and holds the results after the function has run. I've edited the question.
– Jan
Commented Nov 25, 2022 at 12:57

This one is a little difficult to transcribe since it uses the same variable temp for two different purposes. Across calls it's your state variable but during the call it holds a temporary variable. We will rewrite this using input $$x[n]$$, state variable $$u[n]$$, output $$y[n]$$ and the filter coefficient $$b = 1/2^6$$

$$k = u[n-1]+x[n] \tag{1}$$ $$y[n] = b \cdot k = \tag{2} b\cdot (u[n-1] + x[n])$$ $$u[n] = k - y[n] = \frac{y[n]}{b} -y[n] = y[n]\frac{1-b}{b}\tag{3}$$

In the last equation we have used a relation ship from eq. (2) $$k = y[n]/b$$ Now we can pop eq. (3) into equation (2) and we get

$$y[n] = b\cdot (u[n-1] + x[n]) = b\cdot \left(\frac{1-b}{b}y[n-1] + x[n]\right)$$ $$= b\cdot x[n] + (1-b)\cdot y[n-1] \tag{4}$$

And now it's finally in the form of a standard difference equation and we can get the z transform as

$$H(z) = \frac{b}{1-(1-b)z^{-1}}$$

This is indeed the cheapest version of a first order lowpass filter but doesn't have a zero at Nyquist. For that the difference equation would have to be

$$y[n] = b/2\cdot x[n] + b/2\cdot x[n-1]+ (1-b)\cdot y[n-1] \tag{5}$$

and you would have to add a second state variable.

You got it!

Except it's a recursive exponential moving average filter with $$\alpha = 2^{-a}$$

The difference equation is indeed: $$y[n] = \alpha x[n] + (1-\alpha)y[n-1]$$ $$\mathcal{Z}$$-transform to get: $$Y(z) - (1-\alpha)Y(z)z^{-1} = \alpha X(z)$$ $$\Leftrightarrow H(z) = \frac{Y(z)}{X(z)} = \frac{\alpha}{1-(1-\alpha)z^{-1}} = \frac{0.0156}{1 - 0.9844z^{-1}}$$

• Nice, thank you so much. I've just edited my question and added calculation steps. Would you agree that these are correct? How did you calculate the difference equation from the C code?
– Jan
Commented Nov 25, 2022 at 14:42
• @Jan Yes, albeit a little complicated but your derivation looks correct, and you get to the correct result! To be honest, through experience you learn to recognize these things and I didn't need to go through a derivation of the difference equation, but if I had to I probably would have done something similar to what you did, or writing out temp and inOut at timesteps 0,1,2 etc and figuring it out. If you're satisfied with an answer (either mine or @hilmar’s), please hit the accept button so others can benefit as well.
– Jdip
Commented Nov 26, 2022 at 15:49