Backsolving, polynomial division and deconvolution

Recall the filter transformation (4) which converts an input $\bold x$ to an output $\bold y$.It could be backsolved, namely, given $\bold y$, find $\bold x$.I copy that transformation (4) to here, swapping the definitions of $\bold x$ and $\bold y$to maintain the convention that the input is $\bold x$and the output is $\bold y$.Truncate the matrix to be square (to keep the number of unknowns equal the number of equations). Also for convention, divide b₀ out from each column thereby introducing a new filter response a_t = b_t/b₀. We have  

$$\bold A \bold y \eq \left[ \begin{array} {cccccc} 1 & 0 & ... ...\\ x_3 \\ x_4 \\ x_5 \end{array} \right] \eq \bold x$$ (23)

This operator goes under the various names, ``backsolving'', ``polynomial division'', and ``deconvolution''. It is straightfoward to solve the first equation for y₀; taking it to be known, solve the next for y₁, continuing to y₂, etc, until y₅. Mathematically, it is known that the back-solution of the adjoint is the adjoint of the backsolution. This is because the transpose of an inverse matrix is the same as the inverse of a transpose matrix.

The leaky integration transformation (19) is a special case of the backsolving transformation when $a_1=-\rho$ and a₂=0. To see this, you need to notice that the matrices in (23) and (19) are mutually inverse.

The transformation (23) can be derived by multiplying two polynomials, Y(Z) = y₀ + y₁ Z+ y₂ Z² + y₃ Z³ + y₄ Z⁴ + y₅ Z⁵ times A(Z) = 1 + a₁ Z + a₂ Z². Identifying the k-th power of Z in the product X(Z)=A(Z)Y(Z) gives the k-th row of the transformation (23). Thus the operator in (23) is Y(Z)=X(Z)/A(Z). The polynomials in Z are called Z transforms. They may be recognized as Fourier transforms where $Z=e^{i\omega\Delta t}$.

Convolution corresponds to multiplying polynomials (convolving the coefficients) and deconvolution (or backsolving) corresponds to polynomial division.

A causal operator is one that uses its present and past inputs to make its current output. Anticausal operators use the future but not the past. Causal operators are generally associated with lower triangular matrices and positive powers of Z, whereas adjoint operators are associated with upper triangular matrices and negative powers of Z.

Ordinary differential equations lead us to the backsolving operator. For example, the damped pendulum leads to equation (23) with the three nonzero filter coefficients. There is a huge literature on finite-difference solutions of ordinary differential equations, especially in connection with electrical circuits. You might notice the nonphysical meaning of the adjoint. It is like the forward operator except that we must begin at the final time and revert towards the first. The adjoint pendulum damps as we compute it backward in time, but that, of course, means that the adjoint pendulum diverges as it is viewed moving forward in time.

A module to backsolve polydiv is polydiv.  

module polydiv {			                  # Polynomial division
real,    dimension (:), pointer :: a
integer, dimension (:), pointer :: lag
#%  _init ( a, lag)
#%  _lop  ( x, y)
integer  ia, ix, iy
real     t
if( adj)
	do ix= size(x), 1, -1 {
		t = y( ix)
		do ia = 1, size(a) {
			iy = ix + lag( ia);     if (iy > size(x)) exit
			t = t - a( ia) * x( iy)
			} 
		x( ix) = x( ix) + t
		}
else 
	do iy= 1, size(x) {
		t = x( iy)
		do ia = 1, size(a) {
			ix = iy - lag( ia);      if (ix < 1)  exit
			t = t - a( ia) * y( ix)
			} 
		y( iy) =  y( iy) + t
		}
}

Inputs required for module polydiv to backsolve polydiv are nb=2 and the filter coefficients and lags in two tables as follows:

                  i    lag(i)   b(i)
                 ---   ------   ----
                  1      1       a1
                  2      2       a2

We could simplify use of module polydiv by internalizing lag(i)=i however, later in this book, we will encounter examples with many zero valued coefficients where great efficiency is gained by omitting the zero-valued coefficients from the tables.