Basic operators and adjoints |
Ordinary differential equations often lead us to the backsolving operator.
For example, the damped harmonic oscillator leads to
a special case of equation
(22),
where
.
There is a huge literature on finite-difference solutions
of ordinary differential equations
that lead to equations of this type.
Rather than derive such an equation on the basis
of many possible physical arrangements,
we can begin from the filter transformation in equation (4),
but put the top square of the matrix on the other side of the equation
so our transformation can be called one of inversion or backsubstitution.
To link up with applications in later chapters,
I specialize to put 1s on the main diagonal and insert some bands of zeros.
A typical row in equation (22) says:
In the same way that equation (4) can be interpreted as , equation (22) can be interpreted as , which amounts to . Thus, convolution is amounts to polynomial multiplication while the backsubstitution we are doing here is called ``deconvolution,'' and it amounts 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 ; whereas, anticausal operators are associated with upper triangular matrices and negative powers of . A transformation like equation (22) but with the transposed matrix would require us to run the recursive solution the opposite direction in time, as we did with leaky integration.
A module to backsolve equation (22) is recfilt.
for (ix=0; ix < nx; ix++) { tt[ix] = 0.; } if (adj) { for (ix = nx-1; ix >= 0; ix-) { tt[ix] = yy[ix]; for (ia = 0; ia < SF_MIN(na,ny-ix-1); ia++) { iy = ix + ia + 1; tt[ix] -= aa[ia] * tt[iy]; } } for (ix=0; ix < nx; ix++) { xx[ix] += tt[ix]; } } else { for (iy = 0; iy < ny; iy++) { tt[iy] = xx[iy]; for (ia = 0; ia < SF_MIN(na,iy); ia++) { ix = iy - ia - 1; tt[iy] -= aa[ia] * tt[ix]; } } for (iy=0; iy < ny; iy++) { yy[iy] += tt[iy]; } } |
We may wonder why the adjoint of
actually is
.
With the well-known fact that
the inverse of a transpose is the transpose of the inverse
we have:
(25) | |||
(26) | |||
(27) | |||
(28) |
Basic operators and adjoints |