Spectral Analysis of Geophysical Data

The usefulness of a geophysical method for a particular purpose depends, to a large extent, on the characteristics of the target proposition (exploration target) which differentiate it from the surrounding media. For example, in the detection of structures associated with oil and gas (such as faults, anticlines, synclines, salt domes and other large scale structures), we can exploit the elastic properties of the rocks. Depending on the type of minerals sought, we can take advantage of their variations with respect to the host environment, of the electric conductivity, local changes in gravity, magnetic, radioactive or geothermal values to provide information to be analysed and interpreted that will lead to parameter estimation of the deposits. This chapter deals with some tools that can be used to analyse and interpret geophysical data so obtained in the field. We shall be having in mind potential field data (from gravity, magnetic or electrical surveys). For example, the gravity data may be the records of Bouguer gravity anomalies (in milligals), the magnetic data may be the total magnetic field intensity anomaly (in gammas) and data from electrical survey may be records of resistivity measurements (in ohm-metre). There are other thematic ways in which data from these surveys can be expressed; depicting some other attributes of the exploration target and the host environments. Potential fields for now are mostly displayed in 1-D (profile form), 2-D (map form) or 3-D (map form and depth display). Whichever form, the 1-D and 2-D data are usually displayed in magnitudes against space (spatial data). When data express thematic values against space (profile distance) or thematic values against time, they are called time-series data or timedomain data. The changes or variations in the magnitudes (thematic values) with space and/or time may reflect significant changes in the nature and attributes of the causative agents. We therefore use these variations to carefully interpret the nature and the structural features of the causative agents. Most at times the picture of the events painted in the time-domain is poor and undiscernable possibly because of noise effects and other measurement errors. A noise in any set of data is any effects that are not related to the chosen target proposition. Even where such noise and measurements errors are minimized, some features of the data need to be gained/enhanced for proper accentuation. The only recourse to this problem is to make a transformation of the time-domain data to other forms or use some time-domain tools to analyse the data for improve signal to noise ratio; both of which must be


Introduction
The usefulness of a geophysical method for a particular purpose depends, to a large extent, on the characteristics of the target proposition (exploration target) which differentiate it from the surrounding media.For example, in the detection of structures associated with oil and gas (such as faults, anticlines, synclines, salt domes and other large scale structures), we can exploit the elastic properties of the rocks.Depending on the type of minerals sought, we can take advantage of their variations with respect to the host environment, of the electric conductivity, local changes in gravity, magnetic, radioactive or geothermal values to provide information to be analysed and interpreted that will lead to parameter estimation of the deposits.This chapter deals with some tools that can be used to analyse and interpret geophysical data so obtained in the field.We shall be having in mind potential field data (from gravity, magnetic or electrical surveys).For example, the gravity data may be the records of Bouguer gravity anomalies (in milligals), the magnetic data may be the total magnetic field intensity anomaly (in gammas) and data from electrical survey may be records of resistivity measurements (in ohm-metre).There are other thematic ways in which data from these surveys can be expressed; depicting some other attributes of the exploration target and the host environments.Potential fields for now are mostly displayed in 1-D (profile form), 2-D (map form) or 3-D (map form and depth display).Whichever form, the 1-D and 2-D data are usually displayed in magnitudes against space (spatial data).When data express thematic values against space (profile distance) or thematic values against time, they are called time-series data or timedomain data.The changes or variations in the magnitudes (thematic values) with space and/or time may reflect significant changes in the nature and attributes of the causative agents.We therefore use these variations to carefully interpret the nature and the structural features of the causative agents.Most at times the picture of the events painted in the time-domain is poor and undiscernable possibly because of noise effects and other measurement errors.A noise in any set of data is any effects that are not related to the chosen target proposition.Even where such noise and measurements errors are minimized, some features of the data need to be gained/enhanced for proper accentuation.The only recourse to this problem is to make a transformation of the time-domain data to other forms or use some time-domain tools to analyse the data for improve signal to noise ratio; both of which must be accomplished without compromising the quality of the original data.Even where the picture of the time-domain data 'looks good", we can perform further analyses on the signal for correlation, improvement and enhancement purposes.There are many time-domain and frequency-domain tools for these purposes.This is the reason for this chapter.We shall be exploring the uses of some tools for the analyses and interpretations of geophysical potential fields under the banner of Spectral Analysis of Geophysical Data.The first section of the chapter covers the treatment and analysis of periodic and aperiodic functions by means of Fourier methods, the second section develops the concept of spectra and possible applications and the third section covers spectrum of random fields; ending with an application to synthetic and real (field) data.

Periodic and aperiodic functions
A periodic function of time, t can be defined as f(t) = f(t+T), where T is the smallest constant, called the period which satisfies the relation.In general f(t) = f(t+NT), where N is an integer other than zero.As an example, we can find the period of a function such as f(t) = cos t/3 + cos t/4.If this function is periodic, with a period, T, then f(t) = f(t+T).Using the relation cosθ = cos(θ + 2 m), m = 0, ±1, ±2, ..., we can compute the period of this function to be T = 24 .Aperiodic function, on the other hand, is a function that is not periodic in the finite sense of time.We can say that an aperiodic function can be taken to be periodic at some infinite time, where the time period, T = ∞.A function f(t) is even if and only if f(t) = f(-t) and odd if f(t) = -f(-t).

