\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}\parAn important concept in wave imaging is the extrapolationof 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 isan introduction to {\em  practice} with Fourier sums.It assumes you already know something of the theoryand takes you through the theory rather quicklyemphasizing practice by examining examples,and by performing two-dimensionalFourier transformation of data and interpreting the result.For a somewhat more theoretical background,I suggest my previous book PVI athttp://sepwww.stanford.edu/sep/prof/.\parThe second half of this chapteruses Fourier transformationto explain the Hankel waveform we observedin 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 analysisthat\par\boxit{ A product in the Fourier domain	is a convolution in the physical domain}\par\noindentLook 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}\parThe 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 determinedfrom the product of the two numbers $X$ and $B$.Instead of considering all possible numerical valueswe 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}\parA \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 buildinga 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 signalthat 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 frequenciesand more columns for more time points.I have made the matrix square in order to show you nexthow 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} \eqe^{2\pi i/N} \label{eqn:3-4}\end{equation}\parThe 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$; andfor (\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 frequencyin 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.\parA picture of the Fourier transform matrix is shown in Figure \ref{fig:matrix}.Notice the Nyquist frequency is the center rowand 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}\parIn 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 unlessthe argument has no physical dimensions.So where we see $\omega$ without $\Delta t$,we know it is in units of radians per sample.\parIn 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$.\parNow 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.}\parWhat if we want to increase the frequency resolution?Then we need to choose $t_{\rm range}$ larger than required tocover 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}\parTypically, 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.\parThere is no universally correct choiceof \bx{scale factor} in Fourier transform: