CA2276571A1 - System and method for diagnosing and controlling electric machines - Google Patents

Abstract

A system and method of diagnosing and controlling electric machines such as engines is provided. The presence of an instability in a physical characteristic of the device is diagnosed by providing a digital input signal representative of the physical characteristic to a computer system. A wavelet transform is performed on the digital input signal to provide a set of wavelet coefficients, and at least a first probabilistic measure is determined with respect to the coefficients. Preferably, more than one probabilistic measure of the wavelet coefficients are determined and used to identify the precursor. The probabilistic measures may be dispersion or other measures provided by a high rank correlation matrix.
The probabilistic measure(s), and possibly other supplemental device related information, are analyzed to identify a precursor associated with the instability. By identifying the precursor prior to the development of a malfunction in the device, a command execution unit may be used to control the device, in accordance with an instability control approach, to avert the development of a malfunction. The system and method may be applied to devices in real time under time varying, as well as stable or unstable, operating conditions.

Description

BP #11031-001 BERESKIN & PARK CANADA
Title: System and Method for Diagnosing and Controlling Electric Machines Inventors: Igor M. Dremine Victor I. Furletov Oleg V. Ivanov Vladimir A. Nechitailo Vladimir G. Terziev Title: System and Method for Diagnosing and Controlling Electric Machines FIELD OF THE INVENTION
The present invention relates to a system and method of monitoring and diagnosing instabilities indicative of impending malfunctions or failure in electric machinery such as engines, compressors or turbines. The present invention further relates to a system and method for controlling such devices to avoid operational malfunction or failure.
BACKGROUND OF THE INVENTION
Many different systems have been developed for analyzing and monitoring the operation of electric machines such as engines, motors, turbines, compressors, pumps, and other electromagnetic rotating machinery. These systems attempt to diagnose or predict impending failure or malfunction of the electric machine due to instability, and, where possible, to avert that failure by controlling the device accordingly.
For these purposes, diagnostic systems are usually required to operate in an on-line manner to provide continuous monitoring of device operation.
Diagnostic systems are generally based on the reading, recording, and subsequent analysis of physical signals. Sensors may be used to obtain actual measurements of vibrations, deformations, pressure variations, temperature, acoustic noise, and other similar physical phenomenon associated with an operating electric machine. Many failures in electric machines are accompanied and/or preceded by discernable changes in such physical signals. For example, the failure of an electric induction motor may result in, or may be the result of, abnormally high vibrations in the motor.
For instance, prior art systems have been used to measure the vibration levels of an electric machine and then continuously monitor the overall vibration in either the time or the frequency domain. Where required in these prior art systems, the time domain vibration data are converted into the frequency domain by a Fourier transform (typically the fast Fourier transform or FFT is used for computational purposes). The measured vibration signal or the transformed vibration signal, often referred to as the vibration signature of the electric machine, may be continuously or periodically analyzed and compared to a reference signature, reference threshold levels, and/or a reference statistical measure for vibrations of a normally operating machine. See, for example, United States Patent No. 4,184,205 to Morrow, United States Patent No. 4,366,544 to Shima et al., United States Patent No. 5,251,151 to Demjanenko et al., and Max, "Methods and technique of analysis of signals in physical measurements", M., Mir, 1983; v.l, ch.2, pp 18-35.
Some prior art systems have also been used to diagnose electric machines based on analysis of the electric current or power supplied to the machine, since many problems in electric machines may result in harmonic electric currents or other power variations. For example in United States Patent No. 5,629,870 to Farag et al., harmonic current analysis is performed on an operating motor by monitoring the spectral content of the power signature. Spectral components of the power or current supplied to the motor are associated with device operating conditions by referencing a stored knowledge of operational characteristics of the motor.
Similarly, United States Patent No. 5,587,931 to Jones et al.
describes a tool monitoring system which operates in two distinct modes.
In a learning mode, the system gathers statistical data on the power consumption of normally operating tools of a selected type. The power consumption signal of the machine tool is decomposed into time-frequency components by wavelet transform analysis and then reconstructed based on certain selected components to reduce the effects of noise. A statistical power threshold function is then generated based on the selectively reconstructed power consumption signal. Subsequently, in a monitoring mode, the power consumption of a tool is compared to the power threshold. Although the wavelet transformation simultaneously provides information on the power consumption signal in both the time and the frequency domains, the wavelet analysis is for the limited purpose of establishing power level thresholds for a normally operating power signature with reduced noise (i.e. in the learning mode). In the monitoring mode, only a time domain power signature signal is analyzed for diagnostic purposes.
In prior art diagnostic systems for electrical machines, statistical or algorithmic analysis of solely the time domain information or solely the frequency domain information present in a signal generally provides only a very short precursor or advance warning of impending device failure. The very short notice provided by such prior art systems is generally insufficient to provide any opportunity to avoid the predicted failure. Also, the effectiveness of the statistical or algorithmic analysis in prior art systems based on Fourier spectral analysis may deteriorate at high frequencies, because at high frequencies the physical signal is not measured at sufficiently short time intervals.
Furthermore, these prior art systems require that information be compiled for a relatively long time, leading to a delayed response in regulation of the system. For instance, the statistical analysis in some prior art diagnostic systems is often based on variations in the dispersion (standard deviation or variance) of the envelope of the spectrum, the main frequencies of oscillation, and the quasi or harmonic frequencies: see for example Karasev, Maksimov, and Sidorenko, "Vibrational diagnostics of aircraft engines", M., Mashinostroenie, 1978, ch.3, pp 60-83, and Maksimov and Rodov, "Methods and tools for diagnosing unstable flows in compressors", Turbines and jet devices, N12-1280, M., CIAM, 1990, p.
132. This type of analysis is complex and requires the accumulation of information or data over a lengthy period time. This delays the time for regulation of the system. Therefore, for practical purposes these systems are only applicable to the analysis signals for stationary applications which are at rest (e.g. in a laboratory), and they are not suitable for performing on-line or real-time analysis for conditions which are not at rest (e.g. in-flight analysis for an aircraft) and require rapid on-line analysis.
Systems for automatically controlling electric machines, such as engines, to avoid potential breakdown may include a sub-system for diagnosing impending device malfunctions. Ideally, the diagnostic subsystem recognizes an instability or potential malfunction at definite stages of device operation, and subsequently the control system automatically carries out appropriate operations on the device to prevent a malfunction or failure from developing.
For instance, systems have been designed to use an adaptive digital system to automatically control a turbine engine in an attempt to provide full hydrodynamical stability and prevent instabilities in turbine operation. However, the diagnostics of such systems, similar to the prior art diagnostic systems described above, are limited to the analysis of the flight conditions only, and cannot provide sufficient precursors or warnings of engine failure resulting from internal device problems. As a result, regulation by prior art control systems leads in most cases to extremely restricted regimes in which the device or engine is operated far from its full capacity. While optimal algorithms for entering into high capacity regimes have been used, these require very extensive analysis and consideration studies and can, nevertheless, lead to operating decisions which are far from ideal. Other control systems have been proposed that react to abruptly cut-off the fuel supply when a malfunction is detected, for example auto-oscillations arising in an engine (see Waters, "Digital controller applied to the limitation of reheat combustion roughness", Proc.
of AGARD conference, 1974, N15-1). In attempting to prevent fully developed oscillations or other malfunctions from negatively influencing the engine, such abrupt interruptions in fuel supply serve to drastically reduce the engine's capacity and often lead to unpredictable results.
Thus, there is a need for a diagnostic system which accurately and reliably predicts instability or impending malfunction in an electrical device with sufficient warning to allow an associated control system to automatically control the operating regime of the electrical device in a safe, rapid, effective, and efficient manner.
SUMMARY OF THE INVENTION
In a first aspect, the present invention provides a method of diagnosing a device to determine the presence of an instability in a physical characteristic of the device, the method comprising the steps of:
(a) receiving a digital input signal representative of the physical characteristic; (b) performing a wavelet transform on the digital input signal to provide a set of wavelet coefficients; (c) determining at least a first probabilistic measure of at least a portion of the set of wavelet coefficients;
(d) analyzing at least the first probabilistic measure to identify a precursor associated with the instability.
Advantageously, steps (a) to (d) may be performed in real time with the identification of the precursor occurring prior to the development of a malfunction in the device.
The precursor may be characterized by a significant change in the probabilistic measure. Preferably, the first probabilistic measure represents the dispersion of the wavelet coefficients at a particular scale.
In a preferred embodiment, step (c) further comprises determining at least a second probabilistic measure of at least a portion of the set of wavelet coefficients, and step (d) further comprises the step of analyzing at least the first and second probabilistic measures to identify a precursor associated with the instability. In this embodiment, the precursor may be characterized by a significant change in both the first probabilistic measure and the second probabilistic measure. Determining a second probabilistic measure preferably includes the step of determining a high rank correlation matrix for the wavelet coefficients.
Step (b) may comprise performing a discrete wavelet transform or a discretized continuous wavelet transform on the digital input signal. Conveniently, the method also includes the steps of measuring the physical characteristic of the device with a sensor to provide an analog input signal and converting the analog input signal into the digital input signal.
In a second aspect, the present invention provides a method of controlling a device driven by a command execution unit comprising the steps of (i) diagnosing an input signal as described above) and (ii) if a precursor associated with the instability is identified in step (i), directing the command execution unit, in accordance with an instability control approach, to control the device to avert the development of the malfunction in the device. In one embodiment, the device forms part of a mufti-functional system and identifying a precursor associated with the instability further comprises analyzing additional data describing the operation of the mufti-functional system.
In a third aspect, the present invention provides a system for diagnosing a device to determine the presence of an instability in a physical characteristic of the device, comprising a computer system for receiving a digital input signal representative of the physical characteristic, the computer system comprising (a) a wavelet coefficient generation module for performing a wavelet transform on the digital input signal to provide a set of wavelet coefficients, (b) a first measure algorithm module for receiving the wavelet coefficients and determining a first probabilistic measure of at least a portion of the set of wavelet coefficients, and (c) an analysis module for receiving and analyzing at least the first probabilistic measure to identify a precursor associated with the instability prior to the development of a malfunction in the device.
Preferably, the computer system further comprises a second measure algorithm module for receiving the wavelet coefficients and determining a second probabilistic measure of at least a portion of the set of wavelet coefficients, and wherein the analysis module further receives and analyzes the second probabilistic measure to identify the precursor.
Also, the first measure algorithm module may determine the dispersion of the wavelet coefficients at a particular scale, and the second measure algorithm module may determine a high rank correlation matrix of the _7_ wavelet coefficients to enable further determination of other probabilistic measures.
In a fourth aspect, the present invention provides a system for controlling a device driven by a command execution unit comprising a system for diagnosing an input signal as described above, wherein the computer system is connected to the command execution unit and the computer system further comprises an automatic control module for directing the command execution unit, in accordance with an instability control approach, to control the device to avert the development of the malfunction. In one embodiment, the device forms part of a multi-functional system including a data unit containing supplemental information about the multi-functional system, the data unit being connected to the computer system for providing the supplemental information to the automatic control module.
The device may be one of the following: an engine, a motor, a turbine, a compressor, a pump, or other electric machine. The physical characteristic of the device may also be one of the following: vibrations, deformations, pressure, acoustic noise, temperature, or power consumed.
The objects and advantages of the present invention will be more clearly apparent with reference to the remainder of the description and the accompanying drawings.
BRIEF DESCRIPTION OF THE DRAWINGS
In the drawings which illustrate, by way of example, preferred embodiments of the invention:
Figure 1 illustrates the concept of time and frequency resolution for wavelet transform analysis;
Figure 2 shows an exemplary plot of a discrete wavelet transform (DWT) scalogram;
Figure 3 illustrates the pyramidal DWT computation technique;

_8_ Figures 4A-4D show the waveforms of several analyzing or mother wavelet window functions;
Figure 5 is a block diagram operational overview of the diagnostic and control system 100 according to a preferred embodiment of the present invention; and Figures 6-8 show examples of the diagnostic precursors of instability which may be provided by the present invention for an aircraft engine having an axial mufti-stage compressor operating at different regimes.
DETAILED DESCRIPTION OF THE INVENTION
Prior art diagnostic systems for electric machines analyze a time varying signal associated with the machine (such as vibration, pressure, or power consumed) by considering, after spectral conversion, only the frequency domain information present in the signal or by considering only the temporal domain information inherently present in the signal. Thus, prior art systems seek to make conclusions from the frequency spectrum without taking into account simultaneous information about the time properties of the signal, or conversely assess the time content of a signal without concurrently considering frequency information. In the former case, this limits the analysis of the time evolution and of the correlation characteristics of the spectrum to the integral form without any local (in time) information. In the latter case, the analysis is entirely localized in time but does not reveal the frequencies involved.
In contrast, the present invention provide a diagnostic system and method based on mufti-resolution wavelet analysis performed on a physical signal associated with an electric machine. Sensors positioned at or near the electrical device may be used to obtain actual time-varying measurements of vibrations, deformations, pressure variations, temperature, and/or acoustic noise for example. Alternatively, or in addition, analysis can be performed on the current or power provided to the device. In general, any physical phenomenon that is associated with an operating electric machine can be used for analysis.
In accordance with the present invention these signals are digitized and provided to a computer or processor which carries out the signal analysis by using a wavelet transform that simultaneously provides both time and frequency (or scale) information about the physical signal.
Wavelet analysis is discussed in detail in: Rioul et al., "Wavelets and signal processing", IEEE Signal Processing Magazine, October 1991, p. 1-38;
Daubechies, "Ten Lectures on Wavelets", Society for Industrial and Applied Mathematics Press, vol 61, CBMS-NSF Regional Conference Series in Applied Mathematics, 1992; Kaiser, A Friendly Guide to Wavelets (6th), The Virginia Center for Signals and Waves, Birkhauser, Boston 1994; Nievergelt, Wavelets Made Easy, Eastern Washington University, Birkhauser, Boston, 1999; and Polikar, "The Wavelet Tutorial Parts III and IV" available in June 1999 at uniform resource locator http://www.public.iastate.edu /~rpolikar/WAVELETS /WTtutorial.html.
The contents of these references are hereby incorporated into the present description for background purposes.
Generally, wavelet analysis allows data to be analyzed at different scales or resolutions (mufti-resolution analysis). Mufti-resolution wavelet analysis first requires choosing an analyzing or base wavelet function (also referred to as a "mother wavelet"). The wavelet transform deconstructs or decomposes or decomposes or decomposes the original time domain signal into scaled and shifted (or translated) versions or windows of the base wavelet. The original signal can be represented by a linear combination of coefficients of these scaled and shifted wavelet functions. As a result, signal analysis can be carried out on the wavelet coefficients which conveniently represent a correlation between the wavelet and a localized part of the original signal. Unlike, the short time Fourier transform (STFT) or other time-frequency distribution transforms which have a constant resolution at all times and frequencies, wavelets are much more suitable for studying unpredictable (or choppy) signals having both low frequency components and sharp, high frequency spikes.
Both a continuous wavelet transform (CWT) and a discrete wavelet transform (DWT) may be performed. While the CWT
theoretically transforms a continuous input signal into a continuous wavelet coefficient transform, in any practical analysis both the CWT and the DWT are determined based on a discretized input data stream (and the CWT will also be discretized in practice).
Generally, for a time domain input signal function f(t) decomposed with an analyzing or base wavelet function w(t), the wavelet coefficients W(a,b) for location a and scale b (the scale is effectively 1/frequency) are given as follows:
W(a,b) = f w'(a,b,t)~f(t)~dt t The scaled and shifted wavelet functions w(a,b,t) are generated from the analyzing or mother wavelet w(t) as w(a,b,t) - ~ w~t_a~
b b so that the mother wavelet corresponds to the case a = 0 and b = 1. The function w*(a,b,t) is the complex conjugate of w(a,b,t), and these are equal when w(a,b,t) is real. Usually, a >_ 0, so that as a increases the wavelet window is translated along the time axis - the translation a can be considered the time elapsed since t = 0. The scale variable b is > 0, such that when 0 < b < 1 the wavelet window is compressed and when b > 1 the wavelet window is dilated. The factor 1 /~ in w(a,b,t) normalizes the wavelet so that W(a,b) has the same energy at every scale. Wavelet windows with a high scale (low frequency) provide a global (non-detailed or localized) view of the input signal whereas wavelets with a low scale provide a detailed or localized view of the input. In other words, the wavelet transform provides good frequency but poor time resolution at high scales (low frequencies) and good time but poor frequency resolution at low scales (high frequencies).
Figure 1 graphically illustrates this trade-off between time and frequency resolution for wavelet transforms. Each window 10 in the time frequency plane in Figure 1 corresponds to a value of the wavelet transform. It may be noted that, in accordance with the Heisenberg uncertainty principle, the simultaneous mapping of both frequency and time can never be achieved. Thus, a specific value in the time-frequency plane can never be known, illustrated by each of the windows 10 in Figure 1 having a non zero area. Each of the windows 10 represents an equivalent portion of the time-frequency plane (i.e. they have equal areas, but may have different dimensions), however at higher frequencies (lower scales) the resolution in time is better or less ambiguous than the resolution in frequency. This is illustrated by window 12 in Figure 1 which has a narrow width and large height. Similarly for window 14, the resolution in frequency is better or less ambiguous than the resolution in time. (Note that for the short-time Fourier transform, the width and height of the transform windows does not change since a constant window is used.) The wavelet coefficients for the continuous transform W(a,b) can be thought of as a varying surface or landscape over a two dimensional scale-translation plane whose axes are the location variable a and the scale variable b. At each point, the CWT is essentially a measure of the correlation or similarity in frequency between a specific wavelet window function w(a,b,t) and the input signal f(t). The CWT wavelet coefficients reflect this similarity of the input to the wavelet at a specific scale and translation. Therefore, if the input f(t) has a significant frequency component corresponding to a certain scale and occurring within a certain translation interval, the CWT coefficient for the corresponding point in the translation-scale plane will be relatively large.
It should be noted that the wavelet transform can also be considered a constant relative bandwidth or "constant Q" decomposition that employs wideband windows at high frequencies and narrowband windows at low frequencies. (This notion is considered further when discussing the DWT below.) For each scale value, the ratio of the frequency to bandwidth, denoted as Q, remains constant. For the CWT, the analyzing or mother wavelet can be regarded as the impulse response of a reference bandpass filter. A set of parallel bandpass filters with constant Q can therefore be used to realize a wavelet transform.
While W(a,b) is theoretically calculated in a continuous manner for all values of b > 0 and all values of a >_ 0 (i.e. at every point in the translation-scale plane), such a complete transform is generally not required since a signal is usually bandlimited, as would be the case for discretized or sampled input data. Also, if necessary, the input signal f(t) can be reconstructed or synthesized from the wavelet coefficients W(a,b) by summing, in a linear combination, the products of the wavelet coefficients with the corresponding wavelets window functions, or equivalently f(t) = f f W(a,b)~w(a,b,t)~da~db a b Reconstruction of the input is generally possible if the wavelet functions act as a set of orthonormal bases. This requires that f w(a ,b ,t)~w(a ,b ,t)dt - 1 ifa =a andb =b t 0 otherwise By choosing appropriate wavelet functions as discussed below, these conditions can be satisfied. In applications where orthonormal wavelet functions are not available, synthesis may also be possible using biorthogonal bases functions or the concept of frames.
As mentioned, in practice, the continuous wavelet transform must be based on data sampled at discrete times or intervals, and only a sampled or discretized form of the CWT can be computed. Thus W(a,b) is represented by a matrix having a resolution corresponding to the precision of analysis in the computer algorithm. When sampling the CWT, it is possible to reduce the sampling rate as the scale increases (frequency decreases), while still remaining at or above a critical sampling rate (i.e.
the Nyquist rate). This is beneficial since it reduces the number of necessary computations. If synthesis of the input signal f(t) is required, then the rate for all scales should be at least equal to the Nyquist sampling rate (i.e.
twice the frequency at that scale) and preferably somewhat higher to reduce aliasing. However, if synthesis is not required, computations can be further reduced, depending on the analysis.
In discretizing the wavelet window functions, the scale parameter b is generally discretized first on a logarithmic basis. This can conveniently be accomplished using a base 2 logarithm for computational simplicity. If so, only wavelet coefficients for b = 2, 4, 8, and so on are computed (i.e. on a dyadic basis). The translation parameter a is then discretized with respect to the scale b so that a different sampling rate is used at each of the scale levels (for dyadic sampling the sampling rate would be reduced in half each time the scale b jumps to the next level).
The discretized CWT described above, also referred to as a wavelet series transform, is different from the discrete wavelet transform (DWT) which is discussed below. In general, the shifting process in the discretized CWT remains relatively smooth across the sampled data, unlike for the DWT.
In the discrete wavelet transform (DWT), a and b can be replaced with m and n respectively where m, n are integers (m > 0, n >_ 0), so that W(m, n) = 2 = ~ f(i)~w(i~2-~' - n) =o where f[i] are the discretized samples of the signal f(t) and where M is the total number of samples. As can be seen from the above equation, for the DWT the scaling or dilation is stepped exponentially by a power of 2, while shifting or translation occurs by integer steps. Other scaling and shifting bases or step sizes are of course possible, but computationally the DWT as expressed above is relatively simple and therefore preferable. Technically, the DWT W(m,n) becomes a two dimensional coefficient matrix evaluated at the points m,n in the scale location plane; however, as with the discretized CWT downsampling with increasing scale is usually performed. An exemplary plot of a DWT, also referred to as a scalogram, is illustrated in Figure 2.
The scale variable m is also referred to as the level of resolution, and the number of distinct values which m takes on is the number of levels in the DWT. The length of the input sample stream being analyzed determines how many frequency resolutions can possibly be represented, and these frequency resolutions are referred to as the levels in the wavelet transform. To illustrate, for an input signal f[i] which has been discretized into M = 2h samples, there are potentially h levels in the discrete wavelet transform. After downsampling, the number of coefficients at a given level or scale m is 2h-m (except for the last computed resolution level which has 2h-m+1 coefficients).
For the DWT, the analyzing or mother wavelet function w(t) is generally associated with a high pass filter. For some analyzing wavelets, the function is given by an explicit formula in time (e.g. the HAAR
wavelet), whereas for others the function is obtained from the coefficients of the associated high pass filter (e.g. the Daubechies family of wavelets).
In addition to the analyzing wavelet function w(t), the DWT also uses a scaling function s(t) associated with a low pass filter. Generally the scaling function is used to define the averages in a signal decomposition.
The wavelet and scaling functions each satisfy a scale recursive equation: a scaling function is the weighted sum of translated and compressed versions of itself; and a wavelet is the weighted sum of translated and compressed versions of the corresponding scaling function.
The mufti-resolution DWT uses the scaling function and wavelet function to decompose the input signal into different frequency bands.

The DWT analysis of an input signal f[i] at different resolutions (i.e. different frequencies) can be conveniently and efficiently computed with pyramidal coding techniques. The resolution of the signal is changed by filtering steps and subsequently the scale is varied by downsampling the filtered signals. To illustrate, a preferred pyramidal DWT computation technique is shown in Figure 3.
Referring to Figure 3, the DWT analysis begins by subjecting the input f[i] to a digital low pass filter LP[i] and a digital high pass filter HP[i] at 30 and 32 respectively. The filters 30 and 32 are half band filters which behave as a quadrature mirror filter (QMF) pair, so that the original bandwidth of the sequence f[iJ is effectively divided in two. Generally, the high pass filter provides details about the input signal while the low pass filter mainly provides approximation information. The filtered outputs from 30 and 32 are of the same scale (or rate) as f[i] but provide a different resolution. Since the outputs of the filters 30 and 32 are half the bandwidth of f[i], they are downsampled by two at 34 and 36 without any resulting loss of information (assuming f[i] was originally sampled at a rate at or above the Nyquist rate). Downsampling by 2 is conveniently accomplished by simply eliminating every second sample. The downsampled sequences 39 and 38 from 34 and 36 have a scale that is twice the scale of f[i].
The sequence 38 and the sequence 39 are defined mathematically (both as a function of j) as f[i]~HP[2j-i]
i f[i]~LP[2j - i]
i respectively. Also, conveniently, the QMF high pass and low pass filters are computationally linked to one another by HP[L- (i+ 1)] _ (-1)'LP[i]
where L is the length of each filter.

The high-pass filtered and downsampled sequence 38 provides a first level (m=1) of DWT coefficients. Where the input sequence f[i] has M samples, there are M/2 first level DWT coefficients.
The processing of blocks 30, 32, 34, and 36 is then repeated at blocks 40, 42, 44, and 46 with the low-pass filtered and downsampled sequence 39 as the input sequence. Consequently a second level (m=2) of DWT coefficients is provided at 48 (for M original input samples, there would be M/4 coefficients at this level). This is then repeated on the sequence 49 to generate DWT coefficients for successive levels, as needed. If a full multi-resolution analysis is performed, the last stage (m = h) will generate a one element sequence 59 from low pass filter 50 and down-sampler 54 and a one element sequence 58 from high pass filter 52 and down-sampler 56.
Each of these sequences form the DWT coefficients for the last computed resolution level, and so the h'th level will have 2 DWT coefficients. Thus, the DWT will have the same amount of coefficients, M, as in the original input sequence f[i].
It should be noted that the pyramidal analysis of Figure 3 may be generalized to provide wavelet packet analysis by decomposing both the high-pass filtered output sequences as well as the low-pass filtered sequences. The input is thus transformed into several possible wavelet decompositions. The wavelet packages are particular linear combinations of these resulting wavelets usually optimized by recursive algorithms or the like.
Like the CWT, the DWT localizes time information at high frequencies better than at low frequencies. Since the information content of most signals is generally present at higher frequencies, this is usually advantageous. Each successive wavelet transform level reveals coarser frequency or change information about a larger part of the original input -lower scales corresponding to more rapid variations and therefore to higher frequencies. The information containing frequencies of the original input map to high amplitude DWT coefficients for the portion of the DWT that includes those particular frequencies.

Because of the orthogonality of the DWT, information represented at a certain scale or level m is disjoint from (i.e. does not overlap with) information provided by other scales or levels in the DWT.
As for the CWT, reconstruction of the input from the DWT again follows when the wavelet bases are orthonormal. In particular, with DWT analysis as in Figure 3, reconstruction simply requires carrying out the analysis in reverse order, i.e. the signals at every level are upsampled by two, passed through synthesis filters LP'[i] and HP'[i], and then added. Conveniently, the analysis (LP[i] and HP[i]) and synthesis filters (LP'[i] and HP'[i]) are identical except for a time reversal. If the filters are not ideal half-band filters, reconstruction is more difficult, but an appropriate choice of wavelet function, e.g. the Daubechies wavelet family, can provide for good reconstruction.
DWT analysis is computationally much faster than CWT
analysis, including discretized CWT analysis. For many applications, even a discretized CWT will contain significant redundant information. Also the DWT is often more suitable for diagnostics or computer solution of some equations. The CWT, on the other hand, provides for better pattern recognition: see Afanasyeva, Dremin, Kotelnikov, "Pattern Recognition", Modern Physics Letters A12 (1997) 1185. As a result, both the CWT and the DWT can be usefully utilized for many applications.
To complete the discussion on the CWT and DWT, brief consideration of the types of mother wavelets suitable for use in these transform algorithms is merited. Four common wavelet window functions are shown in Figures 4A-4D to illustrate some general principles regarding the wavelet functions. These are the Daubechies-5 wavelet (Figure 4A), the Daubechies-8 wavelet (Figure 4B), the HAAR wavelet (Figure 4C), and the Mexican-Hat wavelet (Figure 4D).
As can be seen from Figures 4A-4D, the mother wavelet function is generally a small window of finite length and is zero outside a certain time interval (for example between 0 and a certain constant To) .
The wavelet should be compact, as well as oscillatory to satisfy the zero mean requirement. The zero mean requirement for the wavelet is also referred to as the admissibility condition and requires that the integral of w(t) over time be zero. When this condition is satisfied, the wavelet transform will be invertible so that the input can be reconstructed.
Different base wavelet functions make various trade-offs between how compactly the base function is localized in time and how smooth the function is (i.e. how well it approximates polynomials).
Referring back to Figure 1, different wavelets can have different window areas for the windows 10 in Figure 1 - although a lower limit for the area is set by the uncertainty principle.
As mentioned previously, the wavelet function may be given by an explicit formula in time (e.g. the HAAR wavelet, the Mexican-Hat wavelet, the Morlet wavelet) or the function may be obtained from the coefficients of the associated high pass filter (e.g. the Daubechies family of wavelets).
Different wavelets have different properties and a specific wavelet must be selected based on the required analysis. For example, the HAAR wavelet is a discontinuous function, while the Daubechies wavelets are specific families of wavelet functions that are particularly suitable for representing polynomial behaviour (the Daubechies wavelets are generally sub-categorized by a number indicating the number of coefficients in the high pass filter associated with the wavelet). The Mexican-Hat wavelet does not provide an orthogonal analysis, although this is typically not a concern where reconstruction of the input is not required. Many other wavelets, including the Morlet, the Coiflets family, and the Symlets family exist. Indeed, the list of wavelets continues to grow and to be refined, and potentially this list of possibilities is limitless.
For the purposes of the present invention, the specific wavelet function used during wavelet analysis, whether continuous or discrete, is not critical, and the inventors have found that similar results are attainable with various different wavelet functions, such as with the HAAR and the Daubechies wavelets.

As will be understood from the above description of wavelet analysis, the wavelet coefficients (CWT or DWT) for a physically occurring time signal f(t) contain a considerable amount of useful temporal and frequency information about that signal. Since the original signal f(t) varies in time, a set of wavelet coefficients at any given scale or level of resolution will also vary in time (although the mean or average of the coefficients will be zero). While scalograms such as the one in Figure 2 can be plotted, compared, and analyzed for diagnostic purposes, statistical or algorithmic analysis which can be carried out quickly and efficiently by a computer or processor is preferred.
For instance, in Thurner, Feurstein, and Teich, "Multiresolution Wavelet Analysis of Heartbeat Intervals Discriminates Healthy Patients from Those with Cardiac Pathology", Phys. Rev. Lett. 80 (1998) 1544-1547, the authors describe a method of diagnosing a human heart by analyzing the sequence of time intervals between heartbeats. The time series was transformed into a discrete wavelet transform (DWT) and the wavelet coefficient standard deviation (or dispersion) as a function of scale was calculated. The authors determined that, at intermediate scales, the wavelet coefficients for heart failure patients exhibited substantially lower variability or dispersion than for normal patients. The analysis thereby provided a clinically significant measure of the presence of heart failure with 100% sensitivity and 100% specificity.
While the dispersion analysis performed by Thurner et al. is effective for diagnosing cardiac pathology under at-rest (non-changing) conditions, it is generally unsuitable for analyzing the real-time signals or conditions associated with a "non-stationary" device such as an in-flight aircraft engine. The method of Thurner et al. is also limited to the analysis of a single dispersion parameter, and so may be insufficiently reliable for more complex physical processes.
In accordance with the present invention, a physical time varying signal associated with the operation of a device such as an engine or another type of electric machine is analyzed with discrete and/or continuous wavelet transforms to determine one or more statistical or probabilistic measures of the resulting wavelet coefficients. These coefficient characteristic measures generally include at least the dispersion of the coefficients (i.e. their standard deviation a or variance a2) at various scales. By also determining high order or high rank correlation matrices for the wavelet coefficients, other statistical or probabilistic measures, such as the two dimensional equivalents of the one dimensional measures of skewness or kurtosis, can be provided. Technically, correlation matrices with infinite rank form a complete set, and various probabilistic measures can be expressed as combinations of these matrices. Generally, when the correlation matrix series is truncated, certain measures become preferable, depending on the level of truncation.
In a preferred embodiment, more than one probabilistic wavelet coefficient measure is used to diagnose and assess the physical signal corresponding to the electric machine so that a more reliable and accurate diagnosis is achieved with the present invention. Different computer algorithms can be used to determine different measures of wavelet coefficient characteristics, and this can also be done at different resolution levels. Furthermore, this type of diagnostic analysis can also be carried out simultaneously on more than one physical signal (for example, for pressure, vibrations, and temperature) to separately provide additional diagnostic indicators with respect to the operation of a device.
The dispersion of wavelet coefficients is essentially the spread of the coefficients about their mean, and is generally given by the standard deviation a. For example, for a set of DWT coefficients W(m,n) the standard deviation at a given level of resolution m is given by N- 1 ~=o where N is the number of wavelet coefficients at the scale m and N=int(Ml l 2 "' Correlation is a ratio of the covariance of the two variables to the product of their standard deviations. (Correlation values are merely measures of covariances of standardised values. Standardizing a data set zeros the mean and sets the standard deviation to one, so that the covariance value is equivalent to the correlation value.) Bivariate correlation provides a single number which summarises the relationship between two variables. The correlation coefficient indicates the degree to which variation in one variable is related to variation in another. The correlation coefficients which are generally obtained from a least squares approximation technique, measure the degree of linear correlation between the two variables.
Typically, a correlation matrix is a variance-covariance matrix of standardized data. The diagonal of the resulting correlation matrix is the correlation of one variable with itself and should have a value of 1. In accordance with the present invention, multi-dimensional correlation matrices may be calculated from wavelet coefficients by means of a generating function technique: see De Wolf, Dremin, and Kittel, "Scaling Laws for Density Correlations and Fluctuations in Multiparticle Dynamics", Phys. Reports 270 (1996) p. 1. For example, a set of DWT
wavelet coefficients W(m,n) can be used as the coefficients in a polynomial sequence of two continuous variables, a and v, to provide a generating function G(u,v). A representative G(u,v) in powers of a and v then provides the elements of the correlation matrix F(q,p) as follows:
ym~n)u~~" _ ~ (u 1)q(v 1)P F(q~P) rn,ri q,P q!pl The correlation matrix F(q,p) is theoretically of infinite dimension, but in practice it may be truncated at some qmaX and p,naX imposing, for example, some threshold value on its elements. Generally, F(q,p) indicate the correlations between the W(m,n) set of wavelet coefficients with the parameters q and p defining the rank or order of a particular correlation Fqp. It should be noted that the second order correlation matrix (q = 1, 2 and p = 1, 2) provides a measure of the dispersion of the W(m,n) coefficients.
Higher order or higher rank correlation matrices (i.e. q,~,aX > 2 and p,.naX ~
2) provide distinct probability measures (such as the two dimensional equivalents of the one dimensional measures of skewness or kurtosis) and can be used to determine other types of important information about the wavelet coefficients and hence the original input signal.
Figure 5 shows of block diagram operational overview of the diagnostic and control system 100 according to a preferred embodiment of the present invention. The system 100 diagnoses and controls an electrical device 102 such as an engine. The device or engine 102, in known manner operates to run a system 108 which may be, for example, an aircraft. The device 102 may be operating under at rest conditions or under "non-stationary" conditions which change with time. Furthermore, the device or engine 102 may be at different operational regimes defined by operation at a certain percentage of the device's nominal limits. While the device 102 is at or within any of these regimes, the operation may become unstable.
As will be appreciated by those skilled in the art, instability increasingly becomes a concern as operation of the device approaches its nominal limits, i.e. full capacity. The operating regime of the device 102 (and any changes thereto) is generally specified via 136 by the command execution unit 134.
A sensor or transducer 104 (preferably several like transducers are used to ensure consistency and accuracy of a physical measurement) attached to the device 102 provides a time varying signal f(t) corresponding to the operating characteristic measured by the transducer 104. The transducer 104 may measure any suitable physical parameters such as vibrations, deformations, pressure variations, acoustic noise, temperature etc. Alternatively the signal f(t) may simply represent the power consumed by the device 102, in which case the sensor 104 is not needed (although the signal may need to be converted to an appropriate voltage form). The signal f(t) which at least partly characterizes the operation of the device is sampled and converted into a digital signal f[i] at 106. The digital stream f[i] is provided to computer system 110.
Note that while only one signal f(t) is processed in the system 100, the system 100 may separately perform similar analysis on more than one time varying characteristic of device 102 operation. This allows for alternative and/or additional diagnostic capability, as will be clear to those skilled in the art.
The computer system 110 includes a wavelet coefficient generation module 120, a first measure (M1) algorithm module 122, a second measure (M2) algorithm module 124, and an automatic control module 130, as shown. While the computer system 110 is shown as having these separate modules, this is merely illustrative of system operation, and it should be noted that software running on computer system 110 may be organized in a number of different manners. For instance, there may only be a single measure algorithm module which computes all of the probabilistic measures, as required. The computer system 110 also has a CPU or processor and memory resident therein (not shown) for actually carrying out the operations in modules 120, 122, 124, and 130. The wavelet coefficient generation module 120 calculates discretized CWT and/or DWT
coefficients based on the input samples f[i]. As the sequence of samples f[i]
is continually updated, so too are the wavelet coefficients calculated by module 120. Thus the system 100 is capable of running in an on-line or real-time manner for a device which is not at rest. The module 120 may be programmed to use a specific analyzing wavelet function (e.g. a Daubechies wavelet) or it may be programmed to select among various possible options based on available processing capabilities, the type of operational characteristic f(t) represents, and any other relevant criteria.
Similarly, the module 120 may be programmed to select between a DWT
and a CWT algorithm.

Once calculated, the wavelet coefficients from module 120 are provided to algorithm modules M1 122 and M2 124 which each determine a probabilistic or statistical measure of the wavelet coefficients, such as dispersion at a specific scale or a measure provided by a high rank correlation matrix at a definite, i.e. truncated, correlation rank or level.
Note that although two separate measure algorithms are shown in Figure 5, the system may only determine one measure algorithm (for example dispersion). Preferably, however, more than one probabilistic wavelet coefficient measure is used to diagnose and assess the physical signal f(t), so that a more reliable and accurate diagnosis is achieved. Consequently, the system 100 may also include more than two measure algorithms.
The computed probabilistic measures from 122 and 124 are provided to the automatic control module 130 as shown in Figure 5.
Optionally, feedback 126 and 128 from the control module 130 to the algorithm modules 122 and 124 respectively may be provided, so that the module 130 can specify, for example, that dispersion be determined at a particular level or scale or that a correlation matrix of a particular rank be provided. Changing the level of multi-resolution analysis (or the correlation matrix rank) in such a manner allows the control module 130 to determine the most appropriate criteria for making a diagnosis and can also provide the automatic control module 130 with supplementary diagnostic information. Generally, the use of several criteria and thresholds in control module 130 provides a reliable means to mutually validate and corroborate any diagnostic results or conclusions.
It should be noted that a determination of appropriate criteria (e.g. the resolution level) for analyzing and diagnosing the operation of a device 102 may initially need to be based, at least to some extent, on experiments and tests of different criteria with the specific device type.
The automatic control module 130 analyzes the criteria measures for significant changes, variations, or disruptions within these measures which vary with time, as the sequence of samples f[i) is continually updated. Upon doing so, the control module 130 determines what feedback response 132, if any, should be supplied to command execution unit 134 to change the operating regime of the device 102.
Significant variations or disruptions in the probabilistic measures used as diagnostic criteria are generally indicative of an instability which (unchecked) will lead to an impending malfunction or failure in the operation of the device 102. What will qualify as a significant variation in, for example, the dispersion measure may vary from application to application (and from one resolution level to another), and an exemplary illustration of this is provided below. Generally, when stable operation is detected, i.e. no instabilities or impending malfunctions are diagnosed, the control module 130 may direct the command execution unit 134 to maintain the device at the current operating regime. Optionally, the control module 130 could also direct the command execution unit 134 to increase the operating regime of the device closer to full capacity if the device is not already operating sufficiently near full capacity. On the other hand, when an instability or impending malfunction is detected, the control module 130 directs the command execution unit 134 to appropriately reduce the operating regime of the device in accordance with the instability control approach of the command execution unit 134 -thereby preventing further development of the instability into a malfunction or failure in the engine's operation.
The reliability of the system 100 can be further guaranteed by providing multiple checks of the validity of commands 132 sent by the control module 130, in addition to the mutual comparison of results obtained from different probability measures, for example a dispersion analysis of wavelet coefficients at different scales or a dispersion measure and a skewness measure from a high rank correlation matrices. As shown in Figure 5, data unit/recorder 140 of the multi-functional system 108 can accumulate and provide other information (e.g., the flight data in case of an aircraft) to the automatic control module 130. Thus, the computer analysis and control in module 130 via the resulting execution commands 132 can be based on the simultaneous processing of many different parameters and on a set of operations prescribed for the engine regulation (including instability control). The versatility and flexibility of the control system of the present invention is clearly apparent from the ability to analyze several physical parameters, flight or other supplemental data, and different probability measures, as well as in the selection of an appropriate approach to engine regulation from among the many available choices.
The specific instability control approach chosen for regulation of device or engine 102 does not form part of the present invention.
Known approaches include active control, passive control, and avoidance control, and generally selection of the most suitable approach will depend on the specific application. See generally, Georgantas, "A Review of Compressor Aerodynamic Instabilities", National Aeronautical Establishment (National Research Council Canada), 1994.
Regardless of the instability control approach selected, the system and method of the present invention can be used to effectively extend the stable operating range and the scope of regulation of the device 102, allowing the device to be safely and reliably operated closer to full capacity. The device is therefore operated more efficiently with an extension of the range of engine regulation not provided by prior art systems. Conclusions about any corrective adjustment in the operation regime and its realization are made by a compute program, based on the time development of these parameters and on the stored algorithms. With the early detection with precursors of possible failures and by appropriate and rapid corrective responses, the present invention provides high reliability and fast feedback both for at rest and in-motion (time varying) situations and for stable and unstable operating regimes.
The very fast feedback response in the present invention enables a possible failure to be averted. This ability stems from the highly reliable and very timely diagnostic conclusions or precursors of instability/malfunction generated by the invention during operation of the device. The timely prediction by way of precursors of failure is primarily due to the wavelet transform properties which have been advantageously exploited by the present invention (as discussed in detail above). Such rapid determination of highly accurate precursors of impending malfunctions was, hitherto, unavailable in prior art diagnostic systems.
Furthermore, the system and method of the present invention differs substantially from prior art systems in that diagnosis of instability or impending malfunctions can be made under time varying conditions whether the device is operating within stable or unstable regimes. The ability to apply the invention not only to systems which are at rest, but also during the "non-stationary" conditions of device operation (such as for an in-flight aircraft engine) provides a significant advantage.
Figures 6-8 illustrate, by way of example, the highly accurate and very timely diagnostic precursors associated with an instability or impending malfunction in a device, as may be provided by the present invention. In the experimental examples of Figures 6-8, the device is an aircraft engine having an axial multi-stage compressor. The rotation of the rotor is at 76% of its nominal limit in Figure 6, 81% of its nominal limit in Figure 7, and 100% of its nominal limit in Figure 8. Signals from eight pressure sensors positioned at various places within the compressor were recorded and digitized as described above. As is known to those skilled in the art, an axial multistage compressor is susceptible to the formation of a rotating stall which may be precipitated by a distorted inlet flow, and this may also lead to the very serious problem of engine surge. In an attempt to avoid these problems, the compressor is usually operated below full capacity at which these instabilities are less likely to occur.
For the experimental test situations in each of Figures 6-8, the aircraft engine was physically operated under at rest conditions, but an instability was introduced by increasing the pressure behind the compressor through a slow injection of extra air into the compressor intake. After a few minutes, this led to a full blown rotating stall instability in the compressor.

Figures 6-8 show the variation of the pressure within the compressor for an interval of about 5 seconds. In each case, this interval includes the occurrence of the fully developed instability, indicated at 200.
Each of Figures 6-8 show the time variation from one of the pressure sensors at 210 (all pressure sensors provided similar results) and the dispersion (standard deviation) of a set of DWT coefficients was computed from the pressure input. The time varying dispersion is shown at 220. In each case, the DWT was calculated with the HAAR wavelet and dispersion was measured at the 4-fold resolution interval (i.e. at the 4th level of DWT
resolution) and at a definite (truncated) correlation level .
With reference to Figure 6, it may be seen that the value of the dispersion during normal operation of the compressor abruptly decreases at 230. This change is considered to be a precursor of the stall and the possible destruction of the compressor, reflected in the abrupt, large increase in the dispersion (at 240) near the end of the time interval. The smaller values of the dispersion are due to rotational instabilities in the compressor. Note that the increase in the dispersion measure 220 prior to the precursor 230 in Figures 6 and 7 is merely a result of a threshold effect due to the finite length of the wavelet resolution at a given scale.
As indicated, the frequency of the rotation of the rotor was 76% of its nominal limits (n/nli,r, = 0.76) in Figure 6. In this case, the precursor 230 provided about a 2.0-2.5 second warning before the beginning of stall 200. This time interval and the significantly large decrease in the dispersion are more than sufficient for very reliable diagnostics.
Referring to Figure 7, a similar plot with the frequency of the rotor rotation equal to 81% of its nominal limit (n/nlim = 0.81) is provided.
Once again the stall precursor 230 provides about 2.2 sec warning before onset of the rotating stall. Finally, in Figure 8 (at the limiting frequency of the rotor rotation (n/num =1)). the precursor 230 starts earlier than in the previous examples and provides about 1.0 seconds of warning. In each of Figures 6-8, the precursor was reflected by a drop of about 30-40% in the dispersion characteristic.
Comparative analysis of Figures 6-8 demonstrates that the more intensive regimes start developing the instability earlier and possess shorter precursor warnings. However, all of these time warning intervals are at least about 1 second. The main delay in an engine control feedback system such as in the present invention results from the rate of the physical processes and from the steps used to transfer and process information about those processes. The inventors' have found that the measurements from the primary sensors' outputs, the wavelet transformation analysis of the data accumulated over a period of time, the mutual comparisons and checks of diagnostic results, and the generation of appropriate control commands takes approximately up to about 0.2 sec.
Thus, the precursor times in Figures 6-8 are all much longer than the time generally necessary for reliable computer analysis to be carried out, and there is enough time for smoothly regulating the operating regime before a rotating stall starts. In comparison, prior art attempts to predict the development of a rotating stall with a velocity measuring probe such as a hot wire anemometer provide a precursor or warning of only about 10 milliseconds. The unambiguous detection of a precursor to instability, its rapid recognition, and the resulting long time available for preventing the engine failure are principle advantages of the present invention not found in the prior art.
In all three cases of Figures 6-8, correlation matrices can also be used for further and more accurate analysis of the precursors.
Furthermore, measure from correlation matrices can be used to help reveal the physical origin of a stall.
It will be apparent from Figures 6-8 that no significant precursors are present in the direct measurements of the highly irregular pressure variations of the compressor. The diagnostic characteristics of the wavelet coefficients are clearly much different in principle from those present in the time-varying signal itself. Furthermore, Fourier spectrograms of the pressure variations were also examined, and they failed to reveal rotational instabilities of sufficiently large amplitude within the necessary time constraints (i.e. more than 0.2 seconds before the onset of stall).
Thus it can be seen that the present invention provides a diagnosis of a possible failure in the engine's operation in a very short time. High-speed performance stems from the locality properties of the wavelets and the use of high-speed processing computers. The use of high rank correlation matrix measures together with the dispersions advantageously leads to a practically error-free diagnosis. As a result, the system and method of the present invention can be used for diagnosing the operating regimes of any regularly (in particular, periodically) operating systems (aircraft compressors, turbines, auto-motors, electrical motors, pumps, power plant turbines, etc.) and for preventing their failure. The invention provides a significant improvement in the diagnostics of the operating regimes of existing engines, which is important for preventing their failure and, consequently, for lowering associated economic losses.
Moreover, in comparison with the prior art, this system and method make it possible to regulate engine operation quickly and in a reliable way by taking into account many parameters without necessitating an abrupt change of the operation regime. This is due to the effectiveness of diagnostics based on wavelet analysis in detecting precursors of engine malfunctions/instability and by rapid feedback. Consequently, the present invention can be used for the automatic regulation of the operating regimes (to prevent malfunctions and possible failure) of any electric machine including engines, motors, turbines, compressors, pumps, auto-motors, other electromagnetic rotating devices, and in general any device which operates regularly (in particular periodically). Furthermore, the system and method of the present invention do not require any human intervention at any stage.

While preferred embodiments of the present invention have been described, the embodiments disclosed are illustrative and not restrictive, and the invention is intended to be defined by the appended claims.

Claims (21)

1. A method of diagnosing a device to determine the presence of an instability in a physical characteristic of the device, the method comprising the steps of:
(a) receiving a digital input signal representative of the physical characteristic;
(b) performing a wavelet transform on the digital input signal to provide a set of wavelet coefficients;
(c) determining at least a first probabilistic measure of at least a portion of the set of wavelet coefficients;
(d) analyzing at least the first probabilistic measure to identify a precursor associated with the instability.

2. A method of diagnosing an input signal according to claim 1 wherein steps (a) to (d) are performed in real time and wherein the identification of the precursor occurs prior to the development of a malfunction in the device.

3. A method of diagnosing an input signal according to claim 1 wherein the precursor is characterized by a significant change in the probabilistic measure.

4. A method of diagnosing an input signal according to claim 1 wherein the first probabilistic measure represents the dispersion of the wavelet coefficients at a particular scale.

5. A method of diagnosing an input signal according to claim 1 wherein step (c) further comprises determining at least a second probabilistic measure of at least a portion of the set of wavelet coefficients, and step (d) further comprises the step of analyzing at least the first and second probabilistic measures to identify a precursor associated with the instability.

6. A method of diagnosing an input signal according to claim 5 wherein the precursor is characterized by a significant change in both the first probabilistic measure and the second probabilistic measure.

7. A method of diagnosing an input signal according to claim 5 wherein determining at least a second probabilistic measure includes the step of determining a high rank correlation matrix for the wavelet coefficients.

8. A method of diagnosing an input signal according to claim 1 wherein step (b) comprises performing a discrete wavelet transform on the digital input signal.

9. A method of diagnosing an input signal according to claim 1 wherein step (b) comprises performing a discretized continuous wavelet transform on the digital input signal.

10. A method of diagnosing an input signal according to claim 1 further comprising the steps of measuring the physical characteristic of the device with a sensor to provide an analog input signal and converting the analog input signal into the digital input signal.

11. A method of controlling a device driven by a command execution unit comprising the steps of:
(i) diagnosing an input signal according to claim 2; and (ii) if a precursor associated with the instability is identified in step (i), directing the command execution unit, in accordance with an instability control approach, to control the device to avert the development of the malfunction in the device.

12. A method of controlling a device according to claim 11 wherein step (i) further comprises diagnosing an input signal according to claim 5.

13. A method of controlling a device as claimed in claim 12 wherein the device forms part of a multi-functional system and wherein identifying a precursor associated with the instability further comprises analyzing additional data describing the operation of the multi-functional system.

14. A system for diagnosing a device to determine the presence of an instability in a physical characteristic of the device, the system for diagnosing comprising a computer system for receiving a digital input signal representative of the physical characteristic, the computer system comprising (a) a wavelet coefficient generation module for performing a wavelet transform on the digital input signal to provide a set of wavelet coefficients, (b) a first measure algorithm module for receiving the wavelet coefficients and determining a first probabilistic measure of at least a portion of the set of wavelet coefficients, and (c) an analysis module for receiving and analyzing at least the first probabilistic measure to identify a precursor associated with the instability prior to the development of a malfunction in the device.

15. A system for diagnosing an input signal according to claim 14 wherein the computer system further comprises a second measure algorithm module for receiving the wavelet coefficients and determining a second probabilistic measure of at least a portion of the set of wavelet coefficients, and wherein the analysis module further receives and analyzes the second probabilistic measure to identify the precursor.

16. A system for diagnosing an input signal according to claim 15 wherein the first measure algorithm module determines the dispersion of the wavelet coefficients at a particular scale, and the second measure algorithm module determines a high rank correlation matrix of the wavelet coefficients.

17. A system for controlling a device driven by a command execution unit comprising a system for diagnosing an input signal according to claim 14, wherein the computer system is connected to the command execution unit and the computer system further comprises an automatic control module for directing the command execution unit, in accordance with an instability control approach, to control the device to avert the development of the malfunction.

18. A system for controlling a device according to claim 17 wherein the system for diagnosing an input signal further comprises a system for diagnosing according to claim 15.

19. A system for controlling a device according to claim 17 wherein the device forms part of a multi-functional system including a data unit containing supplemental information about the multi-functional system, the data unit being connected to the computer system for providing the supplemental information to the automatic control module.

20. A system for controlling a device according to claim 17 wherein the device is one of the following: an engine, a motor, a turbine, a compressor, or a pump.

21. A system for controlling a device according to claim 17 wherein the physical characteristic of the device is one of the following:
vibrations, deformations, pressure, acoustic noise, temperature, or power consumed.