# Patent

## Title
Method to extract target signals of a known type from raw data containing an unknown number of target signals, interference, and noise

## Abstract

A signal analysis method extracts transient target signals of known type from a raw data source signal that contains an unknown number of target signals. The method can enhance the analysis of data obtained from in-line oil-debris sensors. The method comprises steps of: defining signatures of the known target signal, and of at least one of the intrinsic noise and interfering signals; performing a mathematical transform that decomposes the raw data into distinct data sets; using the signal signatures to identify and nullify the data sets containing noise and interfering signal signatures; using the target signal signatures to identify the data sets containing target signal components, or may further use a thresholding rule to remove intrinsic noise from said data sets, and finally applying the inverse transform to the processed data sets in order to reconstruct an enhanced output signal.

## Background
FIELD OF THE TECHNIQUE
The technique generally relates to data processing and more particularly, to the improvement of in-line oil debris sensor capability.
BACKGROUND OF THE TECHNIQUE
Early detection of mechanical damage by continuously monitoring the metallic debris in lubrication oil of various mechanical systems such as aircraft engines, is extremely important for fault diagnosis, maintenance decision-making, and accident prevention. In-line oil debris detection sensors have attracted attention recently because of the advantage of monitoring machine conditions continuously and eliminating some of the laboratory work required in off-line oil sample analysis. However, most oil debris sensors are not only sensitive to metal particles but also susceptible to various noises and vibrations. Therefore, there are two major issues to be addressed in real time oil debris detection. The first issue is how to use an existing sensor to detect the weakest possible debris or particle signals that are concealed in noise, such that quick responses to incipient failures can be achieved. The second is how to eliminate spurious signals caused by vibrations in the sensor working environment such that reliable maintenance decisions can be made.
The size of a particle that can be detected is mainly limited by background noise, while the challenge of eliminating spurious signals lies in their similarity to particle signals. This calls for an effective de-noising approach to address the two issues simultaneously.
Many de-noising techniques have been reported in the literature. Though de-noising has been studied by many researchers in various fields, work on purification of oil debris signal from the collected noisy data with the effects of vibration has not yet been reported in the accessible literature.
A collected oil debris signal is a mixture of several components including intrinsic noise, interfering vibration signals, and possible target oil-debris signals. However, most of the existing de-noising methods focus on the reduction of background Gaussian noise only, and thus cannot be directly applied to purify the collected oil debris signals.
Therefore there is a need for an improved data processing technique, particularly for enhancing in-line oil debris sensor capability.
SUMMARY OF THE DISCLOSURE
In one aspect of the technique, there is provided a method for extracting intermittent, randomly occurring transient target signals of known type from a raw data source signal containing an unknown number of the target signals in addition to interfering signals and/or intrinsic noise signals which comprises (a) defining quantifiable signatures or characteristics that respectively and independently represent the target signals, the interfering signals or the intrinsic noise signals; (b) transforming the raw data source signal with a mathematical transform that effects a decomposition of the source signal into distinct data sets; (c) processing said data sets to identify a first group of data sets which display the signatures or characteristics representing the interfering signals or intrinsic noise signals and a second group of data sets which display the signatures or characteristics representing the target signals, (d) setting the data sets of the first group to zero; and (e) applying an inverse transform to the processed data sets in order to reconstruct a processed output signal.
In another aspect of the technique, there is provided a method of enhancing the capability of a sensor for sensing a particle or bubble in a fluid flow, the sensor being adapted for generating transient target signals when particles or bubbles in the fluid pass through the sensor, which comprises the steps of (a) defining quantifiable signatures of the target signals, having a known frequency range of the target signals and a known kurtosis value range of wavelet coefficients of respective wavelet scales of interest that reflect a transient nature of the target signals; (b) defining quantifiable signatures of interfering signals and intrinsic noise signals, having two known kurtosis value ranges of wavelet coefficients of each of wavelet scales that correspond to the interfering signals and correspond to the intrinsic noise signals respectively; (c) obtaining a raw data sample signal in a time domain from the sensor in use for monitoring the fluid flow; (d) applying a Time-Invariant Wavelet Transform (TIWT) to said raw data sample signal to effect a decomposition of said raw data sample signal into a plurality of data sets in a form of wavelet coefficients of respective wavelet scales; (e) processing the data sets to calculate a kurtosis value for each data set to compare the kurtosis value for each data set with the kurtosis value range defined in step (a) or the kurtosis value ranges defined in step (b), in order to identify a first group of data sets which display the signatures representing the interfering signals and the intrinsic noise signals and a second group of data sets which display the signatures representing the target signals, and then setting the data sets of the first group to zero; (f) applying a thresholding rule to process the data sets of the second group in order to reduce any intrinsic noise signals which remain in the data sets of the second group; and (g) constructing a processed output signal by applying an inverse transform to the data sets processed in step (e) and (f).
Other features and advantages of the technique will be better understood with reference to the preferred embodiments described hereinafter.



