Signal filtering method

Abstract

[Signal Filtering Method] This method for filtering at least one signal, the filtering being determined from representative measurements of the motion of a rotating shaft, comprises the following steps: - definition of at least one indicator characterizing the rotating shaft, the characterization indicator being obtained from a frequency analysis of the rotating shaft's behavior (step 120); - determination of a filtering template based on the characterization indicator (step 130); and - application of the filtering template to the signal to be filtered (step 140). Figure for the abstract: Fig. 3

Description

Title of the invention: Method for filtering a signal Technical field of the invention

[0001] The present invention relates to the characterization of a shaft rotating around an axis of rotation.

[0002] The present invention relates more particularly to the filtering of an indicator for characterizing the gyroscopic effects and deformations of a rotating shaft in order to facilitate its control, modeling, monitoring, diagnosis and maintenance. State of the prior art

[0003] A rotating system, such as a rotating shaft, perfectly balanced rotates around an axis of rotation so that said axis of rotation is stationary relative to its theoretical positioning, and therefore relative to the environment close to the rotating shaft.

[0004] However, defects, deformations, gyroscopic effects, or even effects due to the force of gravity can induce an offset of the axis of rotation relative to its theoretical position during the rotation of said rotating shaft. This is the case, for example, of deformations of the rotating shaft due to flexible natural modes, these modes deforming the shaft as a function of the rotation speed of said shaft.

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

[0006] The properties of a rotating shaft are conventionally understood by 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.

[0007] It is often necessary to characterize the changes in these quantities depending on whether the shaft is rotating or not. However, multiple causes can influence the rotational behavior of a rotating shaft and it is complex to precisely characterize such behavior using reliable indicators. Summary of the invention

[0008] The present invention therefore aims to overcome the aforementioned drawbacks and to propose filtering of a signal to be filtered characterizing a rotating shaft.

[0009] The present invention therefore relates to a method for filtering at least one signal at filtering, the filtering being determined from measurements representative of the movement of a rotating shaft, the method comprising the following steps:

[0010] - definition of at least one indicator for characterizing the rotating shaft, the characterization indicator being obtained from a frequency analysis of the behavior of the rotating shaft;

[0011] - determination of a filtering template as a function of the characterization indicator; And

[0012] - application of the filtering template to said signal to be filtered determined from measurements representative of the movement of the rotating shaft.

[0013] Thus, this method makes it possible to make the characterization indicator more reliable and to isolate certain information facilitating the understanding and rectification of defects on a rotating shaft.

[0014] Advantageously, the characterization indicator is obtained from the following steps:

[0015] - a) acquisition over a predetermined period of at least two input signals, the input signals comprising representative measurements of the movement of the rotating shaft, the measurements, linked to a vector quantity, being non-collinear with each other and optionally coplanar with a plane intersecting the axis of rotation of the shaft;

[0016] - b) determination of two components of the input signals, the components being orthogonal to each other and optionally to the axis of rotation;

[0017] - c) merging the two components into a bivariate signal;

[0018] - d) time-frequency decomposition of the bivariate signal into spectral elements;

[0019] - e) determination of an equivalent ellipse for each spectral element; and

[0020] - f) extraction of the indicator characterizing the rotating shaft from pa parameters of the equivalent ellipse.

[0021] The gyroscopic effects and the rotor deformations, in particular the flexible modes, are simpler to consider under 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 characterized using a frequency decomposition which can highlight a rotating orbit, also called an equivalent ellipse, for each frequency, each orbit rotating at a specific speed corresponding to the frequency considered.

[0022] This characterization of a rotating tree allows for applications such as control, identification, modeling, monitoring or diagnostics.

[0023] 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.

[0024] In various embodiments, the time-frequency decomposition of the signal bivariate on which the filtering template is applied is a time-frequency decomposition obtained by Fourier transform, or by wavelet, or by filtering, or by quaternion Fourier transform.

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

[0026] In a particular embodiment, the method comprises a step of calculating a dual signal of the filtered bivariate signal.

[0027] In one embodiment, the at least one indicator comprises a polarization rate, and / or an indicator of circularity and / or amplitude and / or direction and / or orientation and / or synchronization of an equivalent ellipse.

[0028] Advantageously, the filtering template is defined with respect to a predefined threshold of a polarization rate, and / or of a circularity indicator of an equivalent ellipse.

[0029] Advantageously, the filtering template is defined from at least one standardized indicator between -1 and 1 or between 0 and 1.

[0030] In a particular mode of implementation, the filtering template is a linear bivariate filter and / or is obtained by multiplication of a mask determined as a function of the characterization indicator and a frequency weighting.

[0031] 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 the characterization indicators of the planes.

[0032] The present invention also relates to a computer-implemented method comprising the steps of the method as defined above. The implementation can be carried out on a real-time target, in particular in a magnetic bearing control cabinet. Brief description of the figures

[0033] [Fig.l]

[0034] is a schematic representation of the different steps of a method for obtaining an indicator for the method according to the invention;

[0035] [Fig.2]

[0036] is a representation of an ellipse and its Euler parameters; and

[0037] [Fig.3]

[0038] is a schematic representation of different filtering steps applicable to a signal to be filtered following the steps of [Fig.l]. Detailed description of the invention

[0039] [Fig.l] shows schematically different stages of obtaining a indicator which can, in a particular mode of implementation, be incorporated into the method according to the invention.

[0040] These steps make it possible to characterize the behavior of a rotating shaft. A rotating shaft is understood to mean a movable shaft rotating around an axis of rotation. The rotating shaft is, for example, a shaft of a rotating machine.

[0041] To implement the method, a step 10 of acquiring input signals is first carried out.

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

[0043] The input signals comprise measurements representative of the movement of the rotating shaft. For example, the input signals comprise measurements of position, or speed, or acceleration, or force, or magnetic field, or even magnetic induction of the rotating shaft. The acquisition of the input signals is carried out for example using sensors configured to measure these quantities.

[0044] Generally, a vector quantity is measured according to several components distributed in space sensitive to the rotation of the shaft. The distribution in space of these components makes it possible to evaluate the movement of the rotor.

[0045] A special case is that of measurements that are non-collinear with each other and coplanar with a plane intersecting the axis of rotation of the shaft. From two non-collinear measurements, a characterization of the behavior of the rotating shaft is possible in a plane intersecting the axis of rotation of the shaft. Additional measurements allow greater precision of said characterization.

[0046] Non-coplanar measurements with a plane secant to the axis of rotation of the shaft can nevertheless be brought back by projection into a reference plane secant to the axis of rotation of the shaft.

[0047] For M measurements with M>2, we obtain M input signals u। ( t ) to uM ( t ), t designating the time during the acquisition of the measurements, the set of values ​​of which is designated by T, in other words dates, which t can take.

[0048] A secant plane is understood to mean a plane in which the measurements are made which intersects the axis of rotation of the rotating shaft. The secant plane is, for example, orthogonal to the axis of rotation or non-orthogonal.

[0049] In the case where the secant plane is not orthogonal to the axis of rotation, a step 20 is carried out for determining two components of the input signals, 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 rest of the calculations.

[0050] In the case where the secant plane is orthogonal to the axis of rotation, the previous step is intrinsically carried out.

[0051] In the case where the axis of rotation is not rigorously and precisely defined, the shaft not always being perfectly cylindrical and / or balanced, and / or non-deformable and / or known, said orthogonal secant plane is assimilated to a study plane assumed to be orthogonal to the axis of rotation, or at least sensitive to the rotational movement of the shaft.

[0052] In practice, step 20 of determining the two orthogonal components ux(t) and Uy(t) is carried out by averaging projections of the input signals at UM(t) onto two orthogonal axes x and y, these two axes being ideally orthogonal to the axis of rotation. Thus,

[0053] "A i) = iMi Ut(t) ) ) VteT, uy(t) - moyi=[ [Proj^uM ) )

[0054] The components mv(t) and uy(r) are therefore now in the study plane, ideally orthogonal to the axis of rotation.

[0055] Then, a step 30 is carried out to merge the two orthogonal components ux(t ) and wy( / ) into a bivariate signal r(l) in order to combine these two components in an equivalent manner in the form of a complex signal, such that:

[0056] VîeT. r(t) + U.uy(t)

[0057] It can also be expressed as follows:

[0058] r( / ) - {mv( / ), My(r)}

[0059] Or in vector form:

[0060] z , [M01 L»v(O

[0061] A step 40 of frequency decomposition of the bivariate signal r ( t ) into spectral elements characterized by their time-frequency pair ( / , / ), also called time-frequency cells, is then carried out periodically in order to obtain a decomposition that is both temporal and frequency, called time-frequency decomposition s(t, / ) such that:

[0062] s(t,f) = spectrumt(r(t))

[0063] There are different ways to perform this time-frequency decomposition step 40. These different implementation methods are presented below, with the time-frequency decompositions, or s, or î, or j, being given in continuous, i.e., analog and non-discrete form. A computer-implemented digital implementation then includes a time and frequency sampling step.

[0064] In a first mode of implementation, noted 41 in [Fig.l], the decomposition step 40 time-frequency position is achieved with a filter bank around each frequency f.

[0065] Thus, the time-frequency decomposition f ) is expressed:

[0066] ^feF^ytET, s(t, / )

[0067] with hf the impulse response of a bandpass filter centered around the frequency f, F+ being the set of frequencies greater than or equal to 0.

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

[0069] Basically, the time-frequency decomposition is expressed:

[0070] A .X 1 / X , . > 0 'f = W*' f') =

[0071] In practice, it is preferable to use more efficient estimators based on a division into time intervals, with time weightings, averaging over several time intervals, and overlaps between successive intervals. The formulation is then:

[0072] A / [ / +< 7 X / . y. f SXk( F . / ) = J^- Ux(t ).Wk(r- V and G Syk(Af) = Uy(r).wk(rt)^u2^fydr

[0073] Mtf]-™? VtEFyfEF^ '

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

[0075] The overall time extent of the corresponding cell (t,f) is: dt+ = maxk^p]{dt+k) dt ~ maxk=yj ( dtk ) dt = dt + dt+

[0076] And wk(r) for r G [ - df^ dtk ] for k- [ 1, p] defining a time weighting window for each of these intervals, for example of Hanning type. Classically, wk does not depend on k, identical for each interval.

[0077]

[0078]

[0079]

[0080]

[0081]

[0082]

[0083]

[0084]

[0085]

[0086]

[0087]

[0088]

[0089]

[0090]

[0091]

[0092]

[0093] Depending on the estimated spectral quantity, it is necessary to introduce compensations following: r . = r . V kŒ [ 1, p 1, V k G [ 1, p 1, In a special case, the absence of weighting corresponds to M7=l, in which case bX® = FF = dtk. A. K & In addition, we introduce a subset of dates: = { f-, z = [ 1, fri]} In a third embodiment 43, the time-frequency decomposition step 40 is carried out by Fourier transform, in particular short-term double-sided Fourier transform, also called “Fullspectrum” in English terms, in other words for a set F of frequencies taking values ​​in the set of real numbers. Thus, the basic time-frequency decomposition s(t, f) is expressed as: VfER, S(t, f) = J(T)£- Or with a division into time blocks: V tg Tif V fe F, sk J r(r) .wk (rt) .e-ïi2jTfr£lT V t Tit Y f^F, s(t, f) = avg^) In a fourth embodiment 44, the time-frequency decomposition step 40 is carried out by Fourier transform, in particular short-term quaternion Fourier transform, involving the unitary imaginary numbers verifying li2 = lj2 = 1k2 = li.lj.lk = -1. Thus, the basic time-frequency decomposition s(t, f) is expressed as: V f > 0, s(t,f J (_dtr ( r ) Or with a division into time blocks: d+dt^ VIeT^ / eF* sk(t,f)=]t.dfk r(T).wk(t-1) VtEL, V / eF^. n ,() j +- \ jf • Æ=|i,p]yy The expressions of the Fourier transforms given in implementation modes 41 to 44 above constitute applicable formulas, producing efficient estimators, particularly in their formulation with time weightings, averaging over several time intervals, overlaps between successive intervals, for which a more general expression can be given: n f) = r(r)AV^

[0094] wi being the time weighting, and being defined such that l^2 = - 1.

[0095] Finally, [0°96] § (A y) = mOyk (f) j

[0097] The index k applied to wk (above) or pk (below) means that the weighting can differ from one interval to another.

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

[0099] A more general expression is:

[0100] □ f) = r(r).pk(T, f

[0101] More generally still, it is possible to vary the time division and the weighting as a function of the time-frequency cell (t, / ):

[0102] D Wf) r\r)-Pk^f^rydr

[0103] The trajectory (orbit) is obtained by inverse transformation, of the type:

[0104] rf+dP(t,f ) zn n \ f) (?) = (jU v)^(L f) (r, v) +s(f, -v)4(lf) (r, - v) j Jv

[0105] The frequency range of the cell (t, f) then being df = d f+ + df, the time range of the cell (t, f) being dt = dl' + df.

[0106] From each spectral element t, f ), at a fixed date t and frequency f, we can deduce from the instantaneous orbit an equivalent ellipse. We can describe it as local in the temporal and frequency sense, in other words at date t and frequency f.

[0107] A step 50 is therefore carried out to determine the equivalent ellipse for each spectral element. Whatever 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 Sx», and 0y the instantaneous amplitudes and phases:

[0108] ,,( t, f) = >• / ) )

[0109] s,(t,f) = f).cos(0,(t,f))

[0110] For a fixed frequency, going through time t allows obtaining the trajectory at this frequency.

[0111] This instantaneous orbit takes the form of an equivalent ellipse, provided that Sx(t,f]^y(t,f)), 0y(t, f) vary slowly as a function of time, for example during a complete revolution at frequency / and in all cases during the duration dt of the time interval defining the spectral element.

[0112] In this case, we can distinguish the time t of the spectral element considered from the time T that evolves in this spectral element and. Moreover, by introducing the frequency / in From the expression of the phases assumed to be quasi-linear, we obtain the equivalent ellipse in the form of a time-parametric equation for each component x and y:

[0113] (r, f) (T) = ax (f).cos (2jr.fr + <px ( t, f ) ) [°114] sft,f)(T) = ay(t, f).COS fa.fr+

[0115] Another way to write these parametric equations is, for example: [0H6] sx(t, f) (r) - ax (t,f).( cos (2jr.fr). <py(t, f) ) - sin(2jTfr).sin[<py(t, f) ) )

[0118] Pour obtenir ces équations paramétriques, notamment dans le cas d’une décomposition fréquentielle par banc de filtres, les paramètres ax(L f), <p (. t, f), av( t, f ) , (p {t, f) de l’ellipse équivalente peuvent être obtenus par estimations, tels que moyennage direct ou quadratique ou par valeur efficace, à partir des grandeurs instantanées (t, f ) pour T&[t-df,t + dt+].

[0119] In the case of a frequency decomposition by Fourier transform, we obtain the parametric equation of the equivalent ellipse by inverse Fourier transform of the frequency element:

[0120] sx (t,f)(r)= 2Re{ sx (t, f).e'i2^) Sy (t,f)(r)= 2 Re ( {t, f ) )

[0121] By forming the associated complex quantity — sft, f) + f ), in other words the bivariate representation, we highlight the decomposition of a parametric equation of the equivalent ellipse into two circles of opposite directions:

[0122] s(t, f) (t) = aft, f).eli^W ffe-^f'le^^

[0123] With the following parameters obtained following a calculation step 59: a+ = + a} + 2.ax.ay.sin ( <px - (py ) has. = f + ay + IcpctySin ( (py - tpx j 3+ - arctan2 ( a^sin ( (px ) + ay.cos ( (p ), ax.cos ( <px ) - ay.sin ( ¢7 ) j 3- = arctan2(axsin( <px) -ay.cos(<py), ax.cos(q> x) + avsin( <py) J

[0124] The characterization of the equivalent ellipse is natural with the double-sided Fourier transform (“fullspectrum”) which can be directly written in polar form:

[0125] sUj)=aJy,f).e^^

[0126] Generally, an equivalent ellipse can be determined and plotted by obtaining its Euler parameters.

[0127]

[0128]

[0129]

[0130]

[0131]

[0132]

[0133]

[0134]

[0135]

[0136]

[0137]

[0138]

[0139]

[0140]

[0141]

[0142] Figure 2 shows an example of an ellipse and its Euler parameters a(t, f), which completely define an ellipse. equivalent, with: a>0 Ge [ - V2, V2 ] X e [ - V4, V4 ] (p E [ - TT, ] 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, whose sign gives the direction of rotation of the equivalent ellipse. The parametric equation of the trajectory of the equivalent ellipse is more convenient in complex form, because: - we express the equation of a circle, notably in terms of cos and sin of the same phase, - the axes are dilated to obtain an ellipse, - a rotation is applied to it to obtain the desired direction, like this: .v( t, f ) ( T) ​​= a( t, f ).( cos(X( f, f ) ).cos(ht.fr + <p(Lf ) ) + ILsin) )-sin( + ^0. / ))) The direction is defined as the angle 0 of the equivalent ellipse between the major axis and the x-axis. The orientation is defined as the direction of rotation describing the trajectory of said ellipse. This characterization of the equivalent ellipse is natural with the time-frequency decomposition by quaternion Fourier transform which can be directly put in Euler form: f ) = Rf, f) This parametric equation has the particular advantage of rejecting the time parameterization towards a single exponential term |orSqUC pon applies the transform inverse Fourier transform to the frequency element in order to obtain the time-parametric equation of the trajectory of the equivalent ellipse: with Projç a projection into the set of complex numbers C h- More specifically, step 50 of determining the equivalent ellipse includes a step of calculating Euler parameters of said equivalent ellipse. Starting from the time-frequency decomposition by filter bank, we obtain a calculation step 51 of ax and ay such that:

[0143] Ti being a subset of T.

[0144] Starting from the time-frequency decomposition by short-term single-sided Fourier transform, we obtain a step 52 of calculation of ax, ay, and Vy such that:

[0145] ... ax{t, f) =2.\sx{t, f)\ V t G Ti V f G FL " J (py(t, f) - Arg(sy(t, f))

[0146] Starting from the time-frequency decomposition by short-term double-sided Fourier transform, we obtain a step 53 of calculation of a+, a~, 8+ and such that: [°147] VttT VfeF* = e+(t,f)=Arg(s(Kf)) ' = |S(Af)\ 8-(Lf) = -Arg(s(t, - / ))

[0148] By an additional calculation step 54 or 55, the following Euler parameters are obtained:

[0149] a = n .... W- o- 2 X = arctan2 ( a+ -a^a++a. )

[0150] Step 56 of inverse calculation is also possible, such that:

[0151] a+-^A cos(x) ) a = f.(œ>(z)-sin(x)) 8+= <p +8 e.= <p-e

[0152] The method optionally includes a step 60 for calculating the spectral density of the time-frequency decomposition. The spectral densities of each component and the cross-components allow for estimating the energy distribution between the 8-^,arctan2(2ax.ay.cos(q>x-(pJ, a}-(L) y = iœrcsin( ratio* (2üx£iy, ay+aL}^in{ m- x and y components.

[0153]

[0154] Starting from the time-frequency decomposition by filter bank, we obtain a step 61 of calculation of the spectral density S such that: p+dt+

[0155] Starting from the time-frequency decomposition by short-term single-sided Fourier transform, we obtain a step 62 of calculating the spectral density g such that:

[0156] ho 7 / M t, f) - 2.df.moyk=ï, ............. L -< j \ £ y vte / , vfe fZ t, / ^ Zdfmay... ( ) n.—l | \ h. / *■ \ To' / s„{ t,f) = W»wM w [ —ÿ— )

[0157] Starting from the time-frequency decomposition by short-term double-sided Fourier transform, we obtain a step 63 of calculating the spectral density S such that: 101581 VteT» VfeF', 5(t,f)=dfJnOy^M^^

[0159] We will also use:

[0160] y < czT f (z F* Ç / tf \ fI f mnv ( ^f )^k ) \ Sdjd )-djjnoy^^--

[0161] Starting from the time-frequency decomposition by short-term quater-nionic Fourier transform, we obtain a step 64 of calculating the spectral density $ such that that :

[0162] VJ T' F rc / , JF / if / 4(Af V tg TV fg S(t,f)= dfmoyk=[ I j + df.moyk=[

[0163] From these spectral densities S,$, S, and we carry out a step 70 of calculating the Stokes parameters 50, 5b S2, and <S3 d’au moins une ellipse équivalente.

[0164] The Stokes parameters are four real numbers and are classically used to describe the polarization state of an electromagnetic wave.

[0165] So is commonly the total intensity of an electromagnetic wave and represents here by extension the total intensity of the bivariate signal, S0>0.

[0166] S] and S2 give the linear polarization intensity <2 2 as well as the management of the polarization of the ellipse () - 1 atan2( S2, S} ) •

[0167] S3 gives the circular polarization intensity |53|, its sign giving the direction of rotation.

[0168] We also define the reduced Stokes parameters: [°169] «5o*O, V / e {1,2,3}

[0170] Starting from the spectral density obtained for the time-frequency decomposition by filter bank, we obtain a step 71 of calculating the Stokes parameters of at least one equivalent ellipse such that: 101711 TW / V t^T„ V feF^ , ,

[0172] Starting from the spectral density obtained for the time-frequency decomposition by short-term single-sided Former transform, we obtain a step 72 of calculating the Stokes parameters of at least one equivalent ellipse such that:

[0173]

[0174]

[0175]

[0176]

[0177]

[0178]

[0179]

[0180] s0(fj')=2(Jta / ) Aa / )) 5,(7, / )=2(5,( / , / )-5,( / . / )) 52( / , / )=4.R<"(5,,( / , / )) S3(z, / )=4Jm(i„(z. / )) Starting from the spectral density obtained for the time-frequency decomposition by short-term double-sided Fourier transform, we obtain a step 73 of calculation of the Stokes parameters of at least one equivalent ellipse such that: 5,,( / , / ) = 2(5( / , / )+S(t - / )) V / eTt, V / eF^ 5,(7, / )=4^(5.( / , / )) 5,( / , / ) =47m(57( / , / ) ) S3( ¢, / )=2(5( / , / )-5( / , - / )) Starting from the spectral density obtained for the time-frequency decomposition by short-term quaternion Fourier transform, we obtain a step 74 of calculating the Stokes parameters of at least one equivalent ellipse such that: 50( / , / ) = Re(5( / , / )) 5,( / , / )=^,,(5( / , / )) yf^F+, 5,( / , / ) = / m,((5( / , / )) 53( / , / ) = / » / ,,(5( / , / ) ) Alternatively, the Stokes parameters Sb, S2, and S3 are calculated directly using the time-frequency decomposition by quaternion Fourier transform: W) S3U / ) From these Stokes parameters, we carry out a step 75 of calculating the polarization rate 0 as well as the reduced Stokes parameters previously mentioned, as well as the polarized part of the intensity Sp. [«SU Life {1,2,3) = V te T;, V / e F^ 0(t,f) = f)+sÿ.t, f) F si(t, f)

[0182] The polarization expressed by the polarization rate 0 makes it possible to highlight the cyclo-stationarity of an equivalent ellipse, namely the regularity of the trajectory of said ellipse with respect to a period of revolution. The polarization reflects the structuring of said trajectory relative to defined frequencies. The polarization rate can thus be calculated 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.

[0183] Furthermore, a calculation step 57 can be implemented to determine Euler parameters a, 0, and X from the polarization spectral density Sp and the reduced Stokes parameters:

[0184] V / eF^, e(t, f) = ~£irctatâ^ {.arcsin( ratio*(.¾(t, f ), Sp( t, f ) ) )

[0185] A fourth Euler parameter V giving the phase at the origin of the trajectory can also be deduced by the following calculation step 58:

[0186] vtE^, V / eFZ (p{t. f) = Arg^Proj^e^W

[0187] Following step 50 of determining an equivalent ellipse, in particular by calculating Euler parameters of said equivalent ellipse, a step 80 is carried out of extracting at least one characterization indicator of the rotating shaft from parameters of the equivalent ellipse, for example Euler parameters, Stokes parameters or parameters derived from the latter.

[0188] 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 carried out in the first place. In this embodiment, it is thus possible to carry out a step of comparing the characterization indicators of the planes.

[0189] For example, the at least one characterization indicator comprises the previously described polarization rate 0, and / or a circularity indicator (linked to | / | k and / or an intensity indicator and / or an amplitude indicator, and / or a rotation direction indicator (linked to the sign of X\ and / or a direction indicator (linked to 0), and / or synchronization (linked to the equivalent ellipse.

[0190] Each elementary indicator contributes to a concrete and useful characteristic of the

[0191]

[0192]

[0193]

[0194]

[0195]

[0196]

[0197]

[0198]

[0199]

[0200]

[0201]

[0202]

[0203]

[0204]

[0205]

[0206]

[0207]

[0208]

[0209]

[0210]

[0211] In a rotating system, for example, a high degree of polarization and high circularity are likely due to the rotating shaft. Distinguishing this part of the signal, with the complement being considered disturbances, can be useful for, for example, identifying the shaft (for modeling purposes) or better controlling it in the case of a shaft equipped with magnetic bearings. Optionally, at least one characterization indicator is normalized, in other words, its value is between -1 and 1 or between 0 and 1. The characterization indicator is given for a plan, it is then noted A, or by comparison between two plans, it is then noted I. The characterization indicator is called a criterion, noted C, when it is constructed from any characteristic quantity re-evaluated at each instant t. The characterization indicator is called stationarity, noted S, when it is based on the stationarity over time of a quantity (maximum when the variation over time of said quantity is zero). Thus, we note the indicator of the polarization rate: T^ V f G F+, AC -polarization^ CJ ) = O( t. f ) For a plan, the intensity indicator is: V te V / g F+, ACintensity- ( t, f ) = normalize^, ( S'o ( t, f ) ) With the normalization function: V te TT, V fe FF, normalizes (x(t, f)) = —;.................. TT and FF denote sets of generic dates and frequencies, taking their value from the set of real numbers R. For a plane, the amplitude indicator is: V t G f^F^, ACampUtude normalize^T^^ (a(t,f)) For a plane, the sign indicator is -1 or 1 depending on the sign of X; V t G TV / eF* ACsign (t,f) = sign Mcf)) For a plane, the circularity indicator is: V t GT^, V f G ACcircuhirity For a plane, the indicator of stationarity of sign, in other words of direction of rotation, in other words of orientation, is: V te T^ ASsign(t, / ) = autostat(ACsign( t, f), you = 0) With the stationarity function: for you > 0: autostat(x(t), you) = V t G TT, for you - 0 : autostat0) - (dx(t) - - 0) For a plane, the circularity stationarity indicator is:

[0212] V / EF^ AScimdari^

[0213] For a plane, the indicator of direction stationarity, in other words of geometric inclination of the ellipse, is:

[0214] YtE^, ASdb.ecti(m(t,f)=autostM

[0215] For a plan, the synchronization stationarity indicator is:

[0216] V te V fs FASynchro (F f) = autostat ((p(t, f) -to^ - 1 / 180)

[0217] For two plans, we redefine each parameter or indicator for a plan ad additional, by noting these new parameters using an apostrophe, for example a'

[0218] The comparative characterization indicator of the amplitude is:

[0219] VtEl], V f G F+, ICampiitude(t,f) = normalizesT^T f))

[0220] The comparative characterization indicator of the sign is a boolean between 0 and 1 according to the result of the following test:

[0221] vteTf yfeF. =

[0222] The comparative circularity characterization indicator is: [02231 V te Tt, V ,fe F. = rMipAC„,Mi (t, f), ACt, f)}

[0224] With the relative error function such that: reldif^x, y) = ra / z6>*(2.|xy|, + |_y| )

[0225] And with the division function with exception such as: . «z , x / y if y^0 ratto ( x9 y) — 0 otherwise

[0226] The indicator of comparative characterization of direction, in other words of geometric inclination of the ellipse, is:

[0227] VteTi, IfeF, IC,

[0228] The comparative synchronization characterization indicator is:

[0229] VtET^YfEF» ICxynch^

[0230] The indicator for characterizing stationarity and comparative circularity is:

[0231] YteT» ISdrcularity(t, f) =autostat(ACcircuiarity(j, f)-ACdrM^t,f), tol^A}

[0232] The indicator of stationarity characterization and comparative direction is:

[0233] VfGï;, Vf^F^ ISdirec[ion(t,f)=autostat(e(t,f)-e'(t,f),tol=l0jr / 18())

[0234] The indicator for characterizing stationarity and comparative synchronization is:

[0235] vte TV fe F+, 1Ssvnc}:m (t, f) = autostat ((p(t, f) - tp' (t, f), you = 1 Qjt / 180)

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

[0237] Optionally, a step 90 is carried out for inserting at least one indicator of characterization in a neural network in order to classify said indicator according to the rotation of the tree and link it to a previously listed anomaly. Anomalies similar to the learning base are then detected.

[0238] Optionally, a step 100 is carried out for displaying the at least one extracted characterization indicator to be readable by an operator and allow his understanding of potential anomalies of the rotating shaft.

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

[0240] [Fig. 3] schematically shows the steps of the filtering method according to an example of the invention. In a particular embodiment, the following steps are carried out following step 80 of extracting a characterization indicator.

[0241] The filtering method is intended to filter at least one signal to be filtered, for example a signal determined from measurements representative of the movement of a rotating shaft. The method is carried out in real time or not. Said signal may 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 the bivariate signal determined in step 30.

[0242] In order to carry out this filtering process, a step 120 is first carried out to define at least one indicator for characterizing the rotating shaft, the characterization indicator being obtained from a frequency analysis of the behavior of the rotating shaft.

[0243] In a particular embodiment, the characterization indicator is an indicator extracted during step 80 described above.

[0244] Thus, the at least one characterization indicator comprises the previously described polarization rate 0, and / or a circularity indicator, and / or an intensity indicator and / or an amplitude indicator, and / or a rotation direction indicator, and / or a direction indicator, and / or synchronization of the equivalent ellipse.

[0245] A step 130 is then carried out for determining a filtering template as a function of the at least one characterization indicator.

[0246] Step 130 of determining a filtering template comprises for example a step 131 of calculating a mask M, such that:

[0247] V f^F^, = combine(ACite!.^)

[0248] Where combine is a function that takes as input one or more characterization indicators with values ​​in [0,1] or [-1,1] and returns as output a real value in [0,1].

[0249] Following step 131 of calculating the mask M, a calculation step 132 of the filter template P and its dual P* is performed such that:

[0250]

[0251]

[0252]

[0253]

[0254]

[0255]

[0256]

[0257]

[0258]

[0259] P(t, f) = M(t, ff) vt E 1 i, V / tr . .. With N a frequency weighting with values ​​in [0,1]. N is an arbitrary frequency weighting, independent of M and is similar to a classic frequency filtering. N allows for example to focus on a frequency band around the speed of the shaft. And the dual (complementary) mask: M* = 1 - M The interest of considering M* and consequently P* is that depending on the application, we can be interested in a certain part of the signal, obtained according to the characterization indicators used, given by M and / or the remainder assimilated to the disturbances given by M*. For example, we can imagine a magnetic bearing controller acting 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 with respect to a predefined threshold of a characterization indicator, for example a predefined threshold of a polarization rate, and / or of a circularity indicator of an equivalent ellipse as described previously. The constancy of a characterization indicator also makes it possible to define a filtering template with respect to said characterization 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 from measurements representative of the movement of the rotating shaft. In one implementation mode, 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 dual R* calculated jointly. Step 140 is then carried out according to different implementation modes. Starting from the time-frequency decomposition by short-term single-sided Fourier transform, we obtain a step 141 of calculation of an intermediate filtered bivariate signal R, of its dual jy and of the components of the input signals determined in step 20 filtered JJ vy*, y*, and (jy. Ù,Mr}=2J <e(yePc(^ V î £ F. V r G [t - dt, t + dt |, ^(6 t) = ).y / , Uy(t, t) = 2.Re(^feFc(f).P\t, c(0) = 0.5 with: z, V / V0. = \

[0260] and [026!] R(t,T) = Üx(t,r) + tiÜy(j,T) V / GTVTG 1t - dt, t + dt |$ / \ (t, t) = Ux(t,?) + lLUy(t, r)

[0262] Starting from the time-frequency decomposition by short-term double-sided Fourier transform, we obtain a step 142 of calculation of the intermediate filtered bivariate signal R, of its dual R*, such that:

[0263] ^0,r) = pu,O).s(z,0) +£.( / '( / , f)F^+P(t, -RsU, y mT.VT dt. t^dd\. , R^r) = P'\t,0):sU,()) -fRït,

[0264] Starting from the time-frequency decomposition by short-term quaternionic Fourier transform, we obtain a step 143 of calculation of the intermediate filtered bivariate signal R, of its dual R*, such that:

[0265] _ RUT)=ProjJï^ V / g V tgt - dt, t + dt 1, ... z & ( t, T ) = Pr°jC]i (^feFP* ( L f ) -K 6 f ) £lj27TpT )

[0266]

[0267]

[0268]

[0269]

[0270]

[0271]

[0272]

[0273]

[0274]

[0275]

[0276] Finally, the following equality 144 allows us to determine the filtered bivariate signal R and its dual R* by reconstructing the complete time axis piecewise: R([t-dt, t+ dt+] ) = R(t, [t-dt, t + ) V tg T$ R*( [t-dt, t + dt+ ] ) = R (t, [t-dt, t + ) Optionally, a step 150 is performed to obtain filtered input signals. U} to UM, by projection of the filtered bivariate signal R onto the axes collinear with the measured quantities. The filter template thus allows the selection of specific frequency ranges on the input signals. In a particular implementation method, several filtering templates are applied successively and / or in combination. Furthermore, the filtering process can be applied to measurements taken in different planes or for different quantities. In a particular implementation method, all steps of the process are implemented by computer, with the possibility of real-time execution, particularly in a magnetic bearing control cabinet. Some terms and functions used in the previously detailed calculations are redefined below: The set of strictly positive frequencies is: j' — [ 1 «] ] The set of strictly negative frequencies is: F* — -F We also define: F+ = { 0} U We also define: F_ = F* U (0}

[0277] The complete set of frequencies is: F = F- U F+ = F- U {0} U Fl

[0278] The generic stationarity function is: for you >0: autostat(x(t), you) = V re TT, for you = 0 : autostat( x ( f ), 0 ) = (dx(t) = = 0)

[0279] The maximum function with exception is: sz x max(x) if max(.x)^0 max'(x) = v 7 1 otherwise

[0280] The normalization function is: V te TT, VfE FF, normalizes „„ (x(t, f)) = —;——___

[0281] The division function with exception is: . , x / y if y 0 ratio (.x, y) — 0 otherwise

[0282] The relative error function is: reldif(x, y) = ratio* (l.\x - |x| + |y| )

[0283] Let us also recall the definition of the arctan2 function, a formulation equivalent to the argument of the complex number x + :

[0284] Arg (x + 1 / ty) - arctaril (y, x)

Claims

Claims

1. Method for filtering at least one signal to be filtered, the filtering being determined from measurements representative of the movement of a rotating shaft, characterized in that it comprises the following steps: - definition of at least one characterization indicator of the rotating shaft, the characterization indicator being obtained from a frequency analysis of the behavior of the rotating shaft (step 120); - determination of a filtering template as a function of the characterization indicator (step 130); and - application of the filtering template to said signal to be filtered (step 140).

2. Method according to claim 1, in which the characterization indicator is obtained from the following steps: - a) acquisition over a predetermined period of at least two input signals, the input signals comprising measurements representative of the movement of the rotating shaft, the measurements being non-collinear to each other (step 10); - b) determination of two components of the input signals, the components being orthogonal to each other (step 20); - c) merging of the two components into a bivariate signal (step 30) - d) time-frequency decomposition of the bivariate signal into spectral elements (step 40); - e) determination of an equivalent ellipse for each spectral element (step 50); and - f) extraction of the characterization indicator of the rotating shaft from parameters of the equivalent ellipse (step 80).

3. The method of 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 of claim 3, wherein the time-frequency decomposition of the bivariate signal to which the filtering template is applied is a time-frequency decomposition obtained by Fourier transform, or by wavelet, or by filtering, or by quaternionic Fourier transform.

5. Method according to one of claims 3 and 4, further comprising a step (150) of obtaining filtered input signals by projection of the filtered bivariate signal onto the axes collinear with the measured quantities.

6. A method according to any one of claims 3 to 5, further comprising a step of calculating a dual signal of the filtered bivariate signal.

7. A method according to any one of claims 2 to 6, wherein the at least one indicator comprises a polarization rate, and / or an indicator of circularity and / or amplitude and / or direction and / or orientation and / or synchronization of an equivalent ellipse.

8. A method according to any one of claims 2 to 7, wherein the filtering template is defined relative to a predefined threshold of a polarization rate, and / or a circularity indicator of an equivalent ellipse.

9. Method according to any one of claims 1 to 8, in which the filtering template is defined from at least one normalized indicator between -1 and 1 or between 0 and 1.

10. Method according to any one of claims 1 to 9, in which the filtering template is a linear bivariate filter and / or is obtained by multiplication (step 131) of a mask determined as a function of the characterization indicator, and of a frequency weighting.

Images