This paragraph describes only the false starts I made toward the solution.
It seems that the filter should be something like (1,-2,1),
because that filter interpolates on straight lines
and is not far from the feedback coefficients of a damped sinusoid.
(See equation ().)
So I thought about different ways to force the solution
to move in that direction.
Traditional **linear inverse theory** offers several suggestions;
I puzzled over these before I found the right one.
First, I added the obvious stabilizations
and
,but they simply made the filter and the interpolated data smaller.
I thought about changing the identity matrix in to a diagonal matrix
or
.Using , I could penalize the filter
at even-valued lags, hoping that it
would become nonzero at odd lags, but that did not work.
Then I thought of using
,,
, and
,which would allow freedom of choice of the **mean** and **variance**
of the unknowns.
In that case, I must supply
the mean and variance, however, and doing that seems
as hard as solving the problem itself.
Suddenly,
I realized the answer.
It is simpler than anything above,
yet formally it seems more complicated,
because a full inverse **covariance matrix**
of the unknown filter is implicitly supplied.

I found a promising new approach in the **stabilize**d minimization

(4) |

The next question is,
what are the appropriate definitions for *P _{0}* and

To begin with, I think of *P _{0}* as a low-pass filter,
indicating that data is normally oversampled.
Likewise,

Returning to the one-dimensional signal-interlacing problem,
I take *A _{0}*=0 and choose

Figure 12

To understand the coding implied by
the optimization (4),
it is helpful to write the linearized regression.
The training signal *P _{0}* enters as a matrix of shifted columns
of the training signal, say ;and the high-pass filter

(5) |

The top row restates equation (3).
The middle row says that
,and the bottom row says that
.A program that does the job
is `misfip()` .
It closely resembles `missif()` .

# MISFIP -- find MISsing peF and Input data on 1-axis using Prior data. # subroutine misfip( nt,tt, na,aa, np,pp,known, niter) integer nt, na, ip,np, npa, nta, nx,nr, iter,niter, ax, px, qr, tr real pp(np), known(np), aa(na) # same as in missif() real tt(nt) # input: prior training data set. temporary real x(np+na), g(np+na), s(np+na) temporary real rr(np+na-1 +na+nt-1), gg(np+na-1 +na+nt-1), ss(np+na-1 +na+nt-1) npa= np+na-1; nta= nt+na-1 # lengths of outputs of filtering nx = np+na; nr= npa+nta # length of unknowns and residuals px=1; qr=1; ax=1+np; tr=1+npa # pointers call zero( na, aa); aa(1) = 1. call copy( np, pp, x(px)) call copy( na, aa, x(ax)) do iter= 0, niter { call contran( 0, 0, na,aa, np, pp, rr(qr)) call contran( 0, 0, na,aa, nt, tt, rr(tr)) # extend rr with train call scaleit( -1., nr, rr ) call contran( 1, 0, na,aa, np, g(px), rr(qr)) call contran( 1, 0, np,pp, na, g(ax), rr(qr)) call contran( 1, 1, nt,tt, na, g(ax), rr(tr)) do ip= 1, np { if( known(ip) != 0) { g( ip+(px-1)) = 0. } } g( 1 +(ax-1)) = 0. call contran( 0, 0, na,aa, np, g(px), gg(qr)) call contran( 0, 1, np,pp, na, g(ax), gg(qr)) call contran( 0, 0, nt,tt, na, g(ax), gg(tr)) call cgstep( iter, nx, x, g, s, nr, rr, gg, ss) call copy( np, x(px), pp) call copy( na, x(ax), aa) } return; endThe new computations are the lines containing the training data

You might wonder why we need
another program when we could use the old program
and simply append the training data to the observed data.
We will encounter some applications where
the old program will not be adequate.
These involve the boundaries of the data.
(Recall that, in chapter , when seismic events changed their dip,
we used a two-dimensional wave-killing operator
and were careful not to convolve the operator over the edges.)
Imagine a dataset that changes with time (or space).
Then *P _{0}* might not be training data,
but data from a large interval,
while

10/21/1998