BRIEF DESCRIPTION OF THE DRAWINGS
Reference will now be made to the accompanying drawings showing by way of illustration the preferred embodiments of the technique, in which:
FIG. 1 is a cross-sectional view of an oil debris sensor according to one embodiment employing the technology;
FIG. 2 illustrates shapes of ferrous and non-ferrous particle signatures generated by the sensor shown in FIG. 1;
FIG. 3 is a block diagram showing the length of each approximation and detail of signal x(n) after a Discrete Wavelet Transform (DWT) process;
FIG. 4 is a block diagram showing a comparison of the algorithms of DWT and Time-Invariant Wavelet Transform (TIWT);
FIG. 5 is a graphical illustration of signals for example with different sinusoidal compositions having different kurtosis values;
FIG. 6 is a block diagram showing one embodiment for de-noising a collected sample signal and recovering a desired signature;
FIG. 7 is a graphical illustration of extraction of a non-ferrous particle signature from a collected sample signal without vibration components;
FIG. 8 is a graphical illustration of extraction of a non-ferrous particle signature from a mixture of the collected sample signal and simulated vibration signal components;
FIG. 9 is a graphical illustration of extraction of a non-ferrous particle signature from a mixture of a collected sample signal and a periodically modulated vibration;
FIG. 10 is a graphical illustration of extraction of a ferrous particle signature from a sample signal collected from an in-line oil debris sensor mounted on a shaker; and
FIG. 11 is a block diagram showing another embodiment for de-noising a collected sample signal and recovering a desired signature, which is preferred if memory capacity is limited.



DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS
The definition of the following terms are applicable throughout the specification and the claims of this application.
A target signal is a signal which is sought. According to the described technique, target signals must be of finite duration, and must possess a known characteristic profile. Target signals may however occur at random intervals. In the case of an oil-debris sensor, a target signal is a transient signal created by the oil-debris sensor in response to the passage of a suspended particle passing through the sensor.
An intrinsic noise signal refers to a minimum, unavoidable noise of a sensor and measurement system. The intrinsic noise signal is the signal obtained from the sensor in controlled conditions. These conditions are contrived so as to possess no target signals and no significant interfering signals. In the case of an oil-debris sensor, the intrinsic noise signal is obtained from the sensor when operating in quiet laboratory conditions, without the passage of carrier medium (i.e. no metal particles) through the sensor.
An interference signal is any signal other than the target signal and intrinsic noise signal as defined above. In the case of an oil-debris sensor, interference signals may be caused by vibration of the structure where the sensor is mounted. Such interference is called vibration interference hereafter.
A raw data signal is the data collected in full operating conditions of a sensor system. Raw data signals may include intrinsic noise, possible interference signals, and possible target signals. In the case of an oil-debris sensor, raw data is collected over a time window when the carrier medium (e.g., fluid) that may possess suspended particles, passes through the sensor and when the sensor is subjected to interfering signals arising from vibration.
As illustrated in FIG. 1, an oil debris sensor 20, as an example of sensors used in detecting particles in a fluid flow, generally includes a tubular body 22 connected in a pipeline (not shown), such as a lubricant oil line, to allow the oil flow to pass through the tubular body 22 as shown by arrows F. A coil assembly including one sensing coil 24 and two field coils 26, is mounted around the tubular body 22. The sensing coil 24 is located equidistant between the upstream and downstream field coils 26. When ferromagnetic particles or non-ferromagnetic particles including bubbles in the oil flow, pass the sensor 20, the sensing coil 24 generates a target signal of known characteristics. In real operating conditions, this target signal is usually mixed with interfering and/or noise signals and can be definitively identified only after the raw data output signal from the sensing coil 24 is properly processed. The sensor 20 is known and similar to that described in U.S. Pat. No. 5,315,243 granted to Kempster et al. on May 4, 1994, which is incorporated herein by reference.
FIG. 2 illustrates ferrous and non-ferrous target particle signals 28 and 30 of characteristics of the intermittent, randomly occurring and non-periodical transient signal components generated by the sensor 20 in FIG. 1, representing a detected ferrous particle and non-ferrous particle passing through the sensor 20. The ferrous particle signature 28 and non-ferrous particle signals 30 are similar but are 180° out of phase relative to each other. The shapes of the target signals 28 and 30 may be similar to one full period of a sine wave, depending on conditions under which the sensor 20 operates.
In a practical working environment, the particle target signals 28, 30 in a sample signal obtained from the sensor 20, are mixed with other unwanted signal components caused by background noise, vibration interference, etc. The responsiveness level of the sensor 20 depends on the size of the particles to be detected in the oil flow. The detectable particle size is mainly limited by background noise. Spurious signatures such as vibration-like signals may cause false alarms. The challenge of improving the responsiveness and reducing false alarms lies in the very weak particle signals and their similarity to spurious signals.
Therefore, a method of enhancing the capabilities of a sensor for sensing particles or bubbles in a fluid flow such as an oil flow is developed, generally based on the definition of a target signatures or characteristics of the known target signals generated by sensor 20. The target characteristic or signature is defined to enable the data analysis, as described below.
To provide better de-noising results, target signal signatures or characteristics can be defined using general information relating to the target signal. For example, the target particle signal frequency range may be estimated and this frequency can be used as a target signal signature. At the early stage of a fault developed in a machine, the metal debris is very small and it can be assumed that the particle and oil speeds are identical. Hence the particle signal frequency can be estimated based on oil speed. The flow speed of lubricating oil is assumed to be in range [νmin,νmax]. νmin and νmax represent the minimum speed and the maximum speed, respectively. As shown in FIG. 1, the distance lc between the ends of the two field coils 26 in the sensor 20, can be measured and therefore, the time required for a particle to pass through the two field coils 26 can be calculated. A corresponding frequency component can be obtained as a reciprocal of said time period. Therefore, a target particle signal characteristic is a frequency in range [νmin/lc,νmax/lc]. Obviously, any signal with a frequency outside the ranges should be treated as interfering signals and should be filtered out. However, in practice, not all interfering signals and noises can be removed in this way. Filtered signals may still contain noise and vibration-like signals in the particle frequency range. Therefore, the invention provides for a more complete digital filtering process using the above target signal signature.
In a machine monitoring process, for example, the sensor 20 of FIG. 1 is connected to an oil circuit of the machine for detecting particles in the oil flow passing through the oil circuit. A raw data sample signal in a time domain, as the raw data source signal, can thereby be obtained from the sensor 20. The raw data sample signal obtained from the sensor 20 is then processed in a digital filtering process using defined signatures or characteristics of the target signal and of the noise and/or interfering signals. These are used to separate any target signals from the other raw data signal components.
According to the invention, this process is conducted by processing the data, using a mathematical transform that deconstructs the raw data source signal into distinct, transformed data sets. Each data set is analyzed for signatures of the target signal and of the background noise and interferences. Depending on the specific noise, data sets displaying unwanted signal components are either nullified (set to zero) as in the case of data sets containing interference type signatures or, in the case of data sets containing signatures of the target signal, the intrinsic noise can be removed from the data sets which can then be said to be ‘de-noised’ or purified. The processed data sets are then reconstructed to create a processed purified signal.
In a preferred embodiment, several signal signatures are used, including the frequency range of the target particle signal as described above, and a statistical-type signature that distinguishes target particle signal from vibration interference as described later.
The digitally processed sample signal is re-constructed in the time domain. The re-constructed, processed sample signal contains the target signal component in greater proportion with respect to the unwanted noise and vibration-like components, compared to that found in the original, pre-processed raw data sample signal. The capability and responsiveness of sensors such as sensor 20 is therefore enhanced.
According to one embodiment, the mathematical transform utilized is the joint Time-Invariant Wavelet Transform (referred to as a TIWT hereinafter). In this case, each data set that results from the transform is called wavelet coefficients of a scale.
The specific mother wavelets used in the TIWT are pre-selected so as to possess similar morphology (signature) to the particle target signal.
According to one embodiment, signal signatures are determined by a statistical approach such as a kurtosis analysis. The statistical analysis is used to differentiate the target signals from the interference and noise components. Kurtosis analysis is used to: a) remove background noise by setting the wavelet coefficients that are in the scales with kurtosis value between [3−Δ2, 3+Δ2] to zero; b) identify and remove vibration-like signal components by setting the wavelet coefficients that are in the scales with kurtosis value between [1.5−Δ1, 1.5+Δ1] to zero; and c) identify and extract the non-periodic, intermittent, randomly occurring and transient signal component of known form by retaining wavelet coefficients, and hence the associated signal components, in the scales having sufficiently high kurtosis value.
Wavelet transform is used to: a) decompose the raw signal into frequency bands so that some signal components are de-coupled and each can be more easily handled; b) remove any signal component with frequency that is out of range [νmin/lc,νmax/lc] by setting the associated wavelet coefficients to zero (where vmin and vmax are the minimum and maximum particle passage speeds respectively); and c) remove the Gaussian noise from the high-kurtosis scales (where the kurtosis is high and so these scales contain the particle signal, thus the kurtosis cannot be used for removing noise from these scales and instead a thresholding rule is used).
In summary, both wavelet transform and kurtosis analysis contribute to both noise reduction and spurious signal removal. Such tasks cannot be properly done by separately applying wavelet transform and kurtosis analysis. To be effective, the kurtosis analysis has to be embedded and applied to each wavelet transform scale of interest.
In contrast to many existing de-noising algorithms where discrete wavelet transform (referred to as DWT hereinafter) is employed, TIWT is adopted in this technique because compared with DWT, TIWT has the advantage of translation invariance. In most de-noising methods based on DWT, the size of a signal shrinks in the process of wavelet decomposition, with different shrinkage rates in different frequency bands (or scales). This leads to two problems: reduced time resolution and deteriorated statistical accuracy when a statistical method is jointly applied. With the TIWT, a signal is never down-sampled and the filters are different on each decomposition level. The signal can be decomposed into different scales (frequency bands) without any size shrinkage in any scales (frequency bands). Therefore, adopting TIWT can effectively avoid the two problems. As a result, the same good time resolution can be preserved in all scales (frequency bands) and a statistical method can be applied.
Based on the experimental observation of oil debris and spurious, particularly vibration-like signals, it may be suggested that kurtosis, a statistical indicator, be used to distinguish the two. Machine vibration can be characterized by a series of consecutive waves with specific frequency components while the signals generated by particles appear intermittently but not periodically, it is therefore possible to distinguish vibration signal and particle signal through the analysis of signal statistical characteristics. Kurtosis as a statistical indicator is sensitive to the shape of a signal. Therefore, it will be adopted for in-line oil debris detection. A joint TIWT and kurtosis analysis method may be used to de-noise the collected sample signal data and eliminate spurious signals.
Wavelet transform, if properly applied, is able to reveal useful information of a signal in a time-frequency domain. A wavelet ψ(t) is a function of zero average with limited duration on (−∞,∞). The continuous wavelet transform (CWT) of a finite energy signal x(t) can be defined by




  
    
      
        
          W
          
            a
            ,
            b
          
        
        =
        
          
            
              
              a
              
            
            
              
                -
                1
              
              /
              2
            
          
          ⁢
          
            ∫
            
              
                x
                ⁡
                
                  (
                  t
                  )
                
              
              ⁢
              
                
                  ψ
                  _
                
                (
                
                  
                    t
                    -
                    b
                  
                  a
                
                )
              
              ⁢
              
                d
                t
              
            
          
        
      
    
    
      
        (
        1
        )
      
    
  




where a and b are the dilation and time parameters of the wavelet respectively {a>0, bεR}, Wa,b denotes wavelet transform coefficients, ψa,b represents a wavelet formed from a mother wavelet function ψ(t) through rescaling and shifting, ψ stands for complex conjugate. A different scale a corresponds to a specific frequency component range.

Scale-time parameters (a,b) can be discretized by a=2j and b=k2j for fast computing in practical implementation. Accordingly, the discrete wavelet transform (DWT) is given by

Wj,k=2−j/2∫x(t) ψ(2−jt−k)dt  (2)

where j and k are arbitrary integers.

A fast discrete wavelet transform called multi-resolution pyramidal algorithm is known in the art. With this algorithm, the original signal x(t) passes through two complementary filters and emerges as a signal in a low frequency band and a signal in a high frequency band, respectively. Hence, the original signal can be decomposed into different frequency bands through an iteration process. The pyramidal algorithm is implemented using




  
    
      
        
          <
          
            x
            ⁡
            
              (
              t
              )
            
          
        
        ,
        
          
            φ
            
              j
              ,
              k
            
          
          >=
          
            
              ∑
              k
            
            ⁢
            
              
            
            ⁢
            
              <
              
                x
                ⁡
                
                  (
                  t
                  )
                
              
            
          
        
        ,
        
          
            φ
            
              
                j
                -
                1
              
              ,
              k
            
          
          >
          
            g
            ⁡
            
              (
              
                k
                -
                
                  2
                  ⁢
                  n
                
              
              )
            
          
        
      
    
    
      
        (
        3
        )
      
    
  
  
    
      
        
          <
          
            x
            ⁡
            
              (
              t
              )
            
          
        
        ,
        
          
            ψ
            
              j
              ,
              k
            
          
          >=
          
            
              ∑
              k
            
            ⁢
            
              
            
            ⁢
            
              <
              
                x
                ⁡
                
                  (
                  t
                  )
                
              
            
          
        
        ,
        
          
            φ
            
              
                j
                -
                1
              
              ,
              k
            
          
          >
          
            h
            ⁡
            
              (
              
                k
                -
                
                  2
                  ⁢
                  n
                
              
              )
            
          
        
      
    
    
      
        (
        4
        )
      
    
  




where φj,k is a scaling function, g(k) and h(k) are respectively the impulse responses of low and high pass filters, and <·,·> denotes the inner product operation. For a signal x(n), c0,n=x(n) and n is an arbitrary integer. According to equations (3) and (4), the DWT of signal x(n) can be written as





  
    
      
        
          c
          
            j
            ,
            n
          
        
        =
        
          
            ∑
            k
          
          ⁢
          
            
          
          ⁢
          
            
              c
              
                
                  j
                  -
                  1
                
                ,
                k
              
            
            ⁢
            
              g
              ⁡
              
                (
                
                  k
                  -
                  
                    2
                    ⁢
                    n
                  
                
                )
              
            
          
        
      
    
    
      
        (
        5
        )
      
    
  
  
    
      
        
          w
          
            j
            ,
            n
          
        
        =
        
          
            ∑
            k
          
          ⁢
          
            
          
          ⁢
          
            
              c
              
                
                  j
                  -
                  1
                
                ,
                k
              
            
            ⁢
            
              h
              ⁡
              
                (
                
                  k
                  -
                  
                    2
                    ⁢
                    n
                  
                
                )
              
            
          
        
      
    
    
      
        (
        6
        )
      
    
  




where cj,n and wj,n are discrete scaling coefficients and wavelet coefficients, respectively. They represent the approximations and details of x(n) at resolution 2j, respectively.

In equations (5) and (6), a down-sampling algorithm which keeps one point out of two, is used to perform the transform. The length of signal x(n) will reduce by half after each transform. If original signal x(n) has N samples, discrete coefficients cj,n and wj,n have N/2j samples, respectively, j=1, 2, . . . , J. J is the maximum number of discrete wavelet transform iterations as illustrated in FIG. 3. Therefore, the length of each approximation and detail of x(n) becomes shorter after each transform iteration. With the decrease in sample length, the statistical characteristics of the approximation and detail located in specific frequency bands may no longer be accurate.
To alleviate the above difficulty, TIWT is proposed by introducing an up-sampling algorithm into the DWT in an attempt to keep the time-invariance property of signals, which is preferred in statistical signal processing implementation. For function ƒ(n), the up-sampling algorithm can be defined as ƒ(2n)=ƒ(n) and ƒ(2n+1)=0. The basic idea of TIWT is to modify the quadrature mirror filters through an up-sampling algorithm by




  
    
      
        {
        
          
            
              
                
                  g
                  ⁡
                  
                    (
                    
                      2
                      ⁢
                      n
                    
                    )
                  
                
                =
                
                  g
                  ⁡
                  
                    (
                    n
                    )
                  
                
              
            
          
          
            
              
                
                  g
                  ⁡
                  
                    (
                    
                      
                        2
                        ⁢
                        n
                      
                      -
                      1
                    
                    )
                  
                
                =
                0
              
            
          
        
      
    
    
      
        (
        9
        )
      
    
  
  
    
      
        {
        
          
            
              
                
                  h
                  ⁡
                  
                    (
                    
                      2
                      ⁢
                      n
                    
                    )
                  
                
                =
                
                  h
                  ⁡
                  
                    (
                    n
                    )
                  
                
              
            
          
          
            
              
                
                  h
                  ⁡
                  
                    (
                    
                      
                        2
                        ⁢
                        n
                      
                      -
                      1
                    
                    )
                  
                
                =
                0
              
            
          
        
      
    
    
      
        (
        10
        )
      
    
  




before performing a convolution between the filters and the original signal or previous discrete scaling coefficients obtained by wavelet transform, in order to keep the redundancy of wavelet transforms. The key is to modify filters by inserting zeros in them. Consequently, equations (5) and (6) can be revised as





  
    
      
        
          c
          
            j
            ,
            n
          
        
        =
        
          
            ∑
            l
          
          ⁢
          
            
          
          ⁢
          
            
              c
              
                
                  j
                  -
                  1
                
                ,
                
                  n
                  +
                  
                    2
                    
                      j
                      -
                      
                        1
                        l
                      
                    
                  
                
              
            
            ⁢
            
              g
              ⁡
              
                (
                l
                )
              
            
          
        
      
    
    
      
        (
        11
        )
      
    
  
  
    
      
        
          w
          
            j
            ,
            n
          
        
        =
        
          
            ∑
            l
          
          ⁢
          
            
          
          ⁢
          
            
              c
              
                
                  j
                  -
                  1
                
                ,
                
                  n
                  +
                  
                    2
                    
                      j
                      -
                      
                        1
                        l
                      
                    
                  
                
              
            
            ⁢
            
              h
              ⁡
              
                (
                l
                )
              
            
          
        
      
    
    
      
        (
        12
        )
      
    
  




respectively, where l is an arbitrary integer. Therefore, the nature of time-invariant wavelet transform of the original signal is maintained. All the TIWT coefficients at the Jth decomposition level can be written as {Wj,n|j=1, 2, . . . , J+1} consisting of cj,n, wj,n, wJ−1,n, . . . , and wl,n.

The algorithms of DWT as indicated by numeral 32 and TIWT as indicated by numeral 34, are illustrated in FIG. 4, wherein gj and hj represent the filters used at resolution 2j. In the algorithm of TIWT 34, when high-pass and low-pass filters are applied to a signal, the two new sequences will have the same size as the original one. This effectively preserves the time-invariance property of wavelet coefficients, thereby maintaining the time resolution of the original signal and the statistical properties of the wavelet coefficients. Based on the above analysis, the main differences between algorithms DWT 32 and TIWT 34 can be summarized in Table I.







TABLE I







THE DIFFERENCES BETWEEN THE DWT AND THE TIWT










DWT
TIWT















Size of filters
No changes before
Doubled before



g and h
each DWT
each DWT



Size of wavelet
Half of the prior
Same as the original



transform of
signal size
signal size



signals



Translation
No
Yes



invariance



Time resolution
Low
High



Statistical
Low accuracy
High accuracy



properties



Required size
The size should be
The size should be



of original
large enough to
large enough to



signal
capture signatures
capture signatures and





the size should be





divisible by 2











Frequency
Same




resolution










As mentioned above, the oil debris sensors are not only sensitive to metal particles but also susceptible to noise signals, interferences such as vibration signals. Some of these effects may generate signals similar to the particle signals. The similarity between the two types of signals makes their separation a very challenging issue and improper removal of noise or vibration signals can spoil the outcome in later steps. Therefore, the removal of vibration signal should be handled very carefully. In practice, a mechanical vibration signal is periodic, similar to a sinusoidal wave or a combination of several sinusoidal waves with relatively stable frequencies. On the other hand, particle target signals are of finite duration and they occur randomly and intermittently i.e. neither the target signals nor their rate of occurrence is periodic. Therefore, the collective pattern of a vibration signal is different from that of a particle, though their individual waveforms may look very much alike. It should also be noted that the signal of an individual particle is similar to one full period sine wave in the time domain. As kurtosis is well suited to the detection of such transient signals, it should also be effective in detecting particle signals.
Given a time series x=(x1, x2, . . . , xN), the kurtosis value can be easily computed on line by




  
    
      
        kurtosis
        =
        
          
            
              ∑
              i
            
            ⁢
            
              
            
            ⁢
            
              
                
                  (
                  
                    
                      x
                      i
                    
                    -
                    
                      x
                      _
                    
                  
                  )
                
                4
              
              /
              N
            
          
          
            
              (
              
                
                  ∑
                  i
                
                ⁢
                
                  
                
                ⁢
                
                  
                    
                      (
                      
                        
                          x
                          i
                        
                        -
                        
                          x
                          _
                        
                      
                      )
                    
                    2
                  
                  /
                  N
                
              
              )
            
            2
          
        
      
    
    
      
        (
        13
        )
      
    
  




where i=1, 2, . . . , N, and x is the mean value of x, N is the size of x.

More specifically, the kurtosis of a signal consisting of successive but non-periodic individual wave segments is different from that of a periodic signal of the same wave form. This can be demonstrated by simulated signals. Consider the following four signals as indicated by numerals 36, 38, 40 and 42 in FIG. 5:



  
    
      
        x
        ⁢
        
          
        
        ⁢
        1
        ⁢
        
          (
          t
          )
        
      
      =
      
        sin
        ⁡
        
          (
          
            2
            ⁢
            π
            ⁢
            
              
            
            ⁢
            100
            ⁢
            
              
            
            ⁢
            t
          
          )
        
      
    
    ,
    
      

    
    ⁢
    
      
        x
        ⁢
        
          
        
        ⁢
        2
        ⁢
        
          (
          t
          )
        
      
      =
      
        sin
        ⁡
        
          (
          
            2
            ⁢
            π
            ⁢
            
              
            
            ⁢
            15
            ⁢
            
              
            
            ⁢
            t
          
          )
        
      
    
    ,
    
      

    
    ⁢
    
      
        x
        ⁡
        
          (
          t
          )
        
      
      =
      
        
          x
          ⁢
          
            
          
          ⁢
          1
          ⁢
          
            (
            t
            )
          
        
        +
        
          x
          ⁢
          
            
          
          ⁢
          2
          ⁢
          
            (
            t
            )
          
        
      
    
    ,
    
      

    
    ⁢
    and
  




  
    
      x
      ⁢
      
        
      
      ⁢
      3
      ⁢
      
        (
        t
        )
      
    
    =
    
      {
      
        
          
            
              sin
              ⁡
              
                (
                
                  2
                  ⁢
                  π
                  ⁢
                  
                    
                  
                  ⁢
                  5
                  ⁢
                  
                    
                  
                  ⁢
                  t
                
                )
              
            
          
          
            
              0.432
              <
              t
              <
              0.5
            
          
        
        
          
            0
          
          
            else
          
        
      
    
  



The signals between 0 and 1 second time duration and their kurtosis are shown in FIG. 5. It can be seen that 1) the kurtosis is 1.5 for sinusoidal signals 36 and 38 regardless of their frequencies; 2) the kurtosis of signal 40 formed by a combination of multiple sinusoidal functions is 2.25, higher than that of a single sinusoid function such as signals 36, 38; and 3) a signal 42 containing only one segment of a sine wave leads to a much higher kurtosis value, 22.4999.

