Algebraic Aspects of the Exact Signal Demodulation Problem

Algebraic Aspects of the Exact Signal Demodulation Problem

7th IFAC Symposium on System Structure and Control 7th IFAC Symposium on System Structure and Control Sinaia, Romania, September 9-11,Structure 2019 7...

574KB Sizes 0 Downloads 28 Views

7th IFAC Symposium on System Structure and Control 7th IFAC Symposium on System Structure and Control Sinaia, Romania, September 9-11,Structure 2019 7th IFAC Symposium on System and Control Sinaia, Romania, September 9-11, 2019 Available online at www.sciencedirect.com 7th IFAC Symposium on System and Control Sinaia, Romania, September 9-11,Structure 2019 Sinaia, Romania, September 9-11, 2019

ScienceDirect

IFAC PapersOnLine 52-17 (2019) 82–87

Algebraic Aspects of the Exact Signal Algebraic Aspects of the Exact Signal Algebraic Aspects of the Exact Signal Demodulation Problem Algebraic Aspects of the Exact Signal Demodulation Problem Demodulation Problem Demodulation Problem ∗∗∗ ∗∗∗∗ E. Hubert ∗∗ Y. Bouzidi ∗∗ ∗∗ R. Dagher ∗∗∗ A. Barrau ∗∗∗∗

∗∗∗ A. Barrau ∗∗∗∗ E. Hubert ∗∗ Y. Bouzidi ∗∗ ∗∗ R. Dagher ∗∗∗ † E. Hubert ∗ Y. Bouzidi R. Dagher A. Barrau ∗∗∗∗ A. Quadrat † A. Quadrat ∗∗ † E. Hubert Y. Bouzidi R. Dagher ∗∗∗ A. Barrau ∗∗∗∗ A. Quadrat † A. Quadrat ∗ (e-mail: [email protected]) ∗ Safran Tech, France Tech, France (e-mail: [email protected]) ∗ Safran ∗∗ Safran Tech, France (e-mail: [email protected]) Lille, France (e-mail: [email protected]) ∗∗ Inria Inria Lille, France (e-mail: [email protected]) ∗ ∗∗ ∗∗∗ Safran Tech, France (e-mail: [email protected]) Inria Lille, France (e-mail: [email protected]) Inria Chile, Chili (e-mail: [email protected]) ∗∗∗ ∗∗ Inria Chile, Chili (e-mail: [email protected]) ∗∗∗∗ ∗∗∗ Inria Lille, France (e-mail: [email protected]) Inria Chile, Chili (e-mail: [email protected]) Tech, France (e-mail: [email protected]) ∗∗∗∗ Safran ∗∗∗ Safran Tech, France (e-mail: [email protected]) ∗∗∗∗ † Inria Chile, Chili(e-mail: (e-mail: [email protected]) Safran Tech, France [email protected]) Ouragan project-team, IMJ – † Inria Paris, Inria Paris, Ouragan project-team, IMJ – PRG, PRG, Sorbonne Sorbonne University, University, ∗∗∗∗ † SafranOuragan Tech, France (e-mail: [email protected]) Inria Paris, project-team, IMJ – PRG, Sorbonne University, France (e-mail: [email protected]) France (e-mail: [email protected]) † Inria Paris, Ouragan project-team, IMJ – PRG, Sorbonne University, France (e-mail: [email protected]) France (e-mail: [email protected]) Abstract: Abstract: In In this this paper paper we we introduce introduce aa general general class class of of problems problems originating originating from from gearbox gearbox Abstract: In this paper a general class ofdemodulation problems originating from gearbox vibration analysis. Based we onintroduce a previous work where was formulated as a vibration analysis. Based on a previous work where demodulation was formulated as Abstract: In this paper we general classcase ofdemodulation problems from gearbox vibration analysis. Based onintroduce aweprevious work where formulated as aa matrix approximation problem, studya the specific applicableoriginating towas amplitude and phase matrix approximation problem, study the specific case applicable towas amplitude and phase vibration analysis. Based on can awe work where demodulation formulated as a matrix approximation problem, weprevious study the specific case applicable to amplitude phase demodulation. This be rewritten as system. Based algebraic demodulation. This problem problem can be rewritten as aa polynomial polynomial system. Based on onand algebraic matrix approximation problem, we study the specific case applicable to amplitude phase demodulation. This problem can behomological rewritten asalgebra, a polynomial onand algebraic methods such as linear algebra and we focussystem. on the Based characterization of methods such as linear algebra and algebra, we focussystem. on the characterization of demodulation. This problem behomological rewritten a polynomial on algebraic methods suchand as linear and homological we focus on the Based characterization of the problem solve italgebra in thecan noise-free case. asalgebra, the problem in the noise-free case. algebra, we focus on the characterization of methods suchand as solve linearit and homological the problem and solve italgebra in the noise-free case. Copyright © 2019. The Authors. Published by Elsevier Ltd. All rights reserved. the problem and solve it in the noise-free case. Keywords: Demodulation, polynomial systems, linear systems, homological algebra. Keywords: Demodulation, polynomial systems, linear systems, homological algebra. Keywords: Demodulation, polynomial systems, linear systems, homological algebra. Keywords: Demodulation, polynomial systems, linear systems, homological algebra. 1. of 1. INTRODUCTION INTRODUCTION of modulated modulated signals. signals. The The energy energy separation separation algorithm algorithm 1. INTRODUCTION of separation algorithm hasmodulated also been signals. improvedThe by energy the iterative generalized dehas also been improved by the iterative generalized de1. INTRODUCTION of modulated signals. The energy separation algorithm has also been improved by the iterative generalized demodulation method proposed in Feng et al. (2011). In Modulation and demodulation are basic tools of the theory modulation method proposed initerative Feng etgeneralized al. (2011). deIn also been improved thein Modulation and demodulation are basic tools of the theory has method proposed Feng et al.provides (2011). an In the spectral domain, thebyHilbert transform Modulation and demodulation tools of theory of signal processing which canare bebasic described as the varying a modulation the spectral method domain, proposed the Hilbert transform provides an in Feng et al.provides (2011). an In of signal processing which canare be described as the varying a modulation spectralof the Hilbert transform estimation the of modulated Modulation and demodulation tools of theory of signal processing which can bebasic described asamplitude, varying a the parameter of a sine wave, which can be its estimation ofdomain, the envelope envelope of aa mono-carrier mono-carrier modulated parameter of a sine wave, which can be its amplitude, the spectral domain, the Hilbert transform provides an estimation of the envelope of a mono-carrier modulated signal, along with its phase if needed (Potamianos et al. of signal processing which can be described as varying a parameter of a sinetowave, which can be(modulation) its amplitude, phase or frequency, encode a message or signal, along with its phase of if needed (Potamianos et al. the envelope mono-carrier modulated phase or frequency, towave, encode a message (modulation) or estimation with its phase if aneeded (Potamianos et al. (1994)).along Forofmulti-carrier demodulation, most approaches parameter of a message sineto can be(modulation) its amplitude, phase or frequency, encode a message or signal, to recover the fromwhich such a modified sine wave (1994)). For multi-carrier demodulation, most approaches signal, along with phase if needed (Potamianos et al. to recover the message from such a modified sine wave For multi-carrier demodulation, most into approaches are two-steps: the its signal is first decomposed monophase or frequency, to encode atelecommunications, message (modulation) or (1994)). to recover the message frombysuch a modified sine wave (demodulation). Motivated anaare two-steps: the signal is first decomposed into mono(1994)). For multi-carrier demodulation, most approaches (demodulation). Motivated by telecommunications, anaare two-steps: the signal is first decomposed into monocomponent signals which are then demodulated with one to recover the message fromdeveloped such a modified sine wave (demodulation). Motivated telecommunications, ana- are logical first to component signals which are thendecomposed demodulated with one logical methods methods where where first by developed to perform perform these these two-steps: themethods signal isorfirst into monocomponent signals which are then demodulated with one of the previous any variation of them. In (demodulation). Motivated by telecommunications, analogical methods first developed perform operations in thewhere time domain althoughtothe theorythese sup- component of the previous methods or any variation of them. In signals which the are then demodulated with one operations in thewhere time domain althoughtothe theorythese supof the previous methods or any variation of them. In Friedlander et al. (1995), separation is based on bandlogical methods first developed perform operations in was the based time domain although the theory sup- of porting them on spectral considerations. Friedlander et (1995), is based on previous methods oretseparation any variation of them. In porting them was based on spectral considerations. etInal. al.Santhanam (1995), the the separation isthe based on bandbandpassthe filtering. al. (2000), authors take operations in was the based time domain although the theory sup- Friedlander porting them on spectral considerations. pass filtering. In Santhanam et al. (2000), the authors take Friedlander et al. (1995), the separation is based on bandAlthough these early solutions haveconsiderations. encountered a wide pass authors take filtering. In Santhanam et al. (2000), the advantage of the periodicity of the signal. The separation porting them was based on spectral advantage of the periodicity et of the(2000), signal. The separation Although these early solutions have encountered aa wide filtering. In Santhanam authors take wide pass Although these solutions haveapproach encountered success and are early still the standard to demoduadvantage of the periodicity of al. the signal.the The separation is then performed with a periodicity-based algebraic algosuccess and are still the standard approach to demoduis then performed with a periodicity-based algebraic algoAlthough these early solutions haveapproach encountered a wide advantage of the periodicity of the signal. The separation success and are still the standard to demodulation, improving them is an active field of research due thenand performed with a periodicity-based algebraic algorithm the demodulation with lation, improving them isstandard an activeapproach field of research due is with an an energy-based energy-based algorithm and the demodulation success and are still theis to demoduis thenand performed with a periodicity-based algebraic algolation, improving them an active fieldby of digital research due rithm to both the new possibilities brought signal the demodulation with an energy-based algo(PASED). Another method, called Hilbert-Huang to both the new possibilities brought by digital signal rithm (PASED). Another method, called Hilbert-Huang lation, improving them is an active field of research due rithm and the demodulation with an energy-based algoto both theand new by digital signal processing thepossibilities specificities brought of new applications such rithm (PASED). Another method, called Hilbert-Huang transform, introduced in Huang et al. (1998), analyzes introduced in Huang et al. (1998), analyzes processing and thepossibilities specificities brought of new applications such transform, to both theand new by digital(1986)). signal (PASED). Another method, Hilbert-Huang processing specificities of new applications such rithm as monitoring ofthe mechanical systems (McFadden transform, introduced in Huang et called al. (1998), analyzes the signal by taking it apart into intrinsic mode functions as monitoring of mechanical systems (McFadden (1986)). the signal by taking it apart into intrinsic mode functions processing andofthe specificities of new applications such transform, introduced in mode Huang et al. (1998), analyzes as mechanical (McFadden The latter is the motivation for present signal by taking it apart intodecomposition. intrinsic mode A functions (IMF) using the empirical Hilbert Themonitoring latter work work is the initial initialsystems motivation for the the (1986)). present the (IMF) using the empirical mode decomposition. A Hilbert as monitoring of is mechanical (McFadden (1986)). signal by taking it apart into intrinsic modeIMFs. functions The latter work the initialsystems motivation for the present work. The models encountered in these fields differ from the (IMF) using the empirical mode decomposition. A Hilbert transform is then performed on the obtained For transform is then performed on the obtained IMFs. For work. The models encountered in these fields differ from The latter work the initial motivation for the present using the empirical mode decomposition. Ahave Hilbert in these fields differ from (IMF) work. The models the classical sineis encountered wave modulation/demodulation (also transform is then performed on the obtained IMFs. For non-stationary analysis, wavelet based methods the the classical sine wave modulation/demodulation (also non-stationary analysis, wavelet based methods have the work. The models in these fields from is allowing then performed on based theofobtained IMFs. For the classical sine encountered wave modulation/demodulation (also called “mono-carrier modulation”) because thediffer energy of transform non-stationary analysis, wavelet methods have the advantage of demodulation signals whose timecalled “mono-carrier modulation”) because the energy of advantage of allowing demodulation of signals whose timethe classical sine wave modulation/demodulation (also analysis, wavelet methods the called “mono-carrier because theharmonics energy of non-stationary the signal is several advantage of allowing demodulation of signals timefrequency distributions are curvedbased paths (Yu etwhose al.have (2016)). the modulated modulated signalmodulation”) is spread spread over over several harmonics frequency distributions are curved paths (Yu et al. (2016)). called “mono-carrier modulation”) because the energy of of allowing demodulation of signals whose timethe modulatedmodulation) signal is spread over several (multi-carrier and filtering out allharmonics of them advantage frequency distributions are curved paths (Yu et al. (2016)). In mechanical signal processing, multicomponent demodIn mechanical signal processing, multicomponent demod(multi-carrier modulation) and filtering out all of them frequency the modulated signal is spread over distributions are curved paths to (Yufault et al.diagnosis (2016)). (multi-carrier and out allharmonics of them In but one would modulation) mean loosing mostfiltering of the several available informamechanical signalhave processing, multicomponent demodulation algorithms been adapted but one would mean loosing most of the available informaulation algorithms have been adapted to fault diagnosis (multi-carrier filtering out all informaof them In mechanical signalhave processing, multicomponent demodbut would mean loosing most ofin the available tion.one Note this modulation) situation alsoand exists telecommunications ulation algorithms been adapted to fault diagnosis (Feldman (2011)). tion.one Note this mean situation alsomost existsofin telecommunications (Feldman (2011)). have been adapted to fault diagnosis but would loosing the available informa- ulation algorithms tion. Note this situation also exists in telecommunications although it is not standard (Friedlander et al. (1995)). (Feldman (2011)). although itthis is not standard (Friedlander et al. (1995)). In the work Hubert et al. (2018), the amplitude demodtion. Noteit situation also(Friedlander exists in telecommunications (Feldman (2011)). although is not standard et al. (1995)). In the work Hubert et al. (2018), the amplitude demodRegardingit the new possibilities offered et by al. digital signal In the work Hubert the amplitude ulation was studied etas al. an (2018), optimization problemdemodwhere although is not standard (Friedlander (1995)). Regarding the new possibilities offered by digital signal ulation was studied etas al. an optimization problemdemodwhere In work thethe amplitude the new possibilities offered by digital signal Regarding processing, the first one is to go from the time to the ulation was Hubert studied as couple an (2018), optimization problem where thethe carrier/modulation fitting data at best is processing, the first one is to go from the time to the the carrier/modulation couple fitting the problem data at best is Regarding the possibilities by signal wasshowed studied as an optimization where gooffered fromFourier thedigital time to the ulation processing, the new first one isoftothe spectral domain by means Fast Transform the carrier/modulation couple fitting the data atbetween best is sought. We an unexpected correspondence sought. We showed an unexpected correspondence between spectral domain by means of the Fast Fourier Transform processing, the time first one isoftothe goFast fromFourier thedemodulation time to the the couple matrix fitting the data atbetween best is spectral domain by means Transform (FFT). In the domain, multi-carrier sought. We showed anlow-rank unexpected correspondence this carrier/modulation formulation and approximation, and (FFT). In the time domain, multi-carrier demodulation this formulation and low-rank matrix approximation, and spectral domain by means the FastofFourier Transform We showed an unexpected correspondence between time domain, multi-carrier demodulation (FFT). amountsIn tothe estimating the of envelope a signal. The en- sought. formulation and low-rank matrix approximation, leveraged it to design a new “optimal demodulation”and alamountsIntothe estimating the envelope of a signal. The en- this leveraged it to design a new “optimal demodulation”and al(FFT). time domain, multi-carrier demodulation formulation and low-rank matrix approximation, amounts to estimating theachieving envelope of atask signal. The en- this ergy separation algorithm this uses nonlinleveraged it to on design a new “optimal demodulation” algorithm based eigenvector computation. In this paper, ergy separation algorithm achieving this task uses nonlingorithm based on eigenvector computation. In this paper, amounts to estimating the envelopethis ofetatask signal. en- leveraged it to on design a new “optimal demodulation” alalgorithm achieving usesThe nonlinergy separation ear differential operators (Potamianos al. (1994)) such based eigenvector computation. paper, we the same and the approach we keep keep following following the same route route and extend extendIn thethis approach ear differential operators (Potamianos ettask al. (1994)) such gorithm ergy separation algorithm achieving this uses nonlingorithm based on eigenvector computation. In this paper, ear differential operators (Potamianos al. (1994)) such we as the Teager-Kaiser energy operator et derived in Kaiser keep following the same route and extend the approach to combine amplitude and phase demodulation. Under as the Teager-Kaiser energy operator et derived in Kaiser to combine amplitude and phase demodulation. Under ear differential operators (Potamianos al. (1994)) such we following the same anddemodulation. extend the approach as the Teager-Kaiser energy operator derived in Kaiser (1990), which the instantaneous frequency to keep combine amplitude androute phase Under (1990), which also also extracts extracts the instantaneous frequency as the Teager-Kaiser energy operator derived in Kaiser (1990), which also extracts the instantaneous frequency to combine amplitude and phase demodulation. Under (1990), also© 2019. extracts the instantaneous frequency 2405-8963which Copyright The Authors. Published by Elsevier Ltd. All rights reserved.

Copyright © 2019 IFAC 160 Copyright 2019 IFAC 160 Control. Peer review© under responsibility of International Federation of Automatic Copyright © 2019 IFAC 160 10.1016/j.ifacol.2019.11.031 Copyright © 2019 IFAC 160

2019 IFAC SSSC Sinaia, Romania, September 9-11, 2019

E. Hubert et al. / IFAC PapersOnLine 52-17 (2019) 82–87

the assumption that the phase variations stay small with respect to 2 π, the first contribution of the paper is to show that a matrix reformulation, similar to the one developed in Hubert et al. (2018) can be obtained. It yields a constrained rank 2 approximation problem which, to the authors knowledge, has not be given much attention in the linear algebra literature. Our second contribution is to study and effectively solve this problem by means of algebraic methods.

83

Let us now consider a modified version of (3) which encodes the phase and amplitude modulation signal: s(t) = c(t + φ(t)) m(t), m ∈ D(Nm , Tm ), c ∈ D(Nc , Tc ), fc = Tc−1 ,  (4) / {1, . . . , Nc },  ck = 0, k ∈ L(m, c) = mk = 0, for |k| ≥ fc /2,  φk = 0, for |k| ≥ fc /2.

In Section 2, we first derive a novel formulation of the amplitude and phase demodulation as a constrained lowrank matrix approximation problem. In Section 3, we review standard results on linear algebra that will be of constant use. Finally, Section 4 studies the rank 2 decomposition problem corresponding to the exact phase and amplitude demodulation, and gives explicit necessary and sufficient conditions for the existence of a solution.

A study of this general system is difficult. But under the additional hypothesis that the phase modulation is small φ(t) << 2 π, we can recast it into the general framework of (2). Indeed, we can perform a first-order expansion w.r.t. φ and write: c(t + φ(t)) m(t) = (c(t) + c (t) φ(t)) m(t) = c(t) m(t) + c (t) φ(t) m(t). Setting c1 = c, c2 = c , m1 = m, m2 = m φ, (4) yields:

2. PROBLEM FORMULATION

s(t) = c1 (t) m1 (t) + c2 (t) m2 (t), m1 ∈ D(Nm , Tm ), m2 ∈ D(Nm , Tm ), c1 ∈ D(Nc , Tc ), c2 ∈ D(Nc , Tc ), fc = Tc−1 ,   c 1 = c2 ,    (c ) = 0, k∈ / {1, . . . , Nc }, 1 k L(m, c) =  ) = 0, for |k| ≥ fc /2, (m   1 k (m2 /m1 )k = 0, for |k| ≥ fc /2.

2.1 Modulation-based mechanical models Let T ∈ R>0 be a non-negative real number and D(T ) the set of T -periodic functions of L1 ([0, T ]). According to Fourier analysis, for c ∈ D(T ), we have:  T /2  k k 1 i2π T t ck e , ck = c(t) e−i 2 π T t dt. c(t) = T −T /2 k∈Z

(1) of c is the term defined For n, N ∈ N,n the nth harmonic n by c−n e−i 2 π T t + cn ei 2 π T t and we call the first N harmonics of c the following truncation of (1):  k ck ei 2 π T t . Hn (c)(t) = 0≤|k|≤N

Let D(N, T ) be the set formed by all the functions of D(T ) whose at most first N harmonics are non-zero, i.e.: D(N, T ) = {f ∈ D(T ) | |k| > N : ck = 0}. In this paper, we consider a class of time signals arising in the study of gearbox vibrations that are of the form q  mj (t) cj (t), L(m, c) = 0, (2) s(t) = j=1

where s denotes the output of the sensor, m = (mj )j=1,...,q and c = (cj )j=1,...,q are the unknown components to be estimated which are subjected to certain linear constraints L(m, c) = 0. Let us give a few example of such situations. Amplitude modulation without overlap is the prototypical example of (2) which can be obtained by considering s(t) = c(t) m(t), (3) with the constraints m ∈ D(Nm , Tm ), c ∈ D(Nc , Tc ), Tm Nm = T for a certain Nm ∈ N (resp., Tc Nc = T for a certain Nc ∈ N), and  / {1, . . . , Nc }, ck = 0, k ∈ L(m, c) = mk = 0, for |k| ≥ fc /2,

where fc = Tc−1 . Hence, the value of q is 1 and the linear constraints amount to saying that c has frequency fc and the spectral support of m is in ] − fc /2, fc /2[.

161

(5)

2.2 Separation problem

In the phase and amplitude demodulation problem we are interested in solving, two different situations have to be considered. First, the ideal case described in Section 2.1, referred to here as the noise-free problem. Second, the more realistic case where noise is taken into consideration through what we call an optimal demodulation. In the noise-free situation, finding those couples will be referred to as Problem 1. Problem 1. (Exact demodulation). Given N ∈ N and (1) a discrete signal (s(tn ))n=1,...,N of sampling period Ts , i.e., tn = n Ts , and a duration Ttot = N Ts , (2) a frequency fc = kc · ftot , where ftot = 1/Ttot and kc ∈ N is such that N = kc Nc for a certain Nc ∈ N,

the problem aims to find a sequence of carrier/modulation pairs ((cj (tn ), mj (tn )))n=1,...,N satisfying q  ∀ n = 1, . . . , N, s(tn ) = cj (tn ) mj (tn ), j=1

where (cj (tn ))n=1,...,N is Tc -periodic with Tc = Nc Ts , (mj (tn ))n=1,...,N is Tmj -periodic with Tmj = Nmj Ts , and (cj (tn ), mj (tn )) verifying the constraints of (4).

However, when we consider a more realistic signal, i.e., a noisy signal, we have to use an optimal demodulation. Formally, we seek for a sequence of pairs of signals ((ˆ cj , m ˆ j ))j=1,...,q solving the optimization problem. Problem 2. (Optimal demodulation). Given (1) and (2) as stated in Problem 1, find a sequence carrier/modulation pairs ((cj (tn ), mj (tn )))n=1,...,N which minimizes the following cost function

2019 IFAC SSSC 84 Sinaia, Romania, September 9-11, 2019

E. Hubert et al. / IFAC PapersOnLine 52-17 (2019) 82–87

2   N     q  C(c, m) = (cj (tn ) mj (tn )) − s(tn )   n=1  j=1

over all carrier/modulation pairs verifying the constraints of (4), and where (cj (tn ))n=1...N is Tc -periodic and (mj (tn ))n=1...N is Tm -periodic.

Problem 3. (Exact demodulation). If  cj (resp., m  j ) denotes the vector containing the Fourier coefficients of cj (tn ) (resp., mj (tn )), and H the Hermitian transpose operator, then Problem 2 is equivalent to finding the vectors  cj ’s and the vectors m  j ’s satisfying: q   cj m H Ms = j . j=1

2.3 Matrix formulation In practice, a common specific case allows a very handy reformulation of the problem. First, let us transpose it to the Fourier domain. The Discrete Fourier Transform (DFT) of the obtained discrete signal of total length N (s(tn ))n=1,...,N , denoted here with a tilde sign s, namely, ∀ k = 1, . . . , Nc + Nm ,

s[k] =

N 

k

s(tn ) e−i 2 π N n , (6)

n=1

is then a circular convolution of the carrier and modulation q DFTs, i.e., s = j=1  cj ∗ m  j.

In the case where the patterns centered on two consecutive harmonics do not overlap, the spectrum of the modulation is directly accessible from the spectrum of s. This is the well-known condition making demodulation possible. Hypothesis 1. In what follows, we shall assume that: (7) 2 Nm f m < fc . Note that (cj (tn ))n=1,...,N is Tc -periodic with Tc = Nc Ts , i.e., cj (tn+Nc ) = c(tn ), (mj (tn ))n=1,...,N is Tm -periodic with Tm = Nmj Ts , i.e., mj (tn+Nm ) = mj (tn ). Let (s(tn ))n=1,...,N be a discretized time signal and s its DFT defined by (6). The matrix representation we propose consists in cutting s into buckets centered on each harmonic of the carrier and then stacking them vertically as illustrated by Figure 1. Let us write down

Problem 4. (Optimal demodulation). If  cj (resp., m  j ) denotes the vector containing the Fourier coefficients of cj (tn ) (resp., mj (tn )), and H the Hermitian transpose operator, then Problem 2 is equivalent to finding the vectors  cj ’s and the vectors m  j ’s which minimizes the function  2  q     C( c, m)  =   cj m H − M , s j  j=1  F ro 1/2   n m 2 denotes the sowhere |M |F ro = i=1 j=1 |Mij | called Frobenius norm of M ∈ Kn×m .

Due to length constraint of this paper, we shall only consider here Problem 3. For the study of Problem 4, we refer the reader to Hubert et al. (2019). In the application case we are interested in, i.e., phase and amplitude modulation signals, Problem 3 gives a reformulation of (4) in the following matrix form Ms =  c1 m H c2 m H 1 + 2 ,  where  c2 =  c1 is the derivative of  c1 in the frequency domain. This derivative operation can be expressed as  c1 = D c1 , where D is the diagonal matrix computing the derivative in the frequency domain, i.e.: D = i 2 π fc diag (−Nc , . . . , Nc ) . Thus, we obtain the following general formulation of the main problem considered in this paper: Ms =  c1 m H c1 m H 1 +D 2 . 3. SOLVING INHOMOGENEOUS LINEAR SYSTEMS

3.1 Notation & basic homological algebra

Fig. 1. Matrix spectrum construction from the line spectrum vector a mathematical definition of this matrix Ms . Definition 1. (Matrix representation of a spectrum). Let (s(tn ))n=1,...,N be a discrete time signal of length N , s its DFT defined by (6), and kc an integer dividing N . We call matrix representation of s for kc periods the matrix (8) [Ms ]i+1,: = s [Ikc + kc i] ,

where Ikc + kc i denotes the sequence Ikc with kc i added to each element. This definition is illustrated by Figure 1. Problems 1 and 2 can be both reformulated within a matrix formulation. 162

In what follows, K will denote a field (e.g., K = Q, R, C) and A ∈ Kr×s a r × s matrix with entries in K. Associated with A ∈ Kr×s , we can consider the two K-linear maps: .A : K1×r −→ K1×s A. : Ks×1 −→ Kr×1 (9) λ −→ λ A, η −→ A η. If B ∈ Ks×t is such that A B = 0, then we can consider the so-called complexes of K-vector spaces (Rotman (2009)) (10) K1×r .A  K1×s .B  K1×t , (11) Kr×1  A. Ks×1  B. Kt×1 , i.e., linear maps with zero consecutive compositions, i.e.: imK (.A) := K1×r A ⊆ kerK (.B) := {µ ∈ K1×s | µ B = 0}, imK (B.) := B Kt×1 ⊆ kerK (A.) := {η ∈ Ks×1 | A η = 0}. The defect of exactness of (10) (resp., (11)) is defined by H(K1×s ) := kerK (.B)/imK (.A)   resp., H(Ks×1 ) := kerK (A.)/imK (B.) , where E/F stands for the quotient of a K-vector space E by a K-subvector space F . E/F is the K-vector space

2019 IFAC SSSC Sinaia, Romania, September 9-11, 2019

E. Hubert et al. / IFAC PapersOnLine 52-17 (2019) 82–87

defined by the residue classes π(e) for all e ∈ E, where π(e1 ) = π(e2 ) if e1 − e2 ∈ F , endowed with the operations:  π(e1 ) + π(e2 ) := π(e1 + e2 ), ∀ e1 , e2 ∈ E, ∀ k ∈ K, k π(e1 ) := π(k e1 ).

The complex (10) (resp., (11)) is said to be exact at K1×s (resp., Ks×1 ) if H(K1×s ) = 0 (resp., H(Ks×1 ) = 0), i.e., if kerK (.B) = imK (.A) (resp., kerK (A.) = imK (B.)).

 K1×s .B  K1×t An exact sequence of the form 0 means that kerK (.B) = 0, i.e., that .B is injective, whereas  0 means that the exact sequence K1×r .A  K1×s 1×s imK (.A) = kerK (0) = K , i.e., that .A is surjective. Similar comments hold for linear maps of the form (A.). A standard result on the duality of K-vector spaces and K-linear maps asserts that if (10) (resp., (11)) is an exact sequence of finite-dimensional K-vector spaces, then so is (11) (resp., (10)). In homological algebra, we say that the functor homK ( · , K) is exact, where homK (E, K) denotes the K-vector space formed by all the K-linear maps (forms) from a K-vector space E to K. See, e.g., Rotman (2009). Similarly a standard result in linear algebra asserts that if (10) (resp., (11)) is an exact sequence of finite-dimensional K-vector spaces, then, for any q ∈ N, so is Kq×r



resp., Kr×q 

.A A.

 Kq×s Ks×q 

.B B.

 Kq×t ,  Kt×q ,

(12) (13)

−→ K denotes the natural extension where .A : K of (9) defined by (.A)(Λ) := Λ A for all Λ ∈ Kq×r , and similarly for .B, A. and B. respectively. In homological algebra, we say that the tensor functor Kq×1 ⊗K · (resp., · ⊗K K1×q ) is exact, where E ⊗K F stands for the tensor product of two finite-dimensional K-vector spaces E and F . For more details, see, e.g., Rotman (2009). q×r

q×s

 K1×r .A  K1×s .B  K1×t  0 is a short If 0 exact sequence of finite-dimensional K-vector spaces, then the Euler-Poincar´e characteristic asserts that t−s+r = 0. More generally, the Euler-Poincar´e characteristic of a long exact sequence of finite-dimensional K-vector spaces

85

Let us consider the following K-vector space: kerK (.A) := {λ ∈ K1×r | λ A = 0}. Let B ∈ Kq×r be a matrix whose rows form a basis of kerK (.A). In other words, B is a full row rank matrix (i.e., µ B = 0 yields µ = 0 since the rows of B are K-linearly independent) which satisfies: kerK (.A) = imK (.B) := K1×q B. In particular, we have q = dimK (kerK (.A)). Then, we have λ y = λ A x = 0 for all λ ∈ kerK (.A). Hence, B y = 0 is a necessary condition for the solvability of (14). The next standard theorem shows that it is also sufficient. Theorem 1. For a fixed A ∈ Kr×s and a fixed y ∈ Kr×t , the system (14) is solvable, i.e., (14) admits a solution x ∈ Ks×t , iff the following compatibility condition holds B y = 0, (15) 1×q B. Then, all the where kerK (.A) = imK (.B) := K solutions of (14) are of the form (16) ∀ z ∈ Ku×t , x = E y + C z, s×u is a matrix whose columns form a basis of where C ∈ K kerK (A.) := {η ∈ Ks×1 | A η = 0}, i.e., C is a full column rank matrix (i.e., C θ = 0 yields θ = 0) which satisfies kerK (A.) = imK (C.) := C Ku×1 , s×r is a generalized inverse of A, i.e. A E A = A. and E ∈ K 4. SOLUTION OF THE EXACT PROBLEM 4.1 Reformulation of the problem and a few remarks Motivated by Problem 3, the main problem can be stated. Rank 2 decomposition problem: Given M ∈ Kn×m and D ∈ Kn×n , determine − if they exist − u ∈ Kn×1 and v1 , v2 ∈ K1×m satisfying: (17) M = u v1 + D u v 2 . Problem 3 corresponds to the case K = C, q = 2, i.e., the phase and amplitude modulation problem, u =  c1 , H Du =  c2 , v1 = m H , and v = m  . 2 1 2

Let us note: A(u) := (u D u) ∈ Kn×2 , v := (v1T v2T )T ∈ K2×m . Then, (17) can be rewritten as: 0  K1×rn .An  K1×rn−1 .An−1  . . . .A1  K1×r0 0 A(u) v = M. (18) n asserts that i=0 (−1)i ri = 0. See e.g., Rotman (2009). Problem defined in (18) is bilinear in u and v. Now, since the duality, i.e., the functor homK ( · , K), preRemark 2. If M = 0, then u = 0 or v1 = v2 = 0 solve the serves the exactness of long exact sequences of finiteproblem. Hence, in what follows, we suppose that M = 0. dimensional K-vector spaces, the same result holds for a Remark 3. The K-vector space imK (A(u).) is generated by long exact sequence of the form: the two vectors u and D u, which shows that: A . A . A . 0. 0 Kr0 ×1  1 Kr1 ×1  2 . . .  n Krn ×1  rankK (A(u)) := dimK (imK (A(u).)) ≤ 2. Such an exact sequence always splits (Rotman (2009)), i.e., A necessary condition for the solvability of (17) is then: there exist Bi ∈ Kri−1 ×ri for i = 0, . . . , n, such that: (19) rankK (M ) ≤ 2. B0 = Bn+1 = 0, Bi Ai + Ai+1 Bi+1 = Iri , Bi+1 Bi = 0. Remark 4. If (18) is solvable with a non full row rank matrix v, then there exists α := (α1 α2 ) ∈ K1×2 such 3.2 A standard result of linear algebra that α v = α1 v1 + α2 v2 = 0, which yields:     M = D − α2 α1−1 In u v2 , if α1 = 0, r×s a r×s Let K be field (e.g. K = Q, R, C), A ∈ K    M = In − α1 α2−1 D u v1 , if α2 = 0. matrix with entries in K, and y ∈ Kr×t . We first state a standard result on the existence of solutions x ∈ Ks×t of Hence, M = 0 must satisfy rankK (M ) = 1. A necessary the following K-linear inhomogeneous system: condition for the solvability of (18) for a matrix M A x = y. (14) satisfying rankK (M ) = 2 is then that v has full row rank. 163

2019 IFAC SSSC 86 Sinaia, Romania, September 9-11, 2019

E. Hubert et al. / IFAC PapersOnLine 52-17 (2019) 82–87

Hence, we shall only consider full row rank matrix v. 4.2 Necessary condition for the solvability of Problem (17) Let us solve (18). Let L ∈ Kp×n be a matrix whose rows generate a basis of the K-vector space kerK (.M ), i.e., L is a full row rank matrix satisfying: kerK (.M ) = imK (.L) = K1×p L. Let us explicitly characterize p in terms of the rank of M : p = dimK (kerK (.M )) = n − dimK (imK (.M )) = n − dimK (imK (M.)) = n − rankK (M ).

Notice that the compatibility condition (24) depends on u, and thus we seek for u − if it exists − in the solution space of (22) so that (24) holds. Let us reinterpret (24) and (21) in a more intrinsic mathematical setting. By definition of L, we have the following exact sequence of K-vector spaces:  K1×p .L  K1×n .M  K1×m . 0 Then, L A(u) = 0 iff u satisfies (22). If so, then we get the following complex of K-vector spaces: .A(u)

(20)

Now, (18) yields:

L A(u) v = L M = 0. (21) Since v has full row rank, we get L A(u) = 0, i.e., u must satisfy the following K-linear system:  L u = 0, (22) L D u = 0. Remark 5. Since the p rows of L are K-linearly independent, the dimension of the K-vector solution space of L u = 0, i.e., kerK (L.), is n − p = rankK (M ) by (20). The dimension of the solution space of (22) is then at most rankK (M ) (exactly rankK (M ) if, e.g., D = In or L D = 0). Let us now derive an equivalent characterization of (22). By definition of L, we have the following exact sequence:  K1×p .L  K1×n .M  K1×m . 0 Applying the exact functor homK ( · , K) to it, we obtain the following dual exact sequence of K-vector spaces: 0 Kp×1  L. Kn×1  M. Km×1 . Using kerK (L.) = imK (M.), we get L u = 0 is equivalent to the existence of w ∈ Km×1 such that u = M w. Thus, the second equation of (22) is equivalent to L (D M w) = 0, which in turn is equivalent to the existence of w ∈ Km×1 such that D M w = M w , which can be rewritten as:   w = 0. (M − D M ) w Using (19) and the upper bound of Sylvester’s inequality rankK (D M ) ≤ min{rankK (D), rankK (M )}, then (22) is equivalent to: (23) rankK (M − D M ) ≤ 3. Lemma 6. With the above notations, a necessary condition on u for the existence of a solution of Problem (17) is (22) with p ≥ n − 2, or equivalently (19) and (23). 4.3 Necessary and sufficient conditions Let u be a non-trivial solution solution of (22). We can now form the matrix A(u) = (u D u) and we are then led to the study of the linear inhomogeneous system A(u) v = M .  Let L ∈ Kp ×n be a full row rank matrix whose rows form a basis of kerK (.A(u)), i.e., kerK (.A(u)) = imK (.L ), and p = dimK (kerK (.A(u))). By Theorem 1, there exists v ∈ K2×m which satisfies A(u) v = M iff the following compatibility condition holds: (24) L M = 0. 164

 K1×2 .  K1×p .L  K1×n 0 The defect of exactness of this complex at K1×n is then: H(K1×n , u) := kerK (.A(u))/imK (.L). Now, by definition of L , we have the following exact sequence of K-vector spaces:  K1×p 0 Hence, we obtain:



.L

 K1×n

.A(u)

 K1×2 .

(25)

H(K1×n , u) = kerK (.A(u))/ kerK (.M ) = imK (.L )/imK (.L). 

Since imK (.L) ⊆ imK (.L ), there exists L ∈ Kp×p such that L = L L . Since L has full row rank, we obtain the following isomorphism:  H(K1×n , u) = imK (.L )/imK (.L) ∼ = K1×p /imK (.L ). 

Hence, H(K1×n , u) = 0 iff imK (.L ) = K1×p , i.e., iff there  exists X ∈ Kp ×p such that X L = Ip , i.e., iff L admits a left inverse, which is also equivalent to the injectivity of  the K-linear map L . : Kp ×1 −→ Kp×1 .

Applying the exact functor homK ( · , K) to the exact sequence (25), we obtain the exact sequence:

A(u). 0 Kp×1  L . Kn×1  K2×1 . Then, applying the exact functor · ⊗K K1×m to the last exact sequence, we get the following exact sequence: 



A(u).

K2×m . 0 Kp×m  L . Kn×m  Hence, M belongs to imK1×m (A(u).) = A(u) K2×m , i.e., there exists v ∈ K2×m such that M = A(u) v, iff L M = 0. By definition of L, we have L M = 0, i.e., L (L M ) = 0. Therefore, L M = 0 yields L M = 0 iff L . is injective, i.e., iff L admits a left inverse, i.e., iff H(K1×n , u) ∼ =  K1×p /imK (.L ) = 0, i.e., iff: kerK (.A(u)) = kerK (.M ). (26) Problem (17) is reduced to finding 0 = u satisfying (22) such that the K-vector space H(K1×n , u) is trivial, i.e.: dimK (H(K1×n , u)) = 0. Remark 7. Let us give a direct interpretation of (26). Let us first suppose that (17) or equivalently that (18) is solvable and let us compare the K-vector spaces: kerK (.A(u)) = {λ ∈ K1×n | λ A(u) = 0}, kerK (.M ) = {λ ∈ K1×n | λ M = 0}.

We clearly have kerK (.A(u)) ⊆ kerK (.M ). Now, if we consider λ ∈ kerK (.M ), then λ A(u) v = λ M = 0, which yields λ A(u) = 0 since v is full rank rank. Thus, if a solution exists for Problem (17), then (26) holds.

2019 IFAC SSSC Sinaia, Romania, September 9-11, 2019

E. Hubert et al. / IFAC PapersOnLine 52-17 (2019) 82–87

Conversely, if u is a non-trivial solution of (22) such that  (26) holds, then K1×p L = K1×p L, which first shows  that p = p and there also exist U, V ∈ Kp×p such that L = U L and L = V L , which yields (U V −Ip ) L = 0 and (V U − Ip ) L = 0, i.e., U V = Ip and V U = Ip since both L and L have full row rank. The compatibility condition L M = 0 is thus equivalent to L M = 0, which is satisfied by definition of L. Then, Theorem 1 shows that there exists v ∈ K2×m such that A(u) v = M , which solves (17). Let us state the main result of this paper. Theorem 8. With the above notations, Problem (17) is solvable iff there exists 0 = u ∈ Kn×1 satisfying (22) and such that H(K1×n , u) = kerK (.A(u))/ kerK (.M ) = 0. Let us study the K-vector space H(K1×n , u). Applying the Euler-Poincar´e characteristic to the short exact sequence  kerK (.M ) i  kerK (.A(u)) π  H(K1×n , u)  0, 0 we obtain: dimK (H(K1×n , u)) = dimK (kerK (.A(u)) − dimK (kerK (.M )) = p − p.

(27)

Let us now characterize p . Considering the following two short exact sequences of K-vector spaces  kerK (.A(u))  imK (.A(u))  K1×n  0, 0 0 imK (A(u).)  K2×1  kerK (A(u).)  0, the Euler-Poincar´e characteristic then yields:  dimK (imK (.A(u))) = n − dimK (kerK (.A(u))) = n − p , dimK (imK (A(u).)) = 2 − dimK (kerK (A(u).)).

Thus, we obtain p = n − 2 + dimK (kerK (A(u).)). By (20), we have p = n − rankK (M ). Hence, we obtain: H(K1×n , u) = 0 ⇔ p = p ⇔ dimK (kerK (A(u).)) = 2 − rankK (M ).

Corollary 9. With the above notations, Problem (17) is solvable iff there exists 0 = u ∈ Kn×1 satisfying (22) and one of the following two equivalent conditions holds: (1) p = p, i.e., dimK (kerK (.A(u))) = dimK (kerK (.M )), (2) dimK (kerK (A(u).)) = 2 − rankK (M ). Let us now study kerK (A(u).) = {w ∈ K2×1 | A(u) w = 0}. If w = (w1 w2 ) ∈ kerK (A(u).), i.e., u w1 + D u w2 = 0, then, using u = 0, we have w = 0 if w2 = 0, or D u = −w1 w2−1 u if w2 = 0, i.e., u is an eigenvector of D with the eigenvalue −w1 w2−1 ∈ K. Hence, if u is not an eigenvector of D, then kerK (A(u).) = 0. Now, if u is an eigenvalue of D with eigenvalue λ ∈ K, then kerK (A(u).) = imK (K.), where K = (−λ 1)T , is a K-vector space of dimension 1. Corollary 10. With the above notations, Problem (17) is solvable iff there exists 0 = u ∈ Kn×1 satisfying (22) and: (1) If rankK (M ) = 2, then u is not an eigenvector of D for an eigenvalue λ ∈ K. Then, there exists a unique v ∈ K2×m satisfying A(u) v = M defined by v = E M, 2×n denotes a left inverse of A(u). where E ∈ K

165

87

(2) If rankK (M ) = 1, then u is an eigenvector of D with an eigenvalue λ ∈ K. Then, all the solutions v ∈ K2×m satisfying A(u) v = M are defined by ∀ z ∈ K1×m , v = E M + K z, where K = (−λ 1)T and E ∈ K2×n denotes a generalized inverse of A(u). Remark 11. Let M satisfy rankK (M ) = 1 and D = 0. Then, (20) shows that p = n − 1, and thus there exists a non-trivial solution u of L u = 0 when n ≥ 2, i.e., 0 = u ∈ im(M.). Hence, 2 of Corollary 10 shows that M can always be written as M = u v1 , where u ∈ Kn×1 and v1 = E M ∈ K1×m , where E u = 1. 5. CONCLUSION This paper develops a general class of problems originating from gearbox vibration analysis. Based on Hubert et al. (2018) that made a link between the fields of signal demodulation and low-rank approximation, in this paper, we studied the specific case of phase and amplitude demodulation. Rewritten it as a polynomial problem, a solution has been proposed for the exact problem along with the necessary and sufficient conditions for its solvability. In Hubert et al. (2019), the optimal problem will be studied as well as generalizations stated in a more general framework. REFERENCES Feldman, M. (2011). Hilbert transform in vibration analysis. Mechanical Systems and Signal Processing, 25(3):735–802. Feng, Z., Chu, F. and Zuo, M. J. (2011). Time–frequency analysis of time-varying modulated signals based on improved energy separation by iterative generalized demodulation. Journal of Sound and Vibration, 330(6):1225–1243. Friedlander, B. and Francos, J. M. (1995). Estimation of amplitude and phase parameters of multicomponent signals. IEEE Transactions on Signal Processing, 43(4):917–926. Huang, N. E., Shen, Z., Long, S. R., Wu, M. C., Shih, H. H., Zheng, Q., Yen, N-C., Tung C. C. and Liu, H. H. (1998). The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. In Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 454:903–995. Hubert, E., Barrau, A. and Badaoui, M. El. (2018). New multi-carrier demodulation method applied to gearbox vibration analysis. In IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2141–2145. Hubert, E., Barrau, A., Bouzidi, Y., Dagher R. and Quadrat, A. (2019). Algebraic aspects of the optimal demodulation problem. Kaiser, J. F. (1990). On a simple algorithm to calculate the energy of a signal. In IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 381–384. McFadden, P. D. (1986). Detecting fatigue cracks in gears by amplitude and phase demodulation of the meshing vibration. Journal of Vibration, Acoustics, Stress, and Reliability in Design, 108(2):165–170. Potamianos, A. and Maragos, P. (1994). A comparison of the energy operator and the Hilbert transform approach to signal and speech demodulation. Signal Processing, 37(1):95–120. Rotman, J. J. (2009). Introduction to Homological Algebra. Springer. Santhanam, B. and Maragos, P. (2000). Multicomponent AM-FM demodulation via periodicity-based algebraic separation and energybased demodulation. IEEE Transactions on Communications, 48(3):473–490. Yu, Z., Sun, Y. and Jin, W. (2016). A novel generalized demodulation approach for multi-component signals. Signal Processing, 118:188–202.