paper.tex

国外免费地震资料处理软件包

\def\CAKEDIR{.}

\title{Waves and Fourier sums}
\author{Jon Claerbout}
\maketitle
\label{paper:ft1}
\long\def\HideThis#1{}

%{\em \today . This chapter is owned by JFC.}

\todo{jarring jump from recommending odd FFT, then using even}

\sx{Fourier transform!discrete}
\par
An important concept in wave imaging is the extrapolation
of a wavefield from one depth $z$ to another.
Fourier transforms are an essential basic tool.
There are many books and chapters of books on the
{\em  theory}
of Fourier transformation.
The first half of this chapter is
an introduction to {\em  practice} with Fourier sums.
It assumes you already know something of the theory
and takes you through the theory rather quickly
emphasizing practice by examining examples,
and by performing two-dimensional
Fourier transformation of data and interpreting the result.
For a somewhat more theoretical background,
I suggest my previous book PVI at
http://sepwww.stanford.edu/sep/prof/.


\par
The second half of this chapter
uses Fourier transformation
to explain the Hankel waveform we observed
in chapter~\ref{vela/paper:vela} and chapter~\ref{krch/paper:krch}.
Interestingly,
it is the Fourier transform of $\sqrt{-i\omega}$,
which is half the derivative operator.

\section{FOURIER TRANSFORM}
We first examine the two ways to visualize polynomial multiplication.
The two ways lead us to the most basic principle of Fourier analysis
that
\par
\boxit{ A product in the Fourier domain
	is a convolution in the physical domain}
\par\noindent
Look what happens to the coefficients when we multiply polynomials.
\begin{eqnarray}
X(Z)\, B(Z) &\eq & Y(Z) \label{eqn:1-1-5} \\
(x_0 + x_1 Z + x_2 Z^2 + \cdots )\, (b_0 + b_1 Z + b_2 Z^2) &\eq &
y_0 + y_1 Z + y_2 Z^2 + \cdots  \label{eqn:1-1-6}
\end{eqnarray}
Identifying coefficients of successive powers of $Z$, we get
\begin{eqnarray}
y_0 &\eq & x_0 b_0 \nonumber \\
y_1 &\eq & x_1 b_0  + x_0 b_1 \nonumber \\
y_2 &\eq & x_2 b_0 + x_1 b_1 + x_0 b_2 \label{eqn:1-1-7} \\
y_3 &\eq & x_3 b_0 + x_2 b_1 + x_1 b_2 \nonumber \\
y_4 &\eq & x_4 b_0 + x_3 b_1 + x_2 b_2 \nonumber \\
    &\eq & \cdots\cdots\cdots\cdots\cdots\cdots  \nonumber
\end{eqnarray}
In matrix form this looks like
\begin{equation}
\left[ 
\begin{array}{c}
  y_0 \\ 
  y_1 \\ 
  y_2 \\ 
  y_3 \\ 
  y_4 \\ 
  y_5 \\ 
  y_6
  \end{array} \right] 
\eq
\left[ 
\begin{array}{ccc}
  x_0 & 0   & 0    \\
  x_1 & x_0 & 0    \\
  x_2 & x_1 & x_0  \\
  x_3 & x_2 & x_1  \\
  x_4 & x_3 & x_2  \\
  0   & x_4 & x_3  \\
  0   & 0   & x_4
  \end{array} \right] 
\; \left[ 
\begin{array}{c}
  b_0 \\ 
  b_1 \\ 
  b_2 \end{array} \right]
\label{eqn:contran2}
\end{equation}
The following equation, called the
``convolution equation,''
carries the spirit of the group shown in (\ref{eqn:1-1-7})
\begin{equation}
y_k \eq  \sum_{i = 0} x_{k - i} b_i
\label{eqn:conv}
\end{equation}

\par
The second way to visualize polynomial multiplication is simpler.
Above we did not think of $Z$ as a numerical value.
Instead we thought of it as ``a unit delay operator''.
Now we think of the product $X(Z) B(Z) = Y(Z)$ numerically.
For all possible numerical values of $Z$,
each value $Y$ is determined
from the product of the two numbers $X$ and $B$.
Instead of considering all possible numerical values
we limit ourselves to all values of unit magnitude
$Z=e^{i\omega}$ for all real values of $\omega$.
This is Fourier analysis, a topic we consider next.

