METHOD FOR FILTERING A SIGNAL

Abstract

A method for filtering at least one signal to be filtered. The filtering being determined based on measurements representing the movement of a rotating shaft. The method includes defining at least one characterisation indicator of the rotating shaft. The characterisation indicator being obtained from a frequency analysis of the behaviour of the rotating shaft (step 120). The method includes determining a filtering template as a function of the characterisation indicator (step 130) and applying the filtering template to said signal to be filtered (step 140).

Description

CROSS-REFERENCE TO RELATED APPLICATION

This application claims priority to French Application No. 2310163, filed Sep. 25, 2023, the entirety of which is hereby incorporated by reference.

FIELD

The present disclosure relates to the characterisation of a shaft rotating about an axis of rotation.

More specifically, the present disclosure relates to filtering a characterisation indicator of the gyroscopic effects and deformations of a rotating shaft in order to facilitate the control, modelling, monitoring, diagnosis and maintenance thereof.

BACKGROUND

A perfectly balanced rotating system, such as a rotating shaft, rotates about an axis of rotation such that said axis of rotation is immobile relative to its theoretical position, and therefore relative to the environment close to the rotating shaft.

However, faults, deformations, gyroscopic effects or even effects due to the force of gravity can cause the axis of rotation to shift relative to its theoretical position during the rotation of said rotating shaft. This is the case, for example, for deformations of the rotating shaft due to flexible Eigen-modes, with these modes deforming the shaft as a function of the speed of rotation of said shaft.

When such defects are present, the axis of rotation is not immobile and vibrates. It moves away from its theoretical position along a particular trajectory. More specifically, each point on the axis of rotation can move away from its theoretical position, and two points on the axis of rotation at two different ends of the shaft can have a different trajectory.

The properties of a rotating shaft are conventionally understood through a lateral analysis of parameters. The parameters are values of radial quantities such as the position or the radial speed of the shaft, having a non-zero component in a plane orthogonal to the axis of rotation.

The changes in these quantities often need to be characterised depending on whether or not the shaft is rotating. However, many causes can influence the rotational behaviour of a rotating shaft and it is difficult to precisely characterise such behaviour using reliable indicators.

SUMMARY

Therefore, the aim of the present disclosure is to overcome the aforementioned disadvantages and to propose filtering of a signal to be filtered characterising a rotating shaft.

Therefore, the aim of the present disclosure is a method for filtering at least one signal to be filtered, with the filtering being determined based on measurements representing the movement of a rotating shaft, the method comprising the following steps:

• • defining at least one characterisation indicator of the rotating shaft, with the characterisation indicator being obtained based on a frequency analysis of the behaviour of the rotating shaft;
  • determining a filtering template as a function of the characterisation indicator; and
  • applying the filtering template to said signal to be filtered that is determined based on measurements representing the movement of the rotating shaft.

Thus, this method allows the characterisation indicator to be rendered more reliable and allows some information to be isolated that facilitates the understanding and the rectification of faults on a rotating shaft.

Advantageously, the characterisation indicator is obtained from the following steps:

• • a) acquiring at least two input signals over a predetermined period, the input signals comprising measurements representing the movement of the rotating shaft, the measurements, associated with a vector quantity, being non-collinear to each other and optionally coplanar with a plane intersecting the axis of rotation of the shaft;
  • b) determining two components of the input signals, the components being orthogonal to each other and optionally to the axis of rotation;
  • c) merging the two components into a bivariate signal;
  • d) time-frequency decomposing the bivariate signal into spectral elements;
  • e) determining an equivalent ellipse for each spectral element; and
  • f) extracting the characterisation indicator of the rotating shaft based on parameters of the equivalent ellipse.

The gyroscopic effects and the deformations of the rotor, notably the flexible modes, are simpler to consider subject to a frequency decomposition of the input signals. Indeed, the rotating effects of these deformations are asynchronous, in other words, their frequencies are not directly linked to the rotation frequency of the shaft, and are therefore more easily characterised using a frequency decomposition that can highlight a rotating orbit, also called an equivalent ellipse, for each frequency, with each orbit rotating at a specific speed corresponding to the considered frequency.

This characterisation of a rotating shaft allows applications of the control, identification, modelling, monitoring or diagnostic type.

In one embodiment, the filtering template is applied to the time-frequency decomposition of the bivariate signal so as to obtain a filtered bivariate signal.

In various embodiments, the time-frequency decomposition of the bivariate signal on which the filtering template is applied is a time-frequency decomposition that is obtained using a Fourier, or wavelet, transform, or by filtering, or using a quaternionic Fourier transform.

Advantageously, the method comprises a step of obtaining filtered input signals by projecting the filtered bivariate signal onto the axes collinear with the measured quantities.

In a particular embodiment, the method comprises a step of computing a dual signal of the filtered bivariate signal.

In one embodiment, the at least one indicator includes a polarisation rate, and/or an indicator of circularity and/or amplitude and/or direction and/or orientation and/or synchronisation of an equivalent ellipse.

Advantageously, the filtering template is defined relative to a predefined threshold of a polarisation rate, and/or a circularity indicator of an equivalent ellipse.

Advantageously, the filtering template is defined based on at least one normalised indicator ranging between −1 and 1 or between 0 and 1.

In a particular embodiment, the filtering template is a linear bivariate filter and/or is obtained by multiplying a mask that is determined as a function of the characterisation indicator and of a frequency weighting.

In one embodiment, steps a) to f) are carried out at least a second time for an additional plane intersecting the axis of rotation of the shaft and distinct from said plane, the method further comprising a step of comparing characterisation indicators of the planes.

A further aim of the present disclosure is a computer-implemented method comprising the steps of the method as defined above. The method can be implemented on a real-time target, notably in a magnetic bearing control cabinet.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a schematic representation of the various steps of a method for obtaining an indicator for the method according to the present disclosure;

FIG. 2 is a representation of an ellipse and of its Euler parameters; and

FIG. 3 is a schematic representation of various filtering steps applicable to a signal to be filtered following the steps of FIG. 1.

DETAILED DESCRIPTION

FIG. 1 schematically shows various steps of obtaining an indicator that, in a particular embodiment, can be incorporated into the method according to the present disclosure.

These steps allow the behaviour of a rotating shaft to be characterised. A rotating shaft is understood to be a shaft that is free to rotate about an axis of rotation. The rotating shaft is a shaft of a rotating machine, for example.

In order to implement the method, a step 10 of acquiring input signals is initially carried out.

At least two input signals are acquired over a predetermined period. The predetermined period is several seconds, for example.

The input signals include measurements representing the movement of the rotating shaft. For example, the input signals include measurements of the position, or speed, or acceleration, or force, or magnetic field, or even of the magnetic induction of the rotating shaft. The input signals are acquired, for example, using sensors configured to measure these quantities.

In general, a vector quantity is measured in terms of several spatially distributed components sensitive to the rotation of the shaft. The spatial distribution of these components is used to evaluate the movement of the rotor.

A particular case is that of measurements that are non-collinear to each other and are coplanar with a plane intersecting the axis of rotation of the shaft. Two non-collinear measurements can be used to characterise the behaviour of the rotating shaft in a plane intersecting the axis of rotation of the shaft. Additional measurements allow greater precision for said characterisation.

Measurements that are not coplanar with a plane intersecting the axis of rotation of the shaft nevertheless can be returned by projecting into a reference plane intersecting the axis of rotation of the shaft.

For M measurements with M≥2, M input signals u₁(t) to u_M(t) are obtained, with t denoting the time covering the acquisition of the measurements, the set of values of which is denoted by T, in other words dates, which t can assume.

An intersecting plane is understood to be a plane in which the measurements are taken that intersect the axis of rotation of the rotating shaft. The intersecting plane is, for example, orthogonal to the axis of rotation or is not orthogonal to the axis of rotation.

If the intersecting plane is not orthogonal to the axis of rotation, a step 20 is carried out of determining two components of the input signals, with the components being orthogonal to each other and optionally to the axis of rotation. In other words, this step 20 is a change of basis facilitating the subsequent computations.

In the event that the intersecting plane is orthogonal to the axis of rotation, the previous step is intrinsically carried out.

In the event that the axis of rotation is not strictly and precisely defined, with the shaft not always being perfectly cylindrical and/or balanced, and/or non-deformable and/or known, said orthogonal intersecting plane is treated as a study plane that is assumed to be orthogonal to the axis of rotation, or at least sensitive to the rotational movement of the shaft.

In practice, the step 20 of determining the two orthogonal components u_x(t) and u_y(t) is carried out by averaging the projections of the input signals u₁(t) to u_M(t) over two orthogonal axes x and y, with these two axes ideally being orthogonal to the axis of rotation. Thus,

$\begin{array}{cc}\forall t\in T,& \begin{array}{c}{u}_x\left(t\right)={\mathrm{moy}}_{i=\left[1,M\right]}\left({\mathrm{Proj}}_x\left(u_i\left(t\right)\right)\right)\\ u_y\left(t\right)={\mathrm{moy}}_{i=\left[1,M\right]}\left({\mathrm{Proj}}_y\left(u_i\left(t\right)\right)\right)\end{array}\end{array}$

• • the components u_x(t) and u_y (t) are now in the study plane, ideally orthogonal to the axis of rotation.

A step 30 is then carried out of merging the two orthogonal components u_x(t) and u_y(t) into a bivariate signal r(t) in order to combine these two components in an equivalent manner in the form of a complex signal, such that:

$\begin{array}{cc}\forall t\in T,& r\left(t\right)=u_x\left(t\right)+1i·u_y\left(t\right)\end{array}$

It also can be expressed as follows:

$r\left(t\right)=\left\{u_x\left(t\right),u_y\left(t\right)\right\}$

Or even in vector form:

$r\left(t\right)=\left[\begin{array}{c}{u}_x\left(t\right)\\ u_y\left(t\right)\end{array}\right]$

A step 40 is then carried out of frequency decomposing the bivariate signal r(t) into spectral elements characterised by their time-frequency pair (t, f), also called time-frequency cells, with this step being carried out periodically in order to obtain both a temporal and a frequency decomposition, called time-frequency decomposition s(t, f), such that:

$s\left(t,f\right)={\mathrm{spectre}}_{t,f}\left(r\left(\tau \right)\right)$

This time-frequency decomposition step 40 can be carried out in various ways. These various embodiments are described hereafter, with time-frequency decompositions s, or ŝ, or {tilde over (s)}, or being provided in continuous form, in other words, in analogue and not discrete form. A computer-implemented digital implementation then includes a time and frequency sampling step.

In a first embodiment, denoted 41 in FIG. 1, the time-frequency decomposition step 40 is carried out with a filter bank around each frequency f.

Thus, the time-frequency decomposition s(t, f) is expressed as follows:

$\begin{array}{cc}\forall f\in F_{+},\forall t\in T,& s\left(t,f\right)=r\left(t\right)*h_f\left(t\right)\end{array}$

with h_f being the impulse response of a bandpass filter centred around the frequency f, with F₊ being all the frequencies that are greater than or equal to 0.

In a second embodiment 42, the time-frequency decomposition step 40 is carried out using a Fourier transform, in particular a short-term single-sided Fourier transform.

As standard, the time-frequency decomposition is expressed as:

$$\begin{gathered} \forall f\ge 0,\begin{array}{c}{\hat{s}}_x\left(t,f\right)=\frac{1}{\mathrm{dt}}·{\int }_{t-\mathrm{dt}}^t{u}_x\left(\tau \right)·e^{-1i\text{.2}\pi .f.\tau }·d\tau \\ \begin{array}{cc} \\ & {\hat{s}}_y\left(t,f\right)=\frac{1}{\mathrm{dt}}·{\int }_{t-\mathrm{dt}}^t{u}_y\left(\tau \right)·e^{-1i\text{.2}\pi .f.\tau }·d\tau \end{array}\end{array} \end{gathered}$$

In practice, more efficient estimators are preferably used, based on a division into temporal intervals, with time weighting, averaging over several temporal intervals, and overlaps between successive intervals. The formulation is then:

$\begin{array}{cc}\forall t\in T_{\downarrow },\forall f\in F_{+},& \end{array}\begin{array}{c}{\hat{s}}_{x_k}\left(t,f\right)={\int }_{t-{\mathrm{dt}}_k^{-}}^{t+{\mathrm{dt}}_k^{+}}{u}_x\left(\tau \right)·w_k\left(\tau -t\right)·e^{-1i\text{.2}\pi .f.\tau }·d\tau \\ {\hat{s}}_{y_k}\left(t,f\right)={\int }_{t-{\mathrm{dt}}_k^{-}}^{t+{\mathrm{dt}}_k^{+}}{u}_y\left(\tau \right)·w_k\left(\tau -t\right)·e^{-1i\text{.2}\pi .f.\tau }·d\tau \end{array}$ $\forall t\in T_{\downarrow },\forall f\in F_{+},\begin{array}{c}{\hat{s}}_x\left(t,f\right)={\mathrm{moy}}_{k=\left[1,p\right]}\left(\frac{{\hat{s}}_{x_k}\left(t,f\right)}{w_k^{\left(1\right)}}\right)\\ {\hat{s}}_y\left(t,f\right)={\mathrm{moy}}_{k=\left[1,p\right]}\left(\frac{{\hat{s}}_{y_k}\left(t,f\right)}{w_k^{\left(1\right)}}\right)\end{array}$

With dt_k⁻ and dt_k⁻ for k=[1, p] defining a division into p intervals, with possible overlap, around a given date t. Conventionally, this division does not depend on t, and is identical for each cell. The temporal extent of interval k is dt_k=dt_k⁻+dt_k⁺.

The overall temporal extent of the corresponding cell (t,f) is:

${\mathrm{dt}}^{+}={\max }_{k=\left[1,p\right]}\left({\mathrm{dt}}_k^{+}\right)$ ${\mathrm{dt}}^{-}={\max }_{k=\left[1,p\right]}\left({\mathrm{dt}}_k^{-}\right)$ $\mathrm{dt}={\mathrm{dt}}^{-}+{\mathrm{dt}}^{+}$

With Et w_k(τ) for τ∈[−dt_k⁻, dt_k⁺] for k=[1, p] defining a temporal weighting window for each of these intervals, for example, of the Hanning type. Conventionally, w_k does not depend on k, which is identical for each interval.

Depending on the estimated spectral quantity, the following compensations need to be introduced:

$\forall k\in \left[1,p\right],$ $\begin{array}{c}{w}_k^{\left(1\right)}={\int }_{t-{\mathrm{dt}}_k^{-}}^{+{\mathrm{dt}}_k^{+}}{w}_k\left(\tau \right)·d\tau \\ w_k^{\left(2\right)}={\int }_{t-{\mathrm{dt}}_k^{-}}^{+{\mathrm{dt}}_k^{+}}{w}_k^2\left(\tau \right)·d\tau \end{array}\forall k\in \left[1,p\right],$ $\begin{array}{c}{w}_k^{\left(1\right)}={\int }_{-{\mathrm{dt}}_k^{-}}^{{\mathrm{dt}}_k^{+}}{w}_k\left(\tau \right)·d\tau \\ w_k^{\left(2\right)}={\int }_{-{\mathrm{dt}}_k^{-}}^{{\mathrm{dt}}_k^{+}}{w}_k^2\left(\tau \right)·d\tau \end{array}$

In a particular case, the absence of weighting corresponds to w_k=1, in which case w_k⁽¹⁾=w_k⁽²⁾=dt_k.

In addition, a sub-set of dates is introduced: T_↑={t′_i, i=[1, m′]}

In a third embodiment 43, the time-frequency decomposition step 40 is carried out using a Fourier transform, in particular a short-term two-sided Fourier transform, also called “Full spectrum” transform, in other words, for a set F of frequencies taking values from the set of real numbers.

Thus, as standard the time-frequency decomposition {tilde over (s)}(t, f) is expressed as:

$\begin{array}{cc}\forall f\in R,& \tilde{s}\left(t,f\right)=\frac{1}{\mathrm{dt}}·{\int }_{t-\mathrm{dt}}^{t}r\left(\tau \right)·e^{-1i\text{.2}\pi .f.\tau }·d\tau \end{array}$

Or with a division into temporal blocks:

$\begin{array}{cc}\forall t\in T_{\downarrow },\forall f\in F,& {\tilde{s}}_k\left(t,f\right)={\int }_{t-{\mathrm{dt}}_k^{-}}^{t+{\mathrm{dt}}_k^{+}}r\left(\tau \right)·w_k\left(\tau -t\right)·e^{-1i\text{.2}\pi .f.\tau }·d\tau \end{array}$ $\forall t\in T_{\downarrow },\forall f\in F,\tilde{s}\left(t,f\right)={\mathrm{moy}}_{k=\left[1,p\right]}\left(\frac{{\tilde{s}}_k\left(t,f\right)}{w_k^{\left(1\right)}}\right)$

In a fourth embodiment 44, the time-frequency decomposition step 40 is carried out using a Fourier transform, in particular a short-term quaternionic Fourier transform, using the unitary imaginary numbers verifying 1i²=1j²=1k²=1i.1j.1k=−1.

Thus, as standard the time-frequency decomposition s(t, f) is expressed as:

$\begin{array}{cc}\forall f\ge 0,& {\stackrel{o}{s}}_k\end{array}\left(t,f\right)=\frac{1}{\mathrm{dt}}·{\int }_{t-\mathrm{dt}}^{t}r\left(\tau \right)·e^{-1j\text{.2}\pi .f.\tau }·d\tau$

Or with a division into temporal blocks.

$\begin{array}{cc}\forall t\in T_{\downarrow },\forall f\in F_{+},& {\stackrel{o}{s}}_k\left(t,f\right)={\int }_{t-{\mathrm{dk}}_k^{-}}^{t+{\mathrm{dk}}_k^{+}}\end{array}r\left(\tau \right)·w_k\left(\tau -t\right)·e^{-1i\text{.2}\pi .f.\tau }·d\tau$ $\begin{array}{cc}\forall t\in T_{\downarrow },\forall f\in F_{+},& {\stackrel{o}{s}}_k\left(t,f\right)=\end{array}{\mathrm{moy}}_{k=\left[1,p\right]}\left(\frac{{\stackrel{o}{s}}_k\left(t,f\right)}{w_k^{\left(1\right)}}\right)$

The expressions of the Fourier transforms provided in embodiments 41 to 44 above form applicable formulae, producing efficient estimators, notably in terms of their formulation with temporal weighting, averaging over several temporal intervals, overlaps between successive intervals, a more general expression of which can be provided with:

$\left(t,f\right)=\frac{1}{{\mathrm{dt}}_k}.{\int }_{t-{\mathrm{dt}}_k^{-}}^{t+{\mathrm{dt}}_k^{+}}r\left(\tau \right)·w_k\left(\tau \right)·e^{-1\mu \text{.2}\pi .f.\tau }·d\tau$

with w_i being the temporal weighting, and u being defined such that 1μ²=−1.

Finally,

$\overbrace{s}\left(t,f\right)={\mathrm{moy}}_k\left(\left(t,f\right)\right)$

The index k applied to w_k (above) or p_k (below) means that the weighting can differ from one interval to another.

In addition, other frequency transforms can be used, such as wavelet decompositions.

A more general expression is:

$\left(t,f\right)={\int }_{t-{\mathrm{dt}}_k^{-}}^{t+{\mathrm{dt}}_k^{+}}r\left(\tau \right)·p_k\left(\tau ,f\right)·d\tau$

Even more generally, the time division and the weighting can be varied as a function of the time-frequency cell (t, f):

$\left(t,f\right)={\int }_{t-{\mathrm{dt}}_k^{-}\left(t,f\right)}^{t+{\mathrm{dt}}_k^{+}\left(t,f\right)}r\left(\tau \right)·p_k\left(t,f\right)\left(\tau \right)·d\tau$

The trajectory (orbit) is obtained by inverse transformation, of the following type:

$s\left(t,f\right)\left(\tau \right)={\int }_{f-{\mathrm{df}}^{-}\left(t,f\right)}^{f+{\mathrm{df}}^{+}\left(t,f\right)}\left(\overbrace{s}\left(t,v\right)·q\left(t,f\right)\left(\tau ,v\right)+\overbrace{s}\left(t,-v\right)·q\left(t,f\right)\left(\tau ,-v\right)\right)·\mathrm{dv}$

With the frequency extent of the cell (t, f) then being df=df⁺+df⁻, the temporal extent of the cell (t, f) being dt=dt⁺+dt⁻.

From each spectral element s(t, f), at a fixed date t and frequency f, it is possible to deduce an equivalent ellipse from the instantaneous orbit. This can be referred to as local in the temporal and frequency sense, in other words, at the date t and frequency f.

Therefore, a step 50 is carried out of determining the equivalent ellipse for each spectral element. Irrespective of the time-frequency decomposition method, the instantaneous orbit can be expressed in the form of a parametric equation such that, for each component x and y, and with s_xo, s_yo, and Ø_x, Ø_y being the instantaneous amplitudes and phases:

$s_x\left(t,f\right)=s_{x_o}\left(t,f\right)·\cos \left({\varnothing }_x\left(t,f\right)\right)$ $s_y\left(t,f\right)=s_{y_o}\left(t,f\right)·\cos \left({\varnothing }_y\left(t,f\right)\right)$

For a fixed frequency f, covering the time t allows the trajectory to be obtained at this frequency.

This instantaneous orbit assumes the form of an equivalent ellipse, provided that s_x(t, f), s_y(t, f)), {dot over (Ø)}_x(t, f), {dot over (Ø)}_y(t, f) vary slowly as a function of time, for example, during a complete revolution at the frequency f, and in any case for the duration dt of the temporal interval defining the spectral element.

In this case, the date t of the considered spectral element can be distinguished from the time τ that evolves in this spectral element. In addition, by introducing the frequency f into the expression of the phases that are assumed to be quasi-linear, the equivalent ellipse is obtained in the form of a temporal parametric equation for each component x and y:

${\stackrel{\prime }{s}}_x\left(t,f\right)\left(\tau \right)=a_x\left(t,f\right)·\cos \left(2·\pi ·f·\tau +{\phi }_x\left(t,f\right)\right)$ ${\stackrel{\prime }{s}}_y\left(t,f\right)\left(\tau \right)=a_y\left(t,f\right)·\cos \left(2·\pi ·f·\tau +{\phi }_y\left(t,f\right)\right)$

These parametric equations also can be written as follows, for example:

${\stackrel{\prime }{s}}_x\left(t,f\right)\left(\tau \right)=a_x\left(t,f\right)·\left(\cos \left(2·\pi ·f·\tau \right)·\cos \left({\phi }_x\left(t,f\right)\right)-\sin \left(2·\pi ·f·\tau \right)·\sin \left({\phi }_x\left(t,f\right)\right)\right)$ ${\stackrel{\prime }{s}}_y\left(t,f\right)\left(\tau \right)=a_y\left(t,f\right).\left(\cos \left(2·\pi ·f·\tau \right)·\cos \left({\phi }_y\left(t,f\right)\right)-\sin \left(2·\pi ·f·\tau \right)·\sin \left({\phi }_y\left(t,f\right)\right)\right)$

In order to obtain these parametric equations, notably in the case of a frequency decomposition using a filter bank, the parameters a_x(t, f), φ_x(t, f), a_y(t, f), φ_y(t, f) of the equivalent ellipse can be obtained by estimates, such as direct or quadratic averaging or by RMS value, from the instantaneous quantities (τ, f) pour τ∈[t−dt⁻, t+dt⁺].

In the case of a frequency decomposition using a Fourier transform, the parametric equation of the equivalent ellipse is obtained using an inverse Fourier transform of the frequency element:

${\stackrel{\prime }{s}}_x\left(t,f\right)\left(\tau \right)=2·\mathrm{Re}\left({\hat{s}}_x\left(t,f\right)·e^{1i\text{.2}.\pi .f.\tau }\right)$ ${\stackrel{\prime }{s}}_y\left(t,f\right)\left(\tau \right)=2·\mathrm{Re}\left({\hat{s}}_y\left(t,f\right)·e^{1i\text{.2}\pi .f.\tau }\right)$

By forming the associated complex quantity ś (t, f)=Ś_x(t, f)+1i.Ś_y(t, f), in other words, the bivariate representation, it is possible to highlight the decomposition of a parametric equation of the equivalent ellipse in two circles in opposite directions:

$\stackrel{\prime }{s}\left(t,f\right)\left(\tau \right)=\mathrm{a\_}\left(t,f\right)·e^{1i·{\theta }_{+}\left(t,f\right)}·e^{1i·2·\pi ·f·\tau }+a_{+}\left(t,f\right)·e^{-1i·\mathrm{\theta \_}\left(t,f\right)}·e^{-1i·2·\pi ·f·\tau }$

With the following parameters obtained following a computation step 59:

$a_{+}=\frac{a}{2}·\sqrt{a_x^2+a_y^2+2·a_x·a_y·\sin \left({\phi }_x-{\phi }_y\right)}$ $a_{-}=\frac{a}{2}·\sqrt{a_x^2+a_y^2+2·a_x·a_y·\sin \left({\phi }_y-{\phi }_x\right)}$ ${\theta }_{+}=\mathrm{arc}\tan 2\left(a_x·\sin \left({\phi }_x\right)+a_y·\cos \left({\phi }_y\right),a_x·\cos \left({\phi }_x\right)-a_y·\sin \left({\phi }_y\right)\right)$ ${\theta }_{-}=\mathrm{arc}\tan 2\left(a_x·\sin \left({\phi }_x\right)-a_y·\cos \left({\phi }_y\right),a_x·\cos \left({\phi }_x\right)+a_y·\sin \left({\phi }_y\right)\right)$

The characterisation of the equivalent ellipse is straightforward using the two-sided (“full spectrum”) Fourier transform, which can be directly written in polar form:

$\forall f>0,\begin{array}{c}\tilde{s}\left(t,f\right)=a_{+}\left(t,f\right)·e^{1i·{\theta }_{+}\left(t,f\right)}\\ \tilde{s}\left(t,-f\right)=a_{-}\left(t,f\right)·e^{-1i·{\theta }_{-}\left(t,f\right)}\end{array}$

In general, an equivalent ellipse can be determined and plotted by obtaining its Euler parameters.

FIG. 2 shows an example of an ellipse and its Euler parameters a(t, f), θ(t, f), χ(t,f), φ(t, f), which fully define an equivalent ellipse, with:

$a\ge 0$ $\theta \in \left[-\pi /2,\pi /2\right]$ $\chi \in \left[-\pi /4,\pi /4\right]$ $\phi \in \left[-\pi ,\pi \right]$

Furthermore, b is the semi-major axis, in other words, the major radius, and c is the semi-minor axis, in other words, the minor radius, the sign of which provides the direction of rotation of the equivalent ellipse.

The parametric equation of the trajectory of the equivalent ellipse is more practical in complex form, since:

• • an equation of a circle is notably expressed as cos and sin of the same phase;
  • the axes are dilated in order to obtain an ellipse;
  • a rotation is applied thereto in order to obtain the desired direction, thus:

$\stackrel{\prime }{s}\left(t,f\right)\left(\tau \right)=a\left(t,f\right)·e^{1i·\theta \left(t,f\right)}·\left(\cos \left(\chi \left(t,f\right)\right)·\cos \left(2·\pi ·f·\tau +\phi \left(t,f\right)\right)+1i·\sin \left(\chi \left(t,f\right)\right)·\sin \left(2·\pi ·f·\tau +\phi \left(t,f\right)\right)\right)$

Direction is understood to mean the angle θ of the equivalent ellipse between the major axis and the x-axis. Orientation is understood to mean the direction of rotation describing the trajectory of said ellipse.

This characterisation of the equivalent ellipse is natural with the time-frequency decomposition using a quaternionic Fourier transform, which can be directly placed into Euler form:

$\stackrel{o}{s}\left(t,f\right)=❘\stackrel{o}{s}\left(t,f\right)❘·e^{1i·\theta \left(t,f\right)}·e^{-1k·\chi \left(t,f\right)}·e^{1j·\phi \left(t,f\right)}$

This parametric equation has the particular advantage of rejecting the temporal parametrisation towards a single exponential term e^1j.2.π.f.τ when the inverse Fourier transform is applied to the frequency element in order to obtain the temporal parametric equation of the trajectory of the equivalent ellipse:

$\stackrel{\prime }{s}\left(t,f\right)\left(\tau \right)\approx \left|\stackrel{o}{s}\left(t,f\right)\right|·p_{{\mathrm{roj}}_{C_{1i}}}\left(e^{1i·\theta \left(t,f\right)}·e^{-1k·\chi \left(t,f\right)}·e^{1j·\phi \left(t,f\right)}·e^{1j·2·\pi ·f·\tau }\right)$

with Proj_C1i being a projection into the set of complex numbers C_1i.

More specifically, the step 50 of determining the equivalent ellipse includes a step of computing the Euler parameters of said equivalent ellipse.

Starting from the time-frequency decomposition using a filter bank, a step 51 of computing a_x and a_y is obtained such that:

$\forall \tau \in T_{\downarrow },\forall f\in F_{+}^{*},\begin{array}{c}{a}_x\left(t,f\right)=\sqrt{\frac{2}{dt}·{\int }_{t-{\mathrm{dt}}^{-}}^{t+{\mathrm{dt}}^{+}}{s}_x^2\left(\tau ,f\right)·d\tau }\\ a_y\left(t,f\right)=\sqrt{\frac{2}{dt}·{\int }_{t-{\mathrm{dt}}^{-}}^{t+{\mathrm{dt}}^{+}}{s}_y^2\left(\tau ,f\right)·d\tau }\end{array}$

with T_↑ being a sub-set of T.

Starting from the time-frequency decomposition using a short-term single-sided Fourier transform, a step 52 of computing a_x, a_y, φ_x and φ_y is obtained such that:

$\forall t\in T_{\downarrow },\forall f\in F_{+}^{*},\begin{array}{c}{a}_x\left(t,f\right)=2·❘{\hat{s}}_x\left(t,f\right)❘{\phi }_x\left(t,f\right)=\mathrm{Arg}\left({\hat{s}}_x\left(t,f\right)\right)\\ a_y\left(t,f\right)=2·❘{\hat{s}}_y\left(t,f\right)❘{\phi }_y\left(t,f\right)=\mathrm{Arg}\left({\hat{s}}_y\left(t,f\right)\right)\end{array}$

Starting from the time-frequency decomposition using a short-term two-sided Fourier transform, a step 53 of computing a₊, a₋, θ₊ and θ₋ is obtained such that:

$\forall t\in T_{\downarrow },\forall f\in F_{+}^{*},\begin{array}{c}{a}_{+}\left(t,f\right)=❘\tilde{s}\left(t,f\right)❘{\theta }_{+}\left(t,f\right)=\mathrm{Arg}\left(\tilde{s}\left(t,f\right)\right)\\ a_{-}\left(t,f\right)=❘\tilde{s}\left(t,-f\right)❘{\theta }_{-}\left(t,f\right)=-\mathrm{Arg}\left(\tilde{s}\left(t,-f\right)\right)\end{array}$

An additional computation step 54 or 55 yields the following Euler parameters:

$a=\sqrt{2}·\sqrt{a_{+}^2+a_{-}^2}a=\sqrt{a_x^2+a_y^2}$ $\phi =\frac{{\theta }_{+}+{\theta }_{-}}{2}\mathrm{ou}\phi =\frac{{\phi }_x+{\phi }_y}{2}$ $\theta =\frac{{\theta }_{+}+{\theta }_{-}}{2}\theta =\frac{1}{2}·\arctan 2\left(2·a_x·a_y·\cos \left({\phi }_x-{\phi }_y\right),a_x^2-a_y^2\right)$ $\chi =\arctan 2\left(a_{+}-a_{-},a_{+}+a_{-}\right)\chi =\frac{1}{2}·\arcsin \left(\mathrm{ratio}*\left(2·a_x·a_y,a_x^2+a_y^2\right)·\sin \left({\phi }_x-{\phi }_y\right)\right)$

An inverse computation step 56 is also possible, such that:

$a_{+}=\frac{a}{2}·\left(\cos \left(\chi \right)+\sin \left(\chi \right)\right)$ $a_{-}=\frac{a}{2}·\left(\cos \left(\chi \right)-\sin \left(\chi \right)\right)$ ${\theta }_{+}=\phi +\theta$ ${\theta }_{-}=\phi -\theta$

The method optionally comprises a step 60 of computing the spectral density of the time-frequency decomposition. The spectral densities of each component and the cross-spectral densities between components are used to estimate the energy distribution between the x and y components.

Starting from the time-frequency decomposition using a filter bank, a step 61 of computing the spectral density S is obtained such that:

$\forall t\in T_{\downarrow },\forall f\in F_{+},\begin{array}{c}{S}_x\left(t,f\right)=\frac{1}{dt}·{\int }_{t-{\mathrm{dt}}^{-}}^{t+{\mathrm{dt}}^{+}}{s}_x^2\left(\tau ,f\right)·d\tau \\ S_y\left(t,f\right)=\frac{1}{dt}·{\int }_{t-{\mathrm{dt}}^{-}}^{t+{\mathrm{dt}}^{+}}{s}_y^2\left(\tau ,f\right)·d\tau \\ S_{xy}\left(t,f\right)=\frac{1}{dt}·{\int }_{t-{\mathrm{dt}}^{-}}^{t+{\mathrm{dt}}^{+}}{s}_x\left(\tau ,f\right)·s_y\left(\tau ,f\right)·d\tau \end{array}$

Starting from the time-frequency decomposition using a short-term single-sided Fourier transform, a step 62 of computing the spectral density Ŝ is obtained such that:

$\forall t\in T_{\downarrow },\forall f\in F_{+}^{*},\begin{array}{c}{\hat{S}}_x\left(t,f\right)=2·\mathrm{df}·{\mathrm{moy}}_{k=\left[1,p\right]}\left(\frac{{❘{\hat{s}}_{x_k}\left(t,f\right)❘}^2}{w_k^{\left(2\right)}}\right)\\ {\hat{S}}_y\left(t,f\right)=2·\mathrm{df}·{\mathrm{moy}}_{k=\left[1,p\right]}\left(\frac{{❘{\hat{s}}_{y_k}\left(t,f\right)❘}^2}{w_k^{\left(2\right)}}\right)\\ {\hat{S}}_{xy}\left(t,f\right)=2·\mathrm{df}·{\mathrm{moy}}_{k=\left[1,p\right]}\left(\frac{{\hat{s}}_{x_k}\left(t,f\right)·\stackrel{\_}{{\hat{s}}_{y_k}(t,f})}{w_k^{\left(2\right)}}\right)\end{array}$

Starting from the time-frequency decomposition using a short-term two-sided Fourier transform, a step 63 of computing the spectral density {tilde over (S)} is obtained such that:

$\forall \tau \in T_{\downarrow },\forall f\in F^{*},\tilde{S}\left(t,f\right)=\mathrm{df}·{\mathrm{moy}}_{k=\left[1,p\right]}\left(\frac{{❘{\tilde{s}}_k\left(t,f\right)❘}^2}{w_k^{\left(2\right)}}\right)$

The following also will be used:

$\forall t\in T_{\downarrow },\forall f\in F_{+}^{*},\widetilde{S_{*}}\left(t,f\right)=\mathrm{df}·{\mathrm{moy}}_{k=\left[1,p\right]}\left(\frac{{\tilde{s}}_k\left(t,f\right)·{\tilde{s}}_k\left(t,-f\right)}{w_k^{\left(2\right)}}\right)$

Starting from the time-frequency decomposition using a short-term quaternionic Fourier transform, a step 64 of computing the spectral density is obtained such that:

$\forall t\in T_{\downarrow },\forall f\in F_{+},\stackrel{o}{S}\left(t,f\right)=\mathrm{df}·{\mathrm{moy}}_{k=\left[1,p\right]}\left(\frac{{❘\stackrel{o}{s_k}\left(t,f\right)❘}^2}{w_k^{\left(2\right)}}\right)+\mathrm{df}·{\mathrm{moy}}_{k=\left[1,p\right]}\left(\frac{\stackrel{o}{s_k}\left(t,f\right)·1j·\stackrel{\_}{\stackrel{o}{s_k}\left(t,f\right)}}{w_k^{\left(2\right)}}\right)$

Based on these spectral densities S, Ŝ, {tilde over (S)}, and , a step 70 is carried out of computing the Stokes parameters S₀, S₁, S₂, and S₃ of at least one equivalent ellipse.

The Stokes parameters are four real numbers and are conventionally used to describe the polarisation state of an electromagnetic wave.

S₀ is commonly the total intensity of an electromagnetic wave and in this case by extension represents the total intensity of the bivariate signal, S₀>0.

S₁ and S₂ provide the intensity of rectilinear polarisation √{square root over (S₁²+S₂²)}, as well as the direction of the polarisation of the ellipse

$\theta =\frac{1}{2}.$

a tan 2(S₂, S₁).

S₃ provides the circular polarisation intensity |S₃|, with its sign providing the direction of rotation.

The reduced Stokes parameters are also defined:

$\mathrm{si}{S}_0\ne 0,\forall i\in \left\{1,2,3\right\}{s}_i\stackrel{\mathrm{def}}{=}\frac{S_i}{S_0}$

Starting from the spectral density obtained for the time-frequency decomposition using a filter bank, a step 71 of computing the Stokes parameters of at least one equivalent ellipse is obtained such that:

$\forall t\in T_{\downarrow },\forall f\in F_{+}^{*},\begin{array}{c}{S}_0\left(t,f\right)=2·\left(S_x\left(t,f\right)+S_y\left(t,f\right)\right)\\ S_1\left(t,f\right)=2·\left(S_x\left(t,f\right)-S_y\left(t,f\right)\right)\end{array}$

Starting from the spectral density obtained for the time-frequency decomposition using a short-term single-sided Fourier transform, a step 72 of computing the Stokes parameters of at least one equivalent ellipse is obtained such that:

$\forall t\in T_{\downarrow },\forall f\in F_{+}^{*},\begin{array}{c}\begin{array}{c}\begin{array}{c}{S}_0\left(t,f\right)=2·\left({\hat{S}}_x\left(t,f\right)+{\hat{S}}_y\left(t,f\right)\right)\\ S_1\left(t,f\right)=2·\left({\hat{S}}_x\left(t,f\right)-{\hat{S}}_y\left(t,f\right)\right)\end{array}\\ S_2\left(t,f\right)=4·\mathrm{Re}\left({\hat{S}}_{xy}\left(t,f\right)\right)\end{array}\\ S_3\left(t,f\right)=4·\mathrm{Im}\left({\hat{S}}_{xy}\left(t,f\right)\right)\end{array}$

Starting from the spectral density obtained for the time-frequency decomposition using a short-term two-sided Fourier transform, a step 73 of computing the Stokes parameters of at least one equivalent ellipse is obtained such that:

$\forall t\in T_{\downarrow },\forall f\in F_{+}^{*},\begin{array}{c}\begin{array}{c}\begin{array}{c}{S}_0\left(t,f\right)=2·\left(\tilde{S}\left(t,f\right)+\tilde{S}\left(t,-f\right)\right)\\ S_1\left(t,f\right)=4·\mathrm{Re}\left(\widetilde{S_{*}}\left(t,f\right)\right)\end{array}\\ S_2\left(t,f\right)=4·\mathrm{lm}\left(\widetilde{S_{*}}\left(t,f\right)\right)\end{array}\\ S_3\left(t,f\right)=2·\left(\tilde{S}\left(t,f\right)-\tilde{S}\left(t,-f\right)\right)\end{array}$

Starting from the spectral density obtained for the time-frequency decomposition using a short-term quaternionic Fourier transform, a step 74 of computing the Stokes parameters of at least one equivalent ellipse is obtained such that:

$t\in T_{\downarrow },\forall f\in F_{+},\begin{array}{c}\begin{array}{c}\begin{array}{c}{S}_0\left(t,f\right)=\mathrm{Re}\left(\stackrel{o}{S}\left(t,f\right)\right)\\ S_1\left(t,f\right)={\mathrm{Im}}_{1j}\left(\stackrel{o}{S}\left(t,f\right)\right)\end{array}\\ S_2\left(t,f\right)={\mathrm{Im}}_{1k}\left(\stackrel{o}{S}\left(t,f\right)\right)\end{array}\\ S_3\left(t,f\right)={\mathrm{Im}}_{1i}\left(\stackrel{o}{S}\left(t,f\right)\right)\end{array}$

As an alternative embodiment, the Stokes parameters S₁, S₂, and S₃ are directly computed using the time-frequency decomposition using a quaternionic Fourier transform:

$\forall t\in T_{\downarrow },\forall f\in F_{+},\begin{array}{c}\begin{array}{c}\begin{array}{c}\stackrel{o}{\sigma }\left(t,f\right)=\stackrel{o}{s}\left(t,f\right)·1j·\stackrel{\_}{\stackrel{o}{s}\left(t,f\right)}\\ S_1\left(t,f\right)={\mathrm{Im}}_{1j}\left(\stackrel{o}{\sigma }\left(t,f\right)\right)\end{array}\\ S_2\left(t,f\right)={\mathrm{Im}}_{1k}\left(\stackrel{o}{\sigma }\left(t,f\right)\right)\end{array}\\ S_3\left(t,f\right)={\mathrm{Im}}_{1i}\left(\stackrel{o}{\sigma }\left(t,f\right)\right)\end{array}$

Based on these Stokes parameters, a step 75 is carried out of computing the polarisation rate Ø, as well as the aforementioned reduced Stokes parameters, as well as the polarised part of the intensity S_p.

$\forall t\in T_{\downarrow },\forall f\in F_{+},\begin{array}{c}\begin{array}{c}\forall i\in \left\{1,2,3\right\}{s}_i\left(t,f\right)={\mathrm{ratio}}^{*}\left(S_i\left(t,f\right),S_0\left(t,f\right)\right)\\ \varnothing \left(t,f\right)=\sqrt{s_1^2\left(t,f\right)+s_2^2\left(t,f\right)+s_3^2\left(t,f\right)}\end{array}\\ S_p\left(t,f\right)=\varnothing \left(t,f\right)·S_0\left(t,f\right)\end{array}$

The polarisation expressed by the polarisation rate Ø highlights the cyclo-stationarity of an equivalent ellipse, namely, the regularity of the trajectory of said ellipse relative to a period of revolution. The polarisation reflects the structuring of said trajectory relative to defined frequencies. The polarisation rate thus can be computed for each time-frequency cell resulting from the time-frequency decomposition in order to highlight the cyclo-stationarity of the rotating shaft for a given frequency.

In addition, a computation step 57 can be implemented in order to determine the Euler parameters a, θ, and χ based on the polarisation spectral density S_p and the reduced Stokes parameters:

$\forall t\in T_{\downarrow },\forall f\in F_{+}^{*},\begin{array}{c}\begin{array}{c}a\left(t,f\right)=\sqrt{S_p\left(t,f\right)}\\ \theta \left(t,f\right)=\frac{1}{2}·\mathrm{arc}\tan 2\left(s_2\left(t,f\right),s_1\left(t,f\right)\right)\end{array}\\ \chi \left(t,f\right)=\frac{1}{2}\left({\mathrm{ratio}}^{*}\left(s_3\left(t,f\right),S_p\left(t,f\right)\right)\right)\end{array}$

A fourth Euler parameter φ providing the phase that is the source of the trajectory also can be deduced from the following computation step 58:

$\forall t\in T_{\downarrow },\forall f\in F_{+}^{*},\phi \left(t,f\right)=\mathrm{Arg}\left(\mathrm{Pro}{j}_{C_{1j}}\left(e^{1k·\chi \left(t,f\right)}·e^{-1i·\theta \left(t,f\right)}·\stackrel{o}{s}\left(t,f\right)\right)\right)$

Following the step 50 of determining an equivalent ellipse, notably by computing the Euler parameters of said equivalent ellipse, a step 80 of extracting at least one characterisation indicator of the rotating shaft from parameters of the equivalent ellipse, for example, Euler or Stokes parameters or parameters derived therefrom, is carried out.

In addition, the preceding steps are optionally carried out at least a second time for an additional plane intersecting the axis of rotation of the shaft and distinct from the plane for which these steps were initially carried out. In this embodiment, it is thus possible to carry out a step of comparing the characterisation indicators of the planes.

For example, the at least one characterisation indicator includes the previously described polarisation rate Ø, and/or a circularity indicator (linked to |χ|), and/or an intensity indicator and/or an amplitude indicator, and/or a direction of rotation indicator (linked to the sign of χ), and/or a direction indicator (linked to θ), and/or synchronisation indicator (linked to φ) of the equivalent ellipse.

Each elementary indicator helps to provide a practical and useful characteristic of the rotating system, for example, a high rate of polarisation and a high degree of circularity are probably due to the rotating shaft. Distinguishing this part of the signal, with the complementary part being treated as disturbances, can be useful, for example, for identifying said shaft (for modelling purposes) or for better controlling it in the case of a shaft equipped with magnetic bearings.

Optionally, the at least one characterisation indicator is normalised, in other words, its value ranges between −1 and 1 or between 0 and 1.

The characterisation indicator is provided for one plane, it is then denoted A, or by a comparison between two planes, it is then denoted I.

The characterisation indicator is referred to as a criterion, denoted C, when it is constructed from any characteristic quantity re-evaluated at each instant t. The characterisation indicator is referred to as stationarity, denoted S, when it is based on the stationarity of a quantity over time (maximum when the variation of said variable over time is zero).

Thus, the polarisation rate indicator is denoted:

$\forall t\in T_{\downarrow },\forall f\in F_{+},A{C}_{pol\mathrm{arisation}}\left(t,f\right)=\Phi \left(t,f\right)$

For one plane, the intensity indicator is:

$\forall t\in T_{\downarrow },\forall f\in F_{+},A{C}_{\mathrm{intensité}}\left(t,f\right)={\mathrm{normalise}}_{\tau \in T_{\downarrow },\nu \in F_{+}}\left(S_0\left(t,f\right)\right)$

With the normalisation function:

$\forall t\in TT,\forall f\in \mathrm{FF},{\mathrm{normalise}}_{\tau \in \mathrm{TT},\nu \in \mathrm{FF}}\left(x\left(t,f\right)\right)=\frac{x\left(t,f\right)}{{\max }_{\tau \in \mathrm{TT},\nu \in \mathrm{FF}}^{*}\left(❘x\left(\tau ,\nu \right)❘\right)}$

TT and FF denote sets of generic dates and frequencies, taking their values from the set of real values R.

For one plane, the amplitude indicator is:

$\forall t\in T_{\downarrow },\forall f\in F_{+},{\mathrm{AC}}_{\mathrm{amplitude}}\left(t,f\right)={\mathrm{normalise}}_{\tau \in T_{\downarrow },v\in F_{+}}\left(a\left(t,f\right)\right)$

For one plane, the sign indicator is −1 or 1 depending on the sign of x:

$\forall t\in T_{\downarrow },\forall f\in F_{+},{\mathrm{AC}}_{\mathrm{signe}}\left(t,f\right)=\mathrm{signe}\left(\chi \left(t,f\right)\right)$

For one plane, the circularity indicator is:

$\forall t\in T_{\downarrow },\forall f\in F_{+},{\mathrm{AC}}_{\mathrm{circularit}\stackrel{\prime }{e}}\left(t,f\right)=\frac{❘\chi \left(t,f\right)❘}{\pi /4}$

For one plane, the sign stationarity indicator, in other words, the direction of rotation, in other words, the orientation, is:

$\forall t\in T_{\downarrow },\forall f\in F_{+},{\mathrm{AS}}_{\mathrm{signe}}\left(t,f\right)=\mathrm{autostat}\left({\mathrm{AC}}_{\mathrm{signe}}\left(t,f\right),\mathrm{tol}=0\right)$

With the stationarity function:

$\forall t\in \mathrm{TT},\begin{array}{c}\mathrm{pour}\mathrm{tol}>0:\mathrm{autostat}\left(x\left(t\right),\mathrm{tol}\right)=e^{-\frac{❘\mathrm{dx}\left(t\right)❘}{\mathrm{tol}}}\\ \mathrm{pour}\mathrm{tol}=0:\mathrm{autostat}\left(x\left(t\right),0\right)=\left(\mathrm{dx}\left(t\right)==0\right)\end{array}$

For one plane, the circularity stationarity indicator is:

$\forall t\in T_{\downarrow },\forall f\in F_{+},{\mathrm{AS}}_{\mathrm{circularit}\stackrel{\prime }{e}}\left(t,f\right)=\mathrm{autostat}\left({\mathrm{AC}}_{\mathrm{circularit}\stackrel{\prime }{e}}\left(t,f\right),\mathrm{tol}=0.1\right)$

For one plane, the direction stationarity indicator, in other words, the geometric inclination of the ellipse, is:

$\forall t\in T_{\downarrow },\forall f\in F_{+},{\mathrm{AS}}_{\mathrm{direction}}\left(t,f\right)=\mathrm{autostat}\left(\theta \left(t,f\right),\mathrm{tol}=10·\pi /180\right)$

For one plane, the synchronisation stationarity indicator is:

$\forall t\in T_{\downarrow },\forall f\in F_{+},{\mathrm{AS}}_{\mathrm{synchro}}\left(t,f\right)=\mathrm{autostat}\left(\phi \left(t,f\right),\mathrm{tol}=10·\pi /180\right)$

For two planes, each parameter or indicator is redefined for an additional plane, by denoting these new parameters with an apostrophe, for example, a′.

The amplitude comparative characterisation indicator is:

$\forall t\in T_{\downarrow },\forall f\in F_{+},{\mathrm{IC}}_{\mathrm{amplitude}}\left(t,f\right)={\mathrm{normalise}}_{t\in T_{\downarrow },v\in F_{+}}\left(\sqrt{a\left(t,f\right)·a^{\prime }\left(t,f\right)}\right)$

The sign comparative characterisation indicator is a Boolean between 0 and 1 depending on the result of the following test:

$\forall t\in T_{\downarrow },\forall f\in F_{+},{\mathrm{IC}}_{\mathrm{signe}}\left(t,f\right)=\left({\mathrm{AC}}_{\mathrm{signe}}\left(t,f\right)=={\mathrm{AC}}_{\mathrm{signe}}^{\prime }\left(t,f\right)\right)$

The circularity comparative characterisation indicator is:

$\forall t\in T_{\downarrow },\forall f\in F_{+},{\mathrm{IC}}_{\mathrm{circularit}\stackrel{\prime }{e}}\left(t,f\right)=1-\mathrm{reldif}\left({\mathrm{AC}}_{\mathrm{circularit}\stackrel{\prime }{e}}\left(t,f\right),{\mathrm{AC}}_{\mathrm{circularit}\stackrel{\prime }{e}}^{\prime }\left(t,f\right)\right)$

With the relative error function such that:

$\mathrm{reldif}\left(x,y\right)={\mathrm{ratio}}^{*}\left(2·❘x-y❘,❘x❘+❘y❘\right)$

And with the division function with exception such that:

${\mathrm{ratio}}^{*}\left(x,y\right)=\begin{array}{cc}x/y& \mathrm{si}y\ne 0\\ 0& \mathrm{sinon}\end{array}$

The direction comparative characterisation indicator, in other words, the geometric inclination of the ellipse, is:

$\forall t\in T_{\downarrow },\forall f\in F_{+},{\mathrm{IC}}_{\mathrm{direction}}\left(t,f\right)=1-\mathrm{reldif}\left(\theta \left(t,f\right),{\theta }^{\prime }\left(t,f\right)\right)$

The synchronisation comparative characterisation indicator is:

$\forall t\in T_{\downarrow },\forall f\in F_{+},{\mathrm{IC}}_{\mathrm{synchro}}\left(t,f\right)=1-\mathrm{reldif}\left(\phi \left(t,f\right),{\phi }^{\prime }\left(t,f\right)\right)$

The stationarity and comparative circularity characterisation indicator is:

$\forall t\in T_{\downarrow }\forall f\in F_{+},{\mathrm{IS}}_{\mathrm{circularit}\stackrel{\text{'}}{e}}\left(t,f\right)=\mathrm{autostat}\left({\mathrm{AC}}_{\mathrm{circularit}\stackrel{\text{'}}{e}}\left(t,f\right)-{\mathrm{AC}}_{\mathrm{circularit}\stackrel{\text{'}}{e}}^{\prime }\left(t,f\right),\mathrm{tol}=0.1\right)$

The stationarity and comparative direction characterisation indicator is:

$\forall t\in T_{\downarrow }\forall f\in F_{+},{\mathrm{IS}}_{\mathrm{direction}}\left(t,f\right)=\mathrm{autostat}\left(\theta \left(t,f\right)-{\theta }^{\prime }\left(t,f\right),\mathrm{tol}=10.\pi /180\right)$

The stationarity and comparative synchronisation characterisation indicator is:

$\forall t\in T_{\downarrow }\forall f\in F_{+},{\mathrm{IS}}_{\mathrm{synchro}}\left(t,f\right)=\mathrm{autostat}\left(\phi \left(t,f\right)-{\phi }^{\prime }\left(t,f\right),\mathrm{tol}=10.\pi /180\right)$

Other indicators can be constructed from the quantities constructed in steps 70 and 50.

Optionally, a step 90 is carried out of introducing the at least one characterisation indicator into a neural network in order to classify said indicator as a function of the rotation of the shaft and to link it to a previously listed anomaly. Anomalies similar to the learning base are then detected.

Optionally, a step 100 is carried out of displaying the at least one extracted characterisation indicator so that it can be read by an operator and so that they can understand potential anomalies in the rotating shaft.

A step 110 of displaying an image of the trajectory of at least one equivalent ellipse also can be carried out.

FIG. 3 schematically shows the steps of the filtering method according to one example the present disclosure. In a particular embodiment, the following steps are carried out following the step 80 of extracting a characterisation indicator.

The filtering method is intended to filter at least one signal to be filtered, for example, a signal determined based on measurements representing the movement of a rotating shaft. The method may or may not be carried out in real time. Said signal can be one of the input signals determined in step 10, and/or one of the components of the input signals determined in step 20, and/or even the bivariate signal determined in step 30.

In order to carry out this filtering method, a step 120 of defining at least one characterisation indicator of the rotating shaft is initially carried out, with the characterisation indicator being obtained from a frequency analysis of the behaviour of the rotating shaft.

In a particular embodiment, the characterisation indicator is an indicator extracted during step 80 described above.

Thus, the at least one characterisation indicator includes the previously described polarisation rate Ø, and/or a circularity indicator, and/or an intensity indicator and/or an amplitude indicator, and/or a direction of rotation indicator, and/or a direction indicator, and/or a synchronisation indicator of the equivalent ellipse.

A step 130 is then carried out of determining a filtering template as a function of the at least one characterisation indicator.

The step 130 of determining a filtering template comprises, for example, a step 131 of computing a mask, such that:

$\forall t\in T_{\downarrow }\forall f\in F_{+},M\left(t,f\right)=\underset{i,j,k,l}{\mathrm{combine}}\left({\mathrm{AC}}_{{\mathrm{item}}_i}\left(t,f\right),{\mathrm{AS}}_{{\mathrm{item}}_j}\left(t,f\right),{\mathrm{IC}}_{{\mathrm{item}}_k}\left(t,f\right),{\mathrm{IS}}_{{\mathrm{item}}_l}\left(t,f\right)\right)$

Where combine is a function that as input takes one or more characterisation indicators at values within [0.1] or [−1.1] and as output renders a real value within [0.1].

Following the step 131 of computing the mask M, a step 132 is carried out of computing the filtering template P and its dual P*, such that:

$\begin{array}{cc}\forall t\in T_{\downarrow }\forall f\in F_{+},& \begin{array}{c}P\left(t,f\right)=M\left(t,f\right).N\left(t,f\right)\\ P^{*}\left(t,f\right)=M^{*}\left(t,f\right).N\left(t,f\right)\end{array}\end{array}$

With N being a frequency weighting with values within [0.1]. N is a random frequency weighting, independent of M and is similar to conventional frequency filtering. N allows, for example, a frequency band to be focused upon around the speed of the shaft.

In addition, the dual mask (complementary): M*=1−M

The advantage of considering M* and then P* is that, depending on the application, it is possible to focus on a certain portion of the signal, obtained according to the characterisation indicators used, provided by M and/or the residue equated with disturbances provided by M*. For example, a magnetic bearing controller can be contemplated that acts differently on the rotating part and on the external disturbances.

In one embodiment, the mask M, and by extension the filtering template P, is defined relative to a predefined threshold of a characterisation indicator, for example, a predefined threshold of a polarisation rate, and/or of a circularity indicator of an equivalent ellipse as described above.

The constancy of a characterisation indicator also can be used to define a filtering template relative to said characterisation indicator.

Finally, a step 140 is carried out of applying the filtering template P to the signal to be filtered, for example, a signal determined based on measurements representing the movement of the rotating shaft.

In one embodiment, the filtering template is applied to the time-frequency decomposition of the bivariate signal so as to obtain a filtered bivariate signal R and its jointly computed dual R*. Step 140 is then carried out according to various embodiments.

Starting from the time-frequency decomposition using a short-term single-sided Fourier transform, a step 141 is obtained of computing an intermediate filtered bivariate signal {grave over (R)}, its dual {grave over (R)}*, and components of the filtered input signals determined in step 20 Ù_x, Ù_x*, Ù_y*, and Ù_y.

$\begin{array}{cc}\forall t\in T_{\downarrow }\forall \tau \in \left[t-{\mathrm{dt}}^{-},t+{\mathrm{dt}}^{+}\right],& \begin{array}{c}{\stackrel{\prime }{U}}_x\left(t,\tau \right)=2.\mathrm{Re}\left(\sum _{f\in F_{+}}c\left(f\right).P\left(t,f\right).{\hat{s}}_x\left(t,f\right).e^{1i\text{.2}\pi .f.\tau }\right)\\ {\stackrel{\prime }{U}}_y\left(t,\tau \right)=2.\mathrm{Re}\left(\sum _{f\in F_{+}}c\left(f\right).P\left(t,f\right).{\hat{s}}_y\left(t,f\right).e^{1i\text{.2}\pi .f.\tau }\right)\\ \stackrel{\prime }{U_x^{*}}\left(t,\tau \right)=2.\mathrm{Re}\left(\sum _{f\in F_{+}}c\left(f\right).P^{*}\left(t,f\right).{\hat{s}}_x\left(t,f\right).e^{1i\text{.2}\pi .f.\tau }\right)\\ \stackrel{\prime }{U_y^{*}}\left(t,\tau \right)=2.\mathrm{Re}\left(\sum _{f\in F_{+}}c\left(f\right).P^{*}\left(t,f\right).{\hat{s}}_y\left(t,f\right).e^{1i\text{.2}\pi .f.\tau }\right)\end{array}\end{array}$ $\begin{array}{cc}\mathrm{avec}:& \begin{array}{c}c\left(0\right)=0.5\\ \begin{array}{cc}\forall f\ne 0,& c\left(f\right)=1\end{array}\end{array}\end{array}$ $\mathrm{and}$ $\begin{array}{cc}\forall t\in T_{\downarrow }\forall \tau \in \left[t-{\mathrm{dt}}^{-},t+{\mathrm{dt}}^{+}\right],& \begin{array}{c}\stackrel{\prime }{R}\left(t,\tau \right)={\stackrel{\prime }{U}}_x\left(t,\tau \right)+1i.{\stackrel{\prime }{U}}_y\left(t,\tau \right)\\ \stackrel{\prime }{R^{*}}\left(t,\tau \right)=\stackrel{\prime }{U_x^{*}}\left(t,\tau \right)+1i.\stackrel{\prime }{U_y^{*}}\left(t,\tau \right)\end{array}\end{array}$

Starting from the time-frequency decomposition using a short-term two-sided Fourier transform, a step 142 of computing the intermediate filtered bivariate signal {grave over (R)}, its dual {grave over (R)}*, is obtained such that:

$\begin{array}{cc}\forall t\in T_{\downarrow }\forall \tau \in \left[t-{\mathrm{dt}}^{-},t+{\mathrm{dt}}^{+}\right],& \begin{array}{c}\stackrel{\prime }{R}\left(t,\tau \right)=P\left(t,0\right).\tilde{s}\left(t,0\right)+\sum _{f\in F_{+}^{*}}\left(P\left(t,f\right).\tilde{s}\left(t,f\right).e^{1i\text{.2}\pi .f.\tau }+P\left(t,f\right).\tilde{s}\left(t,f\right).e^{-1i\text{.2}\pi .f.\tau }\right)\\ \stackrel{\prime }{R^{*}}\left(t,\tau \right)=P^{*}\left(t,0\right).\tilde{s}\left(t,0\right)+\sum _{f\in F_{+}^{*}}\left(P^{*}\left(t,f\right).\tilde{s}\left(t,f\right).e^{1i\text{.2}\pi .f.\tau }+P^{*}\left(t,f\right).\tilde{s}\left(t,f\right).e^{-1i\text{.2}\pi .f.\tau }\right)\end{array}\end{array}$

Starting from the time-frequency decomposition using a short-term quaternionic Fourier transform, a step 143 of computing the intermediate filtered bivariate signal {grave over (R)}, its dual {grave over (R)}*, is obtained such that:

$\begin{array}{cc}\forall t\in T_{\downarrow }\forall \tau \in \left[t-{\mathrm{dt}}^{-},t+{\mathrm{dt}}^{+}\right],& \begin{array}{c}\stackrel{\prime }{R}\left(t,\tau \right)={\mathrm{Proj}}_{C_{1i}}\left(\sum _{f\in F_{+}}P\left(t,f\right).\stackrel{o}{s}\left(t,f\right).e^{1i\text{.2}\pi .f.\tau }\right)\\ {\stackrel{\prime }{R}}^{*}\left(t,\tau \right)={\mathrm{Proj}}_{C_{1i}}\left(\sum _{f\in F_{+}}{P}^{*}\left(t,f\right).\stackrel{o}{s}\left(t,f\right).e^{1i\text{.2}\pi .f.\tau }\right)\end{array}\end{array}$

Finally, the following equality 144 can be used to determine the filtered bivariate signal R and its dual R* by reconstructing the complete time axis as parts:

$\forall t\in T_{\downarrow }\begin{array}{c}\stackrel{\prime }{R}\left(\left[t-{\mathrm{dt}}^{-},t+{\mathrm{dt}}^{+}\right]\right)=R\left(t,\left[t-{\mathrm{dt}}^{-},t+{\mathrm{dt}}^{+}\right]\right)\\ R^{*}\left(\left[t-{\mathrm{dt}}^{-},t+{\mathrm{dt}}^{+}\right]\right)=\stackrel{\prime }{R^{*}}\left(t,\left[t-{\mathrm{dt}}^{-},t+{\mathrm{dt}}^{+}\right]\right)\end{array}$

Optionally, a step 150 is carried out of obtaining filtered input signals U₁ to U_M, by projecting the filtered bivariate signal R onto the axes collinear with the measured quantities.

The filtering template allows specific frequency ranges to be selected from the input signals.

In a particular embodiment, several filtering templates are applied successively and/or in combination. In addition, the filtering method can be applied to measurements carried out in different planes or for different quantities.

In a particular embodiment, all the steps of the method are computer-implemented, with the possibility of real time execution, notably in a magnetic bearing control cabinet.

Some terms and functions used in the computations described above are redefined below:

The set of strictly positive frequencies is: F₊*={f_j, j=[1, n]}

The set of strictly negative frequencies is: F₋*=−F₊

The following is also defined: F₊={0}∪F₊*

The following is also defined: F₋=F*∪{0}

The complete set of frequencies is: F=F₋∪F₊=F₋*∪{0}∪F₊*

The generic stationarity function is:

$\begin{array}{cc}\forall t\in \mathrm{TT},& \begin{array}{c}\mathrm{pour}\mathrm{tol}>0:\mathrm{autostat}\left(x\left(t\right),\mathrm{tol}\right)=e^{-\frac{❘\mathrm{dx}\left(t\right)❘}{\mathrm{tol}}}\\ \mathrm{pour}\mathrm{tol}=0:\mathrm{autostat}\left(x\left(t\right),0\right)=\left(\mathrm{dx}\left(t\right)==0\right)\end{array}\end{array}$

The maximum function with exception is:

${\max }^{*}\left(x\right)=\begin{array}{c}\max \left(x\right)\mathrm{si}\max \left(x\right)\ne 0\\ 1\mathrm{sinon}\end{array}$

The normalisation function is:

$\forall t\in \mathrm{TT},\forall f\in \mathrm{FF},{\mathrm{normalise}}_{\tau \in \mathrm{TT},v\in \mathrm{FF}}\left(x\left(t,f\right)\right)=\frac{x\left(t,f\right)}{{\max }_{\tau \in \mathrm{TT},v\in \mathrm{FF}}^{*}\left(❘x\left(\tau ,v\right)❘\right)}$

The division function with exception is:

${\mathrm{ratio}}^{*}\left(x,y\right)=\begin{array}{c}x/y\mathrm{si}y\ne 0\\ 0\mathrm{sinon}\end{array}$

The relative error function is:

$\mathrm{reldif}\left(x,y\right)={\mathrm{ratio}}^{*}\left(2.❘x-y❘,❘x❘+❘y❘\right)$

The definition of the arctan 2 function is also provided by way of a reminder, with a formulation equivalent to the argument of the complex number x+1μ. y:

$\mathrm{Arg}\left(x+1\mu .y\right)=\arctan 2\left(y,x\right)$

Claims

1. A method for filtering at least one signal to be filtered, with the filtering being determined based on measurements representing the movement of a rotating shaft, the method comprising: defining at least one characterisation indicator of the rotating shaft, with the characterisation indicator being obtained from a frequency analysis of the behaviour of the rotating shaft;determining a filtering template as a function of the characterisation indicator; andapplying the filtering template to said signal to be filtered.

2. The method according to claim 1, wherein the characterisation indicator is obtained from the following steps: acquiring at least two input signals over a predetermined period, the input signals comprising measurements representing the movement of the rotating shaft, the measurements being non-collinear to each other;determining two components of the input signals, the components being orthogonal to each other;merging the two components into a bivariate signal;time-frequency decomposing the bivariate signal into spectral elements;determining an equivalent ellipse for each spectral element; andextracting the characterisation indicator of the rotating shaft based on parameters of the equivalent ellipse.

3. The method according to claim 2, wherein the filtering template is applied to the time-frequency decomposition of the bivariate signal so as to obtain a filtered bivariate signal.

4. The method according to claim 3, wherein the time-frequency decomposition of the bivariate signal on which the filtering template is applied is a time-frequency decomposition that is obtained using a Fourier, or wavelet, transform, or by filtering, or using a quaternionic Fourier transform.

5. The method according to claim 3, further comprising a step of obtaining filtered input signals by projecting the filtered bivariate signal onto the axes collinear with the measured quantities.

6. The method according to claim 3, further comprising a step of computing a dual signal of the filtered bivariate signal.

7. The method according to claim 2, wherein the at least one indicator includes a polarisation rate, and/or an indicator of circularity and/or amplitude and/or direction and/or orientation and/or synchronisation of an equivalent ellipse.

8. The method according to claim 2, wherein the filtering template is defined relative to a predefined threshold of a polarisation rate, and/or a circularity indicator of an equivalent ellipse.

9. The method according to claim 1, wherein the filtering template is defined based on at least one normalised indicator ranging between −1 and 1 or between 0 and 1.

10. The method according to claim 1, wherein the filtering template is a linear bivariate filter and/or is obtained by multiplying a mask that is determined as a function of the characterisation indicator and of a frequency weighting.

11. The method according to claim 4, further comprising a step of obtaining filtered input signals by projecting the filtered bivariate signal onto the axes collinear with the measured quantities.

12. The method according to claim 11, further comprising a step of computing a dual signal of the filtered bivariate signal.

13. The method according to claim 12, wherein the at least one indicator includes a polarisation rate, and/or an indicator of circularity and/or amplitude and/or direction and/or orientation and/or synchronisation of an equivalent ellipse.

14. The method according to claim 13, wherein the filtering template is defined relative to a predefined threshold of a polarisation rate, and/or a circularity indicator of an equivalent ellipse.

15. The method according to claim 14, wherein the filtering template is defined based on at least one normalised indicator ranging between −1 and 1 or between 0 and 1.

16. The method according to claim 15, wherein the filtering template is a linear bivariate filter and/or is obtained by multiplying a mask that is determined as a function of the characterisation indicator and of a frequency weighting.