How to mitigate floating point error in recursive DSP calculations

Question

I have implemented a modified linear prediction coding function like this:

def modified_lpc_process(s: NDArray, p: int, fs: int) -> NDArray:
a = lpc_coefficients(s, p)
a_hat = modified_pole_coeffs(a, fs=fs)
# a_hat = a # debug line
s_hat = numpy.zeros(len(s))

for n, sn in enumerate(s):
    as_accum = 0
    for k in range(len(a)):
        if n - 1 - k < 0:
            break
        as_accum += a[k] * s[n - 1 - k]
    un = sn + as_accum

    as_hat_accum = 0
    for k in range(len(a_hat)):
        if n - 1 - k < 0:
            break
        as_hat_accum += a_hat[k] * s_hat[n - 1 - k]

    s_hat[n] = -as_hat_accum + un
    round(s_hat[n], 8)

return s_hat

Which essentially implements the maths here:

\begin{align} u[n] &= s[n] + \sum_{k=1}^p a_k s[n-k]\\ \hat{s}[n] &= - \sum_{k=1}^p \hat{a}_k \hat{s}[n-k] + u[n] \end{align}

I have noticed that the floating point error in this operation quickly blows up and I end up with sample values too large. In an effort to debug, I've set the modified pole coefficients a_hat to be the same with a, expecting to get the same value for s_hat and s; but the same issue is still there. I've also just realised that in floating point, adding a number and subtracting the same number you end up with slightly different results, which in recursive applications like this can accumulate to large differences.

I wonder what the best practice is in this case to mitigate this error accumulation?

Answer 1

I think I fixed it. It is indeed because of having unstable poles, but I was only looking at each individual transfer function instead of the overall transfer function.

In LPC the individual transfer functions are in the form of

$$\frac{1}{1+\sum_{k=1}^{p}a_kz^{-k}}$$

However the overall transfer in my case, combining both lines of equations for $u[n]$ and $\hat{s}[n]$ becomes: $$ \frac{1+\sum_{k=1}^{p}a_kz^{-k}}{1+\sum_{k=1}^{p}\hat{a}_kz^{-k}} $$

So I needed to correct the poles for this denominator, which has a different sign for the sum term than the general LPC transfer function.

I believe this is because I am using the original signal as input, as opposed to traditional LPC which uses a generated excitation, therefore the transfer function has a different sign.