# Spectral Intrinsic Decomposition Method for Adaptive Signal Representation

Authors: Oumar Niang; Abdoulaye Thioune; Éric Deléchelle; Jacques Lemoine
Journal: ISRN Signal Processing
 (2012)

---

## Abstract

We propose a new method called spectral intrinsic decomposition (SID) for the representation of nonlinear signals. This approach is based on the spectral decomposition of partial differential equation- (PDE-) based operators which interpolate the characteristic points of a signal. The SID’s components which are the eigenvectors of these PDE interpolation operators underlie the new signal decomposition-reconstruction method. The usefulness and the efficiency of this method is illustrated, in signal reconstruction or denoising aim, in some examples using artificial and pathological signals.


---

## Body

## 1. Introduction

The signal decomposition into atoms is an popular approach in signal analysis. The Fourier representation technics and other based on wavelets and time-frequency, or time scales analysis methods [1], and recently the Empirical Mode Decomposition [2] are extensively used in signal and image processing. The objective is to understand the contents of the signal by analyzing its components. It is sometimes desirable to have these components well suited to the separation of the noise or data in some scale analysis. Sparse representations of signals have like pursuit methods [3], the Poper Orthogonal Decomposition (POD) [4], or Singular Value Decomposition (SVD) received a great deal of attentions in last recent years. The problem solved by the sparse representation is the most compact representation of a signal in terms of combination of atoms in an overcomplete dictionary. The Empirical Mode Decomposition [2] is a self adaptive decomposition method which is essentially algorithmic and can decompose a nonlinear signal into Amplitude Modulation-Frequency Modulation (AM-FM) component plus an residue. The characteristic points of a signal like local extrema are very useful in signal analysis as its shown in EMD algorithm. The interpolation of the characteristic points provides a low frequency component of a signal whose iterative extraction is the basis of the EMD sifting process.To overcome the lack of a solid theoretical framework of EMD, we have proposed an analytical approach for sifting process based on partial differential equation (PDE) in [5–8]. We give in particular a noniterative scheme to solve the coupled PDEs system for upper and lower envelopes estimation with an adequate definition of the characteristic points of the signal to be decomposed, see [6, 8]. In this following work, we use all the eigenvectors of upper and lower PDE-envelope operators and propose a new Spectral Intrinsic Decomposition (SID) method for nonlinear signal representation. The decomposition obtained with SID acts like a sparse representation and provides relevant results for signal denoising with an competitive rate of reconstruction. In the following, we introduce in Section 2 the PDE-interpolation operator and its spectral numerical resolution. Section 3 describes the Spectral Intrinsic Decomposition method with an discussion on the SID’s component properties in Section 3.2. In Section 4, some comparison of tests are performed between SID-based method for signal reconstruction and the wavelets-based one. Finally, conclusion and perspectives are given in Section 5.


## 2. PDE-Interpolator Decomposition

As proposed in [5–7], the upper (s+) or lower (s-) envelopes of an signal s0 can be computed as the asymptotic solution of a coupled PDEs system as the following:
(1)∂s±(x,t)∂t+g±(x,t)(α∂2s±(x,t)∂x2+(1-α)∂4s±(x,st)∂x4)=0,
where α is the tension parameter which ranges from 0 to 1.The initial value solution of this equation iss±(x,t=0)=s0, and g± are the stopping or diffusivity functions depending on signal derivatives, with 0≤g±≤1. An diffusivity function for Maximum Curvature Points (MCP) interpolation of s0 is given by
(2)g±(x)=19[|sgn(∂3s0(x)∂x3)|±sgn(∂2s0(x)∂x2)+1]2,
where sgn denote the sign function.


### 2.1. Spectral Resolution of the Coupled PDEs System

