On Choosing Fourier Transforms for Practical Geoscience Applications

The variety of definitions of Fourier transforms can create confusion for practical applications. This paper examines the choice of formulas for Fourier transforms and determines the appropriate choices for geoscience applications. One set of Discrete Fourier Transforms can be defined that approximate Fourier integrals and provide transforms between sampled continuous functions in both domains. For applications involving transforms between a continuous function and a discrete function a second set of Discrete Fourier Transforms is needed with different scaling factors. Two classes of application are presented: those where either form of transforms can be used and those where it is necessary to use a particular transform to obtain the correct results.


Introduction
The Fourier transform is a widely used tool for many applications.Its value in physics is best described by Lord Kelvin (in a quote reproduced in the frontpiece of the book by Bracewell [1]): "Fourier's theorem is not only one of the most beautiful results of modern analysis, but it may be said to furnish an indispensable instrument in the treatment of nearly every recondite question in modern physics." This is all the more fulsome praise when it is remembered that Kelvin made his comments long before the introduction of the Fast Fourier Transform by Cooley and Tukey [2].Since then the efficient computation of the Fourier Transform has led to its widespread use in signal processing and incorporation into many software packages.Surprisingly, considering its extensive use, there are no standard definitions for the Fourier transform.Different software packages implement different definitions of the Fourier transform, so that a forward transform in one package can be the same as an inverse transform in another package.This can lead to considerable confusion amongst users and interferes with the efficient use of this valuable tool.
This paper reviews the different forms of the Fourier transform and identifies the particular choice for geoscience (and other) applications.First we consider Fourier's theorem and the definition of the Fourier integral and then examine the constraints imposed by using a discrete Fourier transform and how the discrete Fourier transform can approximate the Fourier integral.This in-cludes consideration of the different Fourier transform formulations for continuous and discrete functions.Several examples are given to illustrate the application of the different Fourier transform formulations.These applications are drawn from studies of electromagnetic induction in the Earth, magnetotellurics and geomagnetically induced current effects on power systems.

Fourier Series and Fourier Integral
The fundamental idea described by Fourier is that any function can be represented as a sum of cosine and sine functions.For the applications considered in this paper we will deal with functions varying in time so the cosine and sine functions represent different frequencies.In the general case, any time varying function can be represented by an infinite sum of cosine and sine waves where [1]: If we compute the partial sum of the Fourier series then we obtain a function that approximates time variation the original How good of an approximation that is achieved depends on the length of the series and the freque tent of the time domain function.
have been removed to av ncy con-For practical applications with sampled data we are dealing with a time domain function for which frequentcies above the Nyquist frequency oid aliasing problems.This "band-limited" signal is only an approximation of the complete true natural variation that occurred, but can be considered a true representation of the natural variation within the frequency range of interest for which the recordings were made.We will denote this "band-limited" signal by s(t).Now the signal s(t) can be represented by the partial Fourier series without any approximation In preparation for introducing the standard forms of the Fourier integral and Discrete Fourier Tra onential form leads to the compact ation (8).The terms involving e and n referred to as the terms for positive and gative frequencies.This may be confusing for anyone who tries to ascribe a physical meaning to term ative frequency".Instead the terms "positive frequency" and "negative frequency" should just be treated as labels for the terms involving positive and negative signs in the exponents of the exponentials making up the cosine and sine terms.
In the limit as the interval between frequencies goes to zero, the Fourier series goes to the Fourier integral.The Fourier integral has a va own by Bracewell [1].The customary formulas for the Fourier transform and the inverse Fourier transform given by Bracewell are where Bracewell's formulas have been rewritten in terms of functions f(t) in the time domain quency domain.If the integrals are written in terms of ω and F(f) in the frethere is 1 2π in the inverse transform.Thus the transform pair is: The forward and inverse transforms are essentially to do the same thing so it could be expecte ward and inverse transforms would be symmetrical.To ac d that the forhieve this symmetry people often write the transform pair as Thus we already have three definitions integral.Brigham [3] examines this diff siders the appropriate factors if the transforms are to co (12b) of the Fourier iculty and conmply with Parseval's theorem that the energy computed in the time domain is equal to the energy computed in the frequency domain and for consistency with the Laplace transform and found that these requirements were in conflict with the various definitions of the Fourier transform pair written in terms of ω.Brigham concludes that a logical way to resolve this conflict is to define the Fourier transform pair in terms of frequency f as done in Equation (10) As long as integration is with respect to f, the scale factor 1 2π never appears.Equation (10) of Fourier integrals also adopted by Bracewell [1].ety of ve a constant 1/N se transform and is the system

Discrete Fourier Transform
The Discrete Fourier Transform is defined in a vari ways.All are basically the same but ha either in the forward transform or inver also differ in the sign of the exponent of the exponential term in each transform.Proakis and Manolakis [4] use the definitions For the forward and inverse transforms.In contrast, Bracewell [1] uses and  transform, whil [5] use these same terms in th priate in the forward e Press et al. e inverse transform.This raises the questions as to which terms are the most appro ones to use.The criteria that we will use for answering these questions are: 1) to provide the best approximation to the Fourier integral and 2) for consistency with standard practice in geoscience applications.
The first choice to make concerns the sign of the exponent in the exponential term.The choices are where n is the index for samples in the time k is the index for samples in the frequency do first choice (Equation (15a)) is the simpler fo domain and main.The rm with a simple addition of the cosine and sine terms.Also ications considered here the forward discrete Fourier transform can be considered to produce one of two classes of outputs in the frequency domain: 1) samples of a continuous function, or 2) a set of discrete frequency components.Obtaining these two results may be considered to be using the discrete Fourier transform either as an approximation to the Fourier integral or as an approximation to a Fourier series.

Continuous-Time, Continuous-Frequency Functions
Considering first the discrete Fo approximation to the Fourier integral pair in Equati (10) we can write Press et al. [5] also use the additional term t  in approximating the Fourier integral, but exclude this from their formal definition of the Fourier case there is a continuous function in the time domain th g transform.In this at has been regularly sampled with a samplin interval t  and there is a continuous function in the frequency domain sampled with a sampling interval f  .Thus the time domain and frequency domain functions have equivalent forms and correspondingly there is a symmein the transforms to go between these domains.

Continuous-Time, Discrete-Frequency Functions
Let us now consider a situation where the time do try main set of transform re appropriately written as a simple function is considered to be the summation of a discrete frequencies.The inverse discrete Fourier is then mo summation of these terms without the f  term: The f  term then needs to be included in the forward transform which, with the relation 1 t f N    allows the forward Fourier transform to be written This transform pair is not symmetrical but the transforms now connect two dissimilar functions: a continuous function in the time domain and a d in the frequency domain, so it is reasonable that symmetry ction in the time domain and a di iscrete function not occur in this case.Equations (17a) and (17b) represent one of the standard discrete Fourier transform pairs commonly used.Thus this is an appropriate choice when transforming between a continuous fun screte function in the frequency domain.However, if the values in both domains are considered as samples of continuous functions then the symmetrical pair of transforms in Equation ( 16) should be used.

Practical Considerations
Bracewell starts with the standard "definitions" of the Fourier transforms, Equations (10a) and (10b), and then mmation over positive s.This does not show tput array in the frequency domain, starting with ze shows they may be written as a su and negative times and frequencie where the discrete Fourier transform comes from.Here we choose to start with the summation of positive and negative time and frequency as an approximation of the Fourier integrals, then show that because the functions are periodic we can write the summations starting at zero.
Discrete Fourier transforms convert between an array of values in the time domain and an array of values in the frequency domain.Discrete Fourier transforms produce an ou ro frequency, then positive frequencies and then negative frequencies, as shown in Figure 1.This placement of "negative frequencies" at the end of the Fourier transform array is explained in a number of books (e.g.[5]); however, the placement of "negative times" is not usually discussed.Normally it is not a consideration because we are dealing with a time series and it is easier to think of the time series starting at t = 0 as in Figure 2(b) rather than starting at 2 t T   as in Figure 2(a).However, knowledge of the "negative t" placement is shown to be necessary in the example considered later in Application 2 where we do an inverse Fourier transform of a transfer function in the fr domain to obtain the impulse response in the time domain.

