METHOD, DEVICE AND COMPUTER PROGRAM FOR MONITORING A ROTATING MACHINE OF AN AIRCRAFT

Abstract

The invention relates to a method (1) for monitoring a rotating machine (100) of an aircraft, wherein a measurement signal is acquired from the rotating machine. According to the invention, instantaneous frequencies (f.sub.K(t)) of sinusoidal components of the measurement signal are estimated, and, using a computing module (12), a plurality of successive iterations are carried out in each of which: complex envelopes of the components are updated (C1), parameters of a model of a noise of the signal are updated (C21) on the basis of the envelopes, whether the model has converged from the preceding iteration to the present iteration is tested (C4), with a view to: o if not, carrying out a new iteration, o if so, performing a computation (D) of the complex envelopes on the basis of the iterations that have been carried out.

Claims

1. A method for monitoring at least one rotating machine of an aircraft, the method comprising the following steps: acquiring at least one measurement signal of the rotating machine is by an acquisition module, estimating instantaneous frequencies of sinusoidal components of the at least one measurement signal by an estimation module, executing by a calculation module several successive iterations, each of which comprising: updating complex envelopes of the sinusoidal components, updating parameters of a noise model of the at least one measurement signal from the complex envelopes having been updated, performing a convergence test as to whether the noise model converges from the preceding iteration to the current iteration, to: in the negative, redo a new iteration, in the affirmative, perform a calculation of the complex envelopes from the iterations carried out, updating at each iteration regularisation parameters of the noise model by the calculation module from the complex envelopes having been updated, and, updating, at each iteration, an auxiliary decoupling variable of the complex envelopes of the components between them in the noise model by the calculation module from the complex envelopes having been updated and from the parameters of the noise model having been updated.

2. The method according to claim 1, wherein the noise model of the at least one measurement signal comprises a structured interfering noise, whereof a distribution of probability a priori of spectrum has a heavy-tailed law calculated by the calculation module at each iteration.

3. The method according to claim 1, wherein the noise model of the at least one measurement signal comprises a non-structured interfering noise, whereof a distribution of probability a priori is given by a third distribution of Gaussian type calculated by the calculation module at each iteration.

4. The method according to claim 1, wherein the noise model of the at least one measurement signal resulting from a structured interfering noise model and from a non-structured interfering noise model is a convolution of a Gaussian law and of a heavy-tailed law whose parameters are calculated by the calculation module at each iteration.

5. The method according to claim 1, wherein a conjugated law of probability a priori is allocated to hyperparameters of laws a priori of the parameters of the noise model.

6. The method according to claim 1, wherein a Gamma law of probability a priori is allocated to at least one of the regularisation parameters and results from a Gamma law a posteriori whose parameters are calculated by the calculation module at each iteration.

7. The method according to claim 1, wherein the auxiliary decoupling variable is defined by a model having a distribution of probability of a type of circularly symmetrical normal complex law, which is calculated by the calculation module at each iteration.

8. The method according to claim 1, wherein the convergence test comprises testing by the calculation module at each current iteration that the difference between an absolute value of the complex envelope calculated from the current iteration and the absolute value of the complex envelope calculated from the preceding iteration is less than a prescribed threshold.

9. The method according to claim 1, comprising analysing by the calculation module the complex envelopes having been calculated, and triggering a communication of an alarm by the calculation module on at least one physical outlet as a function of the complex envelopes having been analysed.

10. A monitoring device of at least one rotating machine of an aircraft, the device comprising a module for acquiring at least one measurement signal of the rotating machine, wherein the device comprises a module for estimating instantaneous frequencies of sinusoidal components of the at least one measurement signal, a calculation module configured, at each of several successive iterations carried out, to: update complex envelopes of the sinusoidal components, update parameters of a noise model of the at least one measurement signal from the complex envelopes having been updated, test whether the noise model converges from the preceding iteration to the current iteration, to: in the negative, redo a new iteration, in the affirmative, perform a calculation of the complex envelopes from the iterations carried out, the calculation module being configured to, at each of the successive iterations carried out: update regularisation parameters of the noise model from the complex envelopes having been updated, update an auxiliary decoupling variable of the complex envelopes of the components between them in the noise model, from the complex envelopes having been updated and from the parameters of the noise model having been updated.

11. The monitoring device according to claim 10, wherein the calculation module is configured to analyse the complex envelopes having been calculated, and is able to trigger a communication of an alarm on at least one physical outlet of the device as a function of the complex envelopes having been analysed.

12. The monitoring device according to claim 10, wherein the auxiliary decoupling variable is defined by a model having a distribution of probability of a type of circularly symmetrical normal complex law, which is calculated by the calculation module at each iteration.

13. A computer program comprising code instructions for executing the monitoring method according to claim 1, when it is executed on a computer.

14. An aircraft comprising at least one rotating machine and a monitoring device of the at least one rotating machine according to claim 10.

Description

[0110] The invention will be better understood from the following description given solely by way of non-limiting example in reference to the figures of the appended drawings, in which:

[0111] FIG. 1 schematically illustrates the evolution of the signal-to-noise ratio of a component estimated according to a method of the prior art,

[0112] FIG. 2 illustrates an explanatory diagram of the device and of the monitoring method according to an embodiment of the invention,

[0113] FIG. 3 illustrates an explanatory flowchart of a part of the monitoring method according to an embodiment of the invention,

[0114] FIG. 4a shows the initial spectrogram of a measured signal,

[0115] FIG. 4b shows an example of an envelope having been calculated by the device and the monitoring method according to an embodiment of the invention,

[0116] FIG. 4c shows the spectrogram of a residual signal obtained from FIGS. 4a and 4b by the device and the monitoring method according to an embodiment of the invention,

[0117] FIG. 4d shows the spectrogram of a residual signal obtained by adjusting a method of the prior art,

[0118] FIG. 4e shows the spectrogram of a residual signal obtained by other adjustment of a method of the prior art,

[0119] FIG. 4f shows the spectrogram of a residual signal obtained by even other adjustment of a method of the prior art,

[0120] FIG. 5 illustrates a modular block diagram of the monitoring device according to an embodiment of the invention.

[0121] The monitoring device 1 of one or more rotating machines 100 embedded on an aircraft is described hereinbelow in reference to the figures. The rotating machine 100 can be a rotating propulsion machine of the aircraft, which can be one or more turbomachines, such as for example one or more turbojets or one or more turboprop. The aircraft can be a plane or a helicopter.

[0122] The monitoring device 1 performs the steps of the monitoring method of the rotating machine 100, described hereinbelow.

[0123] In FIGS. 2 and 5, the monitoring device 1 comprises an acquisition module 10 for acquiring one or more measurement signals y(t) for measuring of one or more magnitudes on the rotating machine 100 while this rotating machine 100 is operating as a function of the time t. The measurement signal y(t) can be a measurement signal y(t) of vibrations or acoustics, for example of a rotor and/or a stator, the signal being due to rotation of the rotor relative to the stator. The measurement signal y(t) of vibrations can be measured on a helicopter on a gas generator and/or on a gearbox and/or on a meshing and/or on a bearing. The measurement signal y(t) of vibrations can be measured on a plane, on a stator of a turbomachine. The acquisition module 10 can comprise one or more sensors, embedded on or near the rotating machine 100 on the aircraft, for measuring the signal or signals y (t), which can comprise for example one or more measurement accelerometers of the vibratory signal y(t) and/or one or more measurement microphones of the vibratory signal y(t) acoustic, and/or one or more sensor (encoder or other) of rotation speed of a shaft or of a rotor of the rotating machine 100, or other.

[0124] According to an embodiment, the acquisition module 10 enables to measure and acquire a vibratory or acoustic signal generated by the rotating machine 100 operating. The sensor enables to obtain an analogue signal representative of the vibratory or acoustic signal which is connected to an acquisition chain 13 configured to provide the digital measurement signal y(t)=y (nTe) with T.sub.e the non-zero sampling period, where T=nT.sub.e represents different sampling instants, where n can be for example natural integer numbers. As a consequence, hereinbelow the variable n is a temporal variable. This acquisition chain 13 in this case comprises a conditioner, an anti-aliasing analogue filter, a sampler blocker and an analogue-digital converter.

[0125] The measurement signal y(t) comprises K sinusoidal components of interest x.sub.k=a.sub.k(n)exp(jΦ.sub.k(n)), the sum of which is equal to x(n) according to the equation hereinbelow:

[00007]$\begin{array}{cc}x\left(n\right)=\underset{k=1}{\overset{K}{.Math.}}{a}_k\left(n\right)\exp \left(j{\varphi }_k\left(n\right)\right)& \left[\mathrm{Math}.29\right]\end{array}$

