Tutorial : Spectral bandwidth extension — Invention versus harmonic extrapolation

Various postprocessing methods can be applied to seismic data to extend the spectral bandwidth and potentially increase the seismic resolution. Frequency invention techniques, including phase acceleration and loop reconvolution, produce spectrally broadened seismic sections but arbitrarily create high frequencies without a physical basis. Tests in extending the bandwidth of lowfrequency synthetics using these methods indicate that the invented frequencies do not tie high-frequency synthetics generated from the same reflectivity series. Furthermore, synthetic wedge models indicate that the invented high-frequency seismic traces do not improve thin-layer resolution. Frequency invention outputs may serve as useful attributes, but they should not be used for quantitative work and do not improve actual resolution. On the other hand, under appropriate circumstances, layer frequency responses can be extrapolated to frequencies outside the band of the original data using spectral periodicities determined from within the original seismic bandwidth. This can be accomplished by harmonic extrapolation. For blocky earth structures, synthetic tests show that such spectral extrapolation can readily double the bandwidth, even in the presence of noise. Wedge models illustrate the resulting resolution improvement. Synthetic tests suggest that the more complicated the earth structure, the less valid the bandwidth extension that harmonic extrapolation can achieve. Tests of the frequency invention methods and harmonic extrapolation on field seismic data demonstrate that (1) the frequency invention methods modify the original seismic band such that the original data cannot be recovered by simple band-pass filtering, whereas harmonic extrapolation can be filtered back to the original band with good fidelity and (2) harmonic extrapolation exhibits acceptable ties between real and synthetic seismic data outside the original seismic band, whereas frequency invention methods have unfavorable well ties in the cases studied.


INTRODUCTION
Resolution of seismic data is a function of data bandwidth and dominant frequency (Widess, 1973;Kallweit and Wood, 1982).Due to the attenuation of high frequencies during wave propagation in an attenuating earth, seismic frequency content may not be adequate for seismic interpretation purposes in specific cases.A seismically thin layer is commonly defined to be a layer thinner than approximately one-fourth wavelength of the dominant frequency of the data.For such layers, the reflection events of opposite sign from the top and base of the layer cannot be independently mapped, and the time difference between them is generally not sensitive to the actual thickness of a thin layer (Widess, 1973).In practice, because many reservoirs or flow units within reservoirs are seismically thin, seismic resolution improvement is often desirable.The best way to accomplish this is to design better data acquisition and to process the data in such a way as to maximize resolution; however, the desired resolution may not be achievable given practical constraints.Given an existing data set of a given resolution, any postprocessing steps that can be applied to improve resolution are, if valid, desirable.
Resolution improves as the signal dominant frequency is increased and as its bandwidth is widened.However, one can conceive of many ways of data manipulation transforming seismic data to higher frequencies that do not actually improve seismic resolution (e.g., Young and Wild, 2005;Stark, 2009).We refer to such methods that arbitrarily create amplitude at frequencies outside the band of the original data as frequency invention techniques because they are not based on the physics of seismic wave propagation.Chávez-Pérez (2015) discusses the fact that these methods are indeed used in the industry and uses canonical examples to show their lack of efficacy.
That bandwidth extension is possible under appropriate circumstances should not be a matter of debate: The only question should be the general applicability of a given method.Spectral extrapolation is recognized in telecommunications as a viable means of recovering high-frequency speech components that are not transmitted (see, e.g., Avendano et al., 1995;Pyke, 1997).In the case of seismic data, there is a strong link between the amplitudes of various frequency components of the reflectivity spectrum for individual layers out to the Nyquist frequency (see, e.g., Partyka et al., 1999;Marfurt and Kirlin, 2001).Consequently, missing frequencies in the data spectrum can potentially be extrapolated by applying predictive operators directly in the frequency domain.Oldenburg et al. (1983) state "the Fourier transform of a reflectivity function for a layered earth can be modeled as an autoregressive (AR) process.The missing high and low frequencies can thus be predicted from the band-limited reflectivity function by standard techniques."Walker and Ulrych (1983) deal with novel extensions to the AR algorithm and demonstrate the efficacy of the method for sparse reflectivity series.In this tutorial, we discuss the use of sparse seismic inversion to achieve the bandwidth extension.
In the special case in which the earth model contributing to a given seismic reflector is caused by a single layer, Puryear and Castagna (2008) show that, with a known wavelet and in the absence of noise, the reflectivity spectrum associated with that layer, which is periodic in frequency, can be determined.Thus, the frequency response of that layer reflectivity can be extrapolated outside of the original band.Other simple layer models can also be readily extrapolated if there is enough bandwidth in the recorded signal to properly characterize the layering.Puryear and Castagna (2008) and Zhang and Castagna (2011) show how sparse inversions can be used to produce bandwidth-extended seismic sections, even when layers are not isolated, by solving for many superposing layer responses simultaneously.Hargreaves et al. (2013) explain that all sparse inversion methods can extend the seismic bandwidth to some extent.Zhang and Castagna (2011) show that sparse layer inversion can resolve thin layers better than sparse spike inversion (Taylor et al., 1979;Oldenburg et al., 1983) and should be more effective for spectral extension.Spectra can also be extended, with some reported success, directly by extrapolating time-frequency analyses produced by spectral decomposition (Smith et al., 2008).That bandwidth extension can be achieved should not be surprising; what may be unexpected is how useful such extrapolated frequencies can be.Bandwidth extension does not increase the information content of the data beyond the original information of the data plus any information resulting from the assumptions made to accomplish that extrapolation (such as assuming a blocky earth model).Bandwidth extension should thus be considered, not a means of recovering lost signal, but as a means of transforming the data to have useful characteristicssuch as possibly improved resolution.
When confined to digital filtering (without invoking sparsity constraints), once a given frequency has been filtered away, no linear digital filter can restore it.Thus, where earth filtering has attenuated away high-frequency information, those frequencies can never be recovered.Deconvolution, for example, will just blow up noise at such missing frequencies.This inverse digital filtering paradigm somehow leads to the misconception that prediction of those missing frequencies cannot be validly performed using any method.However, since the inception of sparse inversion methods, it has been well-known that inverted earth models contain frequencies not present in the original data (e.g., van Riel and Berkhout, 1985).In seismic inversion, it is well-established that a priori information and constraints can restore missing frequencies (Tarantola, 2005).The issue then becomes how adequate, how correct, and how useful the constraints are in the presence of noise.
A more damaging contribution to the confusion associated with bandwidth extension is the practice of frequency invention.Because frequencies outside the band of the seismic data are filtered away by the seismic wavelet, any frequency content can be added to this null space and still be compatible with the original seismic trace.Various schemes can be imagined that create high-frequency data that track seismic reflectors.The result is a high-frequency section that may look geologic, but in which the high frequencies amount to no more than fictitious coherent noise.We will describe two representative frequency invention methods that are common practice, including cosmetic loop reconvolution (Young and Wild, 2005) and phase acceleration (Mallat, 2009;Stark, 2009).These methods invent frequencies outside the band of the data without a physical basis.We will show using synthetic wedge models that such methods do not improve resolution.As pointed out by Wild and Young (2005), high-frequency seismic sections produced by frequency invention may serve as a useful seismic attribute for certain interpretational purposes, but we will show in this tutorial that such outputs should not be wrongly attributed with the properties of those produced by true bandwidth extension techniques.
The purpose of this tutorial paper is to help reduce the confusion relating to these issues to avoid improper use or avoidance of bandwidth extension techniques.We illustrate extension of the data bandwidth using frequency periodic layer responses; we refer to this as harmonic extrapolation and show that it can improve seismic resolution under certain conditions.The broadband reflectivity series for a blocky earth model can be obtained using a variety of approaches (Taylor et al., 1979;Puryear and Castagna, 2008;Zhang and Castagna, 2011) other than the algorithm we will describe here.We do not claim that our harmonic extrapolation method is the optimal means of inverting for the layer model, only that it is a representative way of accomplishing the bandwidth extension.We will demonstrate the applicability of harmonic extrapolation with a series of comparative studies using synthetic data, which is the only circumstance in which the "correct answer" is perfectly known.We start with low-frequency synthetic data, use various bandwidth extension methods, and compare to perfect high-frequency synthetics.Resolution improvement, or lack thereof, will be demonstrated on synthetic wedge models.These synthetic tests will address the theoretical validity of the methods in a perfect world; real-world application in the presence of imperfect data is, of course, another matter.We will test the efficacy of the methods on synthetics with added noise and on real data.We will show (1) that frequency invention does not increase seismic resolution and (2) that under the right circumstances, valid spectral broadening can be achieved to some extent using harmonic extrapolation.

THEORY AND METHODS
We assume that the processed seismic data can be represented by the simple convolutional model, where the seismogram sðtÞ is the convolution of a constant seismic wavelet wðtÞ, with a reflectivity series rðtÞ plus random noise nðtÞ: sðtÞ ¼ wðtÞ Ã rðtÞ þ nðtÞ: (1) This model implicitly assumes that the subsurface geology to be horizontally layered with constant rock properties within the layers, whereas the reflections are generated at the seismic impedance contrast boundaries between those adjacent layers and ignores many wave-propagation effects.These are assumed to be corrected in processing.The convolutional model will be used to generate synthetic traces and to invert seismic traces into band-limited reflectivity.Any real-world deviations from this model are implicitly contained in the noise term.

Frequency invention methods
Some of the inspiration for seismic bandwidth broadening methods stems from basic Fourier theorems (Bracewell, 1965) and related variants commonly seen in the field of electrical engineering.For example, the Fourier shift theorem shows that a given signal can be converted to a higher frequency signal by multiplying the signal with high-frequency sinusoids.The resulting higher frequencies can be added back to the original signal to produce a broader band waveform.Whether any apparently broader spectral band being created in this way, irrespective of the underlying physics, is theoretically valid in accordance with the convolutional model should be questioned.
One such typical method called phase acceleration produces high frequencies using instantaneous frequency (Stark, 2009).Analytical instantaneous frequency is a useful attribute related to rock and fluid properties as well as a measure for evaluating layer thickness (Barnes, 2007).In audio processing, it is a common practice to transpose the signal frequency by scaling the phase to shift the harmonics and modify the sound properties (Mallat, 2009).The phase acceleration technique uses the instantaneous frequency of the original data as a fundamental frequency and increases the seismic frequency by increasing the rate of change of the instantaneous phase.
Based on complex trace analysis, the narrowband seismic signal can be viewed as the real component of a complex number representation (called the analytic trace), which provides a definitive separation for amplitude of the trace envelope and local phase content (Taner et al., 1979).This representation permits explicit calculation and modification of the instantaneous frequency, which is the first derivative of the instantaneous phase in the time domain.The analytical signal SðtÞ, of which the seismogram sðtÞ, is the real part is given by where φðtÞ is the instantaneous phase, AðtÞ is the instantaneous amplitude, and s ⊥ ðtÞ is the quadrature series (Hilbert transform) of the real seismic trace.Barnes (2007) shows that instantaneous frequency at the peak of the envelope is the dominant frequency of the data.Consequently, the instantaneous frequency and the dominant frequency of the original trace can be directly changed by multiplying the phase term by any desired amount.The frequency extrapolated trace s E ðtÞ can then be calculated using s E ðtÞ ¼ AðtÞ X k i¼1 cos m i φðtÞ; (3) where m i can be varied to sum an arbitrary number of invented bands to the original data.For instance, a 2 ms sampling seismic trace with a central frequency of 30 Hz can be broadened to a wider bandwidth with the factors m i (k ¼ 7) taking the integers 0, 1, 2, 3, 4, 5, 6, 7, of which m 1 ¼ 1 denotes the original trace.In the case of an invariant instantaneous phase, the summation will produce beating sinusoids.
Because the beating sinusoids have their own envelope, with maxima spacing determined by the beat frequency, the resulting extrapolated trace may result in events with modified envelope.In the case of application to real data, the duration of events may thus be shorter and appear to have better resolution.Figure 1a-1h shows the application of phase acceleration to a 30 Hz Ricker wavelet.The phase is accelerated by factors of two and three yielding higher frequency wavelets (Figure 1b and 1c).Notice that these wavelets have a higher center frequency, but they have the same envelope (dashed lines) and bandwidth as the original waveform (Figure 1e-1g).These are summed in equation 3 with the original signal to produce a waveform with a broad frequency response, high center frequency, and short envelope width (Figure 1d and 1h).
In the loop reconvolution method (Young and Wild, 2005), a spike representation of the seismic data is created first, and then Figure 1.Phase acceleration on an isolated wavelet (30 Hz Ricker) with envelopes (dashed lines) overlaid.(a) Original wavelet, (b) phase accelerated wavelet with factor of two, (c) phase accelerated wavelet with factor of three, (d) high-frequency wavelet by summation of phase accelerated results, (e) spectrum of the original wavelet; (f) spectrum of the phase accelerated result with factor of two, (g) spectrum of the phase accelerated result with factor of three, and (h) spectrum of the resulting bandwidth-extended wavelet.
the resulting broadband spike series is convolved with a high-frequency wavelet (Figure 2a-2c).Young and Wild (2005) oversample the input data to avoid aliasing and zero out all samples of the seismic trace other than samples representing the original peaks and troughs.This produces a "spike train"; unfortunately, this is sometimes referred to by practitioners as a sparse reflectivity, which incorrectly implies that the spikes represent reflection coefficients.In no case, even with only isolated reflectors, is this spike train representative of the true reflectivity because all wavelet side lobes also appear as spikes, which are always located at waveform peaks and troughs (referred to as the "loops" of the waveform).The "high-frequency" seismic trace is achieved by convolving the spike series with an arbitrary high-frequency wavelet of choice.This method produces output that is correlated with the original peaks and troughs, which can lead to a misleadingly high correlation to synthetic data, if the original well tie is also good.This may occur for isolated or major events, for which peaks and troughs of the high-frequency synthetic remain aligned with the loop reconvolution result if the low-frequency synthetic and input seismic data originally had aligned peaks and troughs.
These frequency invention methods are fast and relatively easy to implement.They may have utility as seismic attributesfor example, faults with throw less than one-half the period of the original data may visually become clearer because peaks may then be juxtaposed against troughs on the higher frequency outputs (Young and Wild, 2005).A purpose of this tutorial is to show that, although producing potentially useful attributes, the value of frequency invention to increase resolution or for quantitative analysis is highly questionable at best.

Harmonic extrapolation method
Bandwidth extension can be achieved by extrapolation of the spectrum outside the original sampled seismic bandwidth.We rely on the fact that any transient signal has an infinite frequency response that is potentially predictable to some extent if a sufficient portion of the spectrum is sampled.Our ability to do this of course depends on how simple and characterizable that spectrum is, how clear the frequency periodicities are, the original bandwidth of the data, and the noise level.
To the extent that the convolutional model adequately simulates seismic wave propagation and data processing, we can view seismic data as a band-limited version of the reflection coefficient time series; the reflectivity spectrum within the data band is shaped by the wavelet spectrum, whereas the frequencies outside the band of the wavelet are zeroed out, causing the ambiguity in resolving seismic reflections.A blocky earth structure provides a physical basis for a valid bandwidth extension as the reflectivity spectrum is a superposition of sinusoidal layer frequency responses (Puryear and Castagna, 2008).
As illustrated by Partyka et al. (1999), and quantitatively studied by Marfurt and Kirlin (2001), any layer with a single reflection coefficient at the top and base can be represented as an impulse pair reflectivity series.As discussed below, one can model different types of layers using weighted sums of the even and odd impulse pair symbols (Bracewell, 1965) defined by and where δðtÞ is the Dirac delta function referred to as a unit impulse at time zero in the context of signal processing.A layer with time thickness Δt and reflection coefficients r of equal magnitude and sign at the top and base, can be represented by an even impulse pair rIIðt∕ΔtÞ, with the reflectivity spectrum given by a real function: which is periodic in frequency.For odd impulse pairs, the spectrum is an imaginary sine function:  I I ðfÞ ¼ i2r sinð πΔtfÞ: (7) Figure 3 shows the reflectivity amplitude spectrum of a 50 ms thick layer consisting of two reflection coefficients with equal unit magnitude, but with opposite sign.The hypothetical seismic bandwidth, bandlimited by the seismic wavelet, is denoted by dotted lines.In practice, because the convolution of the reflectivity with the seismic wavelet in the time domain is equivalent to multiplication with the wavelet spectrum in the frequency domain, the reflectivity spectrum is obtained over the band of the data by dividing the seismic data spectrum by the wavelet spectrum, which is assumed to be known in this discussion and must be well-determined a priori in practice.The amplitudes of the Fourier frequencies of the reflectivity spectrum obtained over the seismic band in this way are indicated by solid circles in Figure 3.In the case of a single layer, these amplitudes would readily reveal the sinusoidal shape of the reflectivity spectrum and could thus be fit with a best least mean-squareerror (LMSE) fit (solid line in Figure 3) to the derived reflectivity spectrum.For this simple case, the frequency spectrum outside the band of the data can be readily predicted by extrapolating the resulting LMSE sinusoidal function outside the band of the original data.
The important reciprocal relationship between the spectral periodicity and the time thickness for the real and imaginary components (equations 6 and 7) is dependent upon the symmetrical location of time zero at the center of the layer.In practice, when composite reflections occur, although it is not practical to define a time zero that can be centered for all such physical layers, any reflectivity series can still be represented as a sum of even and odd impulse pairs referenced to an analysis point.Therefore, all reflectivity spectra are sums of such sinusoidal frequency responses, and they are fundamentally predictable if those impulse pairs are known.To clarify what happens in a general multilayer case, we construct a sparse reflectivity sequence (Figure 4a) composed of three individual components: a predominantly even impulse pair (Figure 4b), a predominately odd impulse pair (Figure 4c), and a single impulse (Figure 4d).Because any pair of reflection coefficients can be uniquely represented by the weighted sum of an even pair and an odd pair (Bracewell, 1965), the complex spectrum of each component can be characterized by a pure cosine function (equation 6) as the real part (solid line) and a pure sine function (equation 7) as the imaginary part (dashed line).In this example, the predominantly even reflector pair has a real spectrum with larger amplitude than the imaginary spectrum (Figure 4f) and vice versa for the predominantly odd pair (Figure 4g), whereas those for the single reflector are of the same amplitude (Figure 4h).The even and odd spectra for each individual event show the same frequency periodicity, but they are shifted by one-quarter of the spectral period.The composite reflectivity spectrum (Figure 4e) is a linear combination of those component spectral sinusoids, producing irregular peaks and troughs by local superposition.The frequencies within the hypothetical data bandwidth (denoted with red) can then be used to determine the individual sinusoidal components of the spectrum (the red curves in Figure 4f-4h), and these can then be extended beyond the original bandwidth (the black curves in Figure 4.A multilayered reflectivity series centered at t ¼ 0 with the corresponding complex spectrum.Real and imaginary parts of the spectrum are denoted by solid line and dashed line, respectively.The red portion hypothesizes the available seismic data band (the wavelet has been divided out).(a) Composite reflectivity sequence, (b) predominantly even impulse pair, (c) predominantly odd impulse pair, (d) single-impulse component, (e) complex spectrum of the composite reflectivity, (f) complex spectrum of the predominantly even impulse pair, (g) complex spectrum of the predominately odd impulse pair, and (h) complex spectrum of the single impulse.The vertical axis is the relative amplitude.Figure 4f-4h).These components are then added to obtain the bandwidth-extended spectrum shown in the black curve in Figure 4e.

Spectral bandwidth extension
Beginning with a temporal analysis window of 2N þ 1 discrete points with a sampling rate of dt, any zero-time centered impulse pair gðt; nÞ can be expressed as where n • dt is the half the time thickness of the impulse pair with n varying from 0 to N. The reflector pair can be decomposed into a unit even pair IIðt∕2n • dtÞ and a unit odd pair I I ðt∕2n • dtÞ weighted by coefficients r e and r o , respectively.The two reflection coefficients for the impulse pair can thus be represented as r e þ r o and r e − r o , either of which ranges from −1 to 1 and could be zero denoting a single impulse.Correspondingly, the frequency spectrum of the impulse pair is given by (equations 6 and 7) Gðf;nÞ 1 nested pairs of reflection coefficients.A time sample can have zero reflectivity if r e and r o are zero, but also if elements of the even and odd reflectivity pairs are equal and opposite at a given time.Thus, K is the number of impulse pairs, including all possible time thicknesses (including a zero time thickness) in the analysis window.The broadband spectrum ranging from zero to the Nyquist frequency of any reflectivity series within the analysis window can thus be explicitly expressed as S r ðfÞ ¼ For a real seismic trace, the narrowband data spectrum associated with a sum of reflector pairs is a sum of frequency sinusoids times the wavelet spectrum.Within a usable spectral band of seismic data with an adequate signal-to-noise ratio (S/N), dividing out the wavelet overprint provides a roughly flattened spectrum that coincides with the portion of the broadband reflectivity spectrum (equation 10) within the wavelet spectral limit (the red portion in Figure 4).The relationship between the data spectral constraint (real part and imaginary part separately) and the sinusoidal elements of varying frequency periodicities within the same spectral limit can be formulated as the forward modeling: where S d is the wavelet-flattened data spectrum, Φ is the kernel matrix consisting of sinusoidal atoms, a is the coefficients vector containing r o ðnÞ and r e ðnÞ, and ε is the prediction residual.The kernel has a dimension of M × K, where M is the number of frequencies determined by the available bandwidth and the sampling rate in frequency, K is the number of layers and thus also the number of cosine (real part) or sine (imaginary part) elements indicated by equation 10, which are varied to include all possible frequency periodicities.
We use the basis-pursuit method (Chen et al., 2001) to thus decompose the original band-limited data spectrum into a sparse superposition of layer responses.The frequency extrapolation is accomplished by decomposing the spectrum into sinusoidal frequency responses (atoms), extending them beyond the band of the seismic data, and then summing them to form the broadband seismic trace.The basis pursuit method solves for the coefficients for all sinusoidal atoms in frequency in a direct manner by L 1 -norm global optimization.Zhang and Castagna (2011) describe how the Chen et al. (2001) basis pursuit method can be used to decompose layer responses in the time domain.We use the same algorithm to do the decomposition in the frequency domain.
Basis pursuit solves for the coefficients of superposed sinusoidal frequency responses in equation 11 by minimizing the L 2 -norm of prediction residual term regularized by the L 1 -norm of the solution scaled by a regularization (trade-off) parameter λ: Increasing this trade-off factor λ produces results with higher sparsity, whereas decreasing it may amplify the inversion noise.In real applications, a proper regularization factor for the inversion is data and model dependent and empirically determined by comparing results with well logs used for validation.
In our implementation, the coefficients r e ðnÞ and r o ðnÞ are the output solutions of the basis pursuit inversion, corresponding to the real part and the imaginary part of the available spectrum, respectively.Any missing frequency outside the original bandwidth can then be computed using equation 10.Because equations 8 and 9 share the same coefficients, the resulting reflectivity series can thus be calculated by s r ðtÞ ¼ þ ðr e ðnÞ − r o ðnÞÞδðt − ndtÞÞ: (13) It is important to note that, theoretically, there can be frequency components of the reflectivity spectrum that are entirely outside the band of the original data that violate our assumptions, and they can never be recovered.The extrapolated frequencies are correct to the extent that the layer model used to generate them is a good representation of the true subsurface reflectivity; any deviation from a sparse blocky earth model assumption may give rise to layers not properly constrained by the original band of the data.For example, layers that are too thin or are of high complexity are less likely to be resolved by the layer inversion.Therefore, a successful extrapolation of the complete frequency responses depends on the reasonableness and accuracy of the decomposition of the complex spectrum into the temporal impulse pair frequency spectra.

Practical implementation of the harmonic extrapolation procedure
These steps are followed in the implementation of harmonic extrapolation: 1) Check the validity of the method assumptions/constraints: Is blocky layering an adequate description of the true earth structure?Is the wavelet well-known?Does convolutional modeling provide a good synthetic tie within the original band of the data?If no well logs are available to ensure that the assumptions are applicable, that the wavelet is correct, and to validate the extended frequencies, the extension will have to be interpreted with great caution.For example, should the actual wavelet spectrum have notches, or be of a narrower band than assumed, the result will be noise magnification at certain frequencies.In addition, uncharacterized periodicities in the wavelet spectrum could result in the inversion producing false layering to reproduce those periodicities.2) Check seismic data quality: Is the seismic data quality sufficient such that there is at least a one-octave passband with a high S/N?Again, this is best checked by synthetic seismogram construction.3) Set parameters: Is the sample rate sufficient for extension objectives?If not, resample the original data by Fourier transform to increase the sampling frequency to an acceptable level relative to the desired output frequency band.Select a regularization parameter.For noise-free cases, we use 0.0001.For noisy cases, a regularization parameter on the order of 0.01 is usually acceptable (basis pursuit denoising).This is empirically determined by a best correlation between bandwidth-extended and synthetic data in blind well validation tests.Select a temporal window length.4) Compute spectra and deconvolve the wavelet: For each seismic trace, slide the analysis window along the trace and compute the Fourier transform at every window location.Divide the spectra by the wavelet spectrum (the wavelet spectra can be time varying if necessary).The real part of the resulting spectrum is the spectrum of the even (symmetrical) part of the windowed trace divided by the wavelet spectrum, and it is the superposition of different cosine functions with different periodicities, each representing an even reflection coefficient impulse pair of a given thickness.The imaginary part of the resulting spectrum is the spectrum of the odd (antisymmetrical) part of the windowed trace divided by the wavelet spectrum, and it is the superposition of different sine functions with different periodicities, each representing an odd reflection coefficient impulse pair of a given thickness.5) Decompose local spectra into a superposition of layer spectra: Use a basis-pursuit atomic decomposition algorithm (Chen et al., 2001; for further discussion, see Zhang and Castagna, 2011) to decompose the real and imaginary parts of the spectrum into summations of cosines and sines.At this point, a regularization parameter is required to weight the sparsity constraint (the selection criteria are discussed above).6) Extend the bandwidth of the data: Extrapolate the decomposed sinusoids to the band of interest, multiply by the desired output wavelet spectrum, and inverse Fourier transform to get a time series associated with each window position.Sum the time series over all window positions to form the bandwidth-extended trace.
The windows are tapered such that summing in the overlap region conserves energy.
We will test the efficacy of the harmonic extrapolation synthetically on blocky earth structures, with and without noise.We will see that resolution enhancement is readily demonstrated on wedge models.However, because the method relies on the simplicity and periodicity of the reflectivity spectrum, we also test the method for a highly interbedded earth model with locally highly complex spec-tra.The ability of harmonic extrapolation to achieve useful bandwidth extension for a blocky earth model, and its superiority over frequency invention methods, will be demonstrated in the following comparative studies (synthetic and real data examples).

SYNTHETIC EXAMPLES
The theoretical validity of the various spectral extension methods in increasing the seismic resolution can be reasonably demonstrated on synthetic tests, whereby a low-frequency synthetic is bandwidth extended and compared with a high-frequency synthetic derived from the same earth-reflectivity model.The resulting bandwidthextended data can be visually compared with known variations in layer thickness to see if additional resolution is achieved.Because a perfect bandwidth extension should predict the high-frequency synthetic perfectly, the residual between calculated and predicted high-frequency seismograms provides a firm test of the efficacy of the various methods and indicates the magnitude of the errors involved.All the traces are plotted with the same relative amplitude scale to provide reasonable comparisons and reported root-meansquare (rms) errors are relative to the rms amplitude of the highfrequency target synthetic.

Wedge models
Wedge models are commonly used to study resolution as a layer is thinned from above tuning to zero thickness.Due to the finite data bandwidth, it will not be possible for all sinusoidal atoms of the decomposition to be completely uncorrelated.The presence of a correlation between these atoms over the limited seismic band is an inherent difficulty in determining the exact solution for all reflection coefficients.Figure 5 illustrates the inherent resolving limitation of harmonic extrapolation by a thinning wedge model with equal and opposite reflection coefficients at the top and base.The wedge layer time thicknesses vary from 50 to 1 ms with an interval of 1 ms, yielding the amplitude spectra for each trace (Figure 5a) showing the gradually variable frequency periodicities.The resulting spectra (Figure 5b) extrapolated from the response of this model to a 30 Hz wavelet indicate that although the method can resolve layers below the original tuning thickness, which is approximately 13 ms in this case (Chung and Lawton, 1995), the inverted and extrapolated spectra deteriorate at thicknesses less than 7 ms.This ambiguity arises from the indistinct frequency periodicities produced by the closely spaced reflectors.Obviously, data distortion due to noise or incorrect wavelet determination also interferes with the decomposition of the contributing frequency periodicity to the total spectrum for a given event within the usable band of the data.These factors limit the resolution that can be achieved by harmonic extrapolation.
Figure 6 shows a wedge model with equal reflection coefficients of the same sign at the top and base of the layer.The time thicknesses of the reflector pairs for each synthetic trace are varied from 25 to 0.5 ms to represent a wedge thinning, with an interval of 0.5 ms, which is the sample rate of the data.A 30 Hz Ricker wavelet is used to model the low-frequency response of the wedge (Figure 6a).The top and bottom reflectors can be visually distinguished at or above a time thickness of 12.5 ms at this wavelet center frequency, using the Ricker criterion (flat topping of the waveform).The phase acceleration, loop reconvolution, and harmonic extrapolation methods are applied to the original low-frequency synthetic to restore the missing high frequencies.In practice, doubling the original seismic center frequency is regarded as a significant improvement in resolution.All the broadband frequency extended outputs are thus spectrally shaped to the spectrum of a 60 Hz standard Ricker wavelet by multiplying by the wavelet spectrum.In this case, the low frequencies in the original signal are reduced in amplitude, but are still preserved, and thus present to be shaped back to the original spectrum.The phase acceleration and loop reconvolution methods produce sharper events (Figure 6b and 6c), but neither improve the temporal resolution (defined as the ability to distinguish separate events from the top and base reflectors) as compared with the original synthetic.Furthermore, the loop reconvolution method produces false events above and below the actual reflectors (Figure 6c).Harmonic extrapolation to 60 Hz, on the other hand, resolves the reflections at approximately 7 ms, not quite a twofold improvement, as well as sharpening the events (Figure 6d).The resulting harmonic extrapolation wedge is very similar to, though not exactly the same as, the true high-frequency synthetic target (Figure 7a), calculated directly with the known reflectivity coefficients and a 60 Hz Ricker wavelet.
It is further revealing to look at the prediction residuals by various methods.The harmonic extrapolation computes the high-frequency section with reasonable fidelity with only a relative rms error of 9.5% (Figure 7d).The fact that harmonic extrapolation passes this simple test is to be expected.However, failing the test, as phase acceleration (an rms error of 87.5% in Figure 7b) and loop reconvolution (an rms error of 80.8% in Figure 7c) do, is significant.This shows that these methods must be used with great caution, and not be used for quantitative work.

Blocky structured sparse reflectivity series
A more realistic, but still blocky, synthetic reflectivity series designed to simulate an earth model with multiple layers (Figure 8a-8l), has a sampling rate of 2 ms, resulting in a Nyquist frequency of 250 Hz (Figure 9a-9l).A low-frequency synthetic seismogram (Figure 8b) is generated by convolving the reflectivity series with a 30 Hz Ricker wavelet.A noisy low-frequency synthetic seismogram (Figure 8h) is obtained by adding a band-limited random noise series with signal power that is 10% of that of the noise-free synthetic seismogram.A higher frequency band of the reflectivity series (Figure 8g) is obtained by convolving with a 60 Hz Ricker wavelet, and it is used as a standard against which predicted high-frequency data will be evaluated.
In the frequency domain, the phase acceleration and loop reconvolution can produce broader spectra that are comparable with the high-frequency target (Figure 9g) for noise-free (Figure 9c and 9e) and noisy (Figure 9i and 9k) data.However, even without noise, the misfits (Figure 8d and 8f) between the synthetic and extrapolated high-frequency traces are of the same order of magnitude as the original data.These large errors would hinder the reliability of subsequent quantitative seismic trace analysis (such as seismic inversion or multiattribute analysis).The frequency invention methods, though not very reliable, are relatively immune to noise because they are driven primarily by the location of events, which is for the most part the same on the original and bandwidth-extended data.If the noise is not of a sufficient level to alter the event location, it does not have a great influence on the residual, which is dominated by  how colocated and similar in shape the events are relative to the synthetic.
Using the same band-limited low-frequency synthetic seismograms (Figure 10a and 10g), the broadband synthetic (Figure 10b) with the 60 Hz Ricker wavelet is recovered by harmonic extrapolation with good fidelity (Figure 10c) and only a 6.3% rms error (Figure 10d) for the noise-free data.This noise-free extrapolation error arises from the nonuniqueness in solving for the frequency periodicities within the limited bandwidth.For the noisy result (Figure 10e), the residual of prediction to 60 Hz increases significantly to a 32.5% rms error (Figure 10f) mainly due to the error in the extrapolation caused by the low-frequency noise in the usable band.However, the increment of a 26.2% rms error suggests that the original noise (a power ratio of 10% corresponds to a relative rms of 31.6%) has been somewhat suppressed by a relatively large regularization parameter in the basis-pursuit inversion.More importantly, the fact that the prediction error is as randomly distributed because the original noise implies that harmonic extrapolation does perform robustly against data distortion when the useful signal is dominant.
Extrapolating to higher frequencies, a broader band synthetic (Figure 10h) using a 90 Hz Ricker wavelet is also well-reproduced (Figure 10i) with good, though somewhat reduced, fidelity that is of 14.3% rms error (Figure 10j).The reason for a fidelity reduction from 60 to 90 Hz in the noise-free case is that for an imperfectly matched frequency periodicity determined from the usable bandwidth, the mismatch increases at frequencies further away from  the calibration frequencies.The prediction residual (36.9% rms error in Figure 10l) for this 90 Hz extrapolated result in the presence of noise (Figure 10k) is thus reasonably larger than that of the corresponding 60 Hz case.However, the percentage fidelity reduction in high-frequency recoveries from 60 to 90 Hz is less in the noisy case than in the noise-free case because the error due to noise is added to the inversion error and dominates the fidelity percentage.The corresponding amplitude spectra (Figure 11a-11l) show that the original bandwidth can be at least doubled with minor residuals.
Figure 9. Corresponding amplitude spectra for the test results on the blocky sparse-reflectivity series by frequency invention methods.(a) Spectrum of original reflectivity, (b) spectrum of low-frequency (30 Hz) synthetic, (c) spectrum of phase acceleration result (60 Hz), (d) spectrum of residual for phase acceleration, (e) spectrum of loop reconvolution result (60 Hz), (f) spectrum of residual for loop reconvolution, (g) spectrum of high-frequency (60 Hz) target synthetic, (h) spectrum of low-frequency (30 Hz) synthetic with 10% noise added, (i) spectrum of phase acceleration result (60 Hz) from the noisy data, (j) spectrum of the residual for phase acceleration (noisy), (k) spectrum of loop reconvolution result (60 Hz) from the noisy data, and (l) spectrum of the residual for loop reconvolution (noisy).
Figure 12 shows what happens if recovering the reflectivity spectrum out to the Nyquist frequency is attempted.The full band spectrum (Figure 12c) extrapolated from the low-frequency noisy synthetic by harmonic extrapolation and the corresponding inverted reflectivity (Figure 12d) show good recovery for the original reflectivity series (Figure 12a and 12b), giving a prediction residual (Figure 12e and 12f) with an rms error of 41.7%.This is reasonably larger than that of the corresponding 60 and 90 Hz broadband cases.The residual spectrum shows a roughly increasing trend from low to high frequency.
These synthetic tests suggest that the harmonic extrapolation method can stably predict the missing frequencies under the circumstance of a blocky earth model providing the noise level is small and the bandwidth extension is limited.

Random numbers series for harmonic extrapolation
These favorable harmonic extrapolation results are not unexpected for a blocky earth model in which we observe readily characterizable frequency periodicities.We construct, as a worst-case Figure 11.Corresponding amplitude spectra for the test results on the blocky sparse reflectivity series by the harmonic extrapolation method.(a) Spectrum of low-frequency (30 Hz) synthetic, (b) spectrum of high-frequency (60 Hz) target synthetic, (c) spectrum of harmonic extrapolation result (60 Hz) in the noise-free case, (d) spectrum of residual for the noise-free case (60 Hz), (e) spectrum of harmonic extrapolation result (60 Hz) in the noisy case, (f) spectrum of residual for the noisy case (60 Hz), (g) spectrum of low-frequency (30 Hz) synthetic with 10% noise added, (h) spectrum of high-frequency (90 Hz) target synthetic, (i) spectrum of harmonic extrapolation result (90 Hz) in the noise-free case, (j) spectrum of residual for the noise-free case (90 Hz), (k) spectrum of harmonic extrapolation result (90 Hz) in the noisy case, and (l) spectrum of residual for noisy case (90 Hz).Spectral bandwidth extension W11 scenario (Figure 13a-13h), an earth model with completely random numbers (sampling rate is 2 ms).The sparse structure representation is then not appropriate because the low-frequency seismogram (30 Hz) is then a superposed response from too many very thin interbedded layers to separate (Figure 13e).The high-frequency recoveries by spectral extrapolation in this case are not acceptable; the error (Figure 13d) in recovery of the 60 Hz result is nearly of the same magnitude (an rms error of 68.4%) as the high-frequency synthetic (Figure 13b) and it greatly increases (an rms error of 110.7%) for the 90 Hz result (Figure 13h).
Inspection of the corresponding error spectra (Figure 14a-14h) reveals an interesting fact: The error is low at frequencies at which there is any amplitude in the original spectrum, and it blows up immediately above the highest frequency in the input wavelet (Figure 14d and 14h).This good recovery of the reflectivity spectrum within the data spectral limit is a great virtue for seismic interpretation because the observed behavior of measured data from real wells has shown the reflectivity spectra to be of "blue" color in a global sense, meaning that spectral amplitudes tend to increase with frequency (Walden and Hosken, 1985).Therefore, bluing is an important process to accentuate high frequencies such that the wavelet overprinted seismic spectrum can be similar in shape to the reflectivity spectra.The random number series test suggests that, unlike seismic deconvolution that boosts noise and useful signal by the same factor at any given frequency, harmonic extrapolation, although unable to extrapolate outside the original band in this case, could be useful in boosting the high-frequency signal more than the high-frequency noise in spectral bluing, thus improving the seismic resolution within the original band of the data.

REAL DATA RESULTS
Unlike the synthetic examples in which the original earth model is known and broadband objective reflections are readily simulated for bandwidth extension evaluation, the complete frequency components of the actual reflectivity series in the real world can only be partially represented by a limited number of available wells, which cannot be assumed to be entirely accurate.In this sense, various band-pass versions of the reflectivity series derived from well logs can be used as criteria, only to a certain extent, to determine the validity of newly generated frequency components by those spectral extension techniques.We require a data set with an initial good well tie because any small error in time-depth conversion will become more significant at higher frequencies.Because the initial time-depth function deemed adequate at low frequency will also be used at high frequency, some correlation reduction for this reason alone is to be expected.The Vietnam data set studied by Ha (2014) met our criterion of a good low-frequency well tie (correlation coefficient = 0.87; see Figure 15 obtained using the Hampson-Russell synthetic modeling package, 2008).The wavelet (Figure 16) was extracted using the well-log reflectivity and assuming a constant phase as a function of frequency.The resulting extracted wavelet was nearly zero phase.Figure 17 shows the results of applying frequency inventions and harmonic extrapolation to a seismic line through the well location.The original band of the seismic section (Figure 17a) is roughly included in the band pass 0-10-65-70 Hz, and the spectrally broadened results (see Figure 17b for phase acceleration, Figure 17c for loop reconvolution, and Figure 17d for harmonic extrapolation) are filtered to the same broadband 0-10-120-130 Hz, which about doubles the frequency content of the original data, resulting in sharper events.Looking at the sections alone, the validity, if any, of the bandwidth extension is not clear because the broadband earth model is not known.Nevertheless, the original narrowband data themselves actually provide a basic criterion to evaluate the effects of the different methods on the seismic data because fidelity should not be greatly degraded in the original data band in a valid spectral extension operation.If the original frequencies are preserved, high-frequency data can be low-cut filtered and spectrally shaped to restore the original spectrum; we can thus loosely test whether the methods are invertible to the original data by simple band-pass filtering back to the original frequency band.The significant residuals with the original data by phase acceleration (a 53.8% rms error in Figure 18a) and loop reconvolution (a 30.8% rms error in Figure 18b) imply that the original frequency components are greatly altered.The frequency invention methods badly fail the test.However, the harmonic extrapolation result (Figure 18c) shows only a minor error (10.0%rms).
Another test is to evaluate the synthetic tie for only frequencies (70-80-120-130 Hz) outside the band of the original data.For the purpose of this high-frequency observation, we choose an analysis window of the original synthetic tie (Figure 19a) at the well location with an initial excellent correlation coefficient of 0.96.The correlation has been reduced to 0.72 for the harmonic extrapolation result (Figure 19d), and it is entirely inadequate for the other methods (0.20 for phase acceleration [Figure 19b] and 0.45 for loop reconvolution [Figure 19c]).Possible causes for the reduced correlation coefficient of harmonic extrapolation relative to the original data include: (1) small time-depth errors, (2) magnification of noise, and (3) deviations from the presumed blocky earth model evident in the original well logs (Figure 15).
We can investigate the statistical significance of the correlations between those high-frequency extrapolated data and the well-log synthetic using an f-test (Snedecor and Cochran, 1989).This F-statistic for validating the prediction model is designed to test a null hypothesis stating there is no relationship between two measured phenomena (Fisher, 1925).The Fisher f coefficient is given by f ¼ ðr 2 ∕KÞ∕½ð1 − r 2 Þ∕ðn − K − 1Þ; (14) where r is the correlation coefficient, K is the number of model parameters, and n is the num-ber of data points.We take the number of free parameters to be the number of samples in the seismic wavelet.From this perspective, each well-log reflection coefficient is assumed to be independent from the others, and based on the convolutional model, the wavelet amplitudes represent the least mean-square-error coefficients needed to multiply the time-shifted coefficients by to produce a weighted sum predicting the seismic amplitude at each time.To the extent that the reflection coefficients are indeed independent, an f value significantly greater than unity would suggest that the model variance dominates over the residual variance, and it thus implies that the extrapolated high frequencies giving the correlation can be regarded as consistent with the underlying model.yields an f ¼ 3.19, suggesting that the 0.72 correlation is unlikely to be spurious.On the other hand, for frequency invention methods, even if their very poor correlations are thought to be acceptable, the f-test (f ¼ 0.12 for phase acceleration and f ¼ 0.75 for loop reconvolution) indicates poor statistical significance.
Table 1 includes all the parameters used in the comparative studies for the various methods.Note that for the loop reconvolution method, all the data are oversampled to avoid aliasing and are readily resampled back to any desired sampling rate with an acceptable Nyquist frequency after the reconvolution.