Applications
To illustrate the application of the different formulations of the Discrete Fourier Transf equency orms three applications are involving a pair of transforms in lation can be used, the second into combine an input signal with a transfer function of a physical system (e.g.induction in the Earth) his can be done presented: the first which either formu volving the "continuous-discrete" forward transform, and the third involving the continuous-continuous inverse transform.

Application 1: Low Pass Filter
In many geoscience applications it is necessary to filter a dataset or which is equivalent to a filter response.T by convolution of the input signal with the filter impulse response in the time domain or using multiplication by the filter transfer function in the frequency domain, as shown in Figure 3.The frequency domain method is often the preferred process because of the computational efficiencies obtained by using the Fast Fourier Transform.This involves taking a Fourier Transform of the input time series to obtain the spectrum of the input signal.Then multiply this spectrum by the filter transfer function to obtain the output spectrum.An inverse Fourier Transform is then used to obtain the output time series.
In this application we consider filtering of geomag- This shows a mixture of rapid and slow variations and, as in many other geoscience applications, it is useful to filter out some of the variations to more clearly show the other components.In this case we wish to apply a low pass filter with transfer function, T(f), defined by where f c = 0.00027 Hz (corresponding to a period of determine the impulse response in the time domain.Consider the frequency response already used in Application 1, defined in Equation (18). Figure 5 shows this frequency response, both in terms of positive and negative frequencies and as it is ordered in the FFT array.The Fourier integral of such a rectangular (boxcar) function in the frequency domain is a sinc function in the time domain.Because we want to approximate the Fourier integral in this case we use the "continuous-continuous" version of the inverse transform (Equation ( 16b)).This gives the results shown in Figure 6.Where, again, the function is shown as ordered by position in the FFT output array and as ordered in terms of positive and negative time.The values for the sinc function look small, but integrating the sinc function (sum of values multiplied by ∆t (=60 sec) gives a value of 1.0 as required.This is confirmation that the scaling factor used in the inverse transform was correct.Use of the other version of the transform would have applied a different scaling factor and given an incorrect result.
hour).Using this filter response with the Fourier Transforms as described above gives the smoother signal shown in Figure 4(b).
In this application using the transform pairs in Equation (16) or using the transform pair in Equation ( 17) gives the same result.This is because using a pair of transforms involves applying the same overall scaling factor, either t  and f  in Equation ( 16), where , or by applying 1 N alone in equation ( 17).Thus in this case either transform pair is satisfactory.However, in the next two applications we will see that a particular choice of transform is necessary in or

Application 3: Spectrum Determination
There are applications that do not involve filtering but where it is useful to know the spectrum of the signal.One such example concerns the partial saturation of transformers produced by a combination of AC and DC currents flowing through transformer windings.This can occur because of geomagnetically induced currents (GIC) in power systems [10].Consider a power transformer with normal magnetising current I AC subjected to a DC current I DC .The extra magnetic field created by the DC current creates an offset in the magnetic field inside the der to obtain the correct results.

Application 2: Impulse Response
In some cases it is more appropriate to use the time domain method shown in Figure 3 for filter calculations.The frequency response of the filter is defined in the frequency domain and an inverse Fourier transform used to   transformer that pushes the transformer into the saturation region of the hysteresis curve for part of each cycle (Figure 7).
The partial satuation of the transformer during each cycle creates the distorted magnetising current waveform shown in Figure 8(a).For analysing the impact on power system operation it is necessary to determine the spectral content of the distorted waveform.This can be done by a discrete Fourier transform of the waveform which shows that the signal is comprised of the fundamental plus harmonics of the 60 Hz AC frequency (Figure 8(b)).Thus the frequency domain function is not continuous and only contains discrete frequencies.In this case it is appropriate to use the discrete Fourier transform in Equation (17a).
The appropriateness of the continuous-discrete transform in such a case can also be seen by taking the discrete Fourier transform of a cosine wave of frequency, f 1 and amplitude, A. This results in a frequency spectrum with spikes of amplitude A/2 at frequencies ±f .Combincted amplitude of the cosine wave.

Discussion
Filtering of a signal can be done by taking the Fourier transform of the input time series, multiplication by the transfer function in the frequency domain, and then taking the inverse Fourier transform to obtain the output in the time domain (Figure 3).An alternative, equivalent procedure is to convolve the input with the filter impulse response to directly obtain the output in the time domain The frequency domain me is often used because of the computational efficiencies provided by the Fast Fourier Transform; however, there are occasions where the time domain method is preferable.The impulse response is obtained by taking the inverse Fourier transform of the frequency domain transfer function (T.F.) and this has b) 1 ing these two, as in Equation ( 5), gives the expe .thod already been explained in Application 2. Now we need to consider the appropriate formulas to use for the convolution calculation.
The impulse response values in the time domain are samples of a contiuous function.This is to be convolved with the time domain signal which is itself a time series of values that are samples of a continuous function.Thus we need to perform a discrete convolution that is an approximation of the convolution integral.
The convolution of two functions f(t) and g(t) is roximation of this integral is then given by

