paper.tex

来自「国外免费地震资料处理软件包」· TEX 代码 · 共 1,940 行 · 第 1/5 页

         .& .& .& .& .& 0
        \end{array} \right] \ 
\left[ \begin{array}{c}
        x_1 \\
        x_2 \\
        x_3 \\
        x_4 \\
        x_5 \\
        x_6
        \end{array} \right]
 \label{eqn:ruff1}
\end{equation}
\par
The \bx{filter impulse response} is seen in any column
in the middle of the matrix, namely $(1,-1)$. 
In the transposed matrix,
the filter-impulse response
is time-reversed to $(-1,1)$.
So, mathematically,
we can say that the adjoint of the time derivative operation
is the negative time derivative.
This corresponds also to the fact that
the complex conjugate of $-i\omega$ is $i\omega$.
We can also speak of the adjoint of the boundary conditions:
we might say that the adjoint of ``no boundary condition''
is a ``specified value'' boundary condition.
The last row in equation~(\ref{eqn:ruff1}) is optional.
%and depends not on the code shown, but the code that invokes it.
It may seem unnatural to append a null row, but it can be a small
convenience (when plotting) to have the input and output be the same size.
\par

Equation~(\ref{eqn:ruff1}) is implemented by the code
in module \texttt{igrad1}
which does the operator itself (the forward operator)
and its adjoint.
\opdex{igrad1}{first difference}{25}{40}{filt/proc}

\par \noindent
The adjoint code may seem strange.
It might seem more natural to code the adjoint
to be the negative of the operator itself
and then make the special adjustments for the boundaries.
The code given, however, is correct and requires no adjustments
at the ends.
To see why, notice for each value of \texttt{i},
the operator itself handles one row of
equation~(\ref{eqn:ruff1})
while for each \texttt{i}
the adjoint handles one column.
That's why coding the adjoint in this way
does not require any special work on the ends.
The present method of coding reminds us that
the adjoint of a sum of $N$ terms is a collection of $N$ assignments.

%A complicated way to think about the adjoint of equation
%(\ref{eqn:ruff1}) is to note that it is the negative of the derivative
%and that something must be done about the ends.
%A simpler way to think about it
%is to apply the idea that the adjoint of a sum of $N$ terms
%is a collection of $N$ assignments.

\begin{comment}
\par
The Ratfor90 dialect of Fortran
allows us to write the inner code of
the \texttt{igrad1} module more simply and symmetrically
using the syntax of C, C++, and Java
where
expressions like {\tt a=a+b} can be written more tersely as {\tt a+=b}.
With this, the heart of module \texttt{igrad1} becomes
\begin{verbatim}
if( adj) {   xx(i+1) += yy(i)
             xx(i)   -= yy(i)
           }
else {       yy(i)   += xx(i+1)
             yy(i)   -= xx(i)
           }
\end{verbatim}
where we see that each component of the matrix is handled both by the
operator and the adjoint.  Think about the forward operator
``pulling'' a sum into {\tt yy(i)}, and think about the adjoint
operator ``pushing'' or ``spraying'' the impulse {\tt yy(i)} back into
{\tt xx()}.
\end{comment}

%Something odd happens at the ends of the adjoint only if we take the
%perspective that the adjoint should have been computed one component
%at a time instead of all together.
%By not taking that view, we avoid that confusion.

\inputdir{igrad1}

\par
Figure \ref{fig:stangrad} illustrates the use of module
\texttt{igrad1} for each north-south line of a topographic map.  We
observe that the gradient gives an impression of illumination from a
low sun angle.

\plot{stangrad}{width=6in,height=8.4in}{ Topography
  near Stanford (top) southward slope (bottom).  }

To apply \texttt{igrad1} along the 1-axis for each point on the 2-axis
of a two-dimensional map, we use the loop
\begin{verbatim}
for (iy=0; iy < ny; iy++) 
     igrad1_lop(adj, add, nx, nx, map[iy], ruf[iy])
\end{verbatim}
\begin{comment}
On the other hand, to see the east-west gradient, we use the loop
\begin{verbatim}
do ix=1,nx 
      stat = igrad1_lop( adj, add, map(ix,:), ruf(ix,:))
\end{verbatim}
\end{comment}

%\par
%The C/C++/Java family of computer languages allows
%an abbreviation that is handy for all operator routines;
%it's so handy that Bob Clapp has installed it in our version of Ratfor90.
%The essential feature of operator calculations is summing
%and spraying and both are rendered in Fortran as ``A=A+B''.
%In the C/C++/Java languages A=A+B may be abbreviated A+=B;
%for example,
%{\tt xx(i+1) = xx(i+1) + yy(i)} reduces to {\tt xx(i+1) += yy(i)}.
%The more complicated the subscripting,
%the more desirable the abbreviation.
%Also, the abbreviation makes more transparent the nature
%of operator-adjoint pairs.


\subsection{Transient convolution}
\sx{convolution}
\sx{convolution ! transient}

\par
The next operator we examine is convolution.
It arises in many applications; and it could be derived in many ways.
A basic derivation is from the multiplication of two polynomials, say
$X(Z) = x_1 + x_2 Z + x_3 Z^2 + x_4 Z^3 + x_5 Z^4 + x_6 Z^5$ times
$B(Z) = b_1 + b_2 Z + b_3 Z^2 + b_4 Z^3$.\footnote{
	This book is more involved with matrices than with Fourier analysis.
	If it were more Fourier analysis we would choose notation
	to begin subscripts from zero like this:
	$B(Z) = b_0 + b_1 Z + b_2 Z^2 + b_3 Z^3$.}
Identifying the $k$-th power of $Z$ in the product
$Y(Z)=B(Z)X(Z)$ gives the $k$-th row of the convolution transformation
(\ref{eqn:contran1}).
%Matrix multiplication and transpose multiplication still
%fit easily in the same computational framework
%when the matrix has a special form, such as

\begin{equation}
\bold y \eq
\left[ 
\begin{array}{c}
  y_1 \\ 
  y_2 \\ 
  y_3 \\ 
  y_4 \\ 
  y_5 \\ 
  y_6 \\ 
  y_7 \\ 
  y_8
  \end{array} \right] 
\eq
\left[ 
\begin{array}{cccccc}
  b_1 & 0   & 0    & 0   & 0   & 0   \\
  b_2 & b_1 & 0    & 0   & 0   & 0   \\
  b_3 & b_2 & b_1  & 0   & 0   & 0   \\
  0   & b_3 & b_2  & b_1 & 0   & 0   \\
  0   & 0   & b_3  & b_2 & b_1 & 0   \\
  0   & 0   & 0    & b_3 & b_2 & b_1 \\
  0   & 0   & 0    & 0   & b_3 & b_2 \\
  0   & 0   & 0    & 0   & 0   & b_3 
  \end{array} \right] 
\; \left[ 
\begin{array}{c}
  x_1 \\ 
  x_2 \\ 
  x_3 \\ 
  x_4 \\ 
  x_5 \\ 
  x_6
  \end{array} \right]
\eq \bold B \bold x
\label{eqn:contran1}
\end{equation}
Notice that columns of
equation (\ref{eqn:contran1})
all contain the same signal,
but with different shifts.
This signal is called the filter's impulse response.


\par
Equation~(\ref{eqn:contran1}) could be rewritten as
\begin{equation}
\bold y \eq
\left[ 
\begin{array}{c}
  y_1 \\ 
  y_2 \\ 
  y_3 \\ 
  y_4 \\ 
  y_5 \\ 
  y_6 \\ 
  y_7 \\ 
  y_8
  \end{array} \right] 
\eq
\left[ 
\begin{array}{ccc}
  x_1 & 0   & 0    \\
  x_2 & x_1 & 0    \\
  x_3 & x_2 & x_1  \\
  x_4 & x_3 & x_2  \\
  x_5 & x_4 & x_3  \\
  x_6 & x_5 & x_4  \\
  0   & x_6 & x_5  \\
  0   & 0   & x_6
  \end{array} \right] 
\; \left[ 
\begin{array}{c}
  b_1 \\ 
  b_2 \\ 
  b_3 \end{array} \right]
  \eq \bold X \bold b
\label{eqn:contran2}
\end{equation}
In applications we can choose between
$\bold y = \bold X \bold b$
and
$\bold y = \bold B \bold x$.
In one case
the output $\bold y$
is dual to the filter $\bold b$,
and in the other case
the output $\bold y$
is dual to the input $\bold x$.
Sometimes we must solve
for $\bold b$ and sometimes for $\bold x$;
so sometimes we use equation~(\ref{eqn:contran2}) and 
sometimes (\ref{eqn:contran1}).
Such solutions begin from the adjoints.
The adjoint of (\ref{eqn:contran1}) is
\begin{equation}
\left[ 
\begin{array}{c}
\hat  x_1 \\ 
\hat  x_2 \\ 
\hat  x_3 \\ 
\hat  x_4 \\ 
\hat  x_5 \\ 
\hat  x_6
  \end{array} \right] 
\eq
\left[ 
\begin{array}{cccccccc}
  b_1 & b_2 & b_3 & 0   & 0   & 0   & 0   & 0  \\
  0   & b_1 & b_2 & b_3 & 0   & 0   & 0   & 0  \\
  0   & 0   & b_1 & b_2 & b_3 & 0   & 0   & 0  \\
  0   & 0   & 0   & b_1 & b_2 & b_3 & 0   & 0  \\
  0   & 0   & 0   & 0   & b_1 & b_2 & b_3 & 0  \\
  0   & 0   & 0   & 0   & 0   & b_1 & b_2 & b_3 
  \end{array} \right] 
\; \left[ 
\begin{array}{c}
  y_1 \\ 
  y_2 \\ 
  y_3 \\ 
  y_4 \\ 
  y_5 \\ 
  y_6 \\ 
  y_7 \\ 
  y_8 \end{array} \right]
\label{eqn:ajtcontran1}
\end{equation}
The adjoint {\em  \bx{crosscorrelate}s} with the filter
instead of convolving with it (because the filter is backwards).
Notice that each row in
equation (\ref{eqn:ajtcontran1}) contains all the filter coefficients
and there are no rows where the filter somehow uses zero values
off the ends of the data as we saw earlier.
In some applications it is important not to assume zero values
beyond the interval where inputs are given.
\par
The adjoint of (\ref{eqn:contran2}) crosscorrelates a fixed portion
of filter input across a variable portion of filter output.
\begin{equation}
\left[ 
\begin{array}{c}
\hat  b_1 \\ 
\hat  b_2 \\ 
\hat  b_3
  \end{array} \right] 
\eq
\left[ 
\begin{array}{cccccccc}
  x_1 & x_2 & x_3 & x_4 & x_5 & x_6 & 0   & 0   \\
  0   & x_1 & x_2 & x_3 & x_4 & x_5 & x_6 & 0   \\
  0   & 0   & x_1 & x_2 & x_3 & x_4 & x_5 & x_6   
  \end{array} \right] 
\; \left[ 
\begin{array}{c}
  y_1 \\ 
  y_2 \\ 
  y_3 \\ 
  y_4 \\ 
  y_5 \\ 
  y_6 \\ 
  y_7 \\ 
  y_8
  \end{array} \right]
\label{eqn:ajtcontran2}
\end{equation}
\par
Module \texttt{tcai1}
is used for
$\bold y = \bold B \bold x$
and module~\texttt{tcaf1}
is used for
$\bold y = \bold X \bold b$.
\opdex{tcai1}{transient convolution}{45}{50}{user/gee}
\opdex{tcaf1}{transient convolution}{45}{50}{user/gee}

\par
The polynomials 
$X(Z)$,
$B(Z)$, and
$Y(Z)$
are called $Z$ transforms.
An important fact in real life
(but not important here)
is that the $Z$ transforms are
Fourier transforms in disguise.
Each polynomial is a sum of terms
and the sum
amounts to a Fourier sum when we take $Z=e^{i\omega\Delta t}$.
The very expression
$Y(Z)=B(Z)X(Z)$ says that a product in the frequency domain
($Z$ has a numerical value) is a convolution in the time domain
(that's how we multipy polynomials, convolve their coefficients).

%The length of the output of \bx{transient convolution}
%\sx{convolution ! transient}
%is one less than the length of the filter plus the length of the input.
%Notice, however,
%that in these modules,
%the input length plus the filter length equals the output length,
%because these routines erase the extra memory cell.
%So you must be sure you have allocated that memory cell.
%I erase the extra cell to reduce clutter
%in the many 2-D and 3-D modules
%that build upon transient convolution.
%It is important that the extra memory cell be erased and not simply ignored,
%because otherwise it may contain an impossible value
%that would break many modules such as plot routines and inversion routines.
 
\subsection{Internal convolution}
%\begin{notforlecture}
\sx{convolution ! internal}
Convolution is the computational equivalent of
ordinary linear differential operators (with constant coefficients).
Applications are vast,
and end effects are important.
Another choice of data handling at ends
is that zero data not be assumed
beyond the interval where the data is given.
This is important in data where the crosscorrelation changes with time.
Then it is sometimes handled as constant in short time windows.
Care must be taken that zero signal values not be presumed
off the ends of those short time windows;
otherwise,
the many ends of the many short segments
can overwhelm the results.
\par
In the sets (\ref{eqn:contran1}) and (\ref{eqn:contran2}),
the top two equations explicitly assume that the input data vanishes
before the interval on which it is given, and likewise at the bottom.
Abandoning the top two and bottom two equations in (\ref{eqn:contran2})
we get:
\begin{equation}
\left[ 
\begin{array}{c}
  y_3 \\ 
  y_4 \\