The above observation gives rise to the use of kurtosis value as a classification criterion for separating vibration and particle signals. However, the collected vibration interference signal may be a mixture of several periodic signals of different frequencies such as signal 40. The kurtosis value of such a data mixture is dependent on the mixture contents and hence cannot be directly used for signal classification. For this reason, the collected data mixture has to be decomposed into several frequency bands. A frequency band may be associated with a unique signal, with or without periodicity. Kurtosis for each frequency band should be calculated separately to make meaningful classification decisions.
As illustrated in FIG. 5, the kurtosis of a sinusoidal signal is 1.5 which can be used as a criterion to remove such signals from the collected data. Given signal x(t), all TIWT coefficients at the Jth decomposition level, {Wj,n|j=1, 2, . . . , J+1} can be obtained. For each j corresponding to a specific frequency band, the associated kurtosis value can be calculated. However, considering the fact that the actual vibration signals are unlikely to be truly sinusoidal, a tolerance range is used and the criterion is relaxed to 1.5±Δ1 (Δ1=1 in this embodiment).
In addition to its aforementioned ability to separate periodic and non-periodic individual wave segments signals, kurtosis is also a good indicator of background noise. The background noise is generally Gaussian noise with a zero mean, the kurtosis of which is three. However, the noise may not be truly Gaussian noise, therefore the criterion is relaxed to 3±Δ2 (Δ2 is set to 0.5 in this embodiment). This kurtosis range will be used as a criterion to remove background noise. Then wavelet coefficients in different scales are either set to zero or kept based on their relevance to the transient particle signal. Such decision may be based on the following rules:

    
    
        Rule 1: If the frequency range of a scale is not in the range of interest, i.e., [νmin/lc,νmax/lc], set the wavelet coefficients in the frequency range to zero; and
        Rule 2: If the kurtosis of the wavelet coefficients of a scale is between 1.5±Δ1 or 3±Δ2, set the wavelet coefficients of this scale to zero.
    
    


The remaining coefficients on a certain scale should contain the information of particles. It should be noted that the remaining coefficients also contain the information of noise as well. For this reason, threshold σ√{square root over (2 log(N))} is used to shrink the remaining wavelet transform coefficients wherein N is the number of samples and σ is the standard deviation of the noise. Though the background noise can be reasonably assumed to be Gaussian noise with a zero mean, its strength, reflected by its standard deviation σ, is unknown. In the prior art, σ is usually estimated by MAD/0.675. MAD is the Median Absolute Deviation value of the appropriately normalized fine-scale wavelet coefficients. To provide more accurate results, it may be suggested to calculate the true strength of noise, σj, corresponding to each scale j directly from the raw data collected without the presence of any particles or vibration signals. Since such raw data contains only background noise, the directly calculated σj represents the true noise level. In practice, the raw data collected for calculation of σj values may be obtained from a reference sample signal obtained from the same sensor 20 of FIG. 1 without fluid passing therethrough and in the absence of structural vibration. The intrinsic noise sample signal is decomposed into the same number of scales as the number of scales of the raw data source signal using TIWT.
The main steps of the proposed embodiment are illustrated in FIG. 6. First, as shown in the left two blocks, the background noise is collected and decomposed by TIWT. Then, the standard deviation σj of the coefficients on each scale (j=1, 2, . . . , J) is calculated. The σj values will be saved in the computer and used to calculate the thresholds for de-noising. Next, as represented on the right part of FIG. 6, the raw data that represents a mixture of particle signal, vibration signal, and background noise, is for example, the sample signal obtained from the sensor 20 of FIG. 1 during the monitoring operation of the sensor 20 connected in the oil circuit of a machine, and is decomposed by TIWT. The obtained coefficients will be processed based on rules 1 and 2 described above. Once the two rules are implemented, the undesirable frequency components are filtered out. The signal components in the remaining scales are then de-noised using the hard thresholding rule of σj√{square root over (2 log(N))}. Finally, the particle signal is reconstructed from the processed coefficients. The reconstructed signal can then be used for particle identification.
It should be noted that the kurtosis value ranges 1.5±1 or 3±0.5 used in the above embodiment are estimated based on experiments for description and demonstration of this embodiment. In practice, the kurtosis value ranges to be used in this process may be different.
It should also be noted that since a wavelet basis is not unique, selecting an appropriate wavelet for a specific problem is usually desirable. For example, Daubechies wavelets are orthogonal and compactly supported so as to provide an ideal mathematical tool for constructing conjugate quadrature filters in signal analysis, and are therefore selected for use in this embodiment in order to avoid information leakage in wavelet transform. It may also be preferable that the mother wavelets match the shape of the particle signatures, such as signatures 28 and 30 shown in FIG. 2. For this reason, db4 wavelet function is selected. In other applications, the wavelet basis to be used in this process may be different.
The number of decomposition levels which determines the number of scales, affects the analysis results. Increasing the number of decomposition levels will reduce the band width of the frequency component of each scale and hence the interfering signal components may be better distinguished from the transient particle signature. However, too many decomposition levels inevitably increase computing time which is undesirable for on line applications. To specify appropriate decomposition levels, it may be useful to know the transient particle signal frequency range. The frequency resolution of the highest scale may preferably be high enough to identify the lowest frequency of interest. In other words, the frequency band width of the highest scale may preferably be smaller than the lowest frequency of interest. Denoting ƒw as the frequency band width of the highest scale and L as the number of decomposition levels, ƒw should satisfy the following condition

