You can view this result as the convolution of the signal with an exponentially decaying sine and cosine.

That is, `y(t') = integral e^kt x(t' - t) dt`, with k complex and negative real part.

If you discretize that using simple integration and t' = i dt, t = j dt you get

    y_i = dt sum_j e^(k j dt) x_{i - j}
    y_{i+1} = dt sum_j e^(k j dt) x_{i+1 - j}
            = (dt e^(k dt) sum_j' e^(k j' dt) x_{i - j'}) + x_i 
            = dt e^(k dt) y_i + x_i

If we then scale this by some value, such that A y_i = z_i we can write this as

    z_{i+1} = dt e^(k dt) z_i + A x_i

Here the `dt e^(k dt)` plays a similar role to (1-alpha) and A is similar to P alpha - the difference being that P changes over time, while A is constant.

We can write `z_i = e^{w dt i} r_i` where w is the imaginary part of k

   e^{w dt (i+1)} r_{i+1} = dt e^(k dt) e^{w dt i} r_i + A x_i
             r_{i+1} = dt e^((k - w) dt) r_i + e^{-w dt (i+1) } A x_i
                     = (1-alpha) r_i + p_i x_i

Where p_i = e^{-w dt (i+1) } A = e^{-w dt ) p_{i-1} Which is exactly the result from the resonate web-page.

The neat thing about recognising this as a convolution integral, is that we can use shaping other than exponential decay - we can implement a box filter using only two states, or a triangular filter (this is a bit trickier and takes more states). While they're tricky to derive, they tend to run really quickly.