choice of scale is a matter of convenience.Equations~(\ref{eqn:3-3}) and (\ref{eqn:3-5}) mimic the $Z$-transform,so their scaling factors areconvenient for the convolution theorem---thata product in the frequency domain is a convolution in the time domain.Obviously, the scaling factors ofequations~(\ref{eqn:3-3}) and (\ref{eqn:3-5})will need to be interchanged for thecomplementary theorem that a convolution in the frequency domainis a product in the time domain.I like to use a scale factor that keeps the sums of squaresthe same in the time domain as in the frequency domain.Since I almost never need the scale factor,it simplifies life to omit it from the subroutine argument list.\begin{comment}When a scaling program is desired,we can use a simple one like \texttt{scale()} \vpageref{/prog:scale}.Complex-valued data can be scaled with {\tt scale()}merely by doubling the value of {\tt n}.\progdex{scale}{scale an array}\subsection{The simple FT code}Subroutine \texttt{simpleft()} \vpageref{/prog:simpleft} exhibits featuresfound in many physics and engineering programs.For example, the time-domain signal (which is denoted ``{\tt tt()}"),has {\tt nt} values subscripted, from {\tt tt(1)} to {\tt tt(nt)}.The first value of this signal {\tt tt(1)} is locatedin real physical time at {\tt t0}.The time interval between values is {\tt dt}.The value of {\tt tt(it)} is at time {\tt t0+(it-1)*dt}.We do not use ``{\tt if}'' as a pointer on the frequency axisbecause {\tt if} is a keyword in most programming languages.Instead, we count along the frequency axis with a variable named {\tt ie}.%\progdex{simpleft}{slow FT}%\newslideThe total frequency band is$2\pi$ radians per sample unitor $1/\Delta t$ Hz.Dividing the total interval by the number of points {\tt nf} gives $\Delta f$.We could choose the frequencies to run from 0 to $2\pi$ radians/sample.That would work well for many applications,but it would be a nuisance for applications such as differentiationin the frequency domain, which require multiplication by $-i\omega$including the \bxbx{negative frequencies}{negative frequency}as well as the positive.So it seems more natural to begin at the most negative frequencyand step forward to the most positive frequency.\end{comment}\section{CORRELATION AND SPECTRA}The spectrum of a signal is a positive function of frequencythat says how much of each tone is present.The Fourier transform of a spectrum yields an interesting functioncalled an ``\bx{autocorrelation},''which measures the similarity of a signal to itself shifted.\subsection{Spectra in terms of Z-transforms}Let us look at spectra in terms of $Z$-transforms.Let a \bx{spectrum} be denoted $S(\omega)$, where\begin{equation}S(\omega) \eq |B(\omega)|^2 \eq \overline{B(\omega)}B(\omega)\label{eqn:1-5-3}\end{equation}Expressing this in terms of a three-point $Z$-transform, we have\begin{eqnarray}S(\omega) & = & (\bar{b}_0+\bar{b}_1 e^{-i\omega} + 		\bar{b}_2 e^{-i2\omega})		(b_0 + b_1e^{i\omega} +b_2 e^{i2\omega})							\\S(Z) & = & \left(\bar{b}_0 +\frac{\bar{b}_1}{Z} +		\frac{\bar{b}_2}{Z^2} \right)		(b_0 + b_1Z + b_2Z^2 )			\\S(Z) & = & \overline{B} \left(\frac{1}{Z}\right) B(Z)\end{eqnarray}It is interesting to multiply outthe polynomial $\bar{B}(1/Z)$ with $B(Z)$ in orderto examine the coefficients of $S(Z)$:\begin{eqnarray}S(Z) &=& \frac{\bar{b}_2b_0}{Z^2}  + 	\frac{(\bar{b}_1b_0 + \bar{b}_2b_1)}{Z} +	(\bar{b}_0b_0 + \bar{b}_1b_1 + \bar{b}_2 b_2)    + (\bar{b}_0 b_1 + \bar{b}_1b_2)Z + \bar{b}_0 b_2 Z^2 \nonumber \\S(Z) &=& \frac{s_{-2}}{Z^2} + \frac{s_{-1}}{Z} + s_0 + s_1Z + s_2 Z^2\label{eqn:1-5-8}\end{eqnarray}The coefficient $s_k$ of $Z^k$ is given by\begin{equation}s_k \eq \sum_{i} \bar{b}_i b_{i+k}\label{eqn:1-5-9}\end{equation}Equation (\ref{eqn:1-5-9}) is the\bx{autocorrelation} formula.The autocorrelationvalue $s_k$ at lag $10$ is $s_{10}$.It is a measure of the similarity of $b_i$with itself shifted $10$ units in time.In the mostfrequently occurring case, $b_i$ is real; then, by inspection of (\ref{eqn:1-5-9}),we see that the autocorrelation coefficients are real,and $s_k=s_{-k}$.\parSpecializing to a real time series gives\begin{eqnarray}S(Z) & = & s_0 + s_1\left(Z+\frac{1}{Z}     \right) +		 s_2\left(Z^2 +\frac{1}{Z^2}\right) 		\\S(Z(\omega )) & = & s_0 + s_1(e^{i\omega} + e^{-i\omega}) +		s_2(e^{i2\omega} + e^{-i2\omega})		\\S(\omega ) & = & s_0 + 2s_1\cos \omega + 2s_2 \cos 2\omega	\\S(\omega ) & = & \sum_{k} s_k \cos k\omega						\label{eqn:1-5-13} \\S(\omega ) & = & \mbox{cosine transform of }\;\; s_k\label{eqn:1-5-14}\end{eqnarray}This proves a classic theorem that for real-valued signalscan be simply stated as follows:\par\boxit{ For any real signal, the cosine transform of the \bx{autocorrelation}	equals the magnitude squared of the Fourier transform.	}\subsection{Two ways to compute a spectrum}There are two computationally distinct methods by which we cancompute a spectrum: (1) compute all the $s_k$ coefficientsfrom (\ref{eqn:1-5-9}) andthen form the cosine sum (\ref{eqn:1-5-13}) for each $\omega$; andalternately, (2) evaluate $B(Z)$ for some value of Z on the unit circle,and multiply the resulting number by its complex conjugate.Repeat for many values of $Z$ on the unit circle.When there are more than about twenty lags,method (2) is cheaper, becausethe fast Fourier transform (coming up soon) can be used.\subsection{Common signals}\inputdir{autocor}Figure~\ref{fig:autocor} shows some common signals and their\bx{autocorrelation}s.Figure~\ref{fig:spectra} shows the cosine transforms ofthe autocorrelations.Cosine transform takes us from time to frequency and it also takesus from frequency to time.Thus, transform pairs in Figure~\ref{fig:spectra}are sometimes more comprehensibleif you interchange time and frequency.The various signals are given names in the figures,and a description of each follows:\plot{autocor}{width=6in}{  Common signals and one side of their autocorrelations.}% \newslide\plot{spectra}{width=6in}{  Autocorrelations and their cosine transforms,  i.e.,~the (energy) spectra of the common signals.}%\newslide\begin{description}\item [cos]The theoretical spectrum of a sinusoid is an impulse,but the sinusoid was truncated (multiplied by a rectangle function).The autocorrelation is a sinusoid under a triangle,and its spectrum is a broadened impulse(which can be shown to be a narrow sinc-squared function).\item [sinc]The \bx{sinc} function is $\sin(\omega_0 t)/(\omega_0 t)$.Its autocorrelation is another sinc function, and its spectrumis a rectangle function.Here the rectangle is  corrupted slightly by``\bx{Gibbs sidelobes},''which result from the time truncation of the original sinc.\item [wide box]A wide\bx{rectangle function}has a wide triangle function foran autocorrelation and a narrow sinc-squared spectrum.\item [narrow box]A narrow rectangle has a wide sinc-squared spectrum.\item [twin] Two pulses.\item [2 boxes]Two separated narrow boxes have the spectrum of one of them,but this spectrum is modulated (multiplied) by a sinusoidal functionof frequency, where the modulation frequency measures thetime separation of the narrow boxes.(An oscillation seen in the frequency domainis sometimes called a ``\bx{quefrency}.'')\item [comb]Fine-toothed-\bx{comb}functions are like rectangle functions with a lower Nyquist frequency.Coarse-toothed-comb functions have a spectrum which is a fine-toothed comb.\item [exponential]The autocorrelation of a transient \bx{exponential} functionis a \bx{double-sided exponential} function.\sx{exponential ! double-sided}The spectrum (energy) is a Cauchy function, $1/(\omega^2+\omega_0^2)$.The curious thing about the\bx{Cauchy function}is that the amplitude spectrumdiminishes inversely with frequency to the {\em  first} power;hence, over an infinite frequency axis, the function has infinite integral.The sharp edge at the onset of the transient exponentialhas much high-frequency energy.\item [Gauss]The autocorrelation of a \bx{Gaussian}function is another Gaussian,and the spectrum is also a Gaussian.\item [random]\bxbx{Random}{random}numbers have an autocorrelation that is an impulsesurrounded by some short grass.The spectrum is positive random numbers.\item [smoothed random]Smoothed random numbers are much the same as random numbers, but their spectral bandwidth is limited.\end{description}\section{SETTING UP THE FAST FOURIER TRANSFORM }Typically we Fourier transform seismograms about a thousand points long.Under these conditions another Fourier summation methodworks about a hundred times faster than those already given.Unfortunately, the faster Fourier transform programis not so transparently clear as the programs given earlier.Also, it is slightly less flexible.The speedup is so overwhelming, however,that the fast program is always used in routine work.\parFlexibility may be lost because the basic fast programworks with complex-valued signals,so we ordinarily convert our real signals to complex ones(by adding a zero imaginary part).More flexibility is lost because typical fast FT programs requirethe data length to be an integral power of 2.Thus geophysical datasets oftenhave zeros appended (a process called ``\bx{zero pad}ding")until the data length is a power of 2.From time to time I notice clumsy computer code written to deducea number that is a power of 2 and is larger than thelength of a dataset.An answer is found by rounding up the logarithm to base 2.%The more obvious and the quicker way to get the desired value,%however, is%with the simple Fortran function {\tt pad2()}.%%\progdex{pad2}{round up to power of two}\parHow fast is the fast Fourier transform method?The answer depends on the size of the data.The matrix times vector operation in~(\ref{eqn:3-3})requires $N^2$ multiplications and additions.That determines the speed of the slow transform.For the fast method the number of adds and multipliesis proportional to $N \,\log_2 N$.Since $2^{10}=1024$, the speed ratio is typically 1024/10 or about 100.In reality, the fast method is not quite that fast,depending on certain details of overhead and implementation.\par\begin{comment}Below is {\tt ftu()}, a version of the \bx{fast Fourier transform} program.\sx{Fourier transform!fast}There are many versions of the program---I have chosenthis one for its simplicity.Considering the complexity of the task,it is remarkable that no auxiliary memory vectors are required;indeed, the output vector lies on top of the input vector.To run this program, your first step might be to copyyour real-valued signal into a complex-valued array.Then append enough zeros to fill in the remaining space.%\progdex{ftu}{unitary FT}\parThe following two lines serve to Fourier transforma vector of 1024 complex-valued points,and then to \bx{inverse Fourier transform} them back to the original data:\sx{Fourier transform!inverse}\begin{verbatim}    call ftu(  1., 1024, cx)    call ftu( -1., 1024, cx) \end{verbatim}\par\end{comment}A reference given at the end of this chaptercontains many %other versions of the FFT program.One version transforms real-valued signals to complex-valuedfrequency functions in the interval $0 \le \omega < \pi$.Others that do not transform data on top of itselfmay be faster with specialized computer architectures.\subsection{Shifted spectrum}Subroutine \texttt{simpleft()} \vpageref{/prog:simpleft}sets things up in a convenient manner:The frequency range runs from minus Nyquistup to (but not including) plus Nyquist.Thus there is no problem with the many (but not all)user programs that have trouble with aliased frequencies.Subroutine \texttt{ftu()} \vpageref{/prog:ftu}, however has a frequency rangefrom zero to double the Nyquist.Let us therefore define a friendlier ``front end'' to {\tt ftu()}which looks more like {\tt simpleft()}.\parRecall that a time shift of $t_0$ can be implemented in the Fourier domain by multiplication by$e^{-i\omega t_0}$.Likewise, in the Fourier domain,the frequency interval, % used by subroutine~\texttt{ftu()} \vpageref{/prog:ftu},namely, $ 0 \le \omega < 2\pi$,can be shifted to the friendlier interval$ -\pi \le \omega < \pi$by a weighting function in the time domain.That weighting function is $e^{-i\omega_0 t}$where $\omega_0$ happens to be the Nyquist frequency,i.e.~alternate points on the time axis are to be multiplied by $-1$.\begin{comment}A subroutine for this purpose is {\tt fth()}.\progdex{fth}{FT, Hale style}\newslide\par\noindentTo Fourier transform a 1024-point complex vector {\tt cx(1024)}and then inverse transform it, we would write\begin{verbatim}call fth( 0, 1., 1, 1024, cx)call fth( 1, 1., 1, 1024, cx)\end{verbatim}\noindentYou might wonder about the apparent redundancy of using boththe argument {\tt adj} and the argument {\tt sign}.Having two arguments instead of one allowsus to define the {\em  forward} transform for a {\em  time} axiswith the opposite sign as             the      forward  transform for a {\em  space} axis.\parThe subroutine {\tt fth()} is somewhat cluttered bythe inclusion of afrequently needed practical feature---namely,the facility to extract vectors from a matrix,transform the vectors, and then restore them into the matrix.\end{comment}