\section{SETTING UP 2-D FT}\begin{comment}The program {\tt fth()} is set up so that the vectors transformedcan be either rows or columns of a two-dimensional array.In any computer language there is a way to extracta vector (column or row) from a matrix.In some languages the vector can be processed directly without extraction.To see how this works in \bx{Fortran},recall a matrix allocated as{\tt (n1,n2)}can be subscripted as a matrix{\tt (i1,i2)}or as a long vector{\tt (i1 + n1*(i2-1),1)},and {\tt call sub(x(i1,i2))} passes the subroutinea pointer to the {\tt (i1,i2)} element.To transform an entire axis, the subroutines{\tt ft1axis()} and{\tt ft2axis()} are given.For a two-dimensional FT,we simply call both{\tt ft1axis()} and{\tt ft2axis()} in either order.%\progdex{ft1axis}{FT 1-axis}\progdex{ft2axis}{FT 2-axis}\end{comment}\par%I confess that there are faster ways to do things%than those I have shown you above.%When we are doing many FTs, for example,%the overhead calculations done the first time%should be saved for use on subsequent FTs,%as in the subroutine {\tt rocca()} included in IEI.%Further, manufacturers of computers for heavy numerical use%generally design special FT codes for their architecture.%Although the basic fast FT used here%ingeniously stores its output on top of its input,%that feature is not compatible with vectorizing architectures.\subsection{Basics of two-dimensional Fourier transform}\inputdir{plane4}\parLet us review some basic facts about\bxbx{two-dimensional Fourier transform}    {Fourier transform!two-dimensional}.A two-dimensional function is represented in a computer as numerical values in a matrix, whereasa one-dimensional Fourier transform in a computeris an operation on a vector.A 2-D Fourier transform can be computed by a sequenceof 1-D Fourier transforms.We can first transform each column vector of the matrix and theneach row vector of the matrix.Alternately, we can first do the rows and later do the columns.This is diagrammed as follows:$$\def\normalbaselines{\baselineskip =\normalbaselineskip                       \lineskip10pt \lineskiplimit10pt}  \matrix{   p(t,\ x) & \smash{\mathop{\longleftrightarrow}}     & P(t, \ k_x)\cr   \lineskip10pt   \Big \updownarrow &               &\Big\updownarrow\cr   \lineskip10pt   P(\omega, \ x)& \smash{\mathop{\longleftrightarrow}}     & P(\omega,\ k_x)\cr  }$$\parThe diagram has thenotational problem that we cannot maintainthe usual convention of using a lower-case letter for the domainof physical space and an upper-case letter for the Fourier domain,because that convention cannot includethe mixed objects  $ P(t ,  k_x )$  and  $ P ( \omega , x )$.Rather than invent some new notation, it seems best to let the readerrely on the context:the arguments of the function must help name the function.\parAn example of\bxbx{two-dimensional Fourier transform}    {Fourier transform!two-dimensional}son typical deep-oceandata is shown in Figure~\ref{fig:plane4}.%\plot{plane4}{width=6in,height=6in}{  A deep-marine dataset  $p(t , x)$  from Alaska  (U.S. Geological Survey)  and the {\em  real} part of various  Fourier transforms of it.  Because of the long traveltime through the water,  the time axis does not begin at  $t=0$.}%\newslideIn the deep ocean, sediments are fine-grained and deposit slowly inflat, regular, horizontal beds.The lack of permeable rocks such as sandstone severely reducesthe potential for petroleum production from the deep ocean.The fine-grained shales overlay irregular, igneous, \bx{basement rock}s.In the plot of  $P(t ,  k_x )$,  the lateral continuity of thesediments is shown by the strong spectrum at low  $k_x$.The igneous rocks show a  $k_x$  spectrumextending to such large  $k_x$  that the deep data may be somewhat\bxbx{spatially aliased}{spatial alias}(sampled too coarsely).The plot of  $P( \omega , x)$  shows that the data contains nolow-frequency energy.% At large  $\omega$  the energy is not dropping off% as fast as one might like, which indicates temporal frequency aliasing.% This aliasing is also apparent in the plot of  $p(t , x)$,% showing up as the steplike shape of the sea-floor arrival.The dip of the sea floor shows up in  $( \omega , k_x )$-space asthe energy crossing the origin at an angle.\parAltogether, the\bxbx{two-dimensional Fourier transform}    {Fourier transform!two-dimensional}of a collection of seismogramsinvolves only twice as much computation as the one-dimensionalFourier transform of each seismogram.This is lucky.Let us write some equations to establish that the asserted proceduredoes indeed do a 2-D Fourier transform.Say first that any function of  $x$  and  $t$  maybe expressed as a superposition of sinusoidal functions:\begin{equation}p(t , x) \ \ =\ \  \int \int \  e^{ -i \omega t + i k_x x } \ P( \omega , k_x ) \  d \omega \  d k_x \label{eqn:2.9}\end{equation}The double integration can be nested to showthat the temporal transforms are done first (inside):\begin{eqnarray*}p(t , x) \ \ &=&\ \  \int \  e^{ i \,k_x x } \  \left[ \int \ e^{ -i \omega t } \  P ( \omega , k_x ) \  d \omega \right]\,dk_x\\ \ \ &=&\ \ \int \  e^{ i \,k_x x } \  P(t ,  k_x ) \  dk_x\end{eqnarray*}The quantity in brackets is a Fourier transform over  $\omega$  donefor each and every  $ k_x $.Alternately, the nesting could be done with the $ k_x $-integralon the inside.That would imply rows first instead of columns (or vice versa).It is the separability of  $\exp (-i \omega t \,+\, i\, k_x x ) $  into aproduct of exponentials that makes the computation easy and cheap.\subsection{Signs in Fourier transforms}\par\label{'sign convention'}In Fourier transforming  $t$-, $x$-, and $z$-coordinates,we must choosea sign convention for each coordinate.Of the two alternative \bx{sign convention}s,electrical engineers have chosen one and physicists another.While both have good reasons for their choices,our circumstances more closely resemble those of physicists,so we will use their convention.For the {\em  inverse} Fourier transform, our choice is\begin{equation}p(t,x,z) \eq \int \int \int \ e^{ -i \omega t \,+\, ik_x x \,+\, ik_z z}\ P ( \omega , k_x , k_z ) \ d  \omega \,  dk_x \, dk_z\label{eqn:6.1}\end{equation}For the {\em  forward} Fourier transform,the space variables carry a {\em  negative} sign, andtime carries a {\em  positive} sign.\parLet us see the reasons why electrical engineers have madethe opposite choice,and why we go with the physicists.Essentially, engineers transform only the time axis,whereas physicists transform both time and space axes.Both are simplifying their lives by their choiceof sign convention,but physicists complicate their time axisin order to simplifytheir many space axes.The engineering choice minimizes the number of minus signsassociated with the time axis, because for engineers,$d/dt$ is associated with $i\omega$ instead of,as is the case for us and for physicists,with $-i\omega$.We confirm this with equation~(\ref{eqn:6.1}).Physicists and geophysicists deal with many more independent variablesthan time.Besides the obvious three space axesare their mutual combinations, such as midpoint and offset.\parYou might ask,why not make {\em  all} the signs positive in equation~(\ref{eqn:6.1})?The reason is that in that casewaves would not move in a positive direction along the space axes.This would be especially unnatural when the space axis was a radius.Atoms, like geophysical sources,always radiate from a point to infinity,not the other way around.Thus, in equation (\ref{eqn:6.1}) the sign of the spatial frequenciesmust be opposite that of the temporal frequency.\parThe only good reason I know to choose theengineering conventionis that we might compute with anarray processor built and microcoded by engineers.Conflict of sign convention is not a problemfor the programs that transform complex-valued time functions to complex-valued frequency functions,because there the sign convention is under the user's control.But sign conflict does make a difference when we use any programthat converts real-time functions to complex frequency functions.The way to live in both worlds is to imaginethat the frequencies produced by such a programdo not range from  $0$  to  $+ \pi$  asthe program description says, but from  $0$  to  $- \pi$.Alternately, we could always take the complex conjugate of the transform,which would swap the sign of the $\omega$-axis.\subsection{Simple examples of 2-D FT}\inputdir{ft2d}\parAn example of a \bxbx{two-dimensional Fourier transform}    {Fourier transform!two-dimensional}of a pulse is shown in Figure~\ref{fig:ft2dofpulse}.%\plot{ft2dofpulse}{width=6in}{  A broadened pulse (left) and the real part of its FT (right).}%\newslideNotice the location of the pulse.It is closer to the time axis than the space axis.This will affect the real part of the FT in a certain way(see exercises).Notice the broadening of the pulse.It was an impulse smoothed over time (vertically) by convolutionwith (1,1) and over space (horizontally) with (1,4,6,4,1).This will affect the real part of the FT in another way.\parAnother example of a two-dimensional Fourier transformis given in Figure~\ref{fig:airwave}.This example simulates an impulsive air wave originating at a pointon the $x$-axis.We see a wave propagating in each directionfrom the location of the source of the wave.In Fourier space there are also two lines, one for each wave.Notice that there are other lines which do not go through the origin;these lines are called ``\bx{spatial alias}es.''Each actually goes through the originof another square plane that is not shown,but which we can imagine alongside the one shown.These other planes are periodic replicas of the one shown.%\plot{airwave}{width=6in}{  A simulated air wave (left) and the amplitude of its FT (right).}%\newslide%\subsection{Dot-product test of real to complex transforms}%My previous book, PVI,%makes the error of mixing up real number pairs with complex numbers.%The difference became a roadblock%when Dave Nichols and I were attempting to write a universal program in C++%for performing dot product tests.%(In Fortran one generally writes a different test program for each operator).%\par%Defining a complex number as a real number pair%works as you expect for addition, subtraction, and length $\overline{z} z$,%but products of complex numbers do not behave like%those of real number pairs.%To illustrate this, let us represent the operator $\Re$ that%takes the real part of a complex number by a simple%$1\times 2$ matrix multiplication%and then examine the dot-product test for this matrix.%\begin{equation}%\left[%   \begin{array}{c}%     \tilde r%     \end{array} \right]%\eq %\left[%   \begin{array}{cc}%     1 & 0%     \end{array} \right]%\left[%   \begin{array}{c}%     x \\%     y \end{array} \right]%\end{equation}%The adjoint operation is%\begin{equation}%\left[%   \begin{array}{c}%     \tilde x \\%     \tilde y \end{array} \right]%\eq%\left[%   \begin{array}{c}%     1 \\%     0 \end{array} \right]%\ %\left[%   \begin{array}{c}%      r%     \end{array} \right]%\end{equation}%%The dot product test assuming that complex numbers are%handled as real number pairs is%\begin{equation}%\tilde r r \eq%\left[%   \begin{array}{cc}%     \tilde x & \tilde y %     \end{array} \right]%\left[%   \begin{array}{c}%     x \\%     y \end{array} \right]%\eq%\tilde x x + \tilde y y%\EQNLABEL{dotpair}%\end{equation}%%%The dot-product test with complex numbers is%%\begin{equation}%%\overline{\tilde z} z \eq%%(\tilde x - i\tilde y)(x+iy)%%\eq%%\tilde x x +%%\tilde y y + i(%%\tilde x y -%%\tilde y x  )%%\EQNLABEL{dotcomplex}%%\end{equation}%%It seems there are two ways%to administer the dot-product test%to a program that inputs real values and outputs complex ones.%The first is to express inputs and outputs as number pairs%and interpret the test as equation~\EQN{dotpair}.%The second is to%%use equation~\EQN{dotcomplex}%express all inputs and outputs as complex numbers.%Then the dot-product test compares two complex numbers, say%\begin{equation}%\tilde z_1' z_1 \eq