Numerical resolutions for coupled PDEs system in (1) are implemented in [5] via classical iterative Crank-Nicolson or Du Fort and Frankel schemes.Equation (1) can be resolved numerically in its discrete implicit unconditionally stable scheme as follows:
(3)Sk+1=Sk+ΔtASk+1,S0=S0,
where S=(s[1],…,s[N])T is the column vector of signal samples for upper or lower envelopes for example, s+ or s-. The time step is denoted by Δt and A is a matrix formed with finite difference approximation coefficients of second and fourth order differential operators (resp., D2 and D4), as
(4)A=G(αD2-(1-α)D4),
with G the diagonal matrix of stopping function values g±[n] constructed with discrete version of stopping function values g(x) as below:
(5)g(δx3s0,δx2s0)=g(D1D2s0,D2s0),
where D2z=D+D-z and D1z=m(D+z,D-z) with D+ and D- forward and backward first difference operators on the x dimension, and where m(a,b) stands for the minmod limiter [6, 7] for derivatives estimation, m(a,b)=0.5(sgna+sgnb)min(|a|,|b|).So the explicit form leads to the following numerical resolution:(6)Sk+1=(I-ΔtA)-1Sk,S0=S0,
with I the identity matrix. Finally (1) can be decomposed into a linear system from implicit numerical scheme (6) by
(7)Sk+1=L-1Sk,S0=S0,k≥0,
where L is the linear operator including stopping function values and differential operator formed by fourth-order and second-order derivative. So, referring to numerical schemes (6), L is given by
(8)L=I-ΔtA.
The operator matrix L, has real-valued eigenvalues that are always greater or equal to 1. Then, eigenvalues, λn, of L-1 are always smaller or equal to 1, (0<λn≤1), see [7]. In Figure 1(d), the sequence of eigenvalues is plotted for an given tested signal.Input chirp signal in (a), some eigenvectors in (b) and all the eigenvalues for the upper envelope operator in (c). A similar result is obtained for lower envelope. Around the intermittency at the coordinate400, each PSMF presents as nonstationarity but contributes everywhere else to the signal composition with a centred and stationary component which are Amplitude Modulation-Frequency Modulation (AM-FM). In (d) all the eigenvalues are in the segment ]0;1] and greater than 1/17.
(a)
(b)
(c)
(d)


### 2.2. The Asymptotic Solution as an Linear Combination of Fixed Vector Point of Upper and Lower Envelope Operators

The iterative scheme (7) can be rewrite in term of initial solution S0 as
(9)Sk=(L-1)kS0,k≥1,
after convergence, the asymptotic solution, S∞, is given by
(10)S∞=(L-1)∞S0.LetV be a matrix of L-1’s sequence of eigenvectors (Vn) and D a diagonal matrix having L-1’s sequence of eigenvalues (λn), at the diagonal. So we have the following decomposition:
(11)L-1=VDV-1.
It is easy to see that
(12)(L-1)k=(VDV-1)k=VDkV-1.
So, the asymptotic solution in (10) is obtained by
(13)S∞=(VD∞V-1)S0.The asymptotic eigenvalue matrixD∞ is a diagonal matrix with eigenvalues λn∞=1 only at loci where matrix G is zeroed, and λn∞=0 where g[n]>0, for example, for λn<1.In the following,E denotes either the upper or lower envelope operator. The upper and lower envelope of the signal are calculated with the eigenvectors associated to eigenvalue λ=1. Hence, as its shown in (13), S∞ is a linear combination of 1-eigenvectors weighted by the signal amplitude. Instead of focusing only on the envelope calculus, we now consider all the set of eigenvalues of the envelope operator E.


## 2.1. Spectral Resolution of the Coupled PDEs System

Numerical resolutions for coupled PDEs system in (1) are implemented in [5] via classical iterative Crank-Nicolson or Du Fort and Frankel schemes.Equation (1) can be resolved numerically in its discrete implicit unconditionally stable scheme as follows:
(3)Sk+1=Sk+ΔtASk+1,S0=S0,
where S=(s[1],…,s[N])T is the column vector of signal samples for upper or lower envelopes for example, s+ or s-. The time step is denoted by Δt and A is a matrix formed with finite difference approximation coefficients of second and fourth order differential operators (resp., D2 and D4), as
(4)A=G(αD2-(1-α)D4),
with G the diagonal matrix of stopping function values g±[n] constructed with discrete version of stopping function values g(x) as below:
(5)g(δx3s0,δx2s0)=g(D1D2s0,D2s0),
where D2z=D+D-z and D1z=m(D+z,D-z) with D+ and D- forward and backward first difference operators on the x dimension, and where m(a,b) stands for the minmod limiter [6, 7] for derivatives estimation, m(a,b)=0.5(sgna+sgnb)min(|a|,|b|).So the explicit form leads to the following numerical resolution:(6)Sk+1=(I-ΔtA)-1Sk,S0=S0,
with I the identity matrix. Finally (1) can be decomposed into a linear system from implicit numerical scheme (6) by
(7)Sk+1=L-1Sk,S0=S0,k≥0,
where L is the linear operator including stopping function values and differential operator formed by fourth-order and second-order derivative. So, referring to numerical schemes (6), L is given by
(8)L=I-ΔtA.
The operator matrix L, has real-valued eigenvalues that are always greater or equal to 1. Then, eigenvalues, λn, of L-1 are always smaller or equal to 1, (0<λn≤1), see [7]. In Figure 1(d), the sequence of eigenvalues is plotted for an given tested signal.Input chirp signal in (a), some eigenvectors in (b) and all the eigenvalues for the upper envelope operator in (c). A similar result is obtained for lower envelope. Around the intermittency at the coordinate400, each PSMF presents as nonstationarity but contributes everywhere else to the signal composition with a centred and stationary component which are Amplitude Modulation-Frequency Modulation (AM-FM). In (d) all the eigenvalues are in the segment ]0;1] and greater than 1/17.
(a)
(b)
(c)
(d)


## 2.2. The Asymptotic Solution as an Linear Combination of Fixed Vector Point of Upper and Lower Envelope Operators

The iterative scheme (7) can be rewrite in term of initial solution S0 as
(9)Sk=(L-1)kS0,k≥1,
after convergence, the asymptotic solution, S∞, is given by
(10)S∞=(L-1)∞S0.LetV be a matrix of L-1’s sequence of eigenvectors (Vn) and D a diagonal matrix having L-1’s sequence of eigenvalues (λn), at the diagonal. So we have the following decomposition:
(11)L-1=VDV-1.
It is easy to see that
(12)(L-1)k=(VDV-1)k=VDkV-1.
So, the asymptotic solution in (10) is obtained by
(13)S∞=(VD∞V-1)S0.The asymptotic eigenvalue matrixD∞ is a diagonal matrix with eigenvalues λn∞=1 only at loci where matrix G is zeroed, and λn∞=0 where g[n]>0, for example, for λn<1.In the following,E denotes either the upper or lower envelope operator. The upper and lower envelope of the signal are calculated with the eigenvectors associated to eigenvalue λ=1. Hence, as its shown in (13), S∞ is a linear combination of 1-eigenvectors weighted by the signal amplitude. Instead of focusing only on the envelope calculus, we now consider all the set of eigenvalues of the envelope operator E.


## 3. The Spectral Intrinsic Decomposition Method

Let us consider all the eigenvalues of the envelope operatorE of a signal s0. The set of eigenvectors of E is an pseudo-dictionary, each of its component called Spectral Proper Mode Function (SPMF), is intrinsic to the signal. The SID decomposition deals to the combination of all SPMF.


### 3.1. On the Properties of SPMF and Spectral Intrinsic Decomposition Principle

The atoms SPMF calculated from the operatorE are adaptive and well localized around the caracteristic points of the signal. In Figure 1, we show the original signal in Figure 1(a), and some eigenvectors (SPMFs V920,V940,V960,V980,V100) associated to lowest eigenvalues for the upper envelope in Figure 1(b). Figure 1(c) presents the SPMFs numbers V20, V40, V60, V80, and V100. Around the intermittency at the coordinate 400, the last SPMFs corresponding to lesser eigenvalues present an nonstationarity but contribute everywhere else to the signal composition, with a centred and stationary component which is Amplitude Modulation-Frequency Modulation (AM-FM).It is interesting to note that the EMD Sifting process [2] allows for tracking of these AM-FM components by searching iteratively around the extrema. Even in SID the role of extrema in the occurrence of nonstationarity is noted. Also, the SPMF contains local frequencies of the signal. Hence localy, the SPMF decomposition (SID) works like the basic EMD’s principle that considers a signal as an superposition of a lower frequency component and a most higher frequency component.When the classical EMDs principle leads to a locally decomposition with two components, the SID decomposition gives a sequence (of number greater than the number of characteristic points) of really localized component as follows:(14)s0=∑k∈{j/λj=1,λjeigenvalueofE}‍VkCk+∑k∉{j/λj=1,λjeigenvalueofE}‍VkCk,
where Vk denotes an eigenvector of E (SPMF) and Ck the decomposition coefficient depending on s0. The first term of (14) corresponds to the envelope of s0. Then it appears that the SID provides a generalization of the EMD’s basic principle because in (14), we have more than the number of maxima or minima components. SPMF participates in the whole dynamic of the signal with a strong localization around the points that generated the eigenvectors.In most cases, an SPMF can be viewed as a (nonlinear) frequency narrow-band waveletφ with Amplitude Modulation by a lower frequency signal a[n]:
(15)SPMFk[n]=Ak[n]φk[n].
The redundancy and the orthogonality of the dictionary of all SPMF depend on the properties of the operator E. When E is symmetric, we can have the orthogonality and SPMF too similar to a wavelet function (see example in Figure 1).


### 3.2. The Spectral Intrinsic Decomposition

The Spectral Intrinsic Decomposition procedure is define as the calculus of all the SPMF for an given signal. Let us take the same notation than in Section3.1 and consider the upper envelope operator E=L-1. The same procedure can be performed with the lower envelope. The eigen decomposition of E gives: [VE,LE]=eig(E), where VE=[V1,…,Vsize(s0)] and LE=[L1,…,Lsize(s0)] (with the possibility zeros to complete the size of the vector) are, respectively, the set of eigenvectors and the set of eigenvalues of E. The coefficient reconstruction of s0 is given by:
(16)C=LEVE-1s0⊥,withs0⊥thetransposeofs0.
Hence s0 is computed by the formula s0=VC.The Spectral Intrinsic Decomposition ofs0 described in Algorithm 1, is given as follows:
(17)s0=∑k=1N‍VkCk.Algorithm 1:Spectral intrinsic decomposition algorithm.
(
1
)  compute  g± from S0, using for example MCP (2)
(
2
)  compute matrix operator L-1=E (8)
(
3
)  perform eigendecomposition of E, [VE,LE]=eig(E)
(
4
)  perform reconstruction coefficients of S0, C=LEVE-1S0-1
(
5
)  set  [Vk], and [Lk] for k=1…N,               ⊳Result
andS0←∑k=1N‍Vk*CkThis decomposition is intrinsic and depends only to the position of caracteristic points ofs0 that define the diffusivity function in the interpolation operator. We notice that the SID versus lower envelope works like the SID with the upper envelope and has the same reconstruction ability.The reason is that the PDE-interpolation operator use all the data ins0 and these SPMFs generate the same functional space. All the SPMF participate locally in the reconstruction of the signal s0. Hence, in the sense of the superposition principle, SID is more general than the EMD’s classical principe.


## 3.1. On the Properties of SPMF and Spectral Intrinsic Decomposition Principle

The atoms SPMF calculated from the operatorE are adaptive and well localized around the caracteristic points of the signal. In Figure 1, we show the original signal in Figure 1(a), and some eigenvectors (SPMFs V920,V940,V960,V980,V100) associated to lowest eigenvalues for the upper envelope in Figure 1(b). Figure 1(c) presents the SPMFs numbers V20, V40, V60, V80, and V100. Around the intermittency at the coordinate 400, the last SPMFs corresponding to lesser eigenvalues present an nonstationarity but contribute everywhere else to the signal composition, with a centred and stationary component which is Amplitude Modulation-Frequency Modulation (AM-FM).It is interesting to note that the EMD Sifting process [2] allows for tracking of these AM-FM components by searching iteratively around the extrema. Even in SID the role of extrema in the occurrence of nonstationarity is noted. Also, the SPMF contains local frequencies of the signal. Hence localy, the SPMF decomposition (SID) works like the basic EMD’s principle that considers a signal as an superposition of a lower frequency component and a most higher frequency component.When the classical EMDs principle leads to a locally decomposition with two components, the SID decomposition gives a sequence (of number greater than the number of characteristic points) of really localized component as follows:(14)s0=∑k∈{j/λj=1,λjeigenvalueofE}‍VkCk+∑k∉{j/λj=1,λjeigenvalueofE}‍VkCk,
where Vk denotes an eigenvector of E (SPMF) and Ck the decomposition coefficient depending on s0. The first term of (14) corresponds to the envelope of s0. Then it appears that the SID provides a generalization of the EMD’s basic principle because in (14), we have more than the number of maxima or minima components. SPMF participates in the whole dynamic of the signal with a strong localization around the points that generated the eigenvectors.In most cases, an SPMF can be viewed as a (nonlinear) frequency narrow-band waveletφ with Amplitude Modulation by a lower frequency signal a[n]:
(15)SPMFk[n]=Ak[n]φk[n].
The redundancy and the orthogonality of the dictionary of all SPMF depend on the properties of the operator E. When E is symmetric, we can have the orthogonality and SPMF too similar to a wavelet function (see example in Figure 1).


## 3.2. The Spectral Intrinsic Decomposition

The Spectral Intrinsic Decomposition procedure is define as the calculus of all the SPMF for an given signal. Let us take the same notation than in Section3.1 and consider the upper envelope operator E=L-1. The same procedure can be performed with the lower envelope. The eigen decomposition of E gives: [VE,LE]=eig(E), where VE=[V1,…,Vsize(s0)] and LE=[L1,…,Lsize(s0)] (with the possibility zeros to complete the size of the vector) are, respectively, the set of eigenvectors and the set of eigenvalues of E. The coefficient reconstruction of s0 is given by:
(16)C=LEVE-1s0⊥,withs0⊥thetransposeofs0.
Hence s0 is computed by the formula s0=VC.The Spectral Intrinsic Decomposition ofs0 described in Algorithm 1, is given as follows:
(17)s0=∑k=1N‍VkCk.Algorithm 1:Spectral intrinsic decomposition algorithm.
(
1
)  compute  g± from S0, using for example MCP (2)
(
2
)  compute matrix operator L-1=E (8)
(
3
)  perform eigendecomposition of E, [VE,LE]=eig(E)
(
4
)  perform reconstruction coefficients of S0, C=LEVE-1S0-1
(
5
)  set  [Vk], and [Lk] for k=1…N,               ⊳Result
andS0←∑k=1N‍Vk*CkThis decomposition is intrinsic and depends only to the position of caracteristic points ofs0 that define the diffusivity function in the interpolation operator. We notice that the SID versus lower envelope works like the SID with the upper envelope and has the same reconstruction ability.The reason is that the PDE-interpolation operator use all the data ins0 and these SPMFs generate the same functional space. All the SPMF participate locally in the reconstruction of the signal s0. Hence, in the sense of the superposition principle, SID is more general than the EMD’s classical principe.


## 4. Application to Signal Reconstruction and Filtering

The SID-based signal denoising principle is derived from the idea that regular signal can be accurately approximated using a small number of approximation coefficients (at a suitably chosen level) or some of the detail coefficients. The SPMF corresponding to smaller eigenvalues contain the noise or the highest frequency component of any given signal. The denoising procedure contains three steps: decomposition-threshold-reconstruction (DTR). The decomposition classically can be a wavelets representation depending on a choice of wavelet and a levelNL. For each level from 1 to NL, a threshold is selected and hard thresholding is applied to the detail coefficients. At last we compute wavelet reconstruction using the original approximation coefficients of level NL and the modified detail coefficients of levels from 1 to N. In following we propose a new reconstruction and denoising method by taking the Spectral Intrinsic instead of the wavelet at the decomposition step of the DTR procedure. There are two denoising approaches that can be explored here for test results. The first consists of taking the wavelet expansion of the signal and keeping the largest absolute value coefficients. In this case, one can set a global threshold and evaluate the denoising performance by the signal-to-noise ratio (SNR), or the relative error of reconstruction. Thus, only a single parameter needs to be selected. The second approach consists of applying visually determined level-dependent thresholds. In Algorithm 2, we have described the SID-based reconstruction method as follows: with the SID of s0, we first compute the set of all reconstruction coefficients by
(18)C=LEVE-1s0⊥,
and retain all the SPMF corresponding to the 1-eigenvalues of E that gives the significant SPMF in SID representation of s0, by calculating the coefficients for 1-eigenvalues denoted by Cm. After, like in the wavelet reconstruction or denoising method, we fix the number of supplementary SPMF (by estimating Cs) to add to the significant set of SPMF retained before. Secondly we reconstruct an approximate and denoised version of s0 by forming the reconstruction coefficient:
(19)CR=[CmCs],
and completing by zeros for equalization of matrix size. In [7] we have proposed an optimal method to compose Cs by regularization technique. The reconstructed signal version is
(20)Sc=VECR.
Finally we compute the SNR and the number of points or number of SPMF retained for the signal reconstruction. In example tests presented in Figures 2(a), 2(b), 3(a), and 3(b) we apply a global thresholding, for a given and unoptimized wavelet choice, to produce a nearly complete square norm recovery for a signal and compare it with the SID version one. For signal 3 reconstruction, the high frequency component is lost with wavelet method while SID-reconstruction retains quite perfectly this essential component. We can clearly see in Figures 4(b), 4(c), 5(a), 5(b), and 6 that the SID-based reconstruction works better than the classical wavelet method with fewest reconstruction error and gives better SNR. Another advantage of the SID based-reconstruction method is its self-adaptability and its unique dependence on signal to be approximate. To compare scores of these methods we compute the retained energy in percentage defined by
(21)-20*log10(∥Sc-s0∥∥s0∥)
and we have compare the percentage of the number of useful points for signal reconstruction, see Figure 4(a).Algorithm 2:Signal reconstruction with spectral intrinsic decomposition approach.
(
1
)  compute matrix operator L-1=E (8)
(
2
)  perform eigendecomposition of E,[VE,LE]=eig(E)
(
3
)  compute reconstruction's coefficients:
C=LEVE-1s0⊥
(
4
)  find retained significant eigenvectors associated to eigenvalue equal to 1,Cm
(
5
)  fix the number of supplementary eigenvectors to add Cs
(
6
)  find retained significant eigenvectors CR=[CmCs]
(
7
)  compute reconstructed signal
Sc=←VECR;
⊳Results
and SNR=←-20*log10(∥Sc-S0∥∥s0∥)Original Signals signal 1 in (a) and signal 2 in (b), their reconstructed version by Wavelets methods and SID-based method. The SID approach works like Wavelets methods with the advantage of the skip of the choice of the mother wavelet.
(a)
(b)Original Signals signal 3 in (a) and signal 4 in (b), their reconstructed version by Wavelets methods and SID-based method. The SID approach works better than wavelets methods with a stronger power of reconstruction and the advantage of the skip of wavelet form choice.
(a)
(b)In (a): The table of comparative results showing the performance of SID. In (b): graphical representation of all the errors. In (c). The SID-based reconstruction method works better than the wavelets one.
(a)
(b)
(c)In (a) the same result with two different scale representations. In (b) the comparative representation of the SNR. SID-based reconstruction method gives satisfactory results compared to wavelet-based method.
(a)
(b)Original Signals. Signal 5, signal 6, and signal 7, respectively, in (a), (b), and (c); their reconstructed version by Wavelets methods and SID-based method. According to the comparison table of SNR and errors in (d), the SID approach works as well as Wavelets methods with the advantage of adaptability.
(a)
(b)
(c)
(d)


## 5. Conclusion

In this paper we have introduced a new decomposition method based on a spectral decomposition of an intrinsic interpolation operator of a signal. The new SID method is self adaptive and works more generally than the EMD basic principle. The SID gives a dictionary in terms of SPMF that are similar to atoms in sparse representation. The test results demonstrate that SID can be used in signal denoising as much as the wavelets technic, with the advantage of self adaptability. The SID is also suitable for signal compression, one issue of our future works.



---

*Source: 457152-2012-12-13.xml*