ractical geoscience
The choice of Discrete Fourier Transform to selection of the signs of the exponent in the exponenttia rse transforms.Here the selections are made based on common practice and to provide the best approximation to Fourier series a grals.

Conclusions
There are a variety of definitions for Fourier integrals and discrete Fourier transforms.This situation is further confused by different software packages using different definitions so that a forward transform in one package can be the same as an inverse transform in another package.This all hinders the production of rigorous reproducile results when using Fourier transforms for p applications.pair reduces l terms and distribution of the scaling factors between the forward and inve nd Fourier inte-The time dependence is chosen as 2π e i ft which is most commonly used in geosciences.This also represents the simple summation of cosine and sine terms, useful at this stage to introduce the exponential forms of th in this equation can be combined into one summation from -N to N 2π e i kn N corresponds to the time dependence of the form 2π e i ft which is most commonly used in magnetotellurics, as shown by the standard papers and textbooks Price [6], Wai 7], Ward and Hohmann [8], and Chave and Wei-].Therefore we will choose the inverse transform with the exponential term 2π e i ft .This automatically requires that the forward transform be written with the exponential term 2π e i ft  .The other consideration our choice of Fourier transform is the placement of any scaling factors (usually 1/N) used with th ummations.For the geoscience appl t

Figure 1 .Figure 2 .
Figure 1.(a) Schematic showing the output array from the Fast Fourier Transform (solid line) and its extension as a periodic function (dashed line); and (b) Allocation of values to positive and negative frequencies.

Figure 3 .
Figure 3. Filtering of a signal via convolution in the time domain or multiplication in the frquency domain.The Fast Fourier Transform (FFT) is used to go between the time domain and frequency domain functions.

Figure 5 .
Figure 5.The frequency response of the boxcar filter.(a) As a function of positive and negative frequencies; (b) As ordered in the array for input to the FFT.

Figure 6 .
Figure 6.Impulse Response (sinc function).(a) As it appears in the output array from the FFT; (b) As a function of positive and negative time.

Figure 7 .Figure 8 .
Figure 7. DC offset in magnetisation of transformer producing a distorted current waveform.
summations and a multiplication with the transfer function or impulse response.

Figure 9
also shows that both calculations involve the scaling factors t   omparing the steps inv calculation (Figure 9(a)) and time domain calculation (Figure 9(b)) we can see that both calcula  and f  .Thus all the same factors are involved in the two calculations showing the equi of the two procedures.valence (a) uency domain; (b) Convolution in the time domain.

Figure 9 .
Figure 9. Computations for filtering of a signal using a filter transfer function (T.F.) by: (a) Multiplication in the freq means that 2π e i ft appears in the inverse transform and consequently the forward transform contains the term 2π e i ft  .he chosen discrete Fourier transform is defined as an approximation to the Fourier integral to transform between s T amples of continuous functions in both the time and frequency domain the (continuous) time domain function is known to be comprised of discrete frequency components instead of a continuous spectrum.In these cases the Discrete Fourier Transforms to go betwe continuous and discrete fun are at involve use of a Fourier transform pair, e.g.filtering by a DFT to the frequency domain, multiplication by the frequency response, followed by inverse DFT, can use either transform pair because the combined scaling factors are the same in each case, of a single transform, either forward or inverse, must use the approisto Pirjola and a referee for useful ERENCES