\subsection{FT as an invertible matrix}
\par
A \bx{Fourier sum} may be written
\begin{equation}
B(\omega) \eq  \sum_t \ b_t \ e^{i\omega t}
	  \eq  \sum_t \ b_t \ Z^t
\end{equation}
where the complex value $Z$
is related to the real frequency
$\omega$ by $Z=e^{i\omega}$.
This Fourier sum is a way of building
a continuous function of $\omega$
from discrete signal values $b_t$ in the time domain.
Here we specify both time and frequency domains by a set of points.
Begin with an example of a signal
that is nonzero at four successive instants,
$( b_0, b_1, b_2, b_3)$.
The transform is
\begin{equation}
B(\omega) \eq  b_0 + b_1 Z + b_2 Z^2 + b_3 Z^3
\end{equation}
The evaluation of this polynomial can be organized as a matrix times a vector,
such as
\begin{equation}
  \left[ \begin{array}{c}
   B_0 \\
   B_1 \\
   B_2 \\
   B_3  \end{array} \right]
\eq   \left[ \begin{array}{cccc}
   1 & 1 & 1 & 1 \\
   1 & W & W^2 & W^3 \\
   1  & W^2 & W^4 & W^6 \\
   1  & W^3 & W^6 & W^9  \end{array} \right]  \;
  \left[ \begin{array}{c}
   b_0 \\
   b_1 \\
   b_2 \\
   b_3  \end{array} \right]   
\label{eqn:3-3}
\end{equation}
Observe that the top row of the matrix evaluates the polynomial at $Z=1$,
a point where also
$\omega=0$.
The second row evaluates $B_1=B(Z=W=e^{i\omega_0})$,
where $\omega_0$ is some base frequency.
The third row evaluates the Fourier transform for $2\omega_0$,
and the bottom row for $3\omega_0$.
The matrix could have more than four rows for more frequencies
and more columns for more time points.
I have made the matrix square in order to show you next
how we can find the inverse matrix.
The size of the matrix in (\ref{eqn:3-3}) is $N=4$.
If we choose the base frequency $\omega_0$ and hence $W$ correctly,
the inverse matrix will be
\begin{equation}
  \left[ \begin{array}{c}
   b_0 \\
   b_1 \\
   b_2 \\
   b_3  \end{array} \right]
\eq   1/N \; 
\left[ \begin{array}{cccc}
   1 & 1 & 1 & 1 \\
   1 & 1/W & 1/W^2 & 1/W^3 \\
   1  & 1/W^2 & 1/W^4 & 1/W^6 \\
   1  & 1/W^3 & 1/W^6 & 1/W^9  \end{array} \right]  \;
  \left[ \begin{array}{c}
   B_0 \\
   B_1 \\
   B_2 \\
   B_3 \end{array} \right]   
\label{eqn:3-5}
\end{equation}
Multiplying the matrix of
(\ref{eqn:3-5}) with that of 
(\ref{eqn:3-3}),
we first see that the diagonals are +1 as desired.
To have the off diagonals vanish,
we need various sums,
such as $1+W  +W^2+W^3$
and     $1+W^2+W^4+W^6$, to vanish.
Every element ($W^6$, for example,
or $1/W^9$) is a unit vector in the complex plane.
In order for the sums of the unit vectors to vanish,
we must ensure that the vectors pull symmetrically away from the origin.
A uniform distribution of directions meets this requirement.
In other words, $W$ should be the $N$-th root of unity, i.e.,
\begin{equation}
W \eq
\sqrt[N]{1} \eq
e^{2\pi i/N} 
\label{eqn:3-4}
\end{equation}
\par
The lowest frequency is zero, corresponding to the top row of
(\ref{eqn:3-3}).
The next-to-the-lowest frequency we find by setting $W$ in
(\ref{eqn:3-4}) to $Z=e^{i\omega_0}$.
So $\omega_0=2\pi /N$; and
for (\ref{eqn:3-5}) to be inverse to (\ref{eqn:3-3}),
the frequencies required are
\begin{equation}
\omega_k \eq { (0, 1, 2, \ldots , N- 1) \, 2\pi \over N} 
\label{eqn:3-2}
\end{equation}



\subsection{The Nyquist frequency}
\inputdir{matrix}
The highest frequency
in equation~(\ref{eqn:3-2}),
$\omega=2\pi (N-1)/N$,
is almost $2\pi$.
This frequency is twice as high as the Nyquist frequency $\omega=\pi$.
The \bx{Nyquist frequency}
is normally thought of as the ``highest possible'' frequency,
because $e^{i\pi t}$, for integer $t$,
plots as $(\cdots ,1,-1,1,-1,1,-1,\cdots)$.
The double Nyquist frequency function,
$e^{i2\pi t}$, for integer $t$,
plots as $(\cdots ,1,1,1,1,1,\cdots)$.
So this frequency above the highest frequency is really zero frequency!
We need to recall that $B(\omega)=B(\omega -2\pi )$.
Thus, all the frequencies near the upper end of the range equation~(\ref{eqn:3-2})
are really small negative frequencies.
Negative frequencies on the interval $(-\pi,0)$
were moved to interval $(\pi,2\pi)$
by the matrix form of Fourier summation.

\par
A picture of the Fourier transform matrix is shown in Figure \ref{fig:matrix}.
Notice the Nyquist frequency is the center row
and center column of each matrix.
\plot{matrix}{width=6in, height=6in}{
  Two different graphical means of showing the
  real and imaginary parts of the Fourier transform
  matrix of size $32\times 32$.
}
%\newslide


\subsection{Laying out a mesh}
\par
In theoretical work and in programs,
the unit delay operator definition $Z=e^{i\omega \Delta t}$
is often simplified to $\Delta t = 1$,
leaving us with $ Z = e^{i\omega} $.
How do we know whether $\omega$
is given in radians per second or radians per sample?
We may not invoke a cosine or an exponential unless
the argument has no physical dimensions.
So where we see $\omega$ without $\Delta t$,
we know it is in units of radians per sample.
\par
In practical work,
frequency is typically given in cycles/sec or \bx{Hertz}, $f$,
rather than radians, $\omega$
(where $\omega = 2\pi f$).
Here we will now switch to $f$.
We will design a computer \bx{mesh} on a physical object
(such as a waveform or a function of space).
\label{'mesh'}
We often take the mesh to begin at $t=0$,
and continue till the end $t_{\rm max}$ of the object,
so the time range $t_{\rm range} = t_{\rm max}$.
Then we decide how many points we want to use.
This will be the $N$ used in the discrete Fourier-transform program.
Dividing the range by the number gives a mesh interval $\Delta t$.
\par
Now let us see what this choice implies in the frequency domain.
We customarily take the maximum frequency to be the Nyquist,
either $f_{\rm max} = .5 /\Delta t$ Hz or
$\omega_{\rm max} = \pi /\Delta t$ radians/sec.
The frequency range $f_{\rm range}$ goes from $-.5/\Delta t$ to $.5/\Delta t$.
In summary:
\begin{itemize}
\item $\Delta t \eq t_{\rm range} / N \quad$ is time \bx{resolution}.
\item $f_{\rm range} \eq 1/\Delta t \eq N /t_{\rm range} \quad$
	is frequency range.
\item $\Delta f \eq f_{\rm range}/N \eq 1/t_{\rm range} \quad$
	is frequency \bx{resolution}.
\end{itemize}
In principle, we can always increase $N$ to refine the calculation.
Notice that increasing $N$ sharpens the time resolution
(makes $\Delta t$ smaller)
but does not sharpen the frequency resolution
$\Delta f$, which remains fixed.
Increasing $N$ increases the frequency
{\em  range,}
but not the frequency
{\em  resolution.}
\par
What if we want to increase the frequency resolution?
Then we need to choose $t_{\rm range}$ larger than required to
cover our object of interest.
Thus we either record data over a larger range,
or we assert that such measurements would be zero.
Three equations summarize the facts:
\begin{eqnarray}
\Delta t \ f_{\rm range} &=& 1			\\
\Delta f \ t_{\rm range} &=& 1			\\
\Delta f \ \Delta t        &=& {1 \over N}
\label{eqn:dtdf}
\end{eqnarray}
\par
\boxit{Increasing {\em  range} in the time domain increases
	{\em  resolution} in the frequency domain and vice versa.
	Increasing \bx{resolution} in one domain
	does not increase \bx{resolution} in the other.
	}

\section{INVERTIBLE SLOW FT PROGRAM}

\par
Typically, signals are real valued.
But the programs in this chapter are for complex-valued signals.
In order to use these programs,
copy the real-valued signal into a complex array,
where the signal goes into the real part of the complex numbers;
the imaginary parts are then automatically set to zero.

\par
There is no universally correct choice
of \bx{scale factor} in Fourier transform: