Abstract:
A method of designing and applying filters for geophysical time series data comprising obtaining a plurality of spatially related geophysical time series and transforming the time series using a wavelet transformation. The wavelet transformation coefficients may be organized into a plurality of subband traces. The method includes modifying one or more transform coefficients within a plurality of the subband traces or within all but one of the subbands of traces and then inverting the subband traces using an inverse transform to produce a filtered version of the transformed portion of the geophysical time series. The method allows for design and analysis of non-stationary filters and filter parameters in untransformed data space. Non-stationary signals may be filtered in all or in windowed portions of the geophysical time series data.

Description:
FIELD OF THE INVENTION  
         [0001]    This invention relates to the field of geophysical prospecting and, more particularly to methods for acquiring and processing seismic data.  
         BACKGROUND OF THE INVENTION  
         [0002]    In the oil and gas industry, geophysical seismic data analysis and processing techniques are commonly used to aid in the search for and evaluation of subterranean hydrocarbon deposits. Generally, a seismic energy source is used to generate a seismic signal that propagates into the earth and is at least partially reflected by subsurface seismic reflectors (i.e., interfaces between underground formations having different acoustic impedances). The reflections are recorded in a geophysical time series by seismic detectors located at or near the surface of the earth, in a body of water, or at known depths in boreholes, and the resulting seismic data may be processed to yield information relating to the location of the subsurface reflectors and the physical properties of the subsurface formations.  
           [0003]    Seismic or acoustic energy may be from an explosive, implosive, swept-frequency (chirp) or random source. A geophysical time series recording of the acoustic reflection and refraction wavefronts that travel from the source to a receiver is used to produce a seismic field record. Variations in the travel times of reflection and refraction events in these field records indicate the position of reflection surfaces within the earth. The analysis and correlation of events in one or more field records in seismic data processing produces an acoustic image that demonstrates subsurface structure. The acoustic images may be used to aid the search for and exploitation of valuable mineral deposits.  
           [0004]    These seismic data are processed to be useful for delineating features in the earth&#39;s subsurface. Geophysical time series data acquired in the field often contain unwanted energy as well as desired energy. Some of the unwanted energy may interfere with, overwhelm or otherwise mask desired seismic energy useful for subsurface investigations. There are many known methods for filtering or separating unwanted energy from seismic signals.  
           [0005]    Prior art processing methods for eliminating unwanted energy include frequency filtering, windowed frequency filtering, and f-k filtering. Filtering may be applied to each trace individually or to multiple traces concurrently. Elimination of unwanted energy traditionally results from the application of filtering to entire traces regardless of whether the unwanted energy is localized in only part of the data trace.  
           [0006]    Traditional filtering is usually implemented in frequency-time space using orthogonal basis functions having perfect localization in frequency but infinite extent in time. Use of the Fourier transform is well known in the art and the most common type of filtering for geophysical data. The underlying assumption of this type of filtering is that the seismic signals are stationary (e.g., do not have time dependent statistical characteristics). However, it is well known in the art that many types of seismic signals recorded in the field are non-stationary. It would be beneficial to filter geophysical time series data containing unwanted non-stationary energy without adversely impacting desired seismic signal.  
           [0007]    An alternative to the Fourier transform is the continuous wavelet transform (CWT). In the continuous wavelet transform, a function ψ used to create a family of wavelets ψ (at+b) where a and b are real numbers, a dilating the function ψ and b translating the function. The CWT turns a signal ƒ(t) into a function with two variables (scale and time), which may be called c(a,b) such that c(a, b)=∫ƒ (t)ψ(at+b)dt. The CWT allows for joint time-frequency localization ideal for non-stationary filtering. The discrete wavelet transform (DWT) is a special case of the CWT.  
           [0008]    Wavelet methods may be designed in transform space rather than data space. However, as is evident from comparing FIG. 2A, data prior to transformation, with FIG. 2B, data from the same seismic record after transformation with wavelet methods, the wavelet transform alters coordinate relationships in the original geophysical time series data ensemble. This altered coordinate relationship is due to an inescapable limitation on joint time-frequency resolution imposed by the transform technique. Though the transformed data contains the required non-stationary time-frequency signal characteristics, signal components are difficult to recognize from their original space-time geometry. To address this dilemma, some practitioners have developed graphing techniques to approximately restore original data-space coordinate relationships to the transformed signals. The details of the technique depend on the particular choice of transform. Deighan and Watts (1997) give an example for a wavelet transform in FIG. 6 of their paper. (GEOPHYSICS, VOL. 62, NO. 6, NOVEMBER-DECEMBER 1997; P. 1896-1903)  
           [0009]    U.S. Pat. No. 5,781,502 to Becquey proposes a method to discriminate elliptical waves propagating in a geologic formation by a combined process that includes measuring and comparing ratios of components along several wave axes. The method calls for the use of multi-component seismic data. The method is directed to removing elliptical waves from seismic data after wavelet analysis is used to determine the surface or tube waves in the data. The detected signals are subtracted from the field records.  
           [0010]    U.S. Pat. No. 5,850,622 to Vassiliou et al. disclose a method to filter seismic data using very short-time Fourier transforms. The method utilizes a very great number of short overlapping Fourier transform windows, together with a Gaussian weight or taper function, to produce a plurality of near single-frequency “sub-band” traces for each seismic trace so analyzed. The transformed data may be used for selective removal of coherent noise events, analysis for seismic attributes related to subsurface features of interest, seismic trace creation by interpolation, and automatic identification and removal of noisy seismic traces.  
         SUMMARY OF THE INVENTION  
         [0011]    A method of designing and applying filters for geophysical time series data is described which comprises obtaining a plurality of spatially related geophysical time series and transforming the time series using a wavelet transformation. The wavelet transformation coefficients may be organized into a plurality of subband traces. The method includes modifying one or more transform coefficients within a plurality of the subband traces or within all but one of the subbands of traces and then inverting the subband traces using an inverse transform to produce a filtered version of the transformed portion of the geophysical time series. The method allows for design and analysis of non-stationary filters and filter parameters in untransformed data space. Non-stationary signals may be filtered in all or in windowed portions of the geophysical time series data. 
       
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS  
       [0012]    The present invention and its advantages will be better understood by referring to the following detailed description and the attached drawings in which:  
         [0013]    [0013]FIGS. 1A and 1B illustrate a  3 D common shot ensemble of geophysical time series seismic signals before and after non-stationary filtering by wavelet transform.  
         [0014]    [0014]FIG. 2A illustrates a shot record containing non-stationary seismic energy.  
         [0015]    [0015]FIG. 2B illustrates the data of FIG. 2A after a wavelet transform has been applied.  
         [0016]    [0016]FIGS. 3 a  through  3   j  illustrates seismic data with varying frequency filters applied.  
         [0017]    [0017]FIG. 4 illustrates filter design methodology in non-transformed data space.  
         [0018]    [0018]FIG. 5 illustrates an embodiment of the wavelet filtering process.  
         [0019]    [0019]FIG. 6 illustrate the application of filter parameterization in transform-data space. 
     
    
       [0020]    While the invention will be described in connection with its preferred embodiments, it will be understood that the invention is not limited thereto. On the contrary, it is intended to cover all alternatives, modifications, and equivalents that may be included within the spirit and scope of the invention, as defined by the appended claims.  
       DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS  
       [0021]    The present invention is a method for designing and applying non-stationary filters using wavelet transforms that overcome the above-described deficiencies of current seismic processing practices. The filtering technique of the present invention provides for very efficient filter design and parameterization as well as full or partial suppression of non-stationary seismic energy in field records. The filters may be applied in a time-localized manner to further reduce impacting desired signal at other record times. Other advantages of the invention will be readily apparent to persons skilled in the art based on the following detailed description. To the extent that the following detailed description is specific to a particular embodiment or a particular use of the invention, this is intended to be illustrative and is not to be construed as limiting the scope of the invention.  
         [0022]    This invention provides effective non-stationary filters for attenuating non-stationary and stationary noise components in geophysical time series data signals. Noise components often interfere with desired signal within time and frequency data windows that vary according to source and receiver (spatial) location. This interference makes the composite signal non-stationary with respect to time and frequency. Groups of such signals with variable spatial location thus form an ensemble that is also non-stationary with respect to space. Ground roll, airwaves, swell noise, and cable noises are some examples of the noise components typically present in seismic data. The parameterization of the present invention prescribes filters adapted to the non-stationary characteristics of such noise.  
         [0023]    Many practitioners have applied non-stationary filters over the time coordinate of seismic signals. However, these filters lack an adequate process for determining time support of the filter within an ensemble of data. For example, the onset and termination times of the strong and unwanted airwave energy vary significantly over a seismic data record. Considering that typical seismic surveys contain hundreds of such ensembles, each one with different unwanted energy time windows, this is a serious shortcoming of previous filter parameterizations. A method that would parameterize analytic trajectories prescribing nominal onset and termination times to eliminate unwanted energy in a time or space localized fashion is highly desirable. Under the new parameterization presented here, analytic trajectories for nominal onset and ending times of noise components (and thus filters) for each subband and spatial coordinate, may be designed and applied for every ensemble of geophysical time series data.  
         [0024]    Many mathematical transforms exist for implementing time-frequency data decompositions necessary for non-stationary filters. Common transformations include wavelet transforms, wave packet transforms, sliding-window Fourier transforms, and matching pursuit decompositions. Each of these has different strengths and weaknesses. They provide different combinations of time-frequency resolution, accuracy and stability, and computational cost.  
         [0025]    The wavelet transform may be applied as a discrete wavelet transform (DWT). The DWT is a convenient way of using the wavelet transform with digital geophysical time series data. In a discrete wavelet transform, a wavelet is translated and dilated only by discrete values. For example, the dilation may be a power of 2 and be of the form ψ(2 j I+k), with j and k whole numbers. Orthogonal wavelets are a special case of discrete wavelets that give a representation without redundancy and lend themselves to fast algorithms.  
         [0026]    Some prior art filtering algorithms, Fourier transform methods for example, mix transform-data coefficients along the space coordinate of geophysical time series ensembles to detect and attenuate noisy coefficients. These approaches assume the desired signal has greater spatial coherence than the noise component. Other approaches zero or mute, to a predetermined level, all coefficients within the time subband support of the filter. This practice can easily harm or delete desirable signal within that same time subband window.  
         [0027]    [0027]FIG. 1A shows an example of non-stationary seismic energy, for instance airwaves  101  or ground roll, in a 3-D common-shot ensemble of geophysical time series data. The time-space distribution of the non-stationary energy component varies considerably, but predictably, throughout the ensemble. This noise concentrates, moreover, within particular subbands of the frequency spectrum. The present invention limits filter action within specified subband(s) to a computed time-space neighborhood of the undesirable signal present in geophysical time series data ensembles. Two analytic trajectories per subband define this neighborhood throughout the ensemble, wherever the noise component resides. The resulting filter passes data outside the prescribed time-space-subband window, unchanged and unharmed.  
         [0028]    Desirable signal may also be present within that time-space-subband window. This parameterization employs a two-stage prescription for the filter action within the window of data to be filtered. This method of design biases filtering towards those data most likely to be noise-saturated.  
         [0029]    This new filter design methodology presented here applies to invertible non-stationary transforms. Given an ensemble of signals, it conveys non-stationary characteristics in data space with the efficacy and cost afforded by that particular transform. One may therefore choose and evaluate filter design parameters directly in data space, incorporating the strengths and limitations imposed by one&#39;s choice of transform.  
         [0030]    The new filter application parameterization of the present invention establishes a data-adaptive threshold with the assumption that only noisy data can exceed it. A complementary parameter defines a second data-dependent level. This two-stage filtering mechanism reduces to the second, lower level only those coefficients that exceed the first threshold, established for noise. This new filter parameterization applies to invertible non-stationary transforms. It prescribes a non-stationary filter consistent with the efficacy and cost afforded by that particular transform.  
         [0031]    [0031]FIG. 1A compared to FIG. 1B gives a particular example of the application of a wavelet transform filter designed with the present invention. FIG. 1A contains non-stationary energy  101 , for example airwave or ground roll energy. After filtering particular subbands of the data using the wavelet transform, the strong non-stationary energy  101  of FIG. 1A has been removed  201  in FIG. 1B demonstrating non-stationary filtering as prescribed by this new parameterization and implemented via wavelet transform.  
         [0032]    This invention processes ensembles of seismic signals to produce output that is optimal for non-stationary filter design. The processed output allows analysis of non-stationary time-frequency characteristics in the non-transformed data space, thus permitting filter specification and testing directly in non-transformed data space. Non-transformed data space naturally preserves space-time coordinate relationships required for effective filter design. Moreover, practitioners ultimately judge effectiveness of the filter from its impact in non-transformed data space.  
         [0033]    Undesirable noise signal components often combine and interfere with the desired signal at variable times and with various frequencies. This interference results in non-stationary composite seismic signals. Non-stationary filters are therefore required to attenuate these undesirable components. Ground roll, airwaves, swell noise, and cable noises are examples of recorded seismic record components usually considered undesirable within the context of data processing objectives. FIG. 2A illustrates a shot record containing non-stationary seismic energy, ground roll and airwave components  101 .  
         [0034]    This invention specifies the support of a non-stationary filter by a nominal beginning time t b  and ending time t e  as, for example, illustrated in FIG. 4. Parameter functions t b   (j)  and t e   (j)  give a pair of times at every spatial coordinate x in the data ensemble, and apply for every subband j of the filter. These functions encompass families of two-parameter trajectories, including the linear and hyperbolic examples depicted in FIG. 4. In FIG. 4 the beginning filter time t b  is an example of a hyperbolic trajectory 401 and the ending filter time t e  trajectory is an example of a linear function  403 . In these two examples, scalars τ b   (j)  and τ e   (j)  take on the interpretation of beginning and ending times, respectively, at spatial coordinate x=0. For these linear and hyperbolic forms, scalars p b   (j)  and p e   (j)  assume the role of reciprocal velocity.  
         [0035]    The horizontal axis of FIG. 4 represents the space coordinate 415. The vertical axis of FIG. 4 represents the time coordinate  413 . Any arbitrary space coordinate position x  417  will have an associated time t b    405  and a t e    407  representing the beginning and end of the filter support region  409  for that space coordinate. The data areas  411  outside of the filter support region  409  are not filtered.  
         [0036]    [0036]FIG. 6 is a transform-data space representation of the filter application parameters. This figure is a generalization of how the parameters are used to apply the filtering. The vertical axis in FIG. 6 is the transform time with an example 16 transform data coefficients, ( 1 ) through ( 16 ). The FIG. 6 horizontal axis represents the transformed data coefficient values with positive values to the right of transform coefficient {circumflex over (d)} j =0 and negative values to the left. In FIG. 6, the example filter parameterization acts on only three of these 16 transform coefficients ( 1 ) through ( 16 ) along the FIG. 6 {circumflex over (t)} transform time vertical axis. This filter parameterization treats the remaining thirteen as desirable signal, preserving them without change. ρ j  prescribes the filter attenuation level as a fraction of the {overscore (s)} j  data statistic. These two parameters act in concert with the filter time-support parameters to detect and attenuate noise components in a data-dependent fashion.  
         [0037]    [0037]FIG. 6 graphs the roles of transform-data scalar parameters σ j  and ρ j  in this invention. σ j  specifies a noise detection threshold as a percentage of {overscore (s)} j , a data statistic computed from the transform coefficients in subband j of the ensemble. Reasonable choices for {overscore (s)} j  include the mean, the median and variance of the data transform coefficients, and many others are possible. When transform data coefficients {circumflex over (d)} j  in the filter region  619  exceed the threshold determined by σ j |{overscore (s)} j |  601  or +σ j |{overscore (s)} j |  607 , the muting function may reduce the data coefficient {circumflex over (d)} j  to the values−ρ j |{overscore (s)} j |  603  or +ρ j |{overscore (s)} j |  605 . For example, data coefficient {circumflex over (d)} j  ( 2 ) has a magnitude  621  greater than +σ j |{overscore (s)} j |. However, data coefficient {circumflex over (d)} j  ( 2 ) is not within window filter region  619 , so no filtering of the data occur. Data coefficient {circumflex over (d)} j  ( 7 ) also has a magnitude  613  which is greater than +σ j |{overscore (s)} j |  611 . Data coefficient {circumflex over (d)} J  ( 7 ) is between the transform time coordinates t b  at the beginning of the window filter region  619  window and the end of the window t e . The filtering operation will attenuate the magnitude of the data coefficient {circumflex over (d)} j  ( 7 ) to+ρ j |{overscore (s)} j |  609  along the line represented by +ρ j |{overscore (s)} j |  605 . Similarly, data coefficient {circumflex over (d)} j  ( 8 ) with transform coefficient  617  with a greater magnitude than−σ j |{overscore (s)} j |  601  would be attenuated to −ρ j |{overscore (s)} j |  603 . Data coefficient {circumflex over (d)} j  ( 10 ) exceeds the threshold magnitude so the data coefficient will be attenuated to +ρ j |{overscore (s)} j |  605 . The other data coefficients within the filter region  619  do not exceed the threshold magnitudes and would not be changed.  
         [0038]    Importantly, one can readily determine this suite of parameters for the seismic noise components enumerated above. One gleans this information from analysis of a few representative data ensembles, selected from the recorded survey. To attenuate the airwaves or other non-stationary seismic energy in FIG. 1A or FIG. 2A we may choose, for example, linear forms for the filter trajectories. That is, t b   (j) (x)=τ b   (j) +p b   (j) x and t e   (j) (x)=p e   (j) +ρ e   (j) x, where x specifies the scalar distance from the source position to every receiver position in this ensemble. The remaining time-support parameters are determined as: t b   (j) −0, τ e   (j)   −0.25, p   b   (J) ={fraction (1/300)}, p e   (j) ={fraction (1/300)}. The airwave speed is about 300 m/s. The time support of the airwave is about 0.25 s. Further analysis shows this airwave to be concentrated in subbands j=5 and j=6 of the wavelet transform. For this example the data statistic, s , may be chosen to be the mean of the transform data. The data coefficient threshold and attenuation parameters may, again for this example, be chosen respectively as σ 5 =σ 6 =10, ρ 5 =ρ 6 ={fraction (1/10)}. This parameter suite produces the filtered ensemble in FIG. 1B where the non-stationary energy present in FIG. 1A has been greatly attenuated or eliminated.  
         [0039]    The present invention enables non-stationary filter design in data space, the natural setting to observe space-time relationships. The present method restores space-time coordinate relationships to the data ensemble via the inverse transform, rather than through graphing in transform space. The problems inherent in designing filters in transform space may be seen by comparing FIG. 2B with FIG. 2A. FIG. 2B contains the data in FIG. 2A after applying a wavelet transform. There is a limitation on joint time-frequency resolution that makes the signal components difficult to recognize from their original space-time geometry. Importantly, the present invention preserves time-frequency signal characteristics required for non-stationary filter specification, which facilitates convenient filter specification and design. The method utilizes a sequence of muting operators in transform space, each one followed by the inverse transform. FIGS. 3 a - 3   j  are used to illustrate an example.  
         [0040]    [0040]FIG. 3 a  through FIG. 3 j  is a common shot ensemble of filter panels from the seismic data of FIG. 2A. FIG. 3 a  is the same record as FIG. 2A. FIG. 3 b  is the wavelet-transformed data of FIG. 3 a  (i.e. the same record as FIG. 2B). FIGS. 3 c  through  3   j  are a progression of individual frequency bands inverse transformed after muting has been applied to wavelet transformed data of FIG. 2A. Each figure panel is produced by selectively muting the wavelet-transformed data followed by inverting the transformed data.  
         [0041]    In the FIG. 3 a - 3   j  examples, the restored, linear geometry of the undesirable signal components primarily identifies them with frequency bands represented by FIG. 3 e  and FIG. 3 f . The practitioner would therefore specify a non-stationary filter that attenuates subbands in FIG. 3 e  and FIG. 3 f  of the wavelet transform. Moreover, filter attenuation may be limited to a neighborhood of linear space-time coordinates occupied by the ground roll and airwave, such that no other data are filtered.  
         [0042]    To use this invention for design of non-stationary filters for ensembles of seismic signals, operator T −1 M j T describes the entire filtering procedure in matrix notation. One applies non-stationary transform T to each input data vector in the ensemble, followed by filter operator M, and inverse transform T −1 . Operator M j  modifies or otherwise filters all but the jth subband of the transform. FIGS. 3 c - 3   j  demonstrate the operator for various subbands of a wavelet transform. For this transform the (i,l)th matrix element of M j  is·m il   (j) =δ ip δ pi  (where p=2 j +k and k indexes the time samples within subband j.  
         [0043]    To use this invention for applying the designed parameterization, operator T −1 {tilde over (M)} j T likewise describes the entire filtering procedure in matrix notation. One applies non-stationary transform T to each input data vector in the ensemble, followed by filter operator {tilde over (M)} j  and inverse transform T −1 . Operator {tilde over (M)} j  modifies or otherwise filters one or more (jth) subband of the transform. The details of this matrix operator depend on the details or characteristics of the signal and noise in T. For a wavelet transform, the (i,l)th element of {tilde over (M)} j  is: {tilde over (m)} il   (j) =δ il +δ ip ({tilde over (m)} pq −δ pq )δ ql . Scalars {tilde over (m)} pq  hold the data-dependent attenuation weights, where p=q=2 j +k and k indexes the time samples within subband j.  
         [0044]    A description of the elements of the non-stationary parameterization follows:  
         [0045]    T −1 {tilde over (M)} j T: Filter operator, where {tilde over (M)} j  is the filter, T is a non-stationary transform and j indexes the subband.  
         [0046]    {overscore (s)} j  : A statistic of the transform data (like mean, median, etc.; units are units of the transform data).  
         [0047]    σ j  : Unitless scalar parameter ≧0.  
         [0048]    σ j |{overscore (s)} j |: Threshold for triggering filter.  
         [0049]    ρ j  : Unitless scalar parameter ≧0.  
         [0050]    ρ j |{overscore (s)} j  : Replacement magnitude of filtered data.  
         [0051]    μ j  : Unitless muting coefficient that affects the magnitudes of data coefficient replacement values.  
           μ   j     =         1                 for                          d   ^     j            &lt;       σ   j                 s   _     j                        and                   μ   j         =         ρ   i                   s   _     j         d   ^     j                          for                          d   ^     j            ≥       σ   j                   s   _     j          .                                              
 
         [0052]    {circumflex over (d)} j  : Transform data in subband j.  
           α   j     =       ρ   j                   s   _     j         d   ^     j                                                 
 
         [0053]    We may express μ j  as: μ j =H(σ j |{overscore (s)} j |−|{circumflex over (d)} j |)+α j H(|{circumflex over (d)} j |−σ j |{overscore (s)} j |) where H is the Heaviside function defined as H(t−t′)=0 for t&lt;t′ and H(t−t′)=1 for t&gt;t′.  
         [0054]    Rewriting μ j : 
         μ j =1 −H (| {circumflex over (d)}   j |−σ j   |{overscore (s)}   j |)+α j   H (| {circumflex over (d)}   j |−σ j   |{overscore (s)}   j |)=1−(1−α j ) H (| {circumflex over (d)}   j |−σ j   |{overscore (s)}   j |). 
         [0055]    Note that μ j =μ j ({circumflex over (t)} k ), that is, μ j  depends on the transform time coordinate {circumflex over (t)} k .  
         [0056]    μ j ({circumflex over (t)} k )=1−(1α j ({circumflex over (t)} k ))H(|{circumflex over (d)} j ({circumflex over (t)} k ) |−σ j |{overscore (s)} j |) shows that muting coefficient μ j  depends on {circumflex over (t)} k  thru {circumflex over (d)} j  (and α j ). Note: k indexes time coordinate {circumflex over (t)} within subband j.  
         [0057]    Travel time trajectories are defined to parameterize regions of time-space coordinates over which muting coefficient μ j ({circumflex over (t)} k ) will be applied.  
         [0058]    t b  : Nominal beginning time of filter.  
         [0059]    t e  : Nominal ending time of filter.  
         [0060]    t b,j  (X;τ b,j , P b,j ): Scale dependent travel time trajectory for beginning times as function of space-coordinate x.  
         [0061]    x: Space coordinate.  
         [0062]    τ b,j  : Beginning time at x=0.  
         [0063]    P b,j  : Parameterizes trajectory shape.  
         [0064]    Line trajectory: t b,j =τ b,j +P b,j X.  
         [0065]    Parabola trajectory: t b,j =τ b,j +P b,j X 2 .  
         [0066]    Hyperbola trajectory: t b,j   2 =τ b,j   2 +P b,j   2 X 2 .  
         [0067]    t e,j (X;τ e,j , P e,j ): Scale dependent travel time trajectory for ending times.  
         [0068]    With nominal beginning and ending times given, respectively, as t b,j (x;τ b,j , P b,j ) and t e,j (x;τ e,j ), we may write that 
           {tilde over (m)}   j ( {circumflex over (t)}   k )=1−(1μ j ( {circumflex over (t)}   k ) H ( {circumflex over (t)}   k   −t   b ) H ( {circumflex over (t)}   e   −{circumflex over (t)}   k )=1−(1−μ j ( {circumflex over (t)}   k )(1 −H ( {circumflex over (t)}   k   −t   e ) 
         [0069]    where {tilde over (m)} j ({circumflex over (t)} k ) is an element of the operator {tilde over (M)} j  that applies transform-time dependent muting to the transform data {circumflex over (d)} j ({circumflex over (t)} k )  
         [0070]    [0070]FIG. 5 is an example of a muting operator for a Discrete Wavelet Transform (DWT). Take N−2 j −16, the number of samples in {circumflex over (d)}, the discrete wavelet transform data vector  503  where J=4. For j=2, the subband to be filtered, FIG. 5 is {tilde over (M)} 2 , an example of {tilde over (M)} j . To generalize the {tilde over (M)} j  DWT example for N=2 j  data samples, vector {tilde over (m)} 2 ({circumflex over (t)} k ) may be mapped into  501  matrix {tilde over (m)} il   (j) , the (i,l) th  element of M j (i,l=1,2 . . . , 2 j )  
         [0071]    Define {tilde over (m)} j ({circumflex over (t)} k )  Δ  {tilde over (m)} pq , with p=q=2 j +k, where j=1,2, . . . , J−1 lists eligible subbands (j=0 is not a wavelet subband but is necessary for reconstruction), and k=1,2, . . . , 2 j  lists eligible time indices within subband j. Therefore, p and q range over all integers from 2 j +1 so (k=1) to 2 J-1  where (k=2 j )  
         [0072]    These definitions give {tilde over (m)} il   (l) =δ il +δ ip ({tilde over (m)} pq −δ pq )δ ql  for i,l=1,2, . . . , 2 j ; 
           p=q= 2 j   +k; p,q= 2 j +1,2 j +2, . . . , 2 j+l   ; j= 1,2, . . . ,  J− 1; and  k =1,2, . . . , 2 j . 
         [0073]    Persons skilled in the art will understand that the methodology and processes for filter design and application described herein may be practiced in various forms including but not limited to the particular embodiments described herein. Further, it should be understood that the invention is not to be unduly limited to the foregoing which has been set forth for illustrative purposes. Various modifications and alternatives will be apparent to those skilled in the art without departing from the true scope of the invention, as defined in the following claims.