[0126] I={1, . . . , K} designates the set of indices of components of interest of the measurement signal y(t).

[0127] Φ.sub.k(n) designates the phase of each sinusoidal component of interest x.sub.k=a.sub.k(n)exp(jΦ.sub.k(n)).

[0128] a.sub.k (n) designates the complex amplitude (or complex envelope mentioned hereinbelow) of each component of interest a.sub.k(n)exp(jΦ.sub.k(n)) or the temporal derivative of any order of this component of interest a.sub.k (n)exp(jΦ.sub.k (n)).

[0129] In FIGS. 2 and 5, the monitoring device 1 comprises an estimation module 11 for estimating several instantaneous frequencies f.sub.k(t) of the sinusoidal components x.sub.k=a.sub.k(n)exp(jΦ.sub.k(n)) of the measurement signal y(t), where k is a natural integer number varying from 1 to K, and where K is a prescribed natural integer number greater than 1. These instantaneous frequencies f.sub.k(t) are frequencies to be monitored in the rotating machine 100 (for example the first order of the gas generator, of the free turbine, of the shaft of the reducer, the meshing frequencies of the gearbox and of the meshings in the starter, or other).

[0130] According to an embodiment of the invention, a rotation frequency f.sub.ref(t) of a reference shaft of the rotating machine 100 (which can be for example the shaft of its rotor) can have been prescribed in advance in the estimation module 11. The rotation frequency f.sub.ref(t) can be calculated from the speed signal measured by a rotation speed sensor forming part of the acquisition module 10 or from the vibratory or acoustic signal y(t) provided by the acquisition module 10. The estimation module 11 can calculate the instantaneous frequency f.sub.k(t) of each sinusoidal component x.sub.k=a.sub.k(n)exp(jΦ.sub.k(n)) as being proportional to the rotation frequency f.sub.ref(t) according to the equation f.sub.k(t)=c.sub.kf.sub.ref(t), where each proportionality coefficient c.sub.k of the instantaneous frequency f.sub.k(t) is real. Since the connection between the different components of the kinematic chain of the rotating machine 100 is rigid, the proportionality coefficients c.sub.k can have been prescribed in advance in the estimation module 11. In FIG. 2, the elements of the machine 100 which are monitored are specified in advance and their kinematic is prescribed in the estimation module 11, without any need for the user to give the frequencies f.sub.k(t) in entries.

[0131] For example hereinbelow n is a natural integer number ranging from 1 to N, where N is a prescribed natural integer number greater than 1.

[0132] It is assumed that the measurement signal y(n) is equal to the sum of the K components of interest a.sub.k(n)exp(jΦ.sub.k(n)) to which the non-structured interfering signal (or noise) l(n) of Gaussian type (also called wideband) is added, and the structured interfering signal (or noise) e(n) (also called narrowband), according to the following equation:

y(n)=x(n)+e(n)+l(n)

According to an embodiment of the invention, the structured interfering signal (or noise) e(n) is defined as being the part of the interfering noise whereof the Fourier transform is discrete, that is this noise e(n) can be written as the sum of sinusoidal signals. In practice, this noise e(n) is constituted by the residual sinusoidal components not included in the signal x(n). The Fourier transform of the structured interfering noise e(n) is constituted by peaks each linked to a sinusoid. It is said that the spectrum of the structured noise e(n) is characterised by a parsimonious look.

[0133] According to an embodiment of the invention, the non-structured interfering signal (or noise) l(n) is defined as being that part of the interfering noise whereof the Fourier transform is not discrete. The non-structured interfering noise l(n) is constituted by the non-sinusoidal random components whereof the spectrum is not discrete. This non-structured interfering noise l(n) can be for example aerodynamic noise and is modelled by a Gaussian law.

[0134] According to the invention, e(n) incorporates a random signal following a law other than a normal law and must favour a parsimonious spectrum.

[0135] As a reminder and by definition in general, a normal or complex Gaussian law f(x)=NC(μ,σ.sup.2) in a real dimension is defined as the distribution of average μ and standard deviation a or variance σ.sup.2 having the density of probability f(x) according to the following equation:

[00008]$\begin{array}{cc}f\left(x\right)=\frac{1}{\sigma \pi }{e}^{-{\left(\frac{x-\mu }{\sigma }\right)}^2}& \left[\mathrm{Math}.30\right]\end{array}$

[0136] The signal (or structured interfering noise) e(n) represents the residual harmonic signal, and can be noted in the following form:

[00009]$\begin{array}{cc}e\left(n\right)=\underset{k\in \overline{I}}{.Math.}{a}_k^{\prime }\left(n\right)\exp \left(j{\varphi }_k\left(n\right)\right)& \left[\mathrm{Math}.31\right]\end{array}$

[0137] where

[Math. 32]

Ī

[0138] is the set of indices of the sinusoidal components of the structured interfering signal e (n), which are each distinct from the indices of the set I.

[0139] The total interfering signal or total noise of the measurement signal y (t) is equal to the sum b (n) of the non-structured Gaussian interfering signal (or noise) l(n) and of the structured interfering signal (or noise) e(n) according to the following equation:

b(n)=e(n)+l(n)

Calculation of the sinusoidal components of interest a.sub.k(n)exp(jΦ.sub.k(n)) is of Bayesian type. Therefore, n, x(n), e(n) and l(n) are the temporal successive realisations of random variables the laws of which will be defined hereinbelow. Also, x (n), e(n) and l(n) are independent of each other.

[0140] Hereinbelow, the exponent R designates the real part and the exponent I designates the imaginary part (distinct from the above set I).

[0141] The vector y.sup.R is defined as being the real observation vector, formed from the real N values y (n).sup.R of y (n) for the natural integer number n ranging from 1 to N, according to the following equation:

[Math. 33]

y.sup.R[y(1).sup.R, . . . ,y(n).sup.R, . . . ,y(N).sup.R].sup.T

[0142] where the exponent T designates the transpose.

[0143] The vector y.sup.I is defined as being the imaginary vector of observation, formed by the N imaginary values y(n).sup.I of y (n) for n ranging from 1 to N, according to the following equation:

[Math. 34]

y.sup.I=[y(1).sup.I, . . . ,y(n).sup.I, . . . ,y(N).sup.I].sup.T.

[0144] The vectors x.sup.R,x.sup.I, x.sub.k.sup.R, x.sub.k.sup.I, I.sup.R, I.sup.I,e.sup.R,e.sup.I are also formed by the real or imaginary N values of x (n) or l(n) and e (n), for n ranging from 1 to N.

Noise Model

[0145] A law of observation of noise b(n) or noise model b(n) is defined hereinbelow. This noise model b(n) is calculated by a calculation module 12 (or monitoring unit 12) in FIG. 2.

[0146] The calculation module 12 is digital (as opposed to analogue circuits). The calculation module 12 is automatic in the sense where it automatically performs the steps described hereinbelow. As shown in FIG. 5, the calculation module 12 can be or comprise one or more calculators, and/or one or more computers, and/or one or more calculating machines, and/or one or more processors and/or one or more microprocessors. The calculation module 12 and/or the acquisition module 10 can comprise a computer program loaded in their memory to automatically execute the steps described hereinbelow. Also, the monitoring device 1 can comprise one or more physical outlets SP for presentation of information to present to a user the complex envelopes a.sub.k(n) having been calculated by the module 12, wherein this physical outlet or these physical outlets SP can comprise a display screen for displaying these complex envelopes a.sub.k(n) having been calculated (or their derivative of any order (amplitude, phase, amplitude as a function of speed, the estimated component or the estimated components, etc.)). The memory is provided for storing the outlets a.sub.k(n) having been calculated by the calculation module 12.

Structured Noise Model e(n) The structured interfering noise model e(n) can be realised by one of the embodiments described hereinbelow.

[0147] The noise model e(n) is defined as being a structured interfering signal as follows. According to an embodiment of the invention, the spectrum of the structured interfering noise e(n) is defined by a model having a parsimonious spectrum, that is, a distribution concentrated at zero for favouring zero values with a large probability and having heavier tails than the normal law to authorise the presence of high values. For example, a good distribution of probability of the spectrum of structured interfering noise e(n) can be defined by a model being a mix of a finite/infinite number of different weighted distributions.

[0148] According to an embodiment of the invention, the spectrum of narrowband interfering noise is defined by a weighted sum of a distribution of probability of Gaussian type and of another distribution of probability of non-Gaussian type.

[0149] Of the distributions of infinite mixes, there is for example the Laplace law and the Student law currently used for modelling parsimonious signals and which can be expressed as infinite mixes of Gaussians with parameters of mixes following gamma laws and inverse gamma laws respectively. Of the distributions of finite mixes, there is for example the mix of two distributions: of a Bernoulli distribution and of a Gaussian distribution known as Bernoulli-Gaussian. In the case of the Bernouilli model, the spectrum of the structured interfering noise e(n) is defined by a model having a distribution of probability given by the barycentre between a distribution of Gaussian type NC(0, 2ζ.sup.−1) weighted by a weight (probability) β of the presence of interfering noise e(n) and a Dirac distribution δ.sub.0 weighted by a weight 1−β, for example as per the following equation:

[Math. 35]

∀n∈{1, . . . ,N}[F.sub.Ne].sub.n|β,ζ˜(1−β)δ.sub.0+βNC(0,2ζ.sup.−1)

[0150] where

[0151] F.sub.N designates the discrete Fourier transform operator of size N, this operator being divided by

[Math. 36]

√{square root over (N)},

[0152] δ.sub.0 is the Dirac distribution positioned at the zero abscissa,

[0153] 0≤β≤1, is the probability of the presence of the structured noise,

[0154] NC(0, 2ζ.sup.−1) designates a normal complex distribution of zero average, whereof the variance is equal to 2ζ.sup.−1. Therefore, the real part of this normal complex distribution NC(0, 2ζ.sup.−1) is a Gaussian distribution of zero average and variance ζ.sup.−1. The imaginary part of this normal complex distribution NC(0, 2ζ.sup.−1) is a Gaussian distribution of zero average and variance ζ.sup.−1.

[0155] According to an embodiment, β is different from zero.

[0156] Hereinabove, [F.sub.Ne].sub.n designates the discrete Fourier transform operator of size N, which is divided by

[Math. 37]

√{square root over (N)}

[0157] and is applied to e (n). Therefore, there is

[00010]$\begin{array}{cc}{F}_N=\frac{1}{\sqrt{N}}{\mathrm{DFT}}_N& \left[\mathrm{Math}.38\right]\end{array}$$\mathrm{And}$$\begin{array}{cc}{F}_N^{\top }{F}_N=F_N{F}_N^{\top }=\mathrm{Id}& \left[\mathrm{Math}.39\right]\end{array}$

[0158] where Id designates the matrix identity.

[0159] Unlike the known Vold-Kalman method, where the noise is assumed as a Gaussian centre, this embodiment easily takes structured noise into account in the model to improve the efficacy of the estimation at the time of coupling without need to estimate the interfering components.

[0160] In variable mode, the spectrum of e(n) can have wider peaks. Fourier transforms S.sub.j, limited to a temporal segment of size N.sub.j for e(n) are defined as follows:

[Math. 40]

∀j∈{1, . . . ,J} ∀n∈{1, . . . ,N.sub.j}[S.sub.je].sub.n|β.sub.j,ζ.sub.j˜(1−β.sub.j)δ.sub.0+β.sub.jNC(0,2ζ.sub.j.sup.−1)

[0161] with

[Math. 41]

S.sub.j=F.sub.N.sub.jP.sub.j

[0162] where P.sub.j is a selection matrix of size

[Math. 42]

N.sub.j×N,

[0163] N.sub.j is a prescribed natural integer number less than N (N.sub.j is the length of the window),

[00011]$\begin{array}{cc}{F}_{N_j}=\frac{1}{\sqrt{N_j}}{\mathrm{DFT}}_{N_j}& \left[\mathrm{Math}.43\right]\end{array}$$\mathrm{with}$$\begin{array}{cc}{\mathrm{DFT}}_{N_j}& \left[\mathrm{Math}.44\right]\end{array}$

[0164] the discrete Fourier transform of size N.sub.j,

[0165] n is a natural integer number ranging from 1 to N.sub.j,

[0166] J is a prescribed natural integer number (number of temporal windows j), j is a natural integer number ranging from 1 to J,

[0167] NC(0, 2ζ.sup.−1) designates a normal complex distribution of zero average, whereof the variance is equal to 2ζ.sup.−1,

[0168] β.sub.j is the weight of the presence of structured noise, allocated to NC(0, 2ζ.sup.−1),

[0169] 0≤β.sub.j≤1,

[0170] 1−β.sub.j is the weight allocated to the Dirac distribution positioned at the zero abscissa. Fourier transforms S.sub.j are therefore of shorter duration than F.sub.N.

[0171] If, also, the matrices P.sub.j have disjointed supports, this gives

[Math. 45]

P.sub.jP.sub.j.sup.T=Id

[0172] and therefore

[Math. 46]

S.sub.jS.sub.j.sup.T=Id

[0173] The matrices P.sub.j have disjointed supports if the sum of all the N, is equal to N:

[00012]$\begin{array}{cc}\underset{j=1}{\overset{J}{.Math.}}{N}_j=N& \left[\mathrm{Math}.47\right]\end{array}$

[0174] The selection matrix P.sub.j is such that it selects N.sub.j values belonging to the window j of a vector of size N.

[0175] P.sub.j is a matrix containing N.sub.j lines and N columns, it is defined by:

[0176] For j=1:

[Math. 48]

∀i.sub.0∈{1, . . . ,N.sub.j}i.sub.1∈{1, . . . ,N}P.sub.j(i.sub.0,i.sub.1)=1 if i.sub.1=i.sub.0 if not P.sub.j(i.sub.0,i.sub.1)=0.

[0177] For j>1

[Math. 49]

∀i.sub.0∈{1, . . . ,N.sub.j}i.sub.1∈{1, . . . ,N}P.sub.j∈(i.sub.0,i.sub.1)=1 if i.sub.1=i.sub.0+Σ.sub.j.sub.0.sub.=0.sup.j−1N.sub.j.sub.0 if not P.sub.j(i.sub.0,i.sub.1)=0.

[0178] Unlike the known Vold-Kalman method where noise is assumed to be stationary, which is rarely the case of aerodynamic noise, this embodiment takes into account the non-stationary aspect of noise, and does this by analysing the signal via short-term Fourier transforms S.sub.j mentioned hereinabove. The time axis is divided into J intervals or windows. The parameters of the model can change from one window to the other. Given that the parameters of these laws are unknown (the level of aerodynamic noise is unknown and the coupling instants between the components and the amplitude of the interfering components are unknown), they are also going to be adapted automatically. Laws a priori are selected to model information a priori on these parameters. In general, if there is no information a priori laws called non-informative can be chosen.

Non-Structured Noise Model l(n)

[0179] The non-structured interfering noise model l(n) can be created by one of the embodiments described hereinbelow.

[0180] According to an embodiment of the invention, the spectrum of the non-structured interfering noise l(n) is defined by a model having a distribution of probability of Gaussian type.

[0181] According to an embodiment of the invention, the spectrum of the non-structured interfering noise l(n) is defined by a model having a distribution of normal complex probability NC(0, 2 σ.sup.−1) of zero average and of variance 2σ.sup.−1, for example according to the following equation:

[Math. 50]

∀n∈{1, . . . ,N}[F.sub.Nl].sub.n|σ˜NC(0,2 σ.sup.−1)

[0182] It can be that the background noise is non-stationary. In this case, the procedure can be via temporal segments j as was done for the interfering signal e(n). Fourier transforms S.sub.j, limited to a temporal segment of size N, for l(n), are defined as follows:

[Math. 51]

∀j∈{1, . . . ,J}∀n∈{1, . . . ,N.sub.j}[S.sub.jl].sub.n|β.sub.j,ζ.sub.j˜NC(0,2τ.sub.j.sup.−1)

[0183] where j, J, n, N.sub.j, S.sub.j, β.sub.j, ζ.sub.j have the definitions mentioned hereinabove,

[0184] NC(0, 2τ.sub.j.sup.−1) designates a normal complex distribution of zero average, whereof the variance is equal to 2τ.sub.j.sup.−1.

[0185] In each temporal window j, the non-structured interfering noise l(n) is Gaussian and has a variance 2τ.sub.j.sup.−1. This is shown as

[Math. 52]

[P.sub.jl].sub.n|τ.sub.j˜NC(0,2τ.sub.j.sup.−1)

Total noise model b(n)

[0186] The total interfering noise model b(n) can be created by one of the embodiments described hereinbelow.

[0187] According to an embodiment of the invention, the spectrum of the interfering noise b(n) results from the models of the structured interfering noise and from the non-structured noise model. In fact, since e(n) and l(n) are independent, the density of probability of the sum b (n)=e (n)+l(n) is the convolution of the density of probability of the structured interfering noise e(n) with that of the non-structured noise l(n).

[0188] According to an embodiment of the invention, for example, in the case of a Bernoulli-Gaussian law for the structured noise, the spectrum of the total interfering noise b(n) is defined by a model having a mix of two distributions of Gaussian type each having an associated weight and different variances.

[0189] According to an embodiment of the invention, the spectrum of the total interfering noise b(n) is defined by a model having a distribution given by the barycentre between the distribution NC(0, 2τ.sub.j.sup.−1 of Gaussian type weighted by the weight 1−β.sub.j and a distribution NC(0, 2τ.sub.j.sup.−1+τ.sub.j.sup.−1)) of Gaussian type weighted by the weighting β.sub.j, for example according to the following equation:

[Math. 53]

∀j∈{1, . . . ,J},∀n∈{1, . . . ,N.sub.j}[S.sub.jb].sub.n|τ.sub.j,ζ.sub.j,β.sub.j ˜(1−β.sub.j)NC(0,2τ.sub.j.sup.−1)+β.sub.jNC(0,2(ζ.sub.j.sup.−1+τ.sub.j.sup.−1))

[0190] where j, J, n, N.sub.j, S.sub.j, β.sub.j, ζ.sub.j, τ.sub.j have the definitions mentioned hereinabove,

[0191] NC(0, 2(ζ.sub.j.sup.−1+τ.sub.j.sup.−1)) designates a normal complex distribution of zero average, whereof the variance is equal to 2 (ζ.sub.j.sup.−1+τ.sub.j.sup.−1).

[0192] As follows,

[Math. 54]

η.sub.j=1/(ζ.sub.j.sup.−1+τ.sub.j.sup.−1)

[0193] Therefore, in this embodiment, the spectrum of the total interfering noise b(n) is defined by a model having a distribution given by the barycentre between the distribution NC(0, 2τ.sub.j.sup.−1) of Gaussian type weighted by the weight 1−β.sub.j and a distribution NC(0,2η.sub.j.sup.−1) of Gaussian type weighted by the weight) β.sub.j, for example according to the following equation:

[Math. 55]

∀j∈{1, . . . ,J}∀n∈{1, . . . ,N.sub.j}[S.sub.jb].sub.n|τ.sub.j,ζ.sub.j,β.sub.j˜(1−β.sub.j)NC(0,2τ.sub.j.sup.−1)+β.sub.jNC(0,2η.sub.j.sup.−1)

[0194] where j, J, n, N.sub.j, S.sub.j, β.sub.j, ζ.sub.j, τ.sub.j, η.sub.j have the definitions mentioned hereinabove,

[0195] NC(0,2η.sub.j.sup.−1) designates a normal complex distribution of zero average, whereof the variance is equal to 2.sub.η.sub.j.sup.−1,

[0196] with in η.sub.j<τ.sub.j.

[0197] τ.sub.j and η.sub.j are also called noise precision parameters.

[0198] The law of observation y.sup.R, y.sup.I of the measurement signal y(t) is given by the Gaussian N according to the following hierarchical model:

[00013]$\begin{array}{cc}\{\begin{array}{c}\begin{array}{c}\forall j\in \left\{1,.Math.,J\right\}\forall n\in \left\{1,.Math.,N\right\}{d}_j\left(n\right).Math.{\tau }_j,\\ {\eta }_j,{\beta }_j~\left(1-\beta \right){\delta }_{{\tau }_j^{-1}}+{\beta }_j{\delta }_{{\eta }_j^{-1}}\end{array}\\ y^R.Math.a,d~N\left({\left[\mathrm{Ca}\right]}^R,2\underset{j=1}{\overset{J}{.Math.}}{S}_j^{\top }\mathrm{diag}\left(d_j\right)S_j\right)\\ y^I.Math.a,d~N\left({\left[\mathrm{Ca}\right]}^I,2\underset{j=1}{\overset{J}{.Math.}}{S}_j^{\top }\mathrm{diag}\left(d_j\right)S_j\right)\end{array}& \left[\mathrm{Math}.56\right]\end{array}$

[0199] where

[0200] the matrix C is the diagonal matrix formed by the terms exp(jΦ.sub.k(n)) for n ranging from 1 to N and k ranging from 1 to K,

[0201] a is the vector of the envelopes,

[0202] diag(d.sub.j) is a matrix the diagonal of which is formed by the coefficients d.sub.j and whereof the other coefficients are zero,

[0203] d.sub.j is a variable which follows a distribution given by the barycentre between a Dirac distribution

[Math. 57]

δ.sub.τ.sub.j.sup.9−1

[0204] positioned in the abscissa τ.sub.j.sup.−1 and weighted by the weight 1−β.sub.j and a Dirac distribution

[Math. 58]

δ.sub.η.sub.j.sup.−1

[0205] positioned in the abscissa η.sub.j.sup.−1 and weighted by the weight β.sub.j.

[0206] According to an embodiment of the invention, the law a posteriori of d.sub.j is a distribution given by the barycentre between a Dirac distribution

[Math. 59]

δ.sub.τ.sub.j.sup.−1

positioned in the abscissa τ.sub.j.sup.−1 and weighted by the weight 1-p.sub.j(n) and a Dirac distribution

[Math. 60]

δ.sub.η.sub.j.sup.−1

[0207] positioned in the abscissa η.sub.j.sup.−1 and weighted by the weight p.sub.j(n), for example according to the following equation: P [Math. 61]

∀j∈{1, . . . ,J}

[Math. 62]

∀n∈{1, . . . ,N.sub.j}

[Math. 63]

d.sub.j(n)|τ.sub.j,η.sub.j,β.sub.ja˜(1−p.sub.j(n))δ.sub.τ.sub.j.sup.−1+p.sub.j(n)δ.sub.η.sub.j.sup.−1+p.sub.j(n)δ.sub.η.sub.j.sup.−1

[0208] In the above,

[00014]$\begin{array}{cc}{p}_j\left(n\right)=\frac{w_j\left(n\right)}{1+w_j\left(n\right)}& \left[\mathrm{Math}.64\right]\end{array}$$\begin{array}{cc}{w}_j\left(n\right)=\frac{{\beta }_j}{1-{\beta }_j}\exp \left(-\frac{\left({\eta }_j-{\tau }_j\right){\eta }_j}{2{\tau }_j}{.Math.{\left[S_j\left(\mathrm{Ca}-y\right)\right]}_n.Math.}^2\right)\mathrm{if}n\in \left\{1,N_j\right\}& \left[\mathrm{Math}.65\right]\end{array}$$\begin{array}{cc}{w}_j\left(n\right)=\frac{{\beta }_j}{1-{\beta }_j}\exp \left(-\frac{\left({\eta }_j-{\tau }_j\right){\eta }_j^2}{2{\tau }_j^2}{.Math.{\left[S_j\left(\mathrm{Ca}-y\right)\right]}_n.Math.}^2\right)\mathrm{if}\mathrm{not}& \left[\mathrm{Math}.66\right]\end{array}$

[0209] Also, diag(d.sub.j) or diag(d.sub.j (n)) is a diagonal covariance matrix of the noise.

Iteration of Calculation of the Complex Envelopes a.sub.k(n) of the Sinusoidal Components

[0210] According to the invention, the monitoring device 1 comprises a calculation module 12 configured, at each respective iteration carried out, to: [0211] during a first step C1 of the respective iteration carried out update the laws (calculate or give a value to the parameters (average and variance) of the law) of the complex envelopes a.sub.k(n) of the sinusoidal components, [0212] during a second step C21 of the respective iteration carried out (after the first step C1 of the respective iteration carried out) update the parameters of the noise model b(n) and/or e(n) and/or l(n) of the measurement signal (y(t)) from the updated complex envelopes a.sub.k(n) or given the updated complex envelopes a.sub.k(n), [0213] during a third step C4 of the respective iteration carried out (after the second step C21 of the respective iteration carried out) test whether the model converges from the preceding iteration (t−1) to the current iteration (t), to: [0214] in the negative (NO in FIG. 2), redo a new iteration of the steps C1, C21, C4 described hereinabove, [0215] and in the affirmative (YES in FIG. 2), during the fourth step D following the third step C4, perform a calculation of the complex envelopes a.sub.k(n) from the iterations carried out, to obtain automatic extraction of the sinusoidal components a.sub.k(n).

[0216] The updating of the noise parameters can be done independently on the j ranging from 1 to J. At each iteration t, given the laws a posteriori of the variables at the preceding iteration t−1 and the observation y, laws a posteriori of the different variables are updated by the calculation module 12, for example according to the diagram shown in FIG. 3.

[0217] The third test step C4 enables to verify whether the laws of the different variables of the problem to be estimated at a given iteration are definitely optimal. According to an embodiment of the invention, the test of the third step C4 comprises automatic calculation of one or more indicators by the calculation module 12 from the envelopes a.sub.k(n) having been calculated (for example their absolute value).

[0218] According to an embodiment of the invention, the test of the third step C4 can be performed by the calculation module 12 which compares the absolute value of the difference between the average of the law of the envelope a.sub.k(n) calculated from the current iteration and the average of the law of the envelope a.sub.k(n) calculated from the preceding iteration to a threshold prescribed in advance, to determine whether this difference is less than this threshold (affirmative convergence hereinabove (YES in FIG. 2) when this difference is less than the threshold; negative convergence hereinabove (NO in FIG. 2) when this difference is not less than the threshold). For example, this threshold can be prescribed between 10.sup.−3 and 10.sup.−6 times the norm of the average of the law of the envelope a.sub.k (n) calculated from the preceding or current iteration. The test of the third step C4 can also be performed by the calculation module 12 on one or more of the described variables of the noise model, such as for example an average, a variance or other, for example of the calculated envelope a.sub.k (n), to test if this variable converges from the preceding iteration to the current iteration. As long as the statistics are not stabilised (NO in FIG. 2), the algorithm continues to update the variables iteratively.

[0219] According to an embodiment of the invention, during the fourth step D, the calculation of the complex envelopes a.sub.k (n) can be performed by the calculation module 12 from the average of the estimated law. Additional statistics can also be given as the confidence interval.

[0220] According to an embodiment of the invention, a maximal number of iterations is prescribed in advance in the calculation module 12.

Parameters of the Noise Model

[0221] According to an embodiment of the invention, these parameters can be one or more of the parameters configuring the Gaussian law of the non-structured noise (average, precision) and those configuring the heavy-tailed law of the structured noise (average, precision, form parameter etc) and this for one or more windows of the J windows.

[0222] According to an embodiment of the invention, in the case of a Bernoulli-Gaussian model these parameters can be one or more or all of the τ.sub.j for j ranging from 1 to J, and/or one or more or all of the β.sub.j for j ranging from 1 to J, and/or one or more or all of the η.sub.j for j ranging from 1 to J, and/or one or more or all of the ζ.sub.j for j ranging from 1 to J.

[0223] This embodiment of the parameters can be realised by one or more of the other embodiments described hereinbelow.

[0224] According to an embodiment of the invention, a law of probability a priori is predefined for these parameters.

[0225] According to an embodiment of the invention, conjugated laws are selected. This facilitates calculations by the calculation module 12.

[0226] According to an embodiment of the invention, in the case of a Bernoulli-Gaussian model for structured noise, a law of uniform probability over the interval ranging from 0 to 1 is allocated to one or more or all the β.sub.j, that is

[Math. 67]

β.sub.j˜U(0,1)∀j∈{1, . . . ,J}

[0227] According to an embodiment of the invention, a law of probability of Beta law type is allocated to one or more or all the β.sub.j.

[0228] By way of reminder, by convention and in general, a Beta law (noted Beta(α, H)) having a third adjustment α (strictly positive real) and a fourth adjustment H (strictly positive real) is defined by the function of density of probability f(x; α, H) according to the equation hereinbelow:

[00015]$\begin{array}{cc}f\left(x;\alpha ,H\right)=\{\begin{array}{c}\frac{{x^{\alpha -1}\left(1-x\right)}^{H-1}}{{\int }_0^1{u^{\alpha -1}\left(1-u\right)}^{H-1}\mathrm{du}}\mathrm{for}x\in \left[0,1\right]\\ 0\mathrm{if}\mathrm{not}\end{array}& \left[\mathrm{Math}.68\right]\end{array}$

[0229] According to an embodiment of the invention, the law a posteriori of one or more or all the β.sub.j resulting from a uniform law a priori and given the observations y(n) is a law Beta(2n.sub.1+1, 2n.sub.2+1) according to the following equation:

[00016]$\begin{array}{cc}\forall j\in \left\{1,.Math.,J\right\}& \left[\mathrm{Math}.69\right]\end{array}$$\begin{array}{cc}{\beta }_j.Math.y,a,d,{\tau }_j,{\eta }_j~\mathrm{Beta}\left(2.n_1+1,2.n_2+1\right)& \left[\mathrm{Math}.70\right]\end{array}$$\mathrm{with}$$\begin{array}{cc}{n}_1=.Math.d_j\left(1\right)={\eta }_j^{-1}.Math.+.Math.d_j\left(N\right)={\eta }_j^{-1}.Math.+\underset{n=2}{\overset{N_j-1}{.Math.}}.Math.d_j\left(n\right)={\eta }_j^{-1}.Math.& \left[\mathrm{Math}.71\right]\end{array}$$\begin{array}{cc}{n}_2=.Math.d_j\left(1\right)={\tau }_j^{-1}.Math.+.Math.d_j\left(N\right)={\tau }_j^{-1}.Math.+2\underset{n=2}{\overset{N_j-1}{.Math.}}.Math.d_j\left(n\right)={\tau }_j^{-1}.Math.& \left[\mathrm{Math}.72\right]\end{array}$

[0230] where the function < > used hereinabove means in general <P>=1 if P is true, or if not <P>=0.

[0231] This law of probability Beta(2n.sub.1+1,2n.sub.2+1) is the law a posteriori which result from the law a priori (uniform law) and from the law of observations (law of y).

[0232] According to an embodiment of the invention, a law of Gamma probability is allocated to one or more or all the τ.sub.j, for example

[Math. 73]

τ.sub.j˜Gamma(a.sub.τ.sub.j,b.sub.τ.sub.j)

[0233] having

[Math. 74]

a.sub.τ.sub.j

[0234] as first adjustment and

[Math. 75]

b.sub.τ.sub.j

[0235] as second adjustment.

[0236] By way of reminder, by convention and in general, a Gamma law (noted Gamma(H, α)) having a first adjustment H (strictly positive real) and a second adjustment α (strictly positive real) is defined by the function of density of probability f(x; H, α) as per the equation hereinbelow:

[00017]$\begin{array}{cc}f\left(x;H,\alpha \right)=\frac{x^{H-1}{e}^{-\alpha x}}{\Gamma \left(H\right){\alpha }^H}& \left[\mathrm{Math}.76\right]\end{array}$

[0237] where Γ(H) is the gamma Euler function according to the equation hereinbelow

[Math. 77]

Γ(H)=∫.sub.0.sup.+∞t.sup.H−1e.sup.−tdt

[0238] According to an embodiment of the invention, a law of Gamma probability is allocated to one or more or all the η.sub.j, for example

[Math. 78]

η.sub.j˜Gamma(a.sub.η.sub.j,b.sub.η.sub.j)

having

[Math. 79]

a.sub.η.sub.j,

[0239] as first adjustment and

[Math. 80]

b.sub.η.sub.j,

[0240] as second adjustment.

[0241] The adjustments

[Math. 81]

a.sub.τ.sub.j,

[Math. 82]

b.sub.τ.sub.j,

[Math. 83]

a.sub.η.sub.j,

[Math. 84]

b.sub.η.sub.j,

[0242] of Gamma laws are each strictly positive.

[0243] According to an embodiment of the invention,

[Math. 85]

a.sub.τ.sub.j,

[Math. 86]

b.sub.τ.sub.j,

[Math. 87]

a.sub.η.sub.j,

[Math. 88]

b.sub.η.sub.j.

[0244] of the Gamma laws are each prescribed in advance in the calculation module 12. Therefore,

[Math. 89]

a.sub.τ.sub.j,

[Math. 90]

b.sub.τ.sub.j,

[Math. 91]

a.sub.η.sub.j,

[Math. 92]

b.sub.η.sub.j.

[0245] of the Gamma laws are each a hyper-parameter quantifying information a priori on the first parameters to be estimated.

[0246] According to an embodiment of the invention, each couple (a, b) of first adjustment and second adjustment, such as for example

[Math. 93]

(a,b)∈{(a.sub.η.sub.j,b.sub.η.sub.j)).sub.j∈{1, . . . ,J},(a.sub.τ.sub.j,b.sub.τ.sub.j).sub.j∈{1, . . . ,J}},

[0247] is prescribed in advance in the calculation module 12 according to information a priori on the distribution of the first parameter B considered with

[Math. 94]

θ∈{η.sub.j∈{1, . . . ,J},τ.sub.j∈{1, . . . ,J}}.

[0248] These adjustments can be prescribed in advance by the user. If the user does not prescribe these adjustments, these adjustments are fixed to very low values close or equal to zero (see the example hereinbelow) to have non-informative laws. This information can be an average of the law considered for θ and/or a variance of the law considered for θ. According to an embodiment of the invention, if the value of θ noted by θ.sub.e>0 is known approximately with uncertainty equal to εθ.sub.e where ε>0, then

[00018]$\begin{array}{cc}a=\frac{{\theta }_e}{\epsilon }\mathrm{and}& \left[\mathrm{Math}.95\right]\end{array}$$\begin{array}{cc}b=\frac{1}{\epsilon }& \left[\mathrm{Math}.96\right]\end{array}$

[0249] are prescribed in advance. If there is no information a priori on the parameter θ, alors a=b=0 are prescribed in advance (the law is non-informative). The presence of information a priori on the parameters most often accelerates convergence.

[0250] According to an embodiment of the invention, the law a posteriori of one or more or all the η.sub.j resulting from its Gamma law a priori of parameters

(

[Math. 97]

a.sub.η.sub.j,

[Math. 98]

b.sub.η.sub.j,

[0251] and the measurement signal y(n) is a law

[Math. 99]

Gamma(ã.sub.η.sub.j,{tilde over (b)}.sub.η.sub.j)

[0252] having

[Math. 100]

ã.sub.η.sub.j

[0253] as first adjustment and

[Math. 101]

{tilde over (b)}.sub.η.sub.j

[0254] as second adjustment according to the following equation:

[Math. 102]

∀j∈{1, . . . ,J}

[Math.103]

η.sub.j|y,a—Gamma(ã.sub.η.sub.j,{tilde over (b)}.sub.η.sub.j)

[0255] According to an embodiment of the invention, the law a posteriori of one or more or all the τ.sub.j resulting from its Gamma law a priori of parameters (

[Math.104]

a.sub.τ.sub.j,

[Math. 105]

b.sub.τ.sub.j)

[0256] and the measurement signal y(n) is a Gamma law

[Math.106]

Gamma(ã.sub.τ.sub.j,{tilde over (b)}.sub.τ.sub.j)

[0257] having

[Math. 107]

ã.sub.τ.sub.j

[0258] as first adjustment and

[Math. 108]

{tilde over (b)}.sub.τ.sub.j,

[0259] as second adjustment according to the following equation:

[00019]$\begin{array}{cc}\forall j\in \left\{1,.Math.,J\right\}& \left[\mathrm{Math}.109\right]\end{array}$$\begin{array}{cc}{\tau }_j.Math.y,a\sim \mathrm{Gamma}\left({\tilde{a}}_{{\tau }_j},{\tilde{b}}_{{\tau }_j}\right)\mathrm{In}\mathrm{the}\mathrm{above},& \left[\mathrm{Math}.110\right]\end{array}$$\begin{array}{cc}{\tilde{a}}_{{\eta }_j}=a_{{\eta }_j}+n_1& \left[\mathrm{Math}.111\right]\end{array}$$\begin{array}{cc}{\tilde{b}}_{{\eta }_j}=b_{{\eta }_j}+\frac{{\eta }_j}{2}\underset{\underset{d_j\left(n\right)={\eta }_j^{-1}}{n}}{.Math.}{.Math.{\left[S_j\left(\mathrm{Ca}-y\right)\right]}_n.Math.}^2& \left[\mathrm{Math}.112\right]\end{array}$$\begin{array}{cc}{\tilde{a}}_{{\tau }_j}=a_{{\tau }_j}+n_2& \left[\mathrm{Math}.113\right]\end{array}$$\begin{array}{cc}{\tilde{b}}_{{\tau }_j}=b_{{\tau }_j}+\frac{{\tau }_j}{2}\underset{\underset{d_j\left(n\right)={\tau }_j^{-1}}{n}}{.Math.}.Math.{\left[S_j\left(\mathrm{Ca}-y\right)\right]}_n{.Math.}^2& \left[\mathrm{Math}.114\right]\end{array}$

Regularisation Parameters γ.SUB.k

[0260] According to an embodiment of the invention, the calculation module 12 is configured, during a fifth step C22 of the respective iteration carried out (following the second step C21 of the respective iteration carried out and prior to the third step C4) to update second regularisation parameters γ.sub.k of the model from the updated complex envelopes a.sub.k(n) or given the updated complex envelopes a.sub.k(n). These second regularisation parameters γ.sub.k are added restrictions imposed to favour regularity of the complex envelopes a.sub.k(n). These second regularisation parameters γ.sub.k characterise suppression of singularities in the components a.sub.k(n).

[0261] This embodiment of the second regularisation parameters γ.sub.k can be realised by one of the other embodiments described hereinbelow.

[0262] Updating of the second parameters can be done independently on the k ranging from 1 to K.

[0263] According to an embodiment of the invention, the density p of probability (a priori) of the complex envelopes a.sub.k(n) in the presence of these second regularisation parameters γ.sub.k is the following, for k ranging from 1 to K:

[00020]$\begin{array}{cc}p\left(a^R,a^I.Math.{\gamma }_1,.Math.,{\gamma }_K\right)=\underset{k=1}{\overset{K}{.Math.}}p\left(a_k^R.Math.{\gamma }_k\right)p\left(a_k^I.Math.{\gamma }_k\right)=C_a\underset{k=1}{\overset{K}{.Math.}}{\gamma }_k^N\exp \left(-{\gamma }_k{.Math.A_p{a}_k.Math.}^2\right)& \left[\mathrm{Math}.115\right]\end{array}$

where C.sub.a is a constant independent from the γ.sub.k and from the complex envelopes a.sub.k(n).

[0264] According to an embodiment of the invention, a law of Gamma probability is allocated to one or more or all the second regularisation parameters γ.sub.k, for example

[Math. 116]

γ.sub.k˜Gamma(a.sub.β.sub.k,b.sub.γ.sub.k)

[0265] having

[Math. 117]

a.sub.γ.sub.k

[0266] as first adjustment and

[Math. 118]

b.sub.γ.sub.k

[0267] as second adjustment. According to an embodiment of the invention, the couple (

[Math. 119]

a.sub.γ.sub.k,

[Math. 120]

b.sub.γ.sub.k)

is prescribed in advance in the calculation module 12 according to information a priori on the distribution of the first parameter γ.sub.k. This information can be an average of the law considered for γ.sub.k and/or a variance of the law considered for γ.sub.k. These adjustments can be prescribed in advance by the user. If the user does not prescribe these adjustments, these adjustments are fixed to very low values close or equal to zero (see the example hereinbelow) to have non-informative laws. According to an embodiment of the invention, if the value of γ.sub.k noted by θ.sub.e>0 is known approximately with uncertainty equal to εθ.sub.e where ε>0, then

[00021]$\begin{array}{cc}{a}_{{\gamma }_k}=\frac{{\theta }_e}{\epsilon }\mathrm{and}& \left[\mathrm{Math}.121\right]\end{array}$$\begin{array}{cc}{b}_{{\gamma }_k}=\frac{1}{\epsilon }& \left[\mathrm{Math}.122\right]\end{array}$

[0268] are prescribed in advance. If there is no information a priori on the parameter θ, then

[Math. 123]

a.sub.γ.sub.k=b.sub.γ.sub.k=0

[0269] are prescribed in advance (the law is non-informative). The presence of information a priori on the parameters most often accelerates convergence.

[0270] According to an embodiment of the invention, a law of probability a posteriori of γ.sub.k (resulting from the Gamma law a priori of adjustment parameters (

[Math. 124]

a.sub.γ.sub.k

[Math. 125]

b.sub.γ.sub.k),

[0271] and from the law a priori of envelopes) is also a law Gamma

[Math. 126]

Gamma(ã.sub.γ.sub.k,{tilde over (b)}.sub.γ.sub.k)

[0272] having

[Math. 127]

ã.sub.γ.sub.k

[0273] as first adjustment and

[Math. 128]

{tilde over (b)}.sub.γ.sub.k

[0274] as second adjustment, for example according to the following equation:

[Math. 129]

[00022]$\forall k\in \left\{1,.Math.,K\right\}$$\begin{array}{cc}{\gamma }_k.Math.a_k\sim \mathrm{Gamma}\left({\tilde{a}}_{{\gamma }_k},{\tilde{b}}_{{\gamma }_k}\right)\mathrm{In}\mathrm{the}\mathrm{above},& \left[\mathrm{Math}.130\right]\end{array}$$\begin{array}{cc}{\tilde{a}}_{{\gamma }_k}=a_{{\gamma }_k}+N& \left[\mathrm{Math}.131\right]\end{array}$$\begin{array}{cc}{\tilde{b}}_{{\gamma }_k}=b_{{\gamma }_k}+\frac{1}{2}{.Math.A_p{a}_k.Math.}^2& \left[\mathrm{Math}.132\right]\end{array}$

[0275] According to an embodiment of the invention, a law a posteriori of the sinusoidal components of interest a.sub.k(n)exp(jΦ.sub.k(n)) is calculated during the first step C1 by the calculation module 12 from of the law of observation of the noise b (n) and from the law of observation of the measurement signal y(t). According to an embodiment of the invention, the law a posteriori of the real part a.sup.R of the sinusoidal components of interest a.sub.k(n)exp(jΦ.sub.k(n)) is defined by a model having a distribution of probability of Gaussian type N(m.sub.x.sup.R,Σ.sub.x) of average m.sub.x.sup.R (real part of m.sub.x) and of variance Σ.sub.x, according to the following equation:

[Math. 133]

a.sup.r|y,d,γ.sub.1, . . . ,γ.sub.k˜n(m.sub.x.sup.R,Σ.sub.x)

[0276] According to an embodiment of the invention, the law a posteriori of the imaginary part a.sup.I of the sinusoidal components of interest a.sub.k(n)exp(jΦ.sub.k(n)) is defined during the first step C1 by a model having a distribution of probability of Gaussian type N(m.sub.x.sup.I,Σ.sub.x) of average m.sub.x.sup.I (imaginary part of m.sub.x) and of variance Σ.sub.x, for example according to the following equation:

[Math. 134]

a.sup.I|y,d,γ.sub.1, . . . ,γ.sub.K˜N(m.sub.x.sup.I,Σ.sub.x)

[0277] The real part a.sup.R of the sinusoidal components of interest a.sub.k(n)exp(jΦ.sub.k(n)) and the imaginary part a.sup.I of the sinusoidal components of interest a.sub.k(n)exp(jΦP.sub.k(n)) can be updated independently during the first step C1. In the case of a very large number of components to be estimated, this first step C1 can be done in parallel because of multi-processor architecture of the module 12.

[0278] In the above,

[00023]$\begin{array}{cc}{m}_x={.Math.}_x{C}^{\top }\underset{j=1}{\overset{J}{.Math.}}{S}_j^{\top }{\mathrm{diag}\left(d_j\right)}^{-1}{S}_{j}y& \left[\mathrm{Math}.135\right]\end{array}$$\begin{array}{cc}{.Math.}_x={\left(C^{\top }\underset{j=1}{\overset{J}{.Math.}}{S}_j^{\top }{\mathrm{diag}\left(d_j\right)}^{-1}{S}_{j}C+A^T\mathrm{diag}\left({\gamma }_k\right)A\right)}^{-1}& \left[\mathrm{Math}.136\right]\end{array}$

[0279] diag(d.sub.j) is a matrix whereof the diagonal is formed by the coefficients d.sub.j and whereof the other coefficients are zero,

[0280] diag(γ.sub.k) is a matrix whereof the diagonal is formed by the coefficients γ.sub.k and whereof the other coefficients are zero,

[0281] A.sub.p is the matrix of the discrete differences ∇.sup.p of order p of diag(γ.sub.k),

[0282] A is the diagonal block matrix formed by the K matrices A.sub.p and is the matrix containing K blocks of A.sub.p.

[0283] the matrix C is the diagonal matrix formed by the terms exp(jΦ.sub.k(n)) for n ranging from 1 to N.

[0284] Auxiliary Decoupling Variable v of the Envelopes

[0285] According to an embodiment of the invention, during the sixth step C3 (following the second step C21 of the respective iteration carried out and/or at the fifth step C22 and prior to the third step C4) an auxiliary decoupling variable v of the complex envelopes a.sub.k(n) of the components between them in the model is calculated by the calculation module 12, from the updated complex envelopes a.sub.k(n), from the first parameters (which can be one or more or all of the τ.sub.j for j ranging from 1 to J, and/or one or more or all of the β.sub.j for j ranging from 1 to J, and/or one or more or all of the η.sub.j for j ranging from 1 to J, and/or one or more or all of the ζ.sub.j for j ranging from 1 to J) and from the second regularisation parameters γ.sub.k. A known problem of optimisation algorithms is the inversion of a full matrix of large dimension, most often badly conditioned. This embodiment resolves this problem by decoupling the components by addition of the auxiliary variable v. The auxiliary variable v takes into account information on the correlations between the components a.sub.k(n). Due to the auxiliary variable v, the components are estimated independently. Coupling occurs implicitly via the auxiliary variable.

[0286] This embodiment resolves the large-dimension digital problem, linked to the above first problem. This embodiment is adapted all the more since a large number K of components is preferred. In fact, the known Vold-Kalman method suffers from the problem of large-dimension matrix inversion above all when the matrix in question is very badly conditioned which is the case for very large values of r.sub.k. This inversion could occur in the known Vold-Kalman method due to an iterative algorithm of conjugated gradient type. This embodiment is not contradicted by this digital problem, given that the different components are decoupled in the algorithm and therefore can be estimated disjointly. The coupling between the different components is considered implicitly by the auxiliary variable v added in the updating of each component.

[0287] This embodiment of the invention on the auxiliary decoupling variable v can be realised by one of the other embodiments described hereinbelow.

[0288] According to an embodiment of the invention, the law of the auxiliary decoupling variable v is defined by a model having a distribution of probability of complex circularly symmetrical normal law type CNC(m.sub.v, Σ.sub.v) having the average m.sub.v and the variance Σ.sub.v. For v there is a conditional law

[Math. 137]

v|a,d˜CNC(m.sub.v,Σ.sub.v)

[0289] The complex circularly symmetrical normal law CNC(m.sub.v, Σ.sub.v) has the average

[00024]$\begin{array}{cc}\left(\begin{array}{c}{m}_v^R\\ m_v^I\end{array}\right)& \left[\mathrm{Math}.138\right]\end{array}$

[0290] and covariance matrix

[00025]$\begin{array}{cc}\frac{1}{2}\left(\begin{array}{cc}{.Math.}_v^R;& -{.Math.}_v^I\\ {.Math.}_v^I;& {.Math.}_v^R\end{array}\right)& \left[\mathrm{Math}.139\right]\end{array}$

[0291] and is defined by

[00026]$\begin{array}{cc}v\sim \mathrm{CNC}\left(m_v,{.Math.}_v\right)\text{<=>}\left(\begin{array}{c}{v}^R\\ v^I\end{array}\right)\sim N\left(\left(\begin{array}{c}{m}_v^R\\ m_v^I\end{array}\right),\frac{1}{2}\left(\begin{array}{cc}{.Math.}_v^R;& -{.Math.}_v^I\\ {.Math.}_v^I;& {.Math.}_v^R\end{array}\right)\right)& \left[\mathrm{Math}.140\right]\end{array}$

[0292] In the above,

[00027]$\begin{array}{cc}{m}_v={.Math.}_{v}a& \left[\mathrm{Math}.141\right]\end{array}$$\begin{array}{cc}{.Math.}_v=\frac{1}{\mu }\mathrm{Id}-C^{\top }\underset{j=1}{\overset{J}{.Math.}}{S}_j^{\top }{\mathrm{diag}\left(d_j\right)}^{-1}{S}_{j}C& \left[\mathrm{Math}.142\right]\end{array}$

[0293] with μ being a real prescribed non-zero verifying

[00028]$\begin{array}{cc}\mu <\frac{1}{K\max \left(d_j\left(n\right)\right)}& \left[\mathrm{Math}.143\right]\end{array}$

[0294] This enables to eliminate problems of inversion of large-dimension matrices in particular for large values of N and K and allows an increase in data, which allows decoupling between the different components a.sub.k (n)exp(jΦ.sub.k (n)).

[0295] Calculation of the Sinusoidal Components a.sub.k of Interest During the First Step C1

[0296] According to an embodiment of the invention, the law a posteriori of the real part a.sub.k.sup.R of the sinusoidal components of interest a.sub.k (n)exp(jΦ.sub.k (n)) is defined by a model having a distribution of probability of normal law type having the average m.sub.k.sup.R and the variance Σ.sub.k, for example according to the following equation:

[Math. 144]

∀k∈{1, . . . ,K}a.sub.k.sup.R|y,d,γ.sub.k,v.sub.k˜N(m.sub.k.sup.R,Σ.sub.k)

According to an embodiment of the invention, the law a posteriori of the imaginary part a.sub.k.sup.I of the sinusoidal components of interest a.sub.k(n)exp(jΦ.sub.k(n)) is defined by a model having a distribution of probability of normal law type having the average m.sub.k.sup.I and the variance Σ.sub.k, for example according to the following equation:

[Math. 145]

∀k∈{1, . . . ,K}a.sub.k.sup.I|y,d,γ.sub.k,v.sub.k˜N(m.sub.k.sup.I,Σ.sub.k)

[0297] The definitions hereinabove of a.sub.k.sup.R and of a.sub.k.sup.I are laws a posteriori of the components.

[0298] In the above,

[00029]$\begin{array}{cc}{m}_k={.Math.}_k\left(C_k^{\top }\underset{j=1}{\overset{J}{.Math.}}{S}_j^{\top }{\mathrm{diag}\left(d_j\right)}^{-1}{S}_{j}y+v_k\right)& \left[\mathrm{Math}.146\right]\end{array}$$\begin{array}{cc}{.Math.}_k={\left(\frac{1}{\mu }\mathrm{Id}+{\gamma }_k{A}_p^{\top }{A}_p\right)}^{-1}& \left[\mathrm{Math}.147\right]\end{array}$

[0299] Therefore, a.sub.k.sup.R and a.sub.k.sup.I can be estimated independently. The coupling is considered only via the auxiliary variable v occurring in the average m.sub.k. It is interesting to note that even if there is a need to make the inversion here, this operation is less expensive. On the one hand, the dimension of the matrix Σ.sub.k is K times smaller than the matrix Σ.sub.x. On the other hand, if the choice is made for conditions of circulating edges (or zero edges by making padding in the zero edges) in the matrix of differences A.sub.p, the matrix Σ.sub.k is then circulating and therefore diagonalisable in the Fourier field.

[0300] According to an embodiment of the invention,

[00030]$\begin{array}{cc}{.Math.}_k={F_N^{\top }\left(\frac{1}{\mu }\mathrm{Id}+{\gamma }_k{Q}_k^{\top }{Q}_k\right)}^{-1}{F}_N& \left[\mathrm{Math}.148\right]\end{array}$

[0301] where Q.sub.k is the diagonal matrix containing the Fourier coefficients of the filter associated with the discrete matrix of difference A.sub.p which is equal to the first column of A.sub.p. Updating of each m.sub.k can also be done by shifting to the Fourier field. Also, there is no need to store Σ.sub.k, or its inverse, but simply the diagonal values of the matrix

[00031]$\begin{array}{cc}\frac{1}{\mu }\mathrm{Id}+{\gamma }_k{Q}_k^{\top }{Q}_k& \left[\mathrm{Math}.149\right]\end{array}$

[0302] According to an embodiment of the invention, the calculation module 12 is configured to perform a seventh decision step E following the fourth step D. According to an embodiment of the invention, the decision during the seventh step E is taken automatically by the calculation module 12 from analyses carried out automatically by the calculation module 12 on the extracted components a.sub.k(n) having been obtained at the fourth step D, and/or from test bench vibratory characterisation (for example identification of eigen modes) carried out by the calculation module 12 on the extracted components a.sub.k(n) having been obtained at the fourth step D of the rotating machine 100. According to an embodiment of the invention, the seventh step E is performed by the calculation module 12 on a monitoring decision and/or on the triggering of an alarm AL or alert AL of defects of the rotating machine 100 (for example: imbalance, meshing, misalignment, mass eccentricity, worn bearings, twisted shafts, misaligned shafts, parallel alignment defect of a rotor shaft of the rotating machine, broken shaft, premature wear of bearings, joints, shafts and couplings of the kinematic chain, resonance). According to an embodiment of the invention, the monitoring device 1 is embedded on the aircraft. According to an embodiment of the invention, the calculation module 12 automatically controls communication of the alert, of the alarm AL or of the decision it has taken at the seventh step E on the physical outlet SP, for example on the display screen, and/or at a fixed ground station SP.

[0303] Updating of the law in the proposed device can be done in two ways.

[0304] On the one hand, it can be approximated by samples. This concerns the Monte Carlo Markov Chain methods (MCMC). Therefore, updating a law means taking a sample of this law. It is noted that all laws are simple to sample. For the parameters of the law of observation and of the regularisation parameter, this is evident. For the law a posteriori of the envelopes, the generation of random variables can be done by shifting to the Fourier field. And given the particular form of Σ.sub.v, the generation of random auxiliary variables can also be done directly based on the fact that CC.sup.T=K. Id and S.sub.jS.sub.j.sup.T=Id. After enough iterations, the random variables generated by the algorithm will follow the preferred laws a posteriori; therefore the statistics (average, variance etc) can be approximated by using empirical estimators with these samples.

[0305] On the other hand, the law a posteriori can be approximated by another simpler one, for example by a Bayesian-Variational (ABV). approximation approach. In this way the statistics of laws of interest can be approximated by those of approximating laws after enough iterations.

[0306] An example of application of the invention on a vibratory signal of a helicopter is given hereinbelow in FIGS. 4a to 4f. The aim is to monitor the component N.sub.2 of the vibratory signal y(t) of the high-pressure shaft of a turbomachine 100, for example of the shaft of the free turbine in the helicopter engine.

[0307] FIG. 4a shows the spectrogram of the measured signal y(t) by indicating, for each instant in abscissae and each frequency on the frequency band [0,700] Hz in ordinates, a pixel whereof the greyscale corresponds according to the right-hand scale of values ECH to a spectral amplitude for this frequency (and similarly for the spectrograms of FIGS. 4c, 4d, 4e and 4f). FIG. 4b shows in ordinates the amplitude of the envelope a.sub.k(n) having been calculated according to the invention as a function of time in abscissae (noise is considered non-stationary and non-Gaussian and the regularisation parameter is estimated automatically). FIG. 4b shows a significant rise in amplitude from the instant 5s, which is the sign of a passage through a resonance. The evaluation of the maximal amplitude during this resonance is highly important practically for vibratory control of the aircraft in transient operation. FIG. 4c shows the spectrogram of the residual signal (the measured signal y(t) from which the component N.sub.2 is subtracted). FIG. 4c shows that the component N.sub.2 of interest has been correctly estimated. FIGS. 4d, 4e and 4f show the spectrograms of the residual signal by using the known method Vold-Kalman with different regularisation parameters adjusted manually (these are values currently used in literature). FIG. 4d shows an underestimation of the amplitude of the component N.sub.2 of interest when passing through resonance for a large value of r.sub.k (the estimated amplitude of the component N.sub.2 of interest is less important than the real one, the estimated envelope is very smoothed excessively, which very aggressively penalises the great variations of the envelope). FIGS. 4e and 4f show that there is overestimation of the amplitude of the component N.sub.2 of interest, for a small value of r.sub.k (the estimated amplitude of the component N.sub.2 of interest is more important than the real amplitude, the estimated envelope is very affected by noise). In fact, the component N.sub.2 of interest having been estimated in FIGS. 4e and 4f is badly filtered and contains much background noise. Also, overestimation of the envelope can result from a badly conditioned matrix and therefore from digital problems (see FIGS. 1 and 4f for example).

[0308] According to another embodiment, β.sub.j is equal to zero. In this case the contribution of the narrowband noise e(n) is omitted. The principal advantage is the reduction in the number of parameters to be estimated (β.sub.j and ζ.sub.j) as well as the calculation cost. Admittedly, this can be a good approximation if no parasite harmonic of the residual signal interferes directly with the signals of interest. This goes back to y|a, τ˜NC (Ca,2τ.sup.−1)

[0309] In this particular case, the advantage relative to the known Vold-Kalman method is the automatic adjusting of the regularisation parameter γ.sub.k. In fact, the Vold-Kalman parameter r.sub.k is linked to the parameters of the invention by r.sub.k=γ.sub.k/τ.

[0310] A disadvantage of this approximation could manifest, for example, in the presence of interception between the harmonics of interest and the residual harmonics. In fact, estimation artefacts can occur at instants of interferences.

[0311] However, if the instants of interferences are known approximately (for example with inspection of the spectrogram), it can be assumed that the noise is Gaussian except for close to coupling, that is at intervals j ∈{1, . . . , J} containing the coupling instants. This means imposing β.sub.j=0 in the other intervals not containing the coupling moments.

[0312] Of course, the embodiments, characteristics, possibilities and examples described hereinabove can be combined with each other or be selected independently of each other.

[0313] If a parsimonious model other than the Bernoulli-Gaussian model is used, the same forms of laws are obtained. In this case, the parameters of the noise models are going to change. The law of the envelopes, of the regularisation parameters and the auxiliary variable will have the same form. The chosen noise model will act only on the values of d.sub.j.