Sampling theorem
A function, f which is defined in some neighbourhood of c is said to be continuous at c provided (1) the function has a definite finite value f(c) at c and (2) as t approaches c, f(t) approaches f(c) as limit, i.e. lim → = .If a function is continuous at all points of an interval a≤t≤b (or a<t<b, etc), then it is said to be continuous on or in that interval.A function f(t) = t 2 is a continuous function that satisfies the two conditions above.A graph of a function that is continuous on an interval a≤t≤b is an unbroken curve over that interval.In practical sense, it is possible to sketch a continuous curve by constructing a table of values (t, f(t)), plotting relatively few points from this table and then sketching a continuous (unbroken) curve through these points.A function that is not continuous at a point t = c is said to be discontinuous or to have a discontinuity at t = c.The function, f(t) = 1/t is discontinuous at t = 0 as the two requirements for continuity above are violated at t = 0.This discontinuity cannot be removed.On the other hand, a function, f(t) = has a removable discontinuity at t = 0, even though the formula is not valid at t = 0. We can therefore extend the domain of f(t) to include the origin, as f(0) = 1.A continuous function, f(t) can be sampled at regular intervals such that the value of the function at each digitized point is f n , with n = 0, ±1, ±2, ±3, ... (Fig. 1).It will appear as if the time domain function f(t) versus t is transformed into a digitized form domain (f n versus n).

Fig. 1. Sampling a continuous signal f n
To recover a function f(t) from its digital form, the sampling points must be sufficiently close to each other.There is actually a maximum sampling period particular to the function concerned, with which the complete recovery may be achieved.
According to the sampling theorem (or Shannon theorem), a function can be fully recovered by the sampling process provided (a) it is a reasonably well-behaved function and (b) bandlimited (Bath, 1974).Condition (a) implies that the function be continuous with no abnormal behaviour such as discontinuities (sharp breaks or singularities).This condition is practically always fulfilled in the case of functions representing natural processes (such as Laplacian fields).The bandlimitedness of a function as a second condition, refers to a function which possesses a Fourier transform of non-zero value outside it.Most functions involving natural processes also fulfill this condition.If a function is digitized with equal sampling interval, τ, the period present in that function which can be recovered by the process is 2τ, since we need a minimum of two sampling intervals to define one period.The equivalent frequency ( = ) has a special significance in the subject of digitization.The frequency, f N is referred to as the folding or Nyquist frequency.The parameter, τ represents the maximum limit for the sampling period with which we can fully specify a function whose lowest period is 2τ.However, the sampling cutoff frequency, f C (= ) and so f C = 2f N .
It has been shown that the continuous function, f(t) which is band-limited can be reconstructed from the digital values f n by using the formula (Papoulis, 1962;Bath, 1974) = ∑ = ∑ Where c = 2 f c .We see that to recover the continuous function, f(t) from the digital version, the n th sample is replaced by the sinc function (i.e.[sin c t]/ c t) which is scaled by the sample value f n and placed at time, nτ.The scaled and shifted sinc functions are then added together to give the original time function, f(t).Complete recovery also means running values from n = -∞ to +∞, which is practically impossible.
Aliasing is a kind of spectrum distortion which is brought about as a result of too coarse sampling.Fine sampling (implying f N >f c ) and critical sampling (implying f N = f c ), produce no aliasing effects.Coarse sampling (implying f N <f c ) is undersampling and there will be considerable overlap between adjacent spectra in the recovered analogue function.Thus to avoid aliasing effect, we must make the sampling frequency sufficiently high to ensure making the Nyquist frequency at least equal to the cutoff frequency of the original signal.

Fourier analysis
Fourier analysis is the theory of the representation of a function of a real variable by means of a series of sines and cosines.The discussion of Fourier analysis starts with a statement of Fourier theorem.Fourier (1768 -1830) stated without proof and used in developing a solution of the heat equation, the so-called Fourier theorem.

Fourier theorem
Any single-valued function, f(x) defined in the interval [-, ] may be represented over this interval by the trigonometric series, i.e.
The expansion coefficients, a n and b n are determined by use of Euler's formulas: Fourier investigated many special cases of the above theorem, but he was unable to develop a logical proof of it.Dirichlet (1805 -1859) in 1829 formulated the restrictions under which the theorem is mathematically valid.These restrictions are normally called Dirichlet conditions, and are summarized as follows: that for the interval [-, ], the function, f(x) must (1) be singlevalued, (2) be bounded, (3) have at least a finite number of maxima and minima, (4) have only a finite number of discontinuities: piece-wise continuous and (5) periodic, i.e. f(x + 2 ) = f(x).For values of x outside of [-, ], the series in equation (1) above converges to f(x) at values of x for which f(x) is continuous and to [ + + − ] at points of discontinuity.
The quantities f(x+0) and f(x-0) refer to the limits from the right and left respectively of the point of discontinuity.The coefficients are still given by the Euler's formulas.