ƒw=ƒs/2L+1<ƒl  (12)

where ƒl=νmin/lc is the lowest frequency of interest, and ƒs denotes the sampling frequency. Then,





  
    
      
        L
        >
        
          
            
              
                ln
                ⁢
                
                  
                
                ⁢
                
                  f
                  s
                
              
              -
              
                ln
                ⁢
                
                  
                
                ⁢
                
                  f
                  l
                
              
            
            
              ln
              ⁢
              
                
              
              ⁢
              2
            
          
          -
          1
        
      
    
    
      
        (
        13
        )
      
    
  




For example, if the sampling frequency is 8000 Hz, the decomposition level should satisfy





  
    
      
        L
        >
        
          
            
              
                ln
                ⁢
                
                  
                
                ⁢
                8000
              
              -
              
                ln
                ⁡
                
                  (
                  
                    
                      v
                      min
                    
                    /
                    
                      l
                      c
                    
                  
                  )
                
              
            
            
              ln
              ⁢
              
                
              
              ⁢
              2
            
          
          -
          1
        
      
    
    
      
        (
        14
        )
      
    
  




FIG. 7 illustrates a sample signal 44 including a particle target signal generated by a non-ferrous particle of size 211 μm mixed with a background noise signal and shows the de-noising result in accordance with the described method, in which the noise signal components are removed and the non-ferrous particle signal 46 is identified.
FIGS. 8 and 9 illustrate sample signals 48 and 52 including particle target signals 50 and 54 generated by a non-ferrous particle of size 211 μm. As seen in the figures the target signals are essentially hidden in the noise and simulated vibration signals. Therefore the signatures 50 and 54 in signals 48 and 52 are not visually observable. It is impossible to successfully process the signals obtained from the oil debris sensors in a conventional manner such as described in U.S. Pat. No. 5,315,243. However, the target signals 50 and 54 of the non-ferrous particles are recovered by removing the noise and vibration signal components in accordance with the above-described technique.
FIG. 10 illustrates a further example of the digital filtering result using the method developed in this technique. A target signal 58 of a ferrous particle of size 125 μm is submerged in the original signal sample 56 and as a result of the filtering according to the invention, is recovered.
FIG. 11 illustrates another de-noising algorithm in accordance with another embodiment. Modification has been made to the de-noising algorithm illustrated by the flow chart of FIG. 6, in order to reduce the microprocessor memory required to accommodate many data samples during the decomposition process for on-line applications. The embodiment shown in FIG. 6 provides better enhancement performance than the embodiment of FIG. 11, if memory capacity permits. However, if the memory capacity is limited, the embodiment of FIG. 11 is preferred.
In the version of TIWT analysis as shown in FIG. 11, following the calculation of kurtosis value, the wavelet coefficients resulting from every decomposition step are down-sampled. Down-sampling is performed according to the level number. For example if the kurtosis value is in the accepted range and it is decided to use the nth decomposition level for the reconstruction purpose, such coefficients will then be down-sampled with the factor of 2n. With this approach the same memory space as the number of acquired samples is needed. The reconstruction would be performed in a way similar to the inverse discrete wavelet transform (IDWT). In comparison to the algorithm shown in FIG. 6, the memory space saving is (L−1)S, where L=the number of decomposition levels and S=the number of samples. The aspects of FIG. 11 similar to those of FIG. 6 will not be redundantly described.
The enhanced ability of oil debris sensors in detecting very small particles provides better responsiveness to machinery fault signals and makes it possible to issue warnings at a much earlier stage. The added capacity in eliminating vibration or vibration-like spurious signals can improve the quality of maintenance decisions and thus reduce false alarms in fault detection.
In the described method, the accuracy of the analysis results depends, to a certain extent, on the calculation of kurtosis. If the size of wavelet coefficients decreases, the accuracy of kurtosis will decrease accordingly. This weakness can be compensated for by jointly using TIWT in the algorithm because the TIWT can effectively maintain constant wavelet coefficient size and hence the accuracy of kurtosis can be preserved.
Modifications and improvements brought about by above-described embodiments may become apparent to those skilled in the art. For example, the general principles of the above-described embodiments are not limited to the processing of transient signals generated by sensors for the detection of particles in fluid flow. The principles of the above-described embodiments may be applicable to data processing in other situations where finite-duration target signals of known characteristics must be identified from raw data input source signal in a time domain containing an unknown number of the target signals and unwanted interfering or noise signals. The invention in this broader form could, in principle, be applied to: acoustic data analysis for detection of finite duration target signals such as animal calls, mechanical acoustic signals relating to machine fatigue or failure, gunshots; optical signals for detection of emboli (gas bubbles) in blood transfusion lines; acoustic or electronic signals for the detection of switching noise in communication lines; electronic signals associated with radar or sonar; cleaning of ‘pop’ noise in musical recordings, and other applications. The foregoing description is intended to be exemplary rather than limiting.

