|Publication number||US5799114 A|
|Application number||US 08/608,381|
|Publication date||Aug 25, 1998|
|Filing date||Feb 28, 1996|
|Priority date||May 5, 1993|
|Also published as||CA2198685A1, EP0793182A1|
|Publication number||08608381, 608381, US 5799114 A, US 5799114A, US-A-5799114, US5799114 A, US5799114A|
|Inventors||Martin Joseph Dowling|
|Original Assignee||Liberty Technologies, Inc.|
|Export Citation||BiBTeX, EndNote, RefMan|
|Patent Citations (4), Non-Patent Citations (10), Referenced by (6), Classifications (7), Legal Events (6)|
|External Links: USPTO, USPTO Assignment, Espacenet|
This application is a continuation-in-part of copending application Ser. No. 08/058,113 filed May 5, 1993, now U.S. Pat. No. 5,517,585. This invention was made with government support under Grant No. DMI-9303205 awarded by the National Science Foundation. The government has certain rights in this invention.
The present invention generally relates to a system and method for stably identifying and characterizing isolated data patterns, or transients, in a signal representing some physical entity, and more specifically to situations where the transients may be arbitrarily aligned with the sample points arising from an analog to digital (A/D) process applied to the signal.
In general, when a sampled digital signal from an A/D process is operated on by a signal processing function, such as correlation, convolution, or transform, an arbitrary shift in the input signal with respect to the sample points causes changes in the output of the signal processing function, as well. This effect is aggravated when the sampling rate is reduced close to the critical sampling rate. FIG. 1 illustrates an arbitrary shift in the sampling of a signal. For many digital signal processes, such shifts produce a significant change in the output of the signal processing function. For example, the discrete Fourier transform computes inner products of time data and a complex exponential. If a signal is sampled and then resampled at the same rate but in such a manner that the second set of samples does not line up with the first, then the real, imaginary, and phase spectra will be different. This variation in the end result arising from an arbitrary shift in the sample points on the signal is referred to herein as shift variance. In the case of the Discrete Fourier Transform, shift variance is not generally a problem because the commonly used magnitude spectrum is not shift variant.
For certain processes, shift variance can be aggravated where sampling is sparse on certain scales. An example is the affine discrete (especially critically sampled) wavelet transform. FIG. 2 shows a transient (FIG. 2a) and its wavelet decomposition (FIGS. 2b, 2c and 2d) before and after a shift (FIGS. 2e, 2f, 2g, 2h). The decomposition is totally different because of a shift of just one sample.
Thus, though the continuous form of the wavelet transform is shift-invariant, the discrete form is not. FIG. 2 illustrates that the discrete wavelet transform is sensitive to the exact positioning of the transient in a sampled time record. Note that dramatically different coefficients result, depending on where the transient happens to occur in the sampled time record. For physical signals, the location of the transient in the sampled time record is usually unknown and completely arbitrary. Thus, for physical signals, the discrete wavelet transform (used in digital signal processing) is unreliable in its ability to detect, characterize and classify transients. This is unfortunate because as a continuous transform (unrealizable in practice), the wavelet is unsurpassed in its ability to identify and characterize transients.
For many processes, such as the discrete wavelet transform, a shift in sampling with respect to the input signal results in a different, not simply shifted, output (e.g., different wavelet transform coefficients). This is especially true for non-integer sample (sub-sample) shifts. The shift-variance problem is aggravated in multiscale analysis, such as the discrete wavelet transform, because even an integer sample shift on the most resolute scale is a non-integer shift on less resolute scales.
Conventional attempts to overcome the shift-variance problem are generally inadequate and typically rely on either oversampling, or resampling the data pattern based on zero-crossings or maxima locations. There also exist specialized transforms designed to be immune to particular types of shifts, but these specialized transforms are flawed due to their lack of generality.
Oversampling reduces the sampling interval so that the extent of arbitrary shifting is confined. For example, if the original sampling rate were 1 KHz, the worst case sub-sample positional shift would be 0.5 msec in the time record. If the sampling rate were increased to 10 KHz, the worse case would be 0.05 msec. However, such improvement is achieved at the expense of the need to process a much greater amount of data. Note that the oversampling approach makes no attempt to solve the underlying shift-variance problem. Rather, the oversampling approach simply bounds the error. oversampling can be implemented in hardware by increasing the rate of the analog-to-digital converter, or in software by interpolation (such as upsampling with zeroes followed by low-pass filtering). In either case, the cost to achieve oversampling is significant. In summary, oversampling is flawed because it does not optimally solve the shift-variance problem, and greatly increases the burden of computation.
The flaws associated with oversampling are aggravated for methods, such as the discrete wavelet transform, which analyze signals on multiple scales. Wavelets are defined as the sets of functions: ##EQU1## The sets are formed by dilations, a, and translations, b, of a function g(x), called the analyzing wavelet. The set of functions maintains the shape of the analyzing wavelet throughout the time-frequency or time-scale plane.
The discrete wavelet transform of a signal s(x) can be represented by: ##EQU2## which is essentially an inner product of the discrete form of the conjugate of the wavelet, g(.), with the discrete form of the signal, s(.), where T is the sampling interval (T=1/fs). By varying the dilation, a, the transform provides decomposition at many scales. By analogy to map-making, "scale" is the level of detail provided at dilation a. Various conditions are usually imposed on wavelets, namely, the analyzing wavelet must have finite energy, and the integral of the wavelet must vanish.
Two currently employed general approaches have been proposed specifically for solving the shift variance problem of the discrete wavelet transform. The first approach estimates the positions of the zero-crossings in the record, then resamples the data in software in order to stabilize the samples. While this approach is an improvement over prior art, it experiences the following difficulties: (a) forming the resampled representation and reconstruction is mathematically intense; (b) it cannot be used with scale-envelopes because envelopes do not have zero-crossings; (c) it is designed for shift-correcting global patterns and thus is of questionable applicability for isolated transients and other isolated patterns; and (d) it does not always produce a unique representation and reconstruction.
The second approach is a variation of the first method. Instead of using zero-crossings, the second approach uses signal maxima as the basis for resampling. The second approach is an improvement upon the first approach, but is still mathematically intensive, of questionable applicability for isolated patterns, and does not guarantee uniqueness.
Recent approaches have attempted to make shift-invariant a more general kind of wavelet transform, known as wavelet packets, by measuring the cost function for each possible indexing combination (odd or even in the dyadic) at each decimation in order to chose the "best" decomposition for invariance. This technique, however, dramatically increases the computational requirements and is not applicable to arbitrary (i.e., sub-sample) shifts and is unsatisfactory for the general case.
Briefly stated, the present invention is directed to a method of stably identifying and characterizing a transient data pattern contained in an electrical signal representing a physical event, a sampled representation of the electrical signal comprising a first data set having a plurality of sample points, wherein alignment of the sample points with the transient data pattern is arbitrary. The method comprises the steps of detecting and isolating the transient data pattern in the first data set; determining a location of the detected and isolated transient data pattern in the first data set; measuring a difference between the location of the detected and isolated transient data pattern and a predetermined location of a predefined mathematical operator at a point of optimal alignment, the predefined mathematical operator including preselected parameters; adjusting the preselected parameters of the predefined mathematical operator in accordance with the measured location difference such that the adjusted mathematical operator is now aligned with the detected and isolated transient data pattern at the point of optimal alignment; analyzing the first data set using the aligned mathematical operator; and characterizing the detected and isolated transient data pattern to identify the physical event causing the transient data pattern.
In a second embodiment, the method comprises the steps of:
(a) analyzing the first data set using a predefined mathematical operator, the predefined mathematical operator including preselected parameters, to produce a second data set in which the transient data pattern can be detected and isolated;
(b) detecting and isolating the transient data pattern in the second data set;
(c) determining a location of the detected and isolated transient data pattern in the second data set;
(d) measuring a difference between the location of the detected and isolated transient data pattern and a predetermined location of the predefined mathematical operator at a point of optimal alignment;
(e) adjusting the location of the detected and isolated transient data pattern in accordance with the measured location difference such that the detected and isolated transient data pattern is now aligned with the predefined mathematical operator;
(f) analyzing the location adjusted transient data pattern using the predefined mathematical operator; and
(g) characterizing the analyzed location adjusted transient data pattern to identify the physical event causing the transient data pattern.
For the purpose of illustrating the invention, there is shown in the drawings embodiments which are presently preferred. It should be understood, however, that the invention is not limited to the precise methods, arrangements, and instrumentalities shown. In the drawings:
FIG. 1 shows two consecutive sample sets of a signal, the second sample set being shifted from the first sample set;
FIGS. 2a-2h show a transient and its wavelet decomposition before and after a shift. (FIG. 2a shows the original signal before the shift. Its 3-scale wavelet transform is shown below in FIGS. 2b-2d. The signal after a one sample shift is shown in FIG. 2e, and the new changed decomposition in FIGS. 2f-2h.);
FIG. 3 is a block diagram of a preferred computer-based system which operates in accordance with the methods of the present invention;
FIG. 4 is a flow chart for analyzing arbitrarily shifted data patterns in accordance with a first embodiment of the present invention;
FIG. 5 is a flow chart for analyzing arbitrarily shifted data patterns in accordance with a second embodiment of the present invention;
FIGS. 6a-6k depict signals for illustrating the operation of the present invention;
FIGS. 7 and 8 show how the proposed method would work in a Gaussian-shaped signal before (FIG. 7) and after (FIG. 8) a shift in sampling; and
FIGS. 9 and 10 show how the proposed method would work on a non-symmetrical triangular shaped signal before (FIG. 9) and after (FIG. 10) a shift in sampling.
The present invention is directed to methods of stably identifying and characterizing a transient contained in a sampled set of physical data where alignment of the sample points with the transient is typically arbitrary. In the following discussion, reference is made to the signal processing of one-dimensional signals. However, it should be understood that the present invention is applicable to other types of processing, such as image processing, and to signals of any number of dimensions.
Referring now to the drawings, in which like numerals indicate like elements throughout, there is shown in FIG. 3 a block diagram of a preferred computer-based system 302 which operates in accordance with the methods of the present invention. The system 302 includes a sensor 304, such as an accelerometer, pressure transducer, microphone, etc. The sensor 304 receives and converts a dynamic physical phenomena (not shown) into an analogous standard instrumentation signal (not shown) using a medium such as electricity, fluid, or light. The signal may represent, for example, an image or element thereof which is being analyzed. More generally, the signal may represent any localized or discrete data pattern.
The system 302 also includes signal conditioning equipment 306, which amplifies, impedance matches, filters, and otherwise improves the signal received from the sensor 304. The system 302 further includes an analog-to-digital (A/D) converter 308, which digitizes the analog signal received from the signal conditioning equipment 306. Signal conditioning equipment and A/D converters are well known to those of ordinary skill in the art and need not be discussed in detail for a complete understanding of the present invention.
A processor 310 receives and analyzes the digital signal from the A/D converter 308 in accordance with the methods of the present invention. Since most modern measurements are based on digitized sensor signals, the present process has wide applicability. Furthermore, it is usually the nature of incipient failures in machinery to appear sporadically at first, and evidence themselves as low-level transient signals in the sample domain. Thus, one important application of the present invention is for detection of incipient failure in machinery monitoring and diagnostics. The steps of the methods of the present invention are preferably embodied in computer software (not shown) which may reside in a random access memory (RAM) 316, or ROM (320) or be otherwise stored and accessible for execution by the processor 310, described below. As will be appreciated, the processor 310 operates in accordance with the computer software. Alternatively, as will be appreciated by those of ordinary skill in the art, the methods of the present invention may be implemented using other means, such as hard wired logic circuits used in combination with or instead of the processor 310 and software (not shown).
The system 302 also includes a keyboard, mouse, or other input device 322, a communication port 318 (for receiving from another location or sending to another location, signals for analysis by the processor 310), a display device 314 and a printer 312. Note that the RAM 316 may store signals for analysis by the processor 310. The structure and operation of the components in the system 302 (other than the computer software which implements the methods of the present invention) are generally well known. Such components of the system 302 may be implemented using any number of well-known, readily available computer-related products.
FIG. 4 is a flow chart in accordance with a first embodiment of the present invention, showing a method of stably identifying and characterizing a transient contained in a first data set, representing physical data, where alignment of the sample data points with the transient may be arbitrary.
In step 408, the processor 310 analyzes the first data set using a predefined mathematical operator to produce a second data set. More particularly, in step 408 the processor 310 analyzes the digital signal (also called herein the raw signal, or first data set) received from the A/D converter 308 (or retrieved from RAM 316, or received from another location via the communication port 318). Such analysis may include the use of an analysis function such as a convolution or correlation, for example. The convolution operation is: ##EQU3## where s(.) is the sequence of signal samples, h(.) is the filter impulse response, and c(.) is the output of the convolution. Alternatively, for cases in which the transient can be detected and isolated sufficiently well in the sample domain (time or space) to establish location, then it is not necessary to analyze the first data set with the predetermined mathematical operator before alignment (i.e. step 408 may be skipped). Thus, in the case where the transient is detected and isolated without applying the predefined mathematical operator to the first data set, the second data set is the same as the first data set.
In step 410, the processor 310 detects and isolates the discrete data pattern (associated with the transient) in the second data set (that is, the output of the analysis function). More particularly, in step 410 the processor 310 detects unmodulated transients in the analysis function output (also called herein the analyzed signal, or second data set). Modulated transients are also preferably detected in step 410 but they can alternatively be detected after amplitude demodulation (step 416). Any conventional detection technique can be used such as comparing the absolute or squared value of every sample against a threshold. Such detection can be enhanced by requiring that several adjacent samples must exceed the threshold before being counted. Another improvement includes ignoring occasional dropouts below the threshold. The detection process in step 410 may be further enhanced by passing demodulated transients through a non-linear filter that triggers only when a high-energy level exists within a certain time restriction (i.e., trigger when a<Δt<b, where a and b are predefined time durations based on some general knowledge of the characteristics of the transients). Note that the upper limit may be long enough to include the possibility that two close transients may superficially appear to be one long signal. In the case of multiscale analysis, such as resulting from application of the affine wavelet transform, detection of an event can be confirmed by looking for simultaneous threshold exceedances on neighboring scales.
After a transient is detected, the detected transient is isolated in step 410 by estimating the boundary of the detected transient. The transient boundary is estimated based on, for example, the region of significant energy above background level.
In steps 412 and 414, the processor 310 tests the analyzed signal to determine if the analyzed signal is modulated, i.e., has the form of a modulated carrier. An auto-regression moving average (ARMA) or fast fourier transform (FFT) approach can be used, for example, to detect modulation. The modulation test can alternatively be done on the raw (unanalyzed) signal, but it is preferably performed after analysis (i.e., after step 408) for better reliability. It should be understood that the sequence of detection (step 410) and testing for modulation (steps 412 and 414) is not critical and should be tailored to the particular application.
If in steps 412 and 414 the transient is found to be essentially a time-limited amplitude-modulated carrier, then better results will often be obtained by calculating the position of the transient (step 418) after first demodulating the analyzed signal. Hence, for modulated transients, the processor 310 in step 416 envelope detects the second data set. That is, the processor 310 calculates the envelope of the transient detected in the analyzed signal. Preferably, amplitude demodulation is performed by creating an analytic signal on each scale. The envelope is then calculated as the magnitude of the analytic signal on each scale. Amplitude demodulation is preferably accomplished by well-known time domain techniques such as using the Hilbert transform, or by using well-known frequency domain techniques which typically include the steps of calculating the FFT of the analyzed signal, zeroing negative components and doubling positive components, and then taking the inverse FFT. While the analytic signal approach to demodulation has been identified as the preferred technique, other amplitude demodulation techniques may also be used in step 416.
In step 418, the processor 310 determines a location of the detected and isolated data pattern in the second data set. That is, the processor 310 determines a location of each transient detected in step 410. The determination (or measurement) conducted in step 418 must be such that small changes in the shape of the analyzed signal should result in small changes in the coefficients, and small errors in the measurement (step 418) should result in small changes in the resulting coefficients. Many techniques satisfying the above criteria exist for measuring the locations of detected transients, including:
a. Start of transient as transient location;
b. Center of transient (as measured by a generic statistic such as the centroid) as transient location;
c. The centroid approach (b) enhanced by weighting regions of low energy with low weights;
d. Minimum distance (such as the well-known Cramer-von Mises technique) fitted to a density function such as normal or the well-known Weibull function (the location of the transient then being taken as that of the density function); and
e. Minimum distance approach (d) used with specific classification models. In this case, classification and location estimation would be accomplished at the same time.
Either the weighted centroid (c) or the minimum distance approach with specific class models (e) is preferred. Oftentimes, the locating techniques are implemented in the presence of noise. In order to overcome the effects of noise, the "tails" of the transient (i.e., the front and back ends of the transient) may be trimmed. The tails of the transient may be eliminated by eliminating the low-energy components of the transient far from the center of the transient. These far points are significantly affected by noise since they are low level and thus, unduly influence the center calculation. The locating techniques can also be enhanced by energy weighting. In energy weighting, the signal is heavily weighted in areas of high energy and lightly weighted in low-energy areas. One method of accomplishing energy weighting is to define the centroid as: ##EQU4##
In step 420, the processor 310 determines the location of the analyzing function. The location of the analyzing function is generally known in advance and is preferably the location parameter used to generate the analyzing function. For example, if a Gaussian-enveloped sinusoid is the analyzing function, then the center of the envelope is preferably used as the location.
In step 422, the processor 310 measures a difference between the location of the detected and isolated data pattern and the location of the analyzing function. That is, the processor 310 measures the distance between the location of each transient (determined in step 418) and the location of the analyzing function (determined in step 420). If the analyzing function has several allowed locations (the usual case), then the nearest location is preferred as the reference. An example of a distance measure is the Euclidean.
In step 424, the processor 310 adjusts preselected parameters of the mathematical operator in accordance with the measured location difference such that the adjusted mathematical operator is aligned with the detected and isolated transient data pattern. That is, in step 424 the processor 310 shifts the analyzing function by the distance that was determined in step 422 to bring the location of the analyzing function into alignment with the location of the transient. Since the analyzing function is a mathematical function, step 424 is easily done by adjusting the location parameter of the analyzing function in accordance with the location difference measured in step 422. For example, if the analyzing function is a Gaussian-modulated cosine, the Gaussian envelope as well as the cosine are both shifted by the same amount.
The realignment of the analyzing function is done by translating the time (or space, in the case of images) reference of the function. The implementation of this process depends on the nature of the analyzing function. If the function is correlational in nature (such as the DFT or cross-correlation), then the function is aligned with the location of the transient by a translation in the same direction. If the function is convolutional in nature, such as an FIR filter, the function is aligned by a translation in the apparent opposite direction. This is because convolution, in effect, reverses the order of the analyzing domain. Thus, "alignment" must be interpreted in accordance with the nature of the analyzing function.
Alternatively, in some cases, it is more accurate to translate the transient into alignment with the analyzing function then vice versa, in which case step 425 would be executed in lieu of step 424. Step 425 would be executed, for instance, in the case of multi-resolution analysis with quadrature mirror filters having an odd number of coefficients. In multi-resolution, only the first scale high and low pass filters need be shifted (subsequent scales use unshifted filters; they are automatically shift-compensated because they are derived from the shifted low-pass filter in the initial dyadic decomposition). However, because of the odd number of coefficients, the high-pass filter has a zero at the folding frequency which interferes with its filtering function. Translating the transient can be accomplished by such means as the following two-step process: (1) Band limit the transient, if required, by low-pass filtering, for example, the cutoff can be placed at 0.9 of folding frequency (for an FIR filter, an odd number of coefficients is used and adjusted for the (integer-sample) delay of the filter, and for an infinite impulse response (IIR) filter, the IIR filter is run in both directions to effect zero phase); and (2) Interpolate the transient by convolving with a windowed (sin(π(n-k)))/(π(n-k)) function, where "n" is the sample value and "k" is the desired sub-sample shift. Even in the example of multi-resolution using odd numbers of coefficients, all scales (except the finest) are stabilized by shifting the analyzing function. Note, the essence of the process remains the same whether the transient or analyzing function is shifted, or whether the analyzing function is correlational or convolutional in nature. In any case, the process causes the sampled transient and analyzing function to be brought into alignment to provide the important feature of shift invariance.
In step 426 the processor 310 reanalyzes the first data set using the aligned mathematical operator. That is, the processor 310 reanalyzes the raw signal using the shifted analysis function. Note that the entire raw signal is not reanalyzed, only those portions containing transients. For signals containing multiple transients, the analysis function is shifted and applied independently for each transient. Note that it is not always possible to separate transients significantly overlapping in both time and frequency (or scale) domains. As a result of step 426, a set of shift-invariant analysis coefficients for each transient is produced.
In the special case where two transients overlap in time but have different frequency (scale) content, the same section of the signal is reanalyzed twice. The analyzing function is shifted to bring it into alignment with the first transient, and then the analyzing function is shifted again but by a different amount to bring it into alignment with the second transient. The above assumes that the transients have been sufficiently isolated and that their centroids or other measures of location have been individually determined.
As is well known, classification models are often used to classify and characterize transients or other data patterns. If classification models are used to classify transients in the present invention, then steps 428-436 are performed. In step 428, the processor 310 determines the nominal location of each classification model. The measure used in step 428 for determining the locations of the classification models should be the same as used in step 418 to determine the locations of the transients.
In step 430, the processor 310 measures the difference between the location (determined in step 418) of the transient which is being classified (called herein the transient of interest) and the location of each classification model (determined in step 428). The distance can be measured by such means as described in step 422.
In step 432, the processor 310 brings each classification model into alignment with the transient. If the classification model is analytic, such alignment is performed by adjusting the location parameter of the classification model to compensate for the difference measured in step 430. For empirical classification models, such alignment is performed by interpolating the classification model, and then synthetically resampling to bring the classification model into alignment with the transient. Alternatively, if desired, the transient can be realigned, in step 433, by such means as that described earlier in which the transient if interpolated by a (sin(π(n-k))/(π(n-k)) function. This may be easier for classification purposes, since only one alignment is required. Otherwise, several classification models must be realigned.
In step 434, the processor 310 measures the difference between each aligned model and the transient of interest. This assumes the transient and the correct model are similar in size. If the transient is a scaled version of one of the models, the models should be scaled to the amplitude and width of the transient, or vice versa, or models and transient can both be normalized in width and amplitude. In step 436, the processor 310 identifies the classification model having the smallest difference with the transient of interest and classifies the transient of interest using the identified classification model. Note that steps 430-436 must be performed for each transient.
As shown in FIG. 4, classification of the transient occurs in parallel to aligning of the analyzing function and reanalysis of the transient. Alternatively, classification may occur after the transient has been reanalyzed to thereby bypass (i.e., eliminate) steps 428 through 434. According to this alternate embodiment, after the data is reanalyzed in step 426, classification (step 436) is directly accomplished on the reanalyzed transient. For particular applications in which the alignment of step 424 cannot be determined with great accuracy, the full classification procedure of FIG. 4, including steps 428-436, are used to provide the best results. However, it is expected that the procedure of the alternate embodiment (described in this paragraph) will typically be both accurate and faster.
The alternate embodiment just described may be enhanced by iterating the procedure of FIG. 4. Specifically, steps 410 and 414-426 are iteratively performed (after a first pass through steps 408-426) until the location difference determined in step 422 is less than a predetermined amount. If the location difference of step 422 is greater than the predetermined amount, then another iteration is performed starting with step 424 and proceeding through steps 426, 410, and 414-422. If the location difference is less than the predetermined amount, then step 436 is performed as described above to classify the transient based on the reanalyzed transient (thereby eliminating steps 428-434). Note that after the first pass is performed, step 414 can be bypassed since it will be known whether the signal is modulated. Also, after the first pass is performed, detection in step 410 can be bypassed since it will be known that the transient exists (although it will still be necessary to perform the isolating procedure of step 410). Further, note that it is expected that in many, if not most, situations, a single pass through the steps of FIG. 4 will be sufficient. However, the iterative procedure just described offers the dual benefits of: (1) increased accuracy, and (2) measurement of convergence.
Probabilistic means can be used to enhance the distance measurement, step 430, in the presence of noise. This includes maximum likelihood and minimum distance estimators. For example, align the transient and candidate model using the centroid. At each sample point evaluate the ordinate of the probability function using the probability function of the background noise added to the model value at that sample. Thus for a normal, zero-mean, noise distribution, the (noiseless) model value serves as the expected value (mean). The dispersion is that of the background noise. In other words, a pdf (probability density function) is created for each model value, using the noise pdf centered around the model value as the expected value, and the pdf ordinate for the corresponding transient sample value within that distribution is evaluated. This can be viewed as the likelihood. The total likelihood is the product of the individual likelihoods and, more useful for very small probabilities, the total log likelihood is the sum of the individual log likelihoods. The alignment between model and transient can then be adjusted back and forth, and the alignment with maximum likelihood or log likelihood selected. This process is then repeated for each model. The model exhibiting the minimum distance to the transient after probabilistic alignment, or, alternatively, the one with the greatest total likelihood or log likelihood, is the classification model.
The alternate embodiments described above with regard to FIG. 4 could be applied to the procedure of FIG. 5, which is a flow chart in accordance with a second embodiment of the present invention of a more specific method of stably identifying and characterizing a transient contained in a first data set, representing physical data, where the analyzing function is the discrete wavelet transform, and the alignment of the sample data points with the transient is arbitrary. The flow chart of FIG. 5 is generally similar to the flow chart of FIG. 4. In particular, step 508 of the second embodiment is generally the same as step 408 of the first embodiment, except that in step 508 the processor 310 analyzes the signal (or first data set) using a discrete wavelet transform.
Step 510 of the second embodiment is generally the same as step 410 of the first embodiment, except that in step 510 the processor 310 must process multiple scales.
Step 512 of the second embodiment is generally the same as step 412 of the first embodiment, except that in step 512, the information needed to determine whether the analyzed signal is modulated may be available directly from the wavelet decomposition. For example, if the energy is much higher on a scale (except for the lowest frequency scale used) relative to its neighboring scales, then a modulated signal is present.
Steps 514 and 516 of the second embodiment are generally the same as steps 414 and 416 of the first embodiment.
Step 518 of the second embodiment is generally the same as step 418 of the first embodiment, for measuring the transient location (and model location), except that in step 518 two-dimensional measures may be used to accommodate the fact that the energy is distributed over two dimensions (time and scale). For example, the bivariate centroid can be calculated using the following formula: ##EQU5## and for real variables: ##EQU6## where t is time, a is the scale, S(.,.) is the wavelet transform value, E(.) is the expected value, and the denominator takes care of normalization. For real, positive variables, such as amplitudes, unsquared values are often used. In the second embodiment, the time location estimate can be enhanced by using only those scales in which the signal predominates. The estimated scale location can be useful for classification.
Steps 520 and 522 in the second embodiment are generally the same as steps 420 and 422 in the first embodiment, except that in steps 520 and 522 the processor 310 must process multiple scales.
Step 524 in the second embodiment is generally the same as step 424 in the first embodiment, except that in step 524 the analyzing wavelet must be shifted on each scale by the same time increment (not necessarily by the same number of samples since the sampling interval changes on different scales in some decompositions). Also, the whole wavelet is shifted as a rigid unit: the real, imaginary, and envelope. This applies also in those cases for which an explicit wavelet is not used, such as multi-resolution analysis. In multi-resolution, low and high pass quadrature mirror filters are utilized in a coordinated fashion to accomplish the same end effect as if bandpass filters (analyzing wavelets) had been used, such that the high- and low-pass filter impulse responses are convolved with the signal. To be equivalent to a shift in the wavelet, the impulse response must be shifted with respect to the samples.
As an alternative to step 524, step 525 may be performed when it is easier to shift the transient data pattern, as opposed to shifting the analyzing function, as previously discussed in reference to step 425.
Step 526 of the second embodiment is generally the same as step 426 of the first embodiment.
Steps 528-536 of the second embodiment are generally the same as steps 428-436 of the first embodiment; however, if classification models (also called prototypes) are developed in the time domain, then they must be decomposed by means of the wavelet transform. Prototypes can also be developed in the time-scale plane if sufficient knowledge is available about their time-scale plane characteristics. If the prototypes are developed in the time-scale plane, then further decomposition is not necessary.
Also, for each prototype, the processor 310 must decide whether to use or not use the envelopes of the scale data. If envelopes are used for the transients (calculated in step 516), then envelopes must also be used for the prototypes. Thus, depending on the nature of the signal, the prototype is characterized by either a regular wavelet decomposition or a decomposition that has been demodulated on each scale.
Additionally, step 534 must take into account multiple scales. Step 534 preferably operates as follows. The difference in value between corresponding samples of each prototype is squared, multiplied by the sampling interval, and then summed over all samples. This is done on each scale for each prototype, and the results summed over all scales. Finally, the square root is taken of the grand total. The operation just described represents the numerator (less the factor 100) in equation (6) below. Low-pass residue may also be used in the distance measure (step 530), although it is not shown in equation (6).
The percentage difference is calculated (in step 534) as the distance (or difference) divided by the sum of the norms of the two functions, times 100, as shown in equation (6): ##EQU7## Wfj and Wgj are the wavelet transforms of functions f and g on scale j, and nΔt is the discretized time. By using equation (6), the percent difference always remains between 0% and 100%, inclusive. The low-pass residue may also be used to calculate the norms, and should be used if the distance measure utilizes it.
The exact order of detection (i.e., step 510) and modulation (i.e., steps 512, 514, and 516) is not crucial, nor are the exact measures used for location (steps 518, 520, and 528) or difference in location (i.e., steps 522 and 530). The optimal approach may vary with particular applications. The above mathematical formulas are presented herein for illustrative purposes only and are not limiting.
The present invention, and in particular the second embodiment of the present invention, shall now be described in greater detail with reference to the waveforms presented in FIGS. 6a-6k.
FIG. 6a illustrates a signal 602 which is to be analyzed by the system 302 in accordance with the second embodiment of the present invention (i.e., FIG. 5). The signal 602 of FIG. 6a contains considerable noise and two barely visible transients.
FIG. 6b illustrates a wavelet transform decomposition of the signal 602 after the processor 310 has performed step 508 (as described above). For simplicity purposes, only three scales are shown. Note in FIG. 6b that the signal to noise ratio has increased due to the decomposition because, in the example of FIG. 6b, the noise is divided among all three scales, each signal resides on only one scale.
As noted above, the order of steps 510, 512, 514, and 516 is not important. In this example, the processor 310 performs steps 512, 514, and 516 before performing step 510. The result of performing steps 512, 514, and 516 is shown in FIG. 6c where the envelopes are indicated by 604.
The processor 310 then performs step 510 to detect the transients in the demodulated signal, as shown in FIG. 6d. Step 510 is performed by comparing the envelopes 604 with thresholds 606. In the simplest implementation, transients are defined as the regions where the envelopes 604 exceed the thresholds 606. Conditions can be associated with this, such that the envelopes 604 must exceed the thresholds 606 for a certain minimum number of samples (in order to avoid irrelevant spikes, for example), or for less than a certain maximum number of samples (in order to discriminate transients from trends, for example). There are many well-known, appropriate detection methods which may be applied.
For illustrative purposes, FIG. 6e depicts the transients themselves (in this case, meaning that part of the envelopes 604, on each scale, that is above the thresholds 606).
Also in step 510, the processor 310 isolates the transients, as depicted in FIG. 6f. Isolation of the transients is relatively simple if the transients are well separated. Overlapping transients on the same scale can be recognized if their combined length is too long for any kind of transient of interest. If such overlapping transients are discovered, they are separated by, for example, splitting them at their minimum point, or by treating them as mixed distributions and using such well-known techniques as expectation maximization.
FIG. 6g depicts the locations 608 of the transients as determined by the processor in step 518. The locations are calculated over multiple scales if the transients cover multiple scales. Preferably, the bivariate centroid is used.
For illustrative purposes, FIG. 6h depicts the real part of the analyzing wavelet, wherein the envelope and the sinusoid peak occur on the same sample. Note that the analyzing wavelet is generated from a mathematical expression and thus can be centered at any location with respect to the sampling. FIG. 6h generally represents the operation of step 520.
FIG. 6i depicts the shifted and aligned analyzing wavelet after operation of steps 522 and 524. In the example of FIG. 6i, the analyzing wavelet has been shifted by the amount that the analyzing wavelet centroid (determined in step 520) is different from the nearest sample in order to bring the center of the analyzing wavelet into alignment with the center of the transient. This is done individually for each transient.
In step 526, the processor 310 reanalyzes the original signal 602 using the shifted analyzing wavelets. An upper transient 610 of the reanalyzed signal is shown in FIG. 6j.
With regard to steps 528-536, FIG. 6j shows the upper transient 610 being compared to a prototype 612 of the wrong class. FIG. 6k shows the upper transient 610 being compared to a prototype 614 of the correct class.
The present invention may be embodied in many forms, for example, as a system to detect rubbing in a journal bearing supporting a steam turbine shaft. The system includes one or more accelerometers mounted on a bearing cap. The contact that occurs during a rub is momentary, and is in the nature of an impulse. The rubbing impact being very brief in time, creates a broad spectrum and excites the natural frequencies in the bearing housing. The system of the present invention detects and analyzes these natural frequency transients. Note that in the earliest stages the rub occurs only occasionally and has low-energy content.
The system of the present invention further distinguishes the rubbing from the following transients: a) electrical transients that occasionally disturb the electrical circuit; b) transients caused by sudden changes in generator load; c) transients caused by sudden motion of injection and extraction valves in the turbine. The latter effects produce a weak reaction at the bearing, but may be in the same order of magnitude as rubbing in incipient stages. In order to automatically distinguish these transients, classification models may be developed in the time-frequency or time-scale plane. In the latter case either an analytic or empirical example of each type of transient is developed or captured and this is analyzed by means of, for instance, the well known affine wavelet transform. The distribution of energy following scale decomposition helps to distinguish the transients (some transients will have more energy on some scales while others may have more energy on other scales). Even after this, however, some transients may be hard to distinguish. Here information regarding the exact shape of the transient may be crucial. To reduce the amount of data without losing essential information, the signal on each scale may be demodulated with a high quality technique such as use of the analytic signal.
Under ideal conditions two transients (the rubbing and electrical interference transients, for example) may be very close in wavelet coefficients but distinguishable. They may take different times to reach full amplitude and to ringdown, and one may have a secondary peak. However when the actual signals are seen under realistic conditions, in which they are arbitrarily positioned with respect to the time sampling function, they may not be distinguishable at all. Because of the shift variance of the discrete wavelet transform, the actual coefficients may differ substantially from the ideal case.
To show how the system of the present invention corrects this problem, consider the case of a Gaussian-shaped demodulated signal on one particular scale of the wavelet transform.
Referring to FIG. 7, assume it is in the first instance centered on a sample so that one of the samples falls at the mean location. Assume further that other samples fall on the standard deviation locations. (None of these details are important to the basic principles disclosed herein.) If only samples having magnitudes of 0.01 or greater are considered, then the sample values are as indicated in the first row of Table 1, which represents coefficients for Gaussian transient, before and after shift.
TABLE 1______________________________________X Abscissa 1 2 3 4 5 6Y Centered 0.054 0.242 0.3989 0.242 0.054Y Offset by 0.5 0.0175 0.1295 0.3521 0.3521 0.1295 0.0175______________________________________
The coefficients in Table 1 are the coefficients after time-scale analysis including demodulation. It can be seen how both the number of coefficients as well as the coefficients themselves have changed because of a simple half-sample shift. If, however, the centroid is calculated (in this case by dividing the sum of the XY products by the sum of Y), it will indicate that the center of the transient is 3 in the first case (ΣY=0.9909, ΣXY=2.9727, ΣXY/ΣY=3) and 3.5 in the second case (ΣY=0.9982, ΣXY=3.4937, ΣXY/ΣY=3.5). Consequently, the results can be stabilized by reanalyzing the data according to the present invention as described above; that is, by shifting the analyzing function by 0.5 for the second case. Since the analyzing function (such as Morlet's wavelet, for example) is expressed analytically as a function of the abscissa, it is usually easy to shift the function by a simple adjustment of the abscissa. When the offset data is reanalyzed, it will have the same coefficients as the centered case and thus will be easily identifiable as belonging to the same class as the centered case. See FIG. 8.
To show that the method of the present invention and as described above works for non-symmetrical transients, a non-symmetrical triangle is used as another specific example. Consider a triangular shaped transient which rises quickly from zero to 4 by means of the rule: Y=2X, 0≦X<2, but then descends slowly to zero by means of: Y=6-X for 2<X≦6. See FIG. 9.
Table 2 tabulates data from FIG. 9. If the transient of FIG. 9 is shifted by, for example, 0.122 units to the left, the Y values at the sample points will be changed, as shown in Table 2 (also see FIG. 10). Yet if the centroid is calculated as described above, it will be found to have shifted by the same amount. For the original signal, (ΣY=12, ΣXY=32, ΣXY/ΣY=2.667), and in the second case (ΣY=12, ΣXY=30.656, ΣXY/ΣY=2.5546). Again, the instability due to the shift can be corrected by reanalyzing the transient according to the present invention as described above, this time with the analyzing function shifted by the amount determined by the shift in centroid.
TABLE 2______________________________________X Abscissa 0 1 2 3 4 5 6Y Original Position 0 2 4 3 2 1 0Y Offset by 0.112 0.224 2.224 3.888 2.888 1.888 0.888 0sample______________________________________
Therefore, the system of the present invention provides a simple, effective way of correcting for the shift sensitivity of the discrete signal processing algorithm. In the case of monitoring the rub in the journal bearing in a steam turbine, the present invention provides stable, repeatable coefficients that can be compared with previous occurrences of the rubbing, thereby allowing for correct classification so that the transient of interest is not confused with other transients. Thus, the present invention provides useful and important diagnostic information, and useful and important maintenance information because it allows for trending by determining if the transients are slowly changing in shape or growing in magnitude. The latter may indicate the rubbing is becoming more of a danger in terms of the potential to cause catastrophic failure. If shift correction is not provided for, different samplings of the same transient can give the misleading appearance of a change in magnitude, or worse, can cause confusion with other types of transients.
The present invention may be embodied in other specific forms without departing from the spirit or essential attributes thereof and, accordingly, reference should be made to the appended claims, rather than to the foregoing specification, as indicating the scope of the invention.
|Cited Patent||Filing date||Publication date||Applicant||Title|
|US5048102 *||Sep 14, 1988||Sep 10, 1991||Commissariat A L'energie Atomique||Multiple interpolation process for image correction|
|US5101441 *||Sep 20, 1989||Mar 31, 1992||Yamaha Hatsudoki K.K.||Method of and device for processing image|
|US5272723 *||Apr 26, 1991||Dec 21, 1993||Fujitsu Limited||Waveform equalizer using a neural network|
|US5517585 *||May 5, 1993||May 14, 1996||Liberty Technologies, Inc.||System and method for stable analysis of sampled transients arbitrarily aligned with their sample points|
|1||I. Daubechies, "Orthonormal Bases of Wavelets With Finite Support-Connection With Discrete Filters", Wavelets: Time-Frequency Methods and Phase Space (Combes, Grossmann and Tchamitchian (Eds.), Springer-Verlag, Second Edition), (1989), pp. 38-66.|
|2||*||I. Daubechies, Orthonormal Bases of Wavelets With Finite Support Connection With Discrete Filters , Wavelets: Time Frequency Methods and Phase Space (Combes, Grossmann and Tchamitchian (Eds.), Springer Verlag, Second Edition), (1989), pp. 38 66.|
|3||Mallat, "Dyadic Wavelets: Energy Zero Crossings", Technical Report MCS-8219196-A02 submitted to IEEE Transactions on Information Theory, (1988), pp. 1-56.|
|4||*||Mallat, Dyadic Wavelets: Energy Zero Crossings , Technical Report MCS 8219196 A02 submitted to IEEE Transactions on Information Theory, (1988), pp. 1 56.|
|5||*||Mallat, Multifrequency Channel Decompositions of Images and Wavelet Models:, IEEE Trans. ASSP, 37 12, (1989) pp. 2091 2110.|
|6||Mallat, Multifrequency Channel Decompositions of Images and Wavelet Models:, IEEE Trans. ASSP, 37-12, (1989) pp. 2091-2110.|
|7||Oppenheim et al., "Chapter 10: Discrete Hilbert Transforms", Discrete-Time Signal Processing, (Oppenheim, Ed., Prentice Hall, Englewood Clifs, NJ), (1989), pp. 676-687.|
|8||*||Oppenheim et al., Chapter 10: Discrete Hilbert Transforms , Discrete Time Signal Processing, (Oppenheim, Ed., Prentice Hall, Englewood Clifs, NJ), (1989), pp. 676 687.|
|9||Tuteur et al., "Chapter 2: Fourier Methods", Modern Electrical Communications: Theory and Systems, (Prentice-Hall, Inc., Englewood Cliffs, NJ), (1979), pp. 66-73.|
|10||*||Tuteur et al., Chapter 2: Fourier Methods , Modern Electrical Communications: Theory and Systems, (Prentice Hall, Inc., Englewood Cliffs, NJ), (1979), pp. 66 73.|
|Citing Patent||Filing date||Publication date||Applicant||Title|
|US6052523 *||Aug 4, 1997||Apr 18, 2000||Mitsubshiki Denki Kabushiki Kaisha||Method of transient analysis considering time steps and a device for implementing the method|
|US6430519 *||Jun 17, 1999||Aug 6, 2002||Philips Electronics No. America Corp.||Method for evaluating an input data signal and circuit system for carrying out said method|
|US6434494 *||Mar 30, 1999||Aug 13, 2002||Simmonds Precision Products, Inc.||Pressure based fluid gauging system|
|US8390513 *||Nov 14, 2008||Mar 5, 2013||Qualcomm Incorporated||GNSS receiver and signal tracking circuit and system|
|US20070285081 *||May 16, 2006||Dec 13, 2007||Carole James A||Method and system for statistical measurement and processing of a repetitive signal|
|US20100210206 *||Nov 14, 2008||Aug 19, 2010||Qual Comm Incorporated||Gnss receiver and signal tracking circuit and system|
|U.S. Classification||382/291, 382/280, 708/403|
|Cooperative Classification||G06K9/00496, G06F17/148|
|Feb 28, 1996||AS||Assignment|
Owner name: LIBERTY TECHNOLOGIES, INC., PENNSYLVANIA
Free format text: ASSIGNMENT OF ASSIGNORS INTEREST;ASSIGNOR:DOWLING, MARTIN JOSEPH;REEL/FRAME:007886/0118
Effective date: 19960227
Owner name: DEERE & COMPANY, ILLINOIS
Free format text: ASSIGNMENT OF ASSIGNORS INTEREST;ASSIGNORS:STEPHENSON, ROGER DALE;JORDAN, MARK ALAN;PARSONS, STEPHEN KENNETH;REEL/FRAME:007900/0851;SIGNING DATES FROM 19960131 TO 19960212
|Jun 8, 1998||AS||Assignment|
Owner name: SILICON VALLEY BANK COMMERCIAL FINANCE DIVISION, C
Free format text: SECURITY AGREEMENT;ASSIGNOR:LIBERTY TECHNOLOGIES, INC.;REEL/FRAME:009207/0750
Effective date: 19980512
|Dec 3, 1999||AS||Assignment|
Owner name: CRANE NUCLEAR, INC., GEORGIA
Free format text: MERGER;ASSIGNOR:LIBERTY TECHNOLOGIES, INC.;REEL/FRAME:010415/0685
Effective date: 19990420
|Mar 12, 2002||REMI||Maintenance fee reminder mailed|
|Aug 26, 2002||LAPS||Lapse for failure to pay maintenance fees|
|Oct 22, 2002||FP||Expired due to failure to pay maintenance fee|
Effective date: 20020825