Filtering of potential field data, although not always referred to as filtering, can be found in the literature as early as 1937. The theoretical methods to calculate the upward and downward continuation and first and second derivatives were first described by Peters (1949) and Henderson and Zietz (1949a and b). The application of these theoretical methods consisted of convolving data sets with a filter in the space domain. The mathematics involved were complex. Owing to the lack of computers to perform these calculations, application of the theoretical or ideal filters was impractical. Many methods to approximate the filtering of data were developed (Griffin, 1949; Peters, 1949; Henderson and Zietz, 1949; Elkins, 1951; Rosenback, 1953; Danes, 1962; Byerly, 1965; Fuller, 1967; Oldham, 1967, Clark, 1969). These methods usually consisted of the creation of coefficient sets to be convolved with select data points to approximate the filter. Instead of having many data values to use in the convolution, only a few values around each data point are needed. The mathematics required less work, but the results were marginal.
Many authors recognized the application of the Fourier integral to the filtering problem. Fourier transforms, initially in the form of the Fourier series, were used by Tsuboi and Fuchida (1937) to estimate subterranean mass distribution from gravity data, Tsuboi (1959) for upward continuation, Danes and Oncley (1962) for the second vertical derivative, and Agarwal and Lal (1972) for the vertical gradient. Bhattacharyya (1965) extended the filtering method to include the reduction to the pole of magnetic data using the Fourier transform. However, these authors used the Fourier transform only to determine supposedly improved sets of coeffiecients for convolution. Calculations were still perfomed in the spatial domain. Nettleton (1954) stated that no single set of coefficients would give good results with all data sets, and Darby and Davies (1967) showed that all coefficient sets are only poor approximations to the ideal filter. Darby and Davies were the first to show that better results can be achieved by filtering directly in the frequency domain. The ideal filter can be used directly if the data are transformed into the frequency domain, filtered, and then transformed back into the space domain. However, Darby and Davies only showed how to perform the transformation using the Fourier integral directly. This took just as much time as convolving and therefore the method was not very useful. They did not make use of the Fast Fourier Transform (FFT) computer algorithm developed in 1965 by Cooley and Tukey. It is the Cooley and Tukey algorithm which makes it advantageous in regards to time to perform filtering directly in the frequency domain with the wavenumber filtering method.
Wavenumber filtering refers to the isolating or enhancing of data in the wavenumber or frequency domain. To perform wavenumber filtering it is necessary to convert anomalies in the gravity field, represented by a data matrix along an X,Y coordinate system, to a two-dimensional set of amplitudes over a range of frequencies or wavenumbers. This is done with the Fourier integral. The Fourier integral can be used to transform a data set in the space domain to the frequency or wavenumber domain. Once in the wavenumber domain, the proper filter can be applied. The filtered data in the wavenumber domain can then be transformed back into the space domain in the same manner using the inverse of the Fourier integral.
It is desirable to filter the data in a way that will isolate or enhance certain features. One may wish to see anomalies of only a certain size or of a certain wavelength. The words frequency, wavelength and wavenumber are often incorrectly used interchangeably. Frequency refers to period, usually dealing with time, as in cycles per second. It may also deal with space, as in cycles per kilometer. Wavelength is the inverse, seconds per cycle or kilometers per cycle. The term wavenumber is in reference to the mathmatics involved in the computer program. After converting a data set from the space domain to the frequency domain, one has a set of values or amplitudes at discrete intervals. The first wavenumber is the frequency associated with the first amplitude. The range of frequencies that can be observed is 0 to 1/(2 times the station spacing). The number of wavenumbers in this range depends on how much data is available. Once the wavenumber spectrum of the data set is known, the wavenumbers that would most likely contain the information of interest can be isolated by cutting out other wavenumbers. One way to achieve this is multiplying the spectrum by a box function. The wavenumbers to be saved are multiplied by 1 and the wavenumbers to be delete are multiplied by 0. When the data set is transformed back into the space domain, only anomalies of a certain size remain. Alternately, certain wavenumbers can be multiplied by factors greater than 1 and other wavenumbers could be multiplied by factors less than 1. This would enhance certain anomalies without actually cutting or deleting any of the spectrum.
From Grant and West (1965), the gravitational potential U at location x,y,z due to masses located in z > 0 where the z-axis is positive downward, is
Equation 1
where 
Since gravity measurements are not of the gravitational potential, but the
z derivative of that potential field, 
. Differentiating equation 1
with respect to z gives
Equation 2
which shows that if 
is known
everywhere on z=0, it can be calculated at any z < 0, that is any level
above the level of the original data. This is known as Upward Continuation.
This equation in its present form cannot be used in real applications.
is not known everywhere on a surface, but only at discrete locations or
intervals over a finite region. In terms of a matrix of gravity values with
Nx columns and Ny rows, the discrete form of equation 2 would be
Equation 3
where 
is the known data at level 0
and 
is the upward
continued data. In this form, a matrix of gravity values at elevation z
can be calculated from a matrix of gravity values at elevation 0.
Equation 3 can be used directly to upward continue a data set, but it is usually not practical to do so. For example, to calculate the upward continuation at one location in a matrix of 128 rows by 256 columns, it would take 32768 summations. For each of the 32768 cells in the matrix, this is over 1 billion summations. Even with today's computers, this is often too time consuming. Performing a transformation on equation 2 from the space domain to the frequency domain using the Fourier Transform simplifies the equation from convolution to simple multiplication. However, as it will be shown, using a direct Fourier transform takes just as much time. Only because of the Fast Fourier Transform computer algorithm can this type of filtering be carried out easily even on small computers.
As shown by Grant and West (1965, p. xx), taking the Fourier transform of both sides of equation 2 simplifies the upward continuation equation. Transforming the data from the time or spacial domain to the frequency domain allows convolution to be replaced by multiplication. An additional benifit of performing Upward Continuation in the frequency domiain is the restriction of z<0 found in equation 2 is no longer necessary. In the frequency domain it is possible to continue upward or downward.
In two dimensions the Fourier integral is
Equation 4
where x and y are space and p and q are the respective wavenumbers.
In the transformation of the upward continuation equation (2), the subscripts o and z are used to denote values of g taken at z=0 and at a depth z. Positive z values refer to downward continuation, while negative z values refer to upward continuation. The Fourier Transform of the left side of equation 2 is
Equation 5
This is denoted as 
. The right hand
side of equation 2 becomes
Equation 6
remembering that
Grant and West (1965, p. 218) show that the second integral in
equation 6 is equal to
Equation 6 becomes
Equation 7
The integral within equation 7 is the Fourier transform of
.
Calling this 
,
the upward continuation equation becomes
Equation 8
Equation 8 shows that in the wavenumber domain, the gravity at level z is
equal to the gravity value at level 0 times a factor
. This is a simple multiplication.
There are no summations to take up computer time. If the value of , the Fourier transform of , is known at wavenumber location p,q, then can be calculated at p,q. To calculate for the entire matrix requires one step for each cell
in the matirx. As in the example above, for a matrix of 128 by 256, this
requires only 32768 steps. Taking the inverse Fourier transform of gives , the gravity
at level z in the space domain. Unlike equations 2 and 3, no singularity
occurs when z is positive. This continuation equation can be used directly
for both negative and positive z, continuation upward and downward.
The method used to filter gravity data using the Fourier Transform can be
seen diagramatically in figure 1. When one begins with the original gravity
data ,
taking Fourier transform of results in
, the gravity data in the wavenumber domain.
Multiplying by a filter gives
, the filtered gravity in the wavenumber domain.
Taking the inverse Fourier transform of gives
the final answer, , the filtered
gravity data.
The Fourier Transforms allows an easy method to calculate the upward or downward continuation of a gravity field to any other level. However to use the results of the continuation filter one must understand what continuation can show. Upward continuation of a field to an elevation z shows what the field would look like if it were to be measured at that elevation z. In regards to determining structure, the gravity field at the surface would be highly affected by small near surface bodies, disguising the gravity effect from deeply buried bodies. The gravity field a great distance above the surface would be less influenced by small bodies, only showing the gravity due to large deep seated structure. As the gravity filed is continued upward, only the most regional trends would remain. In the frequency domain, when the data is multiplied by the filter
Downward continuation allows the user to see what the gravity field would look like if one could measure it closer to the source. Anomalies become sharper as the field is continued downward. Higher wavenumbers are enhanced. From this type of filtering, boundaries between regions of different densities become apparent. However, the highest wavenumber anomalies are usually associated with very shallow bodies or noise. It is often desirable to remove these anomalies so that the effects of the downward continuation on the larger anomalies can be seen. If downward continuation is performed on a potential field to a level Z, one would want there to be no anomalies from bodies above that level Z in the data set. If there is, the field will become unstable. An easy way to avoid this is to first perform band-pass filtering to remove small wavelength and presumably shallow feature anomalies and then continue the data downward.
Derivatives of the gravity field are also useful for structural interpretation. Although the gravity equation can be differentiated in any direction, the derivatives with respect to z are the most useful. Differentiating both sides of equation 8 with respect to z gives
Equation 9
As in the discussion of the upward continuation, is the gravity field at elevation z in the frequency or
wavenumber domain, while is the
gravity field at elevation 0 in the wavenumber domain. Since the
derivative of the original data (the gravity field at elevation 0) is of
interest, z is set to 0.
Equation 10
This equation shows that using the Fourier Transform method, the first
vertical derivative of gravity can be calculated by multiplying each cell in
the wavenumber matrix 
by the filter
. Once
again there are no summations in this equation, so it can be done quickly by a
computer. The steps in performing first derivative filtering are the same as
described in Figure 1. The gravity data is transformed into the wavenumber
domain, multiplied by a filter, and then transformed back into the space
domain.
Equation 10 could have been derived by first differentiating equation 2 and then taking the Fourier Transform, but the mathmatics would have been more difficult.
The second vertical derivative can be obtained by differentiating equation 9 again with respect to z.
Equation 11
As before, setting z=0, the exponential term drops out,
Equation 12
This equation is identical to equation 11 in Darby and Davies (1967). In general, the Nth vertical derivative is
Equation 13
Many authors (Evjen, 1937, Elkins, 1951, Rosenback, 1953, Mescko, 1966) have written about the application of various vertical derivatives of potential fields, but the concept still remains abstract. The first vertical derivative of gravity is also known as the vertical gradient. Evjen (1937) states that the vertical gradient of gravity gives better resolution of structure in certain instances. He gives the example of two cylinders of infinite length buried at depth. The gravity across such a profile will show two distinct anomalies, one for each cylinder, only if the cylinders are not too deep. Below a certain depth the profile will show only one anomaly. If one looks at the vertical gradient for such a profile, the two anomalies will more likely be discernable. The vertical gradient will show certain structure more clearly. In addition, Evjen states that the vertical gradient can be used to estimate the depth of a structure and help define the edges of a structure.
The second derivative of gravity, can also be used to help identify certain features in the structure. Mescko (1966) shows how geological structure can often be delineated better in second derivative maps than maps of the original gravity anomalies. Rosenbach (1953) states that the second derivative can give a clearer picture of the geology because local density anomalies tend to be emphasized over regional features. Elkins (1951) also states the second derivative will show the smaller shallower geologic anomalies at the expense of larger regional fetures.
Since a magnetic field is also a potential field, it can be analyzed using this method as well. Henderson and Zietz (1949a) state that the second vertical derivative of a magnetic anomaly map helps delineate those geologic formations that are magnetic.
Two other types of derivative filters are available in addition to the
vertical derivatives. By differentiating equation 8 N times with respect to z
equation 13, the Nth vertical derivative of , was derived. Differentiating equation 8 N times with respect to x or y
gives the Nth horizontal derivative in the x or y direction. In general
Equation 14
and
Equation 15
where j is the imaginary number 
.
In the horizontal derivative,
when N is odd the filter becomes imaginary and antisymmetric (Bhattacharyya,
1972). This requires special construction of the filter in the computer
program.
The Fourier transform converts data in the space-time domain to the frequency or wavenumber domain. Multiplying the transformed data by a filter allows certain features in the potential field to be enhanced. Equations 8, 13, 14 and 15 show the filters for upward/downward continuation and Nth order vertical and horizontal derivatives. The Fourier method gives the opportunity to filter the data in another way. Features can be isolated by selectively removing wavenumbers. It can be assumed that small bodies as well as near surface bodies will show up as anomalies with small wavelengths or high wavenumbers. Large bodies and deep seated bodies would show up as anomalies with long wavelength or low wavenumbers. If only large or deep seated structures are of interest, it may be useful to remove all anomalies that have wavelengths smaller than a certain value. It also may be useful to cut out large wavelengths, representing the most regional trend, and only analyze the remaining anomalies. This type of filtering is called band-pass filtering. All frequencies or wavenumbers between certain values would be saved or passed, while the rest would be cut. This can be done according to the following equation:
Equation 16
Cutting all wavenumbers below a certain value is called high-pass band-pass filtering. It can be done to remove regional trends. Cutting all wavenumbers greater than a certain value is called low-pass band-pass filtering. It is used to remove small features and noise. Band-pass filters can be constructed to remove both high and low wavenumbers as well.
Another possibility is to isolate certain anomalies with respect to their orientation. It may be desirable to look only at features that have trends between certain orientations. If one wishes to look at all features that have east-west trends, a filter could be constructed that would pass all such features between azimuths of 65 and 115 degrees and cut all others. Strike-pass filtering would follow the equation
Equation 17
It has been shown that gravity data can be isolated or enhanced in various ways by filtering in the frequency or wavenumber domain. Filtering by convolution in the time or space domain requires only multiplication in the frequency or wavenumber domain. Converting from one domain to the other is done with the Fourier Transform. In the next section the method of performing the Fourier Transform on a data set is discussed.
The equations to filter gravity data described in the previous section show how certain aspects of the gravity field can be isolated or enhanced easily and quickly by working in the wavenumber domain. This is done by transforming the data set from the space domain with the Fourier Transform. The Discrete Fourier Transform, (DFT), one that can be carried out on a discrete data set, is
where x and y are functions of n and m respectively. However, it can be seen that the number of summations required to calculate the transformation is just as many as to calculate the upward continuation (equation 2). To perform filtering in the wavenumber domain as described above, both a forward and an inverse transformation must be performed. It seems that taking the Fourier transform of the discrete data set has not gained anything. The computer time required is still prohibitive. To be able to use the simplified filtering equations a method of quickly calculating the Fourier tranform is necessary.
Cooley and Tukey (1965) showed a method of calculating the Fourier transform in a minimum number of steps. As explained by Naidu (1971) this quick method divides
a given sequence into a number of subsequences and to compute the Fourier transform of each, and finally to combine all these transforms into a single transform which correspondes to the tranform of the entire sequence.
If Nx and Ny are highly composite numbers, then the matrix can be factored into smaller matrices before multiplication. Cooley and Tukey showed that for a digital computer, the most efficent calculation of the Fourier transform occurs when the number of rows Ny and the number of columns Nx are powers of 2. The Cooley and Tukey method requires
Ny log`D10`U(Ny) Nx log`D10`U(Nx)
steps to perform the Fourier Transform of a matrix. For a matrix of 128 rows by 256 columns, the total number of steps would be about 150,000. Since 1965, many other methods of calculating the Fourier transform have been published, some faster than others. Winograd (1978) stated that his method requires the same number of additions as Cooley and Tukey, but only about 20% of the number of multiplications. Other methods have been developed to transform data sets that have lengths of powers of 3, 4, 8, etc. The method used in this study is from Reed (1980), modified from Claerbout (1976).
Cooley and Tukey's method and the methods developed by others in later years, are commonly referred as Fast Fourier Transforms, or FFTs. The FFT is simply an algorithm that efficiently computes the discrete Fourier transform. It can transform a set of data from the space or time domain into the wavenumber or frequency domain. It can also perform the inverse: transform a set of data from the wavenumber or frequency domain into the space or time domain. Given a data set of gravity at discrete intervals, the FFT will return the amplitudes and phases at a discrete number of wavelengths that make up the anomalies. The FFT uses every piece of data in the matrix to determine the wavenumbers or frequencies. Its accuracy depends more on the number of significant figures in the data set and the size of a word on the computer (ie., small 8-bit computers may have more truncation errors than 32-bit computers.) A good description of the Fast Fourier Transform and it uses can be found in Brigham (1974).
To use the FFT the spatial data is put into a complex matrix and passed to the FFT routine. Since spatial data is real only, the imaginary portion of the complex matrix is set to zero. The FFT routine returns a set of wavenumbers in complex form. The amplitude spectrum and the phase angle can be calculated from this complex matrix, but for wavenumber filtering, this is not necessary. Wavenumber filtering is performed on the matrix in its complex form. Although attempts have been made to interpret gravity anomalies directly from the frequency spectrum, they have had little success. (e.g., Regan and Hinze, 1976)
If a potential field is known in all locations in every direction, then an infinite number of wavenumbers can be determined. However, data is known only at distinct intervals in a finite region. Therefore the wavenumber content of the data can only be determined over a finite range and at distinct intervals. With the gravity data at equally spaced stations in a matrix, the smallest wavelength that can be determined is called the Nyquist. It has the value of 2 times the station spacing. The Nyquist frequency is the inverse of the Nyquist wavelength. The range of frequencies that can be determined for a data matrix is -Nyquist to +Nyquist. Wavenumbers are integer multiples of the Nyquist. The number of wavenumbers that can be determined from the matrix depends only on the number of rows and columns in the matrix.
Early methods of filtering potential field data were only approximations of the discrete convolution filtering equations. The amount of calculations required for discrete convolution made these equations difficult to use. Approximating the equations allowed calculations to be performed much more quickly, but since the approximation methods did not truly represent the equations, they produced a large number of errors. By using the FFT algorithm, the filtering equations can be used directly. Those errors due to approximating the filtering equations and performing convolution on a limited data set have been edliminated. Any remaining sources of error lie only within the FFT algorithm. By understanding these sources of error, their effect can be limited. Errors due to the FFT method include leakage and Gibb's phenomenon. Leakage refers to the data's amplitude on one side of a matrix affecting the data's amplitude on the other side. Gibb's phenomenon is an oscillation that appears in the data set after the Fourier transform due to breaks or discontinuities in the data set. These errors can highly alter the results of a filtering routine. A discussion of early methods to limit these sources of error can be found in Ku and others (1971).
The filtered gravity equations assume the data set is continuous and infinite along a surface. However the data set is not continuous; it is known only over a small region at discrete locations. It is possible that certain anomalies are not completely represented in the data set. It should be apparent that wavelengths larger than the region of the data set will not be completely defined. Large wavelength anomalies can be caused by very large near surface structure, but also can be due to very strong (high density contrast) deep structure. Studies have indicated that if the depth to the source of an anomaly is more than 1/6 of the length of the profile, errors may be significant (Regan and Hinze, 1976). One must take this into account if one wishes to analyze deep seated structure.
Gibb's phenomenon is one characteristic of the FFT algorithm that causes error. Any large jump in value from one data point to the next causes a discontinuity and will result in oscillations to appear when the data is transformed from one domain to the other. Jumps in value may be real, such as due to an unusually high gradient in the gravity field, but will still cause ripples to appear in the filtered data set. These ripples will be greatest at the point of discontinuity and taper off to both sides. If the discontinuity is small, the oscillations due to Gibb's phenomenon may be tolerable, they may not even be apparent when the data set is contoured, but at times they can cause a great amount of error.
Another characteristic error that occurs due to the FFT is leakage. If a data set is continuous and infinitly long it can be transformed from the space domain to the frequency domain by multiplying it by the Fourier integral. To transform a discrete data set over a finite region one uses the Fast Fourier Transform algorithm. However, the FFT algorithm makes some assumptions that may affect the trueness of the results. The FFT assumes a data set of limited length represents a summation of an infinitely long series of sine and cosine waves. The finite data set is assumed to be one fundamental period of the infinite series. That is, the data set repeats over and over, making it infinitly long. As a result, one end of the data set is placed against the other end. During filtering, this may cause information from each edge of the map to affect data on the opposite edge. This is called leakage (Dobrin, 1976, Reed, 1980).
Leakage of the Fourier signal from one side of the data set to the other can cause great amounts of error. If the two ends of the data set have greatly different amplitudes, the values on one side of the data set may affect values on the other side. Since in potential field studies it is not known what the data should look like, it is not always possible to tell whether leakage has occured. Various methods have been studied to limit leakage. Sato (1954) suggested entering a region of zeros at the ends of the data set. It was believed that any leakage that did occur would appear in the extended regions. After filtering, the extended region would be discarded. However, adding zeros still caused a discontinuity between the data set and the filled region of zeros. Leakage was limited, but Gibb's phenomenon remained. Tsay (1975) improved on Sato's method by suggesting that the extended region be filled with the end values. Again, after filtering, the extended regions would be discarded. Again, a discontinuity still exists, the value at one side of the data set differs from the value at the other, but the discontinuity has been moved away from the real data. Error caused by Gibb's phenomenom remains, but a buffer has been created to absorb the error. (Although the methods presented by Sato and by Tsay were written to be applied to a Fourier series method, the application to the Fast Fourier Transform is the same.) A new method is proposed to limit errors due to leakage and Gibb's phenomenon. In this method the discontinuity from one side of the data set to the other is eliminated by extending the data set and filling the extended region with a reflection of the data set.
To compare the effects of various methods of limiting errors, the results of a filtering routine is compared with data calcuated in another way. If the gravity due to simple shapes is calculated at a level z using a forward modeling routine, the gravity of the same shapes can also be calculated at level z1. By upward continuing the data set at level z0 to level z1 using the wavenumber filtering method, the filtered data can be compared to the actual data at the same level. The difference between the two is a measure of the error in the wavenumber filtering method. This error is assumed to be due to leakage and Gibb's phenomenom.
The results of such a comparision is presented here. A synthetic gravity profile was created by a two-dimensional forward modeling program based on Cady (1980). The program is listed in Appendix A. The profile consists of irregular bodies of various densities (figure 2). All units of distance are in kilometers, densities are in grams per cubic centimeter and the gravity anomaly is in milligals. The gravity of this profile calculated at z=0 is prepared for filtering in various ways to limit error (figure 3). These data sets are upward continued by the wavenumber filtering method to an elevation of z=2 (figure 4). They are then compared to the gravity anomaly calculated directly at elevation z=2 and the error of each method is shown in figure 5.
As was stated above, the FFT assumes that the data set is periodic. The Fourier transform is calculted as if the last point in the data set is followed by a repeat of the first point. This can cause great discontinuities in the data if the value at one end is greatly different than the value at the other. Transforming such a data set with the FFT will cause a significant amount of leakage and Gibb's phenomenon. If by chance, the gravity value at one end of the data set is close to the value at the other end, these errors will be small. However, if there is any type of regional trend across the data set, it is likely that the difference between the two ends will be great and the errors will be severe.
Part (a) of figures 3, 4 and 5 shows the filtering of a data set without using any method to limit errors. Figure 5a shows that the results of filtering is good over the center of the data set, but errors become severe on both ends. A large amount of leakage has occured. Note that for this example, Gibb's phenomenon does not show up well at this scale.
One method of limiting edge effects is to taper the edges of the data set. This is done by multiplying the data set by a window so that each end of the data set tapers to zero. One such window recommended to reduce leakage is the Hanning window, a type of cosine taper (Brigham, 1974, Blackman and Tukey, 1969). The Hanning window tapers both ends of the data set to zero, eliminating any discontinuities from one side of the data to the other. In Part b of figures 3, 4 and 5, filtering has been performed on a data set that has been tapered with a cosine window. Leakage has been limited, but in doing so, the original data set had to be modified. Therefore, the error in the filtered data has been increased.
Another method to limit leakage and Gibb's phenomenon is to add a buffer of data to both ends of the data set so that any errors will happen away from the original data set (Tsay, 1975). If the data set is extended to twice its original length, with the extended region half filled with one end value, and half filled with the other end value, information on one end of the data set will not be able to leak over to the other end. A discontinuity still exists, but it occurs away from the original data set. It is hoped any errors due to this discontinuity will only occur in the extended region and will not affect the original data. Part c of figures 3, 4 and 5 show filtering of the data set extended with its end vaules. In figure 3c note the original data set is the same length as before, but a buffer of another data length has been added. Each end of the data set is extended one half data length and filled with the end values. The discontinuity should now be far enough away from the original data set so that its effect will be minimun. Figure 4c shows the upward continued gravity of this extended data set and 5c shows the relative error. The errors are almost zero all the way across, but the curve does oscillate. Gibb's phenomenon remains. With some data sets, the error due to Gibb's phenomenon is greater than the error due to leakage. One may suppose that extending the data to greater lengths may decrease the amount of error more, but studies indicate that the decrease is not significantly. Computer memory limitations may cause difficulties with this method. Note that for a one-dimensional profile, extending the data set requires twice as much memory. For a two-dimensional matrix, extending requires four times as much memory.
Part (d) of figures 3, 4 and 5 shows a new method to limit leakage and Gibb's phenomenon. Figure 3d shows the data set again extended, but instead of filling the extended region with the last data value, it is filled with a reflection of the data set. In this manner there are no discontinuities from one side of the data set to the other. Figure 4d shows the upward continuation of this reflected data set and figure 5d shows its relative error. The errors are nearly zero across the data set.
It has been shown that the wavenumber filtering method using the Fast Fourier Transform can give good results if one takes steps to limit errors due to leakage and Gibb's phenomenon. One of the best way to limit these error is to extend the data set one data length, filling the extended region with the data sets end vaules, or with the reflection of the data set.
The wavenumber filtering method using the Fast Fourier Transform allows an easy and rapid way of isolating and enhancing certain aspects of a potential field data set. In the previous section the equations to calculate the Upward and Downward Continuation, the Nth Derivative, as well as Band-pass and Strike-pass filters were derived. These equations can be summarized as follows:
This equation states that the Fourier transform of the filtered data is equal to the Fourier transform of the original data times a filter. The shape of the filter F(p,q) determines what type of filtering is to take place. This section describes the application of the wavenumber filtering method to a computer program using the Fast Fourier Transform.
A set of computer programs have been written to carry out wavenumber filtering of potential field data and are listed in Appendix A. These programs are based on programs in Reed (1980), but have the following improvements.
Reed's program created filters by calculating two quadrants of the filter and then taking the complex conjugate to get the other two quadrants. However, since all of the filters except the odd horizontal derivatives are real and have no imaginary part, they are symmetric (Brigham, 1974). That is the imaginary part of the filter is zero. The complex conjugate of a real number is the same real number. Therefore, one quadrant of the filter matrix is only the reflection of another quadrant and not the complex conjugate. In the new computer program, the method of creating real filters has been simplified. One quadrant is now calculated and then reflected into the other three quadrants.
The filtering of potential field data stored in a one or two dimensional array (representing a profile or a matrix) consists of the same steps, no matter which of the filters is being used. These steps are:
Reed's program includes an option to limit errors due to discontinuities in the ideal filter. In the band-pass and strike-pass filter, there are only two values that each cell in the matrix can have: 1 to pass that wavenumber, 0 to remove it. At the boundaries between regions of 1 and regions of 0 there are discontinuities. To eliminate these discontinuities, Reed's option allowed the ideal filter to be tapered. The ideal filter would be transformed from the frequency domain to the space domain, then multiplied by a cosine taper, and transformed back into the frequency domain. This tapered filter would then be multiplied to the transformed data matrix. Such a tapered filter would look like a bell in cross-section (figure 6b), rather than a box for the ideal filter (figure 6a). Since most of the computer processing time is spent on transforming the data to and from the wavenumber domain, performing two extra FFT's will almost double the CPU time. In addition, the tapered filter no longer represented a true band-pass or strike-pass filter. A more efficient method (Brigham, 1974) is to limit discontinuities directly in the frequency domain by separating the regions of 1 and 0 by an intermediate number, such as 1/2 (figure 6c). In this way, there will be less of a discontinuity from the areas that have been cut to the areas that have been passed. This new method has been included in the construction of the filters.
Since it is likely that the user would want to filter the data set in many different ways (upward or downward continue the data to increasing levels, first and second derivatives, band-pass filtering over different band widths, etc.), it would be redundant to transform the gravity data to the wavenumber domian each time. To avoid unnecessary repetition of calculations, steps 1, 2 and 3 have been put into a separate program called SETUP. Steps 4, 5 and 6 are in a program called FILTER. The program SETUP needs only to be performed once per data set. Then the program FILTER can be performed many times for each type of filtering method.
The program SETUP performs the following functions. A one- or two-dimensional array is loaded from a data file into memory. The data set is extended to limit edge effects, and the mean is removed. The program has the option of extending the data set by either filling the extended region with the end values, or reflecting the data set into the extended region. (Reflecting the data set seems to give the best results.) SETUP then transforms the data set into the wavenumber domain, writing the results to a direct or random access output file. (In Fortran, using direct access files is one of the fastest ways to input and output large amounts of data.)
The program FILTER creates the selected ideal filter matrix, loades the transformed data set (the output from SETUP), and multiplies the two. The program then transforms the results back into the space domain and writes the filtered data set to an output file. The user can then see the results by contouring and plotting the output file onto any of the available plotting devices. Listings of SETUP and FILTER can be found in Appendix A.
There are certain characteristics of the FFT that must be considered when the ideal filter matrix is created. Any program using an FFT must know the structure of that particular FFT routine. Various FFT routines may fill their matrices differently. In a one-dimensional array of Nx pieces of data, the FFT will calculate the values of the wavenmbers from -Nyquist to +Nyquist. However, in most FFT routines, the wavenumbers are stored in the array from 0 to +Nyquist, then -Nyquist to 0. In addition, one cell in the array is reserved for the error determination in the FFT. In a two-dimensional array, the two-dimensional Fourier transform is determined by running the FFT subroutine on each row, storing the results in the original matrix, then performing the FFT on each column (Naidu, 1971). Figure 7 shows the method of storage of the FFT routine used in the programs listed in Appendix A for a two-dimensional array.
All of the filters are functions of the wavenumbers, p and q. Note that for all but the odd horizontal derivatives, the filter equation F(p,q) has no imaginary part. Therefore the filter is symmetric. F(p,q) is equal to F(-p,-q), as well as F(p,-q) and F(-p,q). Each quadrant in the filter matrix is a reflection of the others.
The longest period possible in a continuous and infinite data set is
infinity, corresponding to the wavenumber 0. The longest period in a discrete
finite data set is the length of the data set, corresponding to the wavenumber
1. The shortest period that can be seen in a discrete data set is the Nyquist,
defined as 
where 
is the station spacing. The
Nyquist frequency is the inverse, 
.
Figure 8 shows two-dimensional examples of each type of filter in the
wavenumber domain. Figure 8a is the first
vertical derivative filter. The filter ranges from zero in the regions of low
wavenumbers to 
in regions of high
wavenumbers. Figure 8b is of the second
vertical derivative. It ranges from zero to 
.
When the original data is multiplied by one of these filters, low wavenumbers
are surpressed and high wavenumbers are enhanced. Thus when converted back into
the space domain, the anomaly field will enhance smaller or shallower features
at the expense of larger or deeper features. Figure
8c and 8d are of the first horizontal
derivative in the X and Y directions.
Figure 8e shows the ideal filter for upward
continuation to a level z, where in this example,
, the grid
spacing. The filter ranges from 0 in the region of highest wavenumbers to 1
in the region of lowest wavenumbers. This filter enhances the larger
wavelengths and suppresses the smaller. This corresponds with the concept that
if one were to measure the gravity field a distance z above the surface, the
smaller bodies would be less important and the larger bodies more important.
The ideal filter of upward continuation to a different z would have the same
shape as figure 8e, but the contours would
represent different values. An ideal filter for continuation downward is shown
in figure 8f. This figure looks similar to the
upward continuation filter, but the filter ranges from infinity in the region
of high wavenumbers to 1 in the region of low wavenumbers.
The two remaining type of filters mentioned earlier are band-pass and strike-pass filters. In these type of filters, all wavenumbers that the user wants to keep are set to 1, while all wavenumbers that the user wants to cut are set to 0. Figure 8g shows a band-pass filter where all wavenumber between 1/4 Nyquist and 1/2 Nyquist will be passed; all others will be cut. Figure 8h shows an ideal strike-pass filter where all frequencies between the azimuths of 60 and 120 will be passed. Note that to limit errors due to the discontinuity between the regions of 1 and regions of 0, there is an extra contour representing the region set to 1/2.
In addition to the four filters described here, the original programs of Reed (1980) included another filter routine that performs reduction to the pole for magnetic data. Since the present study has been on gravity data only, this filter has not been included in program FILTER listed in Appendix A.
The computer programs SETUP and FILTER allow quick and efficient implimentation of the of the wavenumber filtering method to virtually any size data set. The programs have been written to take full advantage of the present generation of computers. Once a data matrix has been created, it can be filtered in dozens of ways and the results plotted in a matter of hours. This can be a powerful tool in analyzing a geological region.