## Claims

1. A method for extracting intermittent, randomly occurring transient target signals of known type from a raw data source signal containing an unknown number of the target signals in addition to interfering signals or intrinsic noise signals, the method comprising:
a) obtaining the raw data source signal from a sensor for detecting the passage of particles or bubbles in a fluid flow, wherein a target signal is generated by the sensor when detecting a particle or bubble;
b) defining quantifiable signatures or characteristics that respectively and independently represent the target signals, the interfering signals or the intrinsic noise signals, wherein the signatures comprise a frequency range of the target signals and a kurtosis value range of wavelet coefficients of respective wavelet scales of interest that reflect a transient nature of the target signals;
c) applying a Time-Invariant Wavelet Transform (TIWT) to decompose the raw data source signal into distinct data sets in a form of wavelet coefficients of the respective wavelet scales;
d) processing said data sets to identify a first group of data sets which display the signatures or characteristics representing the interfering signals or intrinsic noise signals and a second group of data sets which display the signatures or characteristics representing the target signals;
e) setting the data sets of the first group to zero; and
f) applying an inverse transform to the processed data sets in order to reconstruct a processed output signal.

2. The method as defined in claim 1 wherein step (d) further comprises removing any intrinsic noise signals which remain in the data sets of the second group.

3. The method as defined in claim 1 wherein the signatures comprise two kurtosis value ranges of the wavelet coefficients of each wavelet scale that correspond to the interfering signals and correspond to the intrinsic noise signals, respectively.

4. The method as defined in claim 2 comprising employing a thresholding rule to process the data sets of the second group in order to reduce the intrinsic noise signals.

5. The method as defined in claim 1 wherein an Inverse Time-Invariant Wavelet Transform (ITIWT) is applied to the processed data sets for reconstructing the processed output signal.

6. The method as defined in claim 1 wherein the raw data source signal comprises acoustic signals.

7. The method as defined in claim 1 wherein the raw data source signal comprises optical signals.

8. The method of claim 1 in which the raw data source signal is obtained from sonar or radar apparatus.

9. The method as defined in claim 1 wherein the raw data source signals comprise electronic signals.

10. The method as defined in claim 1 wherein the frequency range of the target signals is determined from information relating to a physical dimension of the sensor and a velocity of the fluid flow.

11. The method as defined in claim 1 wherein the raw data source signal obtained from the sensor is decomposed into a number of scales so that a frequency band width of a highest scale is smaller than a lowest frequency of the target signals.

12. The method as defined in claim 1 wherein the identifying process of the data sets comprises calculating a kurtosis value for each wavelet scale of the raw data source signal and then comparing the same with the kurtosis value range of the wavelet coefficients associated with the defined signatures of the target signals, thereby identifying the respective first and second group of data sets.

13. The method as defined in claim 12 wherein wavelet coefficients in respective wavelet scales are down-sampled prior to an inverse transform according to a decomposition level, following a kurtosis value calculation in the respective scales of said raw data source signal, in order to reduce memory requirement in hardware implementation.

14. The method as defined in claim 13 wherein the inverse transform is an Inverse Discrete Wavelet Transform (IDWT).

15. A method of enhancing the capability of a sensor for sensing a particle or bubble in a fluid flow, the sensor being adapted for generating transient target signals when particles or bubbles in the fluid pass through the sensor, the method comprising steps of:
a) defining quantifiable signatures of the target signals, having a known frequency range of the target signals and a known kurtosis value range of wavelet coefficients of respective wavelet scales of interest that reflect a transient nature of the target signals;
b) defining quantifiable signatures of interfering signals and intrinsic noise signals, having two known kurtosis value ranges of wavelet coefficients of each wavelet scale that correspond to the interfering signals and correspond to the intrinsic noise signals, respectively;
c) obtaining a raw data sample signal in a time domain from the sensor in use for monitoring the fluid flow;
d) applying a Time-Invariant Wavelet Transform (TIWT) to said raw data sample signal to effect a decomposition of said raw data sample signal into a plurality of data sets in a form of wavelet coefficients of respective wavelet scales;
e) processing the data sets to calculate a kurtosis value for each data set to compare the kurtosis value for each data set with the kurtosis value range defined in step (a) or the kurtosis value ranges defined in step (b), in order to identify a first group of data sets which display the signatures representing the interfering signals and the intrinsic noise signals and a second group of data sets which display the signatures representing the target signals, and then setting the data sets of the first group to zero;
f) applying a thresholding rule to process the data sets of the second group in order to reduce any intrinsic noise signals which remain in the data sets of the second group; and
g) constructing a processed output signal by applying an inverse transform to the data sets processed in steps (e) and (f).