Different forms of Fourier series
According to the Fourier theorem, a function such as f(t), which satisfies Dirichlet's conditions can be represented by the following infinite series, the Fourier series, as Where the constants a 0 , a n and b n are given by These coefficients (a 0 , a n and b n ) are called Fourier coefficients and the determination of their values is called Fourier or harmonic analysis.
The three expressions in equation ( 3) are respectively obtained by multiplying both sides of equation ( 2) by 1 (i.e.cos 0t), cos nt and sin nt and integrating with respect to t over the period length 2 .We then use the orthogonality properties of sine and cosine functions to eliminate some of the expressions.The Fourier series in equation (2) will exactly represent the function f(t) only when an infinite number of terms are included in the summation, i.e. when n runs from 1 to infinity (∞).If a finite number of terms is taken, the sum will shoot beyond the value of f(t) in the neighbourhood of a discontinuity (boundaries of the function).The overshoot oscillates about the value with a decreasing amplitude as we move away from the discontinuity.
Increasing the number of terms does not influence the error magnitude at the discontinuity, but only leads to a better approximation for the continuous part of f(t).The overshoot and oscillatory behaviour of f(t) at the vicinity of a discontinuity is known as Gibb's phenomenon.
If f(t) is an even function, then equation (2) becomes Since b n = 0 for even f(t), Similarly for odd f(t), a 0 , a n = 0 and so To solve many physical problems, it is necessary to develop a Fourier series that will be valid over a wider interval.To obtain an expansion valid in the interval [-T, T], we let = + ∑ cos + sin and determine such that f(t) = f(t + 2T).In this case, = , hence we obtain: If f(t) satisfies the Dirichlet conditions in this interval, then the Fourier coefficients a 0 , a n and b n can be computed in a similar fashion as in equation (3).
The complex form of the Fourier series in equation ( 6) can be obtained by expressing cos and sin in the exponential using the Euler identity cos θ + sin θ = e iθ , where i = √-1 is a complex number.Thus the complex form of the Fourier series can be written as Where C 0 = a 0 /2 = F(0), C n = (a n -ib n )/2 = F(n) and C -n = (a n + ib n )/2 = F(-n) Equation ( 7) is the complex form of the Fourier series.On multiplying both sides of this equation by e -in t/T and integrating with respect to t, we obtain indicating that the phase spectrum, is antisymmetrical.Equation ( 7) is the time domain Fourier series and equation ( 8) is the frequency domain Fourier series.The two equations are an example of a Fourier transform pair, or we can say f(t) is the inverse Fourier transform of F(n) or C n .

Application of Fourier series
Fourier series, as we have already seen, is used as an effective tool in the analysis of periodic functions.We have also noted that any periodic function having a period, T (and satisfying Dirichlet conditions) can be represented by the infinite series of trigonometric functions.Thus f(t) is represented by the addition of sinusoidal and cosinusoidal waves whose frequencies are integral multiples of some fundamental unit of frequency in the signal, f(t).The frequency, /T is called the fundamental and its integral multiples are called the harmonics.
If we plot the Fourier coefficients a n and b n as functions of frequency, we obtain a number of discrete spectral lines located at fixed spacing of /T.As T→∞, the spacing of the spectral lines approaches zero.Analysis of the spectrum of a field gives the energy content of the field and may represent phenomenal changes in the attributes of the causative agents of the field.This allows systems to be analyzed in the frequency domain so as to find out the frequency response of the system from the impulse response and vice-versa.It can be used as an intermediate step in more elaborate signal processing techniques (e.g. the fast Fourier transform).Fourier integrals (next section) is useful in the analysis of limited-duration (transient) signals.

Fourier integrals
If we look back at equation ( 6) and decide to put the expressions for a 0 , a n and b n that followed inside it and use the dummy variable, λ, we can write equation ( 6) in a more compact form as We let n = n /T (the n th angular frequency), n-1 = (n-1) /T and therefore ∆ = nn-1 = /T and substituting all these in equation ( 9), we obtain If we let T→∞, the following changes will occur in equation ( 10) a.The first part of the expression in the RHS will vanish to zero.b.The increment ∆ becomes very small and in the limit ∆ →d .Thus the digitally increasing n becomes a continuous variable, .c.The summation can be replaced by an equivalent integral with appropriate limits.Thus equation (10) becomes This is the Fourier integral.We can further analyze equation (11a) by using the fact that cos (t-λ) = cos t cos λ + sin t sin λ and so Where = cos and = sin Equation ( 11b) is the Fourier integral: a Riemann sum of integral consisting of the first of the RHS of the equation, called the Fourier cosine integral and the second part called the Fourier sine integral.This integral equation will converge to f(t) when f(t) is continuous and converges even at points of discontinuity just like a Fourier series.
Note that the function is suppose to be defined from -∞ to +∞, but because of the parity of the function, we only need the function from 0 to ∞.This also means that if we only are interested in the range of 0 to ∞, we can define the function from -∞ to any where we want, then we can have either cosine integral or sine integral by extending the function into negative range either in an even or odd form.Thus Fourier cosine and sine integrals are equivalent to the half-range expansion of Fourier series.

