Spectral Intrinsic Decomposition Method for Adaptive Signal Representation

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.

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 (𝑠⁺) or lower (𝑠⁻) envelopes of an signal 𝑠₀ can be computed as the asymptotic solution of a coupled PDEs system as the following: 𝜕𝑠^±(𝑥,𝑡)𝜕𝑡+𝑔^±()𝛼𝜕𝑥,𝑡²𝑠^±(𝑥,𝑡)𝜕𝑥²+()𝜕1−𝛼⁴𝑠^±(𝑥,s𝑡)𝜕𝑥⁴=0, (1) where 𝛼 is the tension parameter which ranges from 0 to 1.

The initial value solution of this equation is 𝑠^±(𝑥,𝑡=0)=𝑠₀, and 𝑔^± are the stopping or diffusivity functions depending on signal derivatives, with 0≤𝑔^±≤1. An diffusivity function for Maximum Curvature Points (MCP) interpolation of 𝑠₀ is given by 𝑔^±(𝑥)=19𝜕sgn³𝑠₀(𝑥)𝜕𝑥³𝜕±sgn²𝑠₀(𝑥)𝜕𝑥²+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: 𝐒^𝑘+1=𝐒^𝑘+Δ𝑡𝐀𝐒^𝑘+1,𝐒⁰=𝐒₀, (3) where 𝐒=(𝑠[1],…,𝑠[𝑁])^T is the column vector of signal samples for upper or lower envelopes for example, 𝑠⁺ or 𝑠⁻. The time step is denoted by Δ𝑡 and 𝐀 is a matrix formed with finite difference approximation coefficients of second and fourth order differential operators (resp., 𝐃₂ and 𝐃₄), as 𝐀=𝐆𝛼𝐃₂−(1−𝛼)𝐃₄, (4) with 𝐆 the diagonal matrix of stopping function values 𝑔^±[𝑛] constructed with discrete version of stopping function values 𝑔(𝑥) as below: 𝑔𝛿³_𝑥𝑠₀,𝛿²_𝑥𝑠₀=𝑔𝐷₁𝐷₂𝑠₀,𝐷₂𝑠₀, (5) where 𝐷₂𝑧=𝐷⁺𝐷⁻𝑧 and 𝐷₁𝑧=𝐦(𝐷⁺𝑧,𝐷⁻𝑧) with 𝐷⁺ and 𝐷⁻ forward and backward first difference operators on the 𝑥 dimension, and where 𝐦(𝑎,𝑏) stands for the minmod limiter [6, 7] for derivatives estimation, 𝐦(𝑎,𝑏)=0.5(sgn𝑎+sgn𝑏)min(|𝑎|,|𝑏|).

So the explicit form leads to the following numerical resolution: 𝐒^𝑘+1=(𝐈−Δ𝑡𝐀)⁻¹𝐒^𝑘,𝐒⁰=𝐒₀, (6) with 𝐈 the identity matrix. Finally (1) can be decomposed into a linear system from implicit numerical scheme (6) by 𝐒^𝑘+1=𝐋⁻¹𝐒^𝑘,𝐒⁰=𝐒₀,𝑘≥0, (7) where 𝐋 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), 𝐋 is given by 𝐋=𝐈−Δ𝑡𝐀. (8) The operator matrix 𝐋, has real-valued eigenvalues that are always greater or equal to 1. Then, eigenvalues, 𝜆_𝑛, of 𝐋⁻¹ are always smaller or equal to 1, (0<𝜆_𝑛≤1), see [7]. In Figure 1(d), the sequence of eigenvalues is plotted for an given tested signal.

(a)

(b)

(c)

(d)

Figure 1 

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 coordinate 400, 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.

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 𝐒₀ as 𝐒^𝑘=𝐋^−1𝑘𝐒₀,𝑘≥1, (9) after convergence, the asymptotic solution, 𝐒_∞, is given by 𝐒_∞=𝐋^−1∞𝐒₀. (10)

Let 𝐕 be a matrix of 𝐋⁻¹’s sequence of eigenvectors (𝐕_𝑛) and 𝓓 a diagonal matrix having 𝐋⁻¹’s sequence of eigenvalues (𝜆_𝑛), at the diagonal. So we have the following decomposition: 𝐋⁻¹=𝐕𝓓𝐕⁻¹. (11) It is easy to see that 𝐋^−1𝑘=𝐕𝓓𝐕^−1𝑘=𝐕𝓓^𝑘𝐕⁻¹. (12) So, the asymptotic solution in (10) is obtained by 𝐒_∞=𝐕𝓓^∞𝐕⁻¹𝐒₀. (13)

The asymptotic eigenvalue matrix 𝓓^∞ is a diagonal matrix with eigenvalues 𝜆^∞_𝑛=1 only at loci where matrix 𝐆 is zeroed, and 𝜆^∞_𝑛=0 where 𝑔[𝑛]>0, for example, for 𝜆_𝑛<1.

In the following, 𝐄 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), 𝐒_∞ 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 𝐄.

3. The Spectral Intrinsic Decomposition Method