CONCLUSIONS
The subject of bandwidth extension is controversial.We know from digital signal processing that once a given frequency within a signal has been zeroed out, it is not recoverable by inverse linear filtering.Earth attenuation and noise limit the frequency band over which we can obtain positive S/N.Thus, if linear deconvolution operators are pushed too hard to boost high frequencies, the result is primarily noise magnification.This fact can easily lead one to the conclusion that bandwidth extension is theoretically not possible a false and overstated conclusion in the opinion of these authors.
We also know that there are a variety of ways that a seismic trace can be manipulated to produce high-frequency traces.Two of these methods are phase acceleration and loop reconvolution.These methods are not based on physics and essentially invent frequency content outside the seismic band.Wedge models, however, demonstrate that neither of these methods, although sharpening existing events, improve the resolution.The frequency invention methods also exhibit other deleterious consequences such as creating false events (in the case of loop reconvolution) and not filtering back to the input data (in both cases); indicating that the convolutional relationship among the data, the reflectivity, and the wavelet, has been destroyed within the original frequency band.Synthetic tests show poor correlation of the invented frequencies to high-frequency synthetics.Such methods may produce useful attributes for visualization, but we urge extreme caution if quantitative analysis, such as impedance inversion, is to be used on such bandwidth-extended outputs.
Harmonic extrapolation can be a valid bandwidth extension technique, based on physics, that is possibly applicable when the earth conforms well to the a priori assumption of a blocky impedance structure.A sparsely layered earth structure results in correlations in the reflectivity spectrum between frequencies across the digital spectrum out to the Nyquist frequency.In the simplest case of a single layer with two reflection coefficients, the reflectivity spectrum is sinusoidal in frequency, and generally a superposition of a sine and cosine related to the odd and even parts of the reflection coefficient pair.In the more general case of many superposed layer responses, the periodicity may or may not be adequately characterized within the seismic band.If the layers are thick and few enough, and the data are sufficiently noise free, it may be possible to decompose the spectrum to a summation of a limited number of sinusoids if the seismic wavelet spectrum is known.In this situation, frequencies outside the band are readily extrapolated.However, even in a noise-free case, if the reflectivity series is too complicated, as in the case of a random series of closely spaced reflection coefficients, bandwidth extension may not be achievable without additional constraints.Of note, even in this worst-case scenario, harmonic extrapolation can be applied to blue the spectrum within the band of the data, without boosting the high-frequency noise to the extent that would result from a simple spectral shaping or deconvolution.
A real data example shows that harmonic extrapolation produces potentially useful bandwidth-extended data that filter back to the original data with good fidelity.For the real data case studied, outside the seismic band, harmonic extrapolation produces a statistically significant correlation coefficient with a synthetic computed over the same bandwidth, whereas frequency-invention methods exhibit very poor and unacceptable correlation over the same band.

Figure 3 .
Figure3.Reflectivity amplitude spectrum for a simple layer of 50 ms thickness, with equal and opposite reflection coefficients of unit amplitude at the top and base.Two frequency cycles are apparent within a hypothetical seismic bandwidth (dashed lines).Assuming the wavelet is known, a band-limited reflectivity spectrum is obtained by dividing the data spectrum by the wavelet spectrum (solid circles), and a best-fit sinusoid (solid line) is readily fit to these observed amplitudes.It is then a simple matter to extrapolate the layer reflectivity spectrum outside the original seismic band by extending the sinusoid.

Figure 5 .
Figure 5. Amplitude spectra comparison for (a) a wedge model and (b) the frequency extrapolated result by harmonic extrapolation.The original model consists of layers with thicknesses varying from 50 to 1 ms by the 1 ms sampling rate.The reflection coefficients for each trace are equal and opposite at the top and base.(b) The inverted and extrapolated spectra are predicted from a narrow band of the original model by a 30 Hz Ricker wavelet, which gives a tuning thickness approximately 13 ms.

Figure 7 .
Figure 7. Residuals between the high-frequency (60 Hz) synthetic target section (a) and the extrapolated results by various methods: (b) phase acceleration with a 87.5% rms error, (c) loop reconvolution with a 80.8% rms error, and (d) harmonic extrapolation with a 9.5% rms error.

Figure 12 .
Figure 12.Spectra comparison between the original blocky reflectivity model and the full-band spectral extension by harmonic extrapolation from the 30 Hz low-frequency noisy synthetic data.(a) Amplitude spectrum of the original reflectivity, (b) original sparse reflectivity, (c) full-band frequency extrapolated spectrum, (d) inverted reflectivity series, (e) full-band prediction residual spectrum (41.7% rms error), and (f) prediction residual in the time domain.

Figure 14 .
Figure 14.Corresponding amplitude spectra for the test results on the random number series by the harmonic extrapolation method.(a) Spectrum of a random number series, (b) spectrum of the high-frequency (60 Hz) target synthetic seismogram, (c) spectrum of the high-frequency (60 Hz) recovery, (d) spectrum of residual for the 60 Hz result, (e) spectrum of the low-frequency (30 Hz) synthetic seismogram; (f) spectrum of the high-frequency (90 Hz) target synthetic seismogram, (g) spectrum of the high-frequency (90 Hz) recovery; and (h) spectrum of the residual for the 90 Hz result.
Figure15.Original tie between well-log synthetic and seismic data.The panel includes gamma ray, density, and sonic logs, along with the synthetic trace (blue), the composite trace at the well location (red) and the neighboring traces (black).The correlation coefficient is 0.87 for the entire log ranging from 1730 to 2430 ms.

Figure 16 .
Figure 16.Extracted seismic wavelet at the well location assuming the wavelet phase to be constant.The entire well log is used to determine the solution.(a) The wavelet in the time domain and (b) the amplitude and phase spectrum of the wavelet.The narrow usable bandwidth is taken between 10 and 55 Hz.

Figure 18 .
Figure 18.Spectrally extended seismic sections by various methods are filtered back to the original bandwidth and compared with the original data.(a) Residual for phase acceleration result with a 53.8% rms error, (b) residual for loop reconvolution result with a 30.8% rms error, and (c) residual for harmonic extrapolation result with a 10.0% rms error.The apparent error could be further reduced by spectral shaping.

Figure 19 .
Figure 19.Comparison for well-log synthetic ties at (a) the original frequency, and (b-d) with only extended high frequencies by various methods and the original seismic band excluded.The reference synthetic (blue) is generated by convolving the well-log reflectivity with the (a) seismic wavelet or (b-d) an analytical band-pass filter (70-80-120-130 Hz).The composite traces (red) are extracted from the traces (black) near the well location.The correlation window is specified within 2050-2250 ms.(a) Original well tie.The correlation coefficient is 0.96; (b) phase acceleration result.The correlation coefficient is 0.20; (c) loop reconvolution result.The correlation coefficient is 0.45; (d) harmonic extrapolation result.The correlation coefficient is 0.72.

Table 1 .
Comparison of parameters used for all comparative studies.All the methods are implemented on the entire trace.The original traces are oversampled in the loop reconvolution and are then resampled back to the original sample rate.