Fourier transforms and theorems
From the complex form of Fourier series given in equation ( 7) and expression for the coefficients given in equation ( 8), we make a transition T→∞ and introduce again the variable, = n /T, with ∆ = , since ∆n = 1.Hence equations ( 7) and ( 8) would be written as = ∑ ∆ and = , where C T ( ) = TC n / .If we let T→∞, C T ( )→C( ), and so There are other ways of expressing equations ( 12) and ( 13) which are the Fourier transforms of each other.If we let F( ) = 2 C(-), then Equations ( 14) and ( 15) are called the Fourier transform pair.If f(t) satisfies the Dirichlet conditions and the integral | | ∞ is finite, then F( ) exists for all and is called the Fourier transform of f(t).The function f(t) in equation ( 15) is called the Fourier transform (inverse transform) of F( ) and the pair may be represented by a notation f(t)↔F( ).F( ) is called the amplitude of the time domain field, f(t).S o f a r w e h a v e u s e d a s v a r i a b l e s t and ω, representing time and angular frequency, respectively.Mathematics will, of course, be the same if we change the names of these variables.In describing the spatial variations of a wave, it is more natural to use either r or x, y, and z to represent distances.In a function of time, the period T is the time interval after which the function repeats itself.In a function of distance, the corresponding quantity is called the wavelength λ, which is the increase in distance that the function will repeat itself.Thus, if t is replaced by r, then the angular frequency ω, which is equal to 2π/T, should be replaced by a quantity equal to 2π/λ, which is known as the wave number, k.In three dimensions, we can define a Fourier transform pair as Where and are radius and wave number vectors respectively.If the wave number, is along the z-axis of the coordinate space, then .= kr cosθ and d 3 r = r 2 sinθdθdrdϕ and so the Fourier transform of f(r) becomes, = ∞ as ϕ runs from 0 to 2 , r runs from 0 to ∞ and θ runs from 0 to .
Again, how to split 1/(2π) 3 between the Fourier transform and its inverse is somewhat arbitrary.Here we split them equally to conform to most standard expressions.
There are number of useful theorems that connect the Fourier transform pair in equations ( 14) and ( 15) and these can be found in standard texts on the subject.We shall only mention a few ones that are easily provable as follows: a.The Convolution Theorem -If f 1 (t)↔F 1 ( ) and f 2 (t)↔F 2 ( ) are Fourier pairs, then This a very important theorem which has a wide field application.It simply says that the spectrum of a product of two time functions is the convolution of their individual spectra.The theorem can be extended to any number of functions as Note that in general a convolution is expressed as * = − and * = − .
The mathematical operation of convolution consists of the following steps: 1. Take the mirror image of f 2 (τ ) about the coordinate axis to create f 2 (−τ ) from f 2 (τ).
2. Shift f 2 (−τ ) by an amount t to get f 2 (t − τ ).If t is positive, the shift is to the right, if it is negative, to the left.
3. Multiply the shifted function f 2 (t − τ) by f 1 (τ ). 4. The area under the product of f 1 (τ) and f 2 (t−τ ) is the value of convolution at t.The convolution in the time domain is multiplication in the frequency domain and vice versa.Convolution operation has commutative, associative and distributive properties, like most linear systems.We encounter convolution operations in filtering processes, truncating lengthy data using window functions and sampling a signal with a comb function (digitization).b.The Multiplication Theorem -If f 1 (t)↔F 1 ( ) and f 2 (t)↔F 2 ( ) are Fourier pairs, then Where * is the complex conjugate of F 1 ( ).The product * is called the cross power spectrum.If we put then the previous equation reduces to This is called Parseval's theorem and the real quantity [ ] represents the power spectrum (or energy spectrum) of the function f(t).Thus, it is interesting to note that if the amplitude spectrum F( ) of a given is known, it is possible to compute its power spectrum [ ] and its total energy, c.The Correlation Theorem -The Fourier transform of the cross correlation function ϕ 12 (τ) is the cross power spectrum, E 12 ( ) and that of the auto-correlation function, ϕ 11 (τ) is the power spectrum, E 11 ( ).Thus and so the correlation theorem says that The cross correlation function behaves in accordance with the degree of similarity between the two correlated functions.It grows larger when the two functions are similar and diminishes otherwise.This function becomes zero in the case of completely random data.We see that the time domain cross correlation and auto-correlation functions are reduced to multiplication of the amplitude spectrum.Auto-correlation of two wavelets is equal to the convolution of the first wavelet with the time-reverse of the second wavelet.

Fast Fourier transform
It is important to find the discrete forms of the continuous Fourier transform pair in equations ( 12) and ( 13).These can be written as (Smith, 1999) Both G(n), the discrete amplitude spectrum and the digitized time domain signal, g(k) are periodic as they repeat at every N points.The amplitude spectrum is symmetric while the phase spectrum is antisymmetric.
For real applications, each of equations ( 16) and ( 17) will have to be decomposed into the real and imaginary parts with trigonometric argument of = .Thus an N point time domain signal is contained in arrays of N real parts and N imaginary parts for each equation.Calculating the discrete Fourier transform (DFT) equations takes a considerable time even with high speed computers because of the cycles of computations that must be run.The raw computations of DFT can either be done by use of simultaneous equations (very inefficient for practical use) or by correlation method in which signals can be decomposed into orthogonal basis functions using correlation (not too useful a method).The third and the most efficient method for calculating the DFT is by Fast Fourier Transform (FFT).This is an ingenious algorithm that decomposes a DFT with N points into N DFTs each with a single point.The FFT is typically hundreds of times faster than the other methods.In actual practice, correlation is the preferred techniques if the DFT has less than 32 points, otherwise the FFT is used.Cooley & Tukey (1965) have the credit for bringing the FFT to the scientific world.The FFT, a complicated algorithm is based on the complex DFT; a more sophisticated version of the real DFT and operates by first decomposing an N -point time-domain signal into N time domain signals, each composed of a single point.The second step is to calculate the N frequency spectra corresponding to these N time domain signals.Lastly, the N spectra are synthesized into a single frequency spectrum.An interlaced decomposition is used each time a signal is broken into two, that is, the signal is separated into its even and odd numbered samples.There will be Log 2 N stages required in the decomposition and Nlog 2 N multiplications instead of NxN multiplications without the FFT use.For instance, a 16-point signal (2 4 ) requires 4 stages, 512-point signal (2 7 ) requires 7 stages, a 4096-point signal (2 12 ) requires 12 stages of signal decomposition and so on.The FFT algorithm works on a data length of 2 M , where M is a positive integer (≥5).If the data length is not up to the requirement for FFT operation, then "zeros" are sufficiently added.The decomposition is nothing more than a reordering of the samples in the signal.After the decomposition, the FFT algorithm finds the frequency spectra of a 1-point time domain signal (easy!) and then combine the N frequency spectra in the exact reverse order that the time domain decomposition took place.
The computation of the inverse FFT is very similar to the forward FFT because of the identical nature of the two.Estimation of power spectrum of a signal can be done by means of FFT.

The concept of spectra
Though we have mentioned spectra or spectrum in our previous discussions, we shall formally explain it here.The word spectrum (plural: spectra) is used to describe the variation of certain quantities such as energy or amplitude as a function of some parameter, normally frequency or wavelength.Optical spectrum of white light (colour spectrum) dispersed by a glass prism or some other refractive bodies (such as water) is a good example.When a signal is expressed as a function of frequency, it is said to have been transformed into a frequency spectrum.Thus, mathematically, a time-domain signal, f(t) can be expressed by F( ), where represents angular frequency ( = 2 f; f being the linear frequency).The function F( ) is in general complex and may be represented by 1.The sum of real and imaginary parts: 2. The product of real and complex parts: The modulus | | is normally called the amplitude spectrum and the argument is called the phase spectrum.

Spectral analysis of periodic functions
In modern analysis, the time function may be expressed through certain mathematical transformation into a function of frequency.The discrete Fourier transform views both the time domain and the frequency domain as periodic.This can be confusing and inconvenient since most of the real time signals are not periodic.The most serious consequence of time domain periodicity is time domain aliasing.Naturally, if we take time domain signal and pass it through DFT, we find the frequency spectrum.If we could immediately pass this frequency spectrum through the inverse DFT to reconstruct the original time domain signal, we are expected to recover this signal, save for spill over from one period into several periods -a problem of circular convolution.Periodicity in the frequency domain behaves in much the same way (as frequency aliasing), but is more complicated.

Spectral analysis of aperiodic functions
Equations ( 14) and ( 15) or any of their similar versions give the Fourier transform pair for a periodic function.For non-repetitive (or aperiodic) signal, the period T→∞ and the Fourier transform pair are expressed as Note again that f(t) is real time domain signal, while F( ), the amplitude spectrum is a complex function.

Fourier spectrum
A time function, f(t) such as gravity field may be transform into another function, F( ), where the amplitude of all frequency components present in f(t) and their corresponding phases are expressed as function of frequency.The two transform relations are already given in equations ( 18) and ( 19).The complex function, F(ω), is called the Fourier spectrum and its modulus and arguments as earlier explained are called amplitude and phase spectra respectively.The cosine transform part of F(ω)[ = a(ω) + ib(ω)], a(ω) is called the co-spectrum and the sine transform part, b(ω) is called the quadrature spectrum.

Power spectrum
If is the mean power of a real function, f(t) whose period is T, then (Thompson 1982) Where (f(t)) 2 is termed the instantaneous energy and the complete integration in equation ( 20) is the total (mean) energy of the function.
We have already noted that for two Fourier pairs, f 1 (t)↔F 1 ( ) and f 2 (t)↔F 2 ( ), then f 1 (t).f 2 (t)↔ * * and that The power spectrum | | and its total energy E T are then related by where the power spectrum | | is a real quantity.

Spectral windows and their uses
When we were discussing the convolution theorem, we noted that we might run into convolution operations in truncating lengthy data (signal) by use of window functions.Data windowing can be viewed as the truncation of an infinitely long function, f(t).A boxcar function, w(t) = 1, -T<t<T and w(t) = 0 elsewhere, has a value 1 over the required length (2T) and zero elsewhere.The function, w(t), a time window can be used to truncate f(t) and the truncated time function, f tr (t) = f(t).w(t).Using the convolution theorem, the Fourier transform, F tr ( ) of the truncated function, f tr (t) is given by = * , where F( ) and W( ) are the Fourier transforms of f(t) and w(t) respectively.We can compute the W( ) of the rectangular pulse, w(t) as .Thus = * .This shows that truncating a signal brings about a spectrum modification expressed by the convolution operation between the two spectra, F( ) and W( ).The truncation of a signal, therefore, introduces a smoothing effect whose severity depends on the window length.The shorter the window length, the greater the degree of smoothing and vice-versa.The truncated Fourier transform F tr ( ) is often called the average or weighted spectrum (Blackman & Tukey, 1959).Since all observational data or signals have finite length, truncation effect can never be avoided.
In order to minimize spectral distortion from the signal truncation, other types of time windows may be applied.In general, a window which tapers off gradually towards both ends of the signal introduces less distortion than a window which has near-vertical sides (like the box-car function).A least distortive time window should have the following properties: a.The time interval must be as long as possible.This implies that its corresponding Fourier transform or the spectral window has its energy concentrated to its main lobe.b.The shape must be as smooth as possible and free of sharp corners.The more smooth the time window is, the smaller the side lobes of the corresponding spectral window become (the box-car function is a dirty window!).At this juncture, we shall mention some popular time windows.These include the box-car (rectangular), Bartlett (triangular), Blackman, Daniell, Hamming, Hanning (raised cosine), Parzan, Welch and tapered (rectangular windows).Excellent treatise on spectral windows can be consulted.In general, the Fourier transforms of time domain windows have central main lobe and side lobes in each transform and the magnitudes of the side lobes emphasize the differences between them.Ideally, the main lobe width should be narrow, and the side lobe amplitude should be small.Windows are also extensively used in designing filters and the window parameters (side lobe amplitude, transition width and stopband attenuation) must be used for the design.