Let us consider all the eigenvalues of the envelope operator 𝐄 of a signal 𝑠₀. The set of eigenvectors of 𝐄 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 operator 𝐄 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 𝐕₉₂₀,𝐕₉₄₀,𝐕₉₆₀,𝐕₉₈₀,𝐕₁₀₀) associated to lowest eigenvalues for the upper envelope in Figure 1(b). Figure 1(c) presents the SPMFs numbers 𝐕₂₀, 𝐕₄₀, 𝐕₆₀, 𝐕₈₀, and 𝐕₁₀₀. 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: 𝐬₀=_{𝑘∈𝑗/𝜆𝑗=1,𝜆𝑗eigenvalueof𝐄}𝐕_𝐤𝐂_𝑘+_{𝑘∉𝑗/𝜆𝑗=1,𝜆𝑗eigenvalueof𝐄}𝐕_𝐤𝐂_𝑘, (14) where 𝐕_𝐤 denotes an eigenvector of 𝐄 (SPMF) and 𝐂_𝑘 the decomposition coefficient depending on 𝐬₀. The first term of (14) corresponds to the envelope of 𝐬₀. 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 𝑎[𝑛]: 𝐒𝐏𝐌𝐅_𝑘[𝑛]=𝐀_𝑘[𝑛]𝜑_𝑘[𝑛]. (15) The redundancy and the orthogonality of the dictionary of all SPMF depend on the properties of the operator 𝐄. When 𝐄 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 Section 3.1 and consider the upper envelope operator 𝐄=𝐋⁻¹. The same procedure can be performed with the lower envelope. The eigen decomposition of 𝐄 gives: [𝐕_𝐄,𝐋_𝐄]=eig(𝐄), where 𝐕_𝐄=[𝐕_𝟏,…,𝐕_{size(𝐬𝟎)}] and 𝐋_𝐄=[𝐋_𝟏,…,𝐋_{size(𝐬𝟎)}] (with the possibility zeros to complete the size of the vector) are, respectively, the set of eigenvectors and the set of eigenvalues of 𝐄. The coefficient reconstruction of 𝐬₀ is given by: 𝐂=𝐋_𝐄𝐕_𝐄⁻¹𝐬^⟂_𝟎,with𝐬^⟂_𝟎thetransposeof𝐬₀. (16) Hence 𝐬₀ is computed by the formula 𝐬₀=𝐕𝐂.

The Spectral Intrinsic Decomposition of 𝐬₀ described in Algorithm 1, is given as follows: 𝐬₀=^𝑁_𝑘=1𝐕_𝐤𝐂_𝑘. (17)

(1)  compute  𝑔± from 𝐒0, using for example MCP (2) (2)  compute matrix operator 𝐋−1=𝐄 (8) (3)  perform eigendecomposition of 𝐄, [𝐕𝐄,𝐋𝐄]=eig(𝐄) (4)  perform reconstruction coefficients of 𝑆0, 𝐂=𝐋𝐄𝐕𝐄−1𝑆0−1 (5)  set  [𝐕𝐤], and [𝐋𝐤] for 𝑘=1…𝑁,               ⊳𝑅𝑒𝑠𝑢𝑙𝑡    and 𝐒𝟎←∑𝑁𝑘=1𝐕𝐤∗𝐂𝐤

Algorithm 1 

Spectral intrinsic decomposition algorithm.

This decomposition is intrinsic and depends only to the position of caracteristic points of 𝐬₀ 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 in 𝐬₀ and these SPMFs generate the same functional space. All the SPMF participate locally in the reconstruction of the signal 𝐬₀. 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 level 𝐍_𝐿. For each level from 1 to 𝐍_𝐿, 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 𝐍_𝐿 and the modified detail coefficients of levels from 1 to 𝑁. 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 𝐬₀, we first compute the set of all reconstruction coefficients by 𝐂=𝐋_𝐄𝐕_𝐄⁻¹𝐬^⟂_𝟎, (18) and retain all the SPMF corresponding to the 1-eigenvalues of 𝐄 that gives the significant SPMF in SID representation of 𝐬₀, by calculating the coefficients for 1-eigenvalues denoted by 𝐂_𝐦. After, like in the wavelet reconstruction or denoising method, we fix the number of supplementary SPMF (by estimating 𝐂_𝐬) to add to the significant set of SPMF retained before. Secondly we reconstruct an approximate and denoised version of 𝐬₀ by forming the reconstruction coefficient: 𝐂_𝐑=𝐂_𝐦𝐂_𝐬, (19) and completing by zeros for equalization of matrix size. In [7] we have proposed an optimal method to compose 𝐂_𝐬 by regularization technique. The reconstructed signal version is 𝐒𝐜=𝐕_𝐄𝐂_𝐑. (20) 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 −20∗log₁₀𝐒𝑐−𝐬₀𝐬₀ (21) and we have compare the percentage of the number of useful points for signal reconstruction, see Figure 4(a).

(1)  compute matrix operator 𝐋−𝟏=𝐄 (8) (2)  perform eigendecomposition of 𝐄,[𝐕𝐸,𝐋𝐸]=eig(𝐄) (3)  compute reconstruction's coefficients:    𝐂=𝐋𝐸𝐕𝐸−1𝐬⟂0 (4)  find retained significant eigenvectors associated to eigenvalue equal to 1,𝐂𝐦 (5)  fix the number of supplementary eigenvectors to add 𝐂𝑠 (6)  find retained significant eigenvectors 𝐂𝑅=[𝐂𝐦𝐂𝐬]   (7)  compute reconstructed signal    𝐒𝐜=←𝐕𝐄𝐂𝐑;                                     ⊳𝑅𝑒𝑠𝑢𝑙𝑡𝑠     and 𝐒𝐍𝐑=←−20∗log10‖𝐒𝑐−𝐒0‖‖𝐬0‖

Algorithm 2 

Signal reconstruction with spectral intrinsic decomposition approach.

(a)

(b)

Figure 2 

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)

Figure 3 

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)

(c)

Figure 4 

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)

Figure 5 

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)

(c)

(d)

Figure 6 

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.

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.