Stationary convolution and inverse convolution

The convolution of a vector, ${\bf x}=(x_0 \;\; x_1 \;\; x_2 \;\; ... \;\; x_{N-1})^T$, with a causal filter, ${\bf a}=(a_0 \;\; a_1 \;\; a_2 \;\; ... \;\; a_{N_a-1})^T$,whose first element, a₀=1, and whose length, N_a<N, onto an output vector, ${\bf y}=(y_0 \;\; y_1 \;\; y_2 \;\; ... \;\; y_{N-1})^T$, can be defined by the set of equations:  

$$y_k = x_k + \sum_{i=1}^{\min(N_a-1, k-1)} a_{i} \; x_{k-i}.$$ (1)

This can be rewritten in linear operator notation, as ${\bf y}={\bf A}{\bf x}$, where ${\bf A}$ is a lower-triangular Toeplitz matrix representing convolution with the filter, ${\bf a}$.

The adjoint operator, ${\bf A}'$, which describes time-reversed filtering with filter ${\bf a}$, can similarly be expressed by considering the rows of the matrix-vector equation, ${\bf x}={\bf A}' {\bf y}$, as follows,  

$$x_k = y_k + \sum_{i=1}^{\min(N_a-1, N-1-k)} a'_{i} \; y_{k+i}.$$ (2)

The helicon Fortran90 module Claerbout (1998a) exactly implements the linear operator (and adjoint) pair described by equations (1) and (2).

Equation (1) explicitly prescribes internal boundary conditions near k=0; however, since ${\bf a}$ is causal, no particular care is needed near k=N-1. On the other hand, equation (2) explicitly imposes internal boundary conditions near k=N-1, and no care is needed near k=0. It is possible to rewrite equations (1) and (2) in a more symmetric form; however, as written, the equations lead naturally to recursive inverses for operators, ${\bf A}$ and ${\bf A}'$.

Rearranging equation (1), we obtain  

$$x_k = y_k - \sum_{i=1}^{\min(N_a-1, k-1)} a_{i} \; x_{k-i}.$$ (3)

This provides a recursive algorithm, starting from x₀=y₀ for solving the system of equations, ${\bf y}={\bf A}{\bf x}$.Equation (3) describes the exact, analytic inverse of causal filtering with equation (1). In principle, given a filter, ${\bf a}$, and a filtered trace, ${\bf y}$, the above equation can recover the unfiltered trace, ${\bf x}$ exactly; although in practice, with numerical errors, the division may become unstable if ${\bf a}$ is not minimum phase. Similarly, if we inverse filter a trace with equation (3), we can recover the original by causal filtering with equation (1) subject to the stability of the inverse filtering process.

Equation (3) appears very similar to polynomial division. However, the output of polynomial division is an infinite series, while equation (3) is defined only in the range, $0 \leq k \leq N-1$. As such, equation (3) describes polynomial division followed by truncation.

Equation (2) can also be rewritten as  

$$y_k = x_k - \sum_{i=1}^{\min(N_a-1, N-1-k)} a'_{i} \; y_{k+i},$$ (4)

which provides an exact recursive inverse to adjoint operator, ${\bf A}'$, that can be computed starting from y_N-1=x_N-1, and decrementing k.