Fourier spectrum of observational data
To compute the spectrum of a function f(t) which obeys Dirichlet conditions, Fourier transform is applied to it directly, particularly the use of Fourier integral equations.We can use the basic theorems already presented to evaluate the transforms.The common types of functions which are usually subjected to Fourier analysis are those obtained by some kind of physical measurements.Functions which represent observational data are normally converted into digital form (if presented as a continuous plot in profile or map forms) so that their spectra can be computed by numerical Fourier transformation.Observational functions are usually not continuous and not infinite as the theory of Fourier transformation demands.For this reason, observational spectra suffer from two types of distortions: (1) truncation effect (by a window function) and (2) digitization effect.When a signal is digitized, its spectrum becomes periodic and so the original spectrum (scaled by the inverse of sampling period) becomes repetitive with the same frequency as that used in sampling the signal.Coarse digitization results in distorted spectrum.The extent of digitization effect depends on the sampling frequency as well as on the cut-off frequency of the signal (see Section 3).

Random functions
A random variable is a real-valued function defined on the events of probability system.A random variable, f(t) emanates from a random or stochastic process: a process developing in time and controlled by probabilistic laws.A random (or stochastic) process is an ensemble or set of functions of some parameter (usually time, t) together with a probability measure by which we can determine that any member, or a group of members has certain statistical properties.Like any other functions, random processes can either be discrete or continuous.At any point in a medium, a unit mass or a unit charge or a unit magnetic pole experiences a certain force.This force will be a force of attraction in the case of gravitational field.It will be a force of attraction or repulsion when two unit charges or two magnetic monopoles of opposite or same polarity are brought close to each other.Every mass in space is associated with a gravitational force of attraction.This force has both magnitude and direction.For gravitational field, the force of attraction will be between two masses along a line joining the bodies.For electrostatic, magnetostatic and direct current flow fields, the direction of the field will be tangential to any point of observation.These forces produce force fields.These fields, either global or man-made local fields are used to quantitatively estimate some physical properties at every point in a medium.Most geophysical potential fields, in particular gravity and magnetic fields are caused by an ensemble of sources distributed in some complex manner, which may be best described in a stochastic or random framework.We shall examine some characteristics of these fields.

Geophysical potential fields
The potential field, ϕ(x, y, z) in free space (i.e.without sources) satisfies the Laplace equation When sources are present, the potential fields satisfy the so-called Poisson equation Where (x,y,z) stands for the density, magnetization or conductivity depending opun whether stands for gravity, magnetic or electric potential respectively.It is important to know that both global or local fields are subject to inverse square law attenuation of the signal strengths.It is at its simple peak in gravity work where the field due to a point mass is inversely proportional to the square of the distance from the mass, and the constant of proportionality (the gravitational constant, G) is invariant.Magnetic fields, though complex, also obey the inverse square law.The fact that their strength is, in principle, modified by the permeability of the medium, is irrelevant in most geophysical work, where measurements are made in either air or water.Magnetic sources are, however, essentially bipolar and the modifications to the simple inverse-square law due to this fact are important.A dipole here consists of equal-strength positive and negative point sources: a very small distance apart.
Field strength here decreases as the inverse cube of distance and both strength and direction change with "latitude" (inclination) of the Earth's magnetic field.The intensity of the field at a point on a dipole axis is double the intensity at a point, the same distance away on the dipole "equator", and in the opposite direction.
Electric current flowing from an isolated point electrode embedded in a continuous homogeneous ground provides a physical illustration of the significance of the inverse square law.All the current leaving the electrode must cross any closed surface that surrounds it.Usually the surface is spherical, concentric with the electrode and the same function of the total current will cross each unit area on the surface of the sphere.The current per unit area will be inversely proportional to the total surface (half-space) of 2 r 2 .Current flow in the earth, however, is modified drastically by conductivity variation.
The potential fields of either gravity, magnetic or electrical fields are the ones given by either the Laplace or Poisson equations.Some of the useful properties of (x, y, z) are (i) given this potential field (scalar) over any plane, we can compute the primary or force field (vector) at almost all points in the space by analytic continuation and (ii) the points where the force field cannot be computed are the so-called singular points.A closed surface enclosing all such singular points also encloses the sources which give rise to the potential field.Thus the singularities of the potential field are confined to the region filled with sources.All these properties are best described and accentuated in the Fourier domain.We shall therefore express the Fourier transformation of the potential field in two or three dimensions (see Section 4.5).In two dimensions, the Fourier transform pair,  (u, v) and its inverse (x, y) are given by and Where here, u and v are coordinates of the Fourier plane.Equation ( 24) is also known as the Fourier integral representation of (x, y).Equation ( 23) exists only if and only if a condition generally not satisfied in most geophysical situations except for an isolated anomaly (Roy, 2008).However, the Fourier transform of a real function in two dimensions possesses the following symmetry: If (x, y, z) is the potential field on a plane z, satisfying the Laplace equation (equation ( 21)), its Fourier integral representation is given by (Naidu 1987) Where H(u, v, z) is to be determined by requiring that it also satisfies the Laplace equation and is only true if it satisfies the differential equation whose solution is H(u, v, z) = e  s z for all values of z, where =√ + .For z ≥ 0 equation (26) becomes Equation ( 27) is a useful representation of the potential field and forms the basis for derivation of many commonly encountered results.

A random interface
Potential fields that are caused by stochastic sources are of two types: a random interface separating two homogeneous media (e.g.sedimentary rocks overlying a granitic basement) and a horizontal layer of finite thickness within the density or magnetization varying randomly.We can relate the stochastic character of the potential fields to this interface or layer with the aim of determining some gross features of the source model (e.g.depth to interface). Figure 2 shows a homogeneous random interface separating two media.Let f(x, y) be a homogeneous (stationary) stochastic field.The spectral representation (or Cramer representation) of a homogeneous random field is given by (Yaglom, 1962) Where and F(u, v) is the generalized Fourier transform of f(x, y).Some of the properties of dF( Where S f (u, v) is the spectrum (power spectrum) of the stochastic field, f(x, y).We can further obtain some basic properties of the random potential field in free space.If (x, y, z) is a random potential field in free space, then similar to equation ( 26), Where H(u, v, z) is selected so that (x, y, z) satisfies Laplace equation and d(u, v) is a random function.When equation ( 29) is substituted in the Laplace equation, we obtain and so equation (29) becomes Fig. 2. A homogeneous random interface separating two media of uniform density/magnetization The spectral representation of the three components (x, y, z) of the potential field is easily obtained by differentiating equation ( 30) with respect to the coordinate axis , , =− , , , =− , and the spectrum of the potential field and its components may also be obtained from equation ( 30).
The cross-spectrums between different components can also be computed from the field equations resulting from the coordinate components obtained from equation (30).
We again look at the random interface (Fig. 2) for the cases of gravity and magnetic fields.
Here we see that the potential field is observed on a plane, h unit above the interface.Medium I can be replaced with a vacuum (i.e.where density/magnetization = 0) and the lower medium (Medium II) is characterized by density or magnetization equal to the difference in density or magnetization of Medium I and Medium II.It is easy to see that the observed field will consist of two components: a constant part representing the field due to a semi-finite medium with its upper surface as a horizontal surface, and a random component representing the field due to the random interface.
For the gravity aspect, the vertical gravity field due to a random interface on the observation plane may be expressed as (Telford et al., 1990;Roy, 2008) Where ∆ is density contrast between Medium II and Medium I with the mass element located at point (x 0 , y 0 , z 0 ) and ∆z(x 0 , y 0 ) is a homogeneous random function.Equation ( 31) can be further expressed using a generalized Fourier transform, ∆Z(u, v) of the random function, ∆z(x, y) as (Naidu & Mathew, 1998) Where ∆Z n (u, v) is define from the following relation When ∆z(x 0 ,y 0 )<<h The magnetic field due to a random interface may be obtained either by repeating the mathematical procedure for a magnetic field or by making use of the relationship between the magnetic and gravity potential (Poisson's relation).The magnetic potentials and fields can be estimated from gravitational potential using Poisson's relation expressed as = , where φ is the magnetic potential, is the gravitational potential, I, the magnetic polarization, , the volume density of the medium and G is the gravitational constant.The horizontal and vertical components of the magnetic field are respectively, =− and =− .
By using the latter approach, the magnetic potential, f T (x, y, h) is expressed as (Naidu & Mathew, 1998)

Discrete sampling of potential fields
Potential fields are continuous functions of space.They have to be sampled for processing on a digital computer.How do we sample them?A homogeneous random field, f(x) is sampled at x = n∆x, n = 0, 1, 2, ..., where ∆x is the sampling interval.The accuracy of the sampling is judged by whether by using the samples, we can recover the original function with as small an error as possible (i.e. with the mean square error identically zero).This can be achieved as we have earlier noted that (i) the spectrum of the field be band-limited and (ii) sampling rate (= 1/∆x) is at least twice the highest frequency present.Band-limitedness of the spectrum of the field means = | | ≥ , where u 0 = 2 /λ 0 ; λ 0 being the smallest wavelength corresponding to the highest frequency and condition (ii) implies that ∆x = λ 0 /2.A random field in two dimensions can be sampled in more than one way: we may use a square or rectangular grid, polar grid or hexagonal grid.Such a choice of sampling patterns is not available for one dimensional field.The Fourier transform of discrete data turns out to 6.5 Angular and radial spectra of the field Naidu (1969) and Mishra & Naidu (1974) gave the spectrum of a two-dimensional magnetic survey as Where S (u, v) is the power spectrum of a 2-D aeromagnetic field, X (u, v) is the Fourier transform of the field, L x and L y are length dimensions, and u and v are frequency in the x and y directions respectively given as u = 2 /L x and v = 2 /L y .A two-dimensional spectrum can be expressed in a condensed form as one-dimensional spectra; radial and angular (Spector & Grant, 1970;Mishra & Naidu, 1974;Naidu, 1980;Naidu & Mishra, 1980;Naidu & Mathew, 1998).The radial spectrum is defined as Where again =√ + is the magnitude of the frequency vector and  = tan -1 (v/u) is the direction of the frequency vector in the spatial frequency plane v and u (in radian/km) in the x-and y-directions respectively.The angular spectrum is defined as Where, s is the radial frequency band starting from s 0 to s 0 + s, over which the averaging is carried out.It is useful as a rule to look at power spectra in one-dimensional or profile form rather than in two-dimensional or map form.This is because S is a somewhat bumpy function of  (Fedi et al., 1997) when the width of the model is moderately large, and the bumpiness imparts a certain irregularity to the contours (Spector & Grant, 1970).The power spectra in one-dimension also enables ensemble of magnetic block parameters (average magnetic moment/unit depth, average depths to top, average thickness and average widths) to be factored out completely for effective analysis.Usually the angular spectrum is normalized with respect to the radial spectrum so as to free this spectrum from any radial variation.In this case (Naidu & Mathew, 1998) The computations of spectra in equations ( 40) and ( 41) may require the use of template.The radial spectrum is computed by averaging the 2D spectrum over a series of annular regions while the angular spectrum is computed by averaging over angular sectors in the template.Naidu (1980) and Naidu & Mathew (1998) have shown that the angular spectrum of the total field of a uniform magnetized layer having an uncorrelated random magnetization is given by www.intechopen.com

Spectral Analysis of Geophysical Data 47
Where h is depth to the magnetic layer, I 0 is the direction of the current earth's magnetic field, ,  and  are directional cosines while  0 is the declination of the earth's magnetic field.The non-exponent component of equation ( 44) is exactly the same as the square of Γ(u, v) expressed in equation ( 34).The angular spectrum, A norm (θ) may now be expressed as Where Ω is a constant.Thus the angular spectrum is a maximum in the direction of polarization vector.Naidu (1970) had earlier shown that the shape of the angular spectrum is a product of a number of factors such as rock type, strike and polarization vector, but the latter two factors influence the shape of the spectrum significantly.Thus the presence of peaks in the angular spectrum gives an indication of linear features in the map.
On the other hand, the radial spectrum gives a measure of the rate of decay with respect to radial frequency of the spectral power, which may represent a deep-seated phenomenon (Bhattacharyya, 1966;Naidu, 1970, Spector & Grant, 1970).

Estimation of radial and angular spectra of aeromagnetic data
We shall apply the concept of radial and angular spectra to synthetic and real data.For the synthetic data, we have taken the field over a magnetic dipole buried at a depth of 4 km.Other parameters of the dipole include: inclination and declination of the inducing field (and remanent field of the dipole) on the dipole are respectively 6 o and 8 o (i.e. at low magnetic latitudes), where the field strength of the inducing field is 33510 nT and the magnetization intensity is 0.01 A/m, while its susceptibility is assumed as 9.5 x 10 -5 cgs (0.0012 SI).For space, the anomaly map is not shown.In computing the angular spectrum of this dipole anomaly, the 50x50 data grid was padded with sufficient zeros and cosine tapered to make data matrix amenable for FFT computations.This resulted in a data matrix size of 64x64.The angular spectrum is then calculated for three frequency sub-bands (1-10, 10-20 and 20-30 frequency numbers) for angular interval of 180 o .The highest frequency in the data is 32. Figure 3 shows the curves representing the low (1-10), middle (10-20) and high (20-30) frequency numbers.Observation of these curves shows that 8 o and 90 o spectral peaks appear in the three frequency subbands followed by a peak at 156 o in the mid to high frequency sub-bands though with reduced magnitude.Thus we see the emergence in the three frequency sub-bands, of 8 o : the declination of both the inducing and remanents fields.The other angular feature (90 o ) at subdued level is due to the tapered window used in the analysis.The other peak (156 o ) corresponds to the direction of the polarization in the direction of the inducing field at this magnetic latitude.The real data used for the computation of radial and angular spectra were obtained from aeromagnetic total field intensity of the Middle Benue Trough, Nigeria(MBT) collected from 1974 to 1976.The MBT is the central part of the main Benue Trough of Nigeria.The Benue trough is linked genetically to the oil/gas bearing rocks of the inland Nigerian Niger Delta area.With an upbeat in petroleum efforts in the inland basins, attention is now focused on the Benue trough with more emphasis laid on the structural setting of the basin.
The composite total field intensity data were corrected for the main field using the IGRF 1975 model and this resulted in the residual field map (Fig. 4) used for the present analysis.The angular and radial spectra of the residual total field magnetic map (Fig. 4) were computed.First, for detailed delineation of angular features, the map was divided into two sub-areas: north/south or upper/lower.The upper part (enclosing nearly 90% of the sediments) has a size of 128x128 and the lower portion has a size of 57x120.It is computationally efficient to use the Fast Fourier Transform (FFT) algorithm to accomplish the analysis of the angular and radial spectra.The upper part of the map (128x128) is already amenable for FFT application since the matrix size is already 2 N , where N is an integer.The lower part however, is cosine tapered and then padded with sufficient zeros to make data matrix amenable for FFT application.This resulted into another data of 128x128.The 4 o peak appearing in the three sub-bands corresponds to the value of the inclination of the Earth's magnetic field in the area.The peak at 56 o corresponds to the general trend of the Benue Trough (N56E).Other angular values have been fully mapped in the area (Benkhelil, 1982(Benkhelil, , 1989;;Likkason, 2005).

)
Note that the amplitude spectrum, | | = | − | = is symmetrical.From the identity in complex number representation: a + ib = re iθ , with = √ + and tan θ = b the two points overlap [i.e.(u, v) = (u', v' )Where Γ(u, v) = (im x u + im y v -m z s)(i u + i v -s) and (m x , m y , m z ) are components of the magnetization.Note that m x = I x ∆k, m y = I y ∆k and m z = I z ∆k, where ∆k = (k medII -k medI ) is the magnetization contrast between the medium II and medium I and ( , , ) are the direction cosines of the gradient direction and (I x , I y , I z ) are the three components of the inducing magnetic field.When |∆ ,

Fig. 5 .
Fig. 5. Radial spectrum plotted against frequency number of the total field magnetic anomaly over the Middle Benue Trough, Nigeria.Five linear segments are marked (A, B, C, D, E).Estimates of depths from slopes are indicated.