Three-Mode Biomedical Signal Denoising in the Local Phase Space based on a Tensor Approach

The task of noise reduction is a central theme in a wide variety of fields. As the biomedical signals and the random noise often have overlapping bandwidths, the conventional methods based on the spectrum analysis did not work well for this data. Due to its simplicity in implementation and efficiency in computation the nonlinear phase-space projection technique together with singular value decomposition or approximate joint diagonalization a set of time-delayed covariance matrices procedure is an essential tool in noise reduction [1–5], signal detection [7, 8] and biomedical signal processing [9] algorithms. Several phase space projection methods, based on subspace decomposition, were proposed for application to the problem of additive noise reduction in the context of phase space analysis – the global projections method [2], [5] and the local (nearest neighborhoods) phase spaces method [1–4], [6].The local projection approach project the data in the neighborhood onto the hyper-plane and bring the flow pattern deviated back to the real dynamics system [6]. This approach has achieved nice noise reduction effects and has been applied to the speech signals, biomedical signals, mechanical vibration signals, etc. The local models process data in the vicinity of local neighborhoods leading to more detailed and possibly more complex models in comparison with the global projection methods. On the other hand, the local projection noise reduction approach is influenced by the neighborhood selection greatly and usually, neighborhoods merge if all data are contaminated by large amounts of noise, therefore these methods performed well with small or moderate amounts of noise. The objective of this paper is to investigate the denoising performance of an improved local projection noise reduction approach based on 2D model of neighbors. The neighborhood tensor of 2D neighbors – windows with several consecutive vectors of reconstructed phase-space – is computed rather than neighborhood matrix of each vector. Tensor approach is compared with its matrixvalued counterpart, which requires stacking the 2D neighbors into one highly structured neighborhood matrix.


Introduction
The task of noise reduction is a central theme in a wide variety of fields.As the biomedical signals and the random noise often have overlapping bandwidths, the conventional methods based on the spectrum analysis did not work well for this data.Due to its simplicity in implementation and efficiency in computation the nonlinear phase-space projection technique together with singular value decomposition or approximate joint diagonalization a set of time-delayed covariance matrices procedure is an essential tool in noise reduction [1][2][3][4][5], signal detection [7,8] and biomedical signal processing [9] algorithms.Several phase space projection methods, based on subspace decomposition, were proposed for application to the problem of additive noise reduction in the context of phase space analysis -the global projections method [2], [5] and the local (nearest neighborhoods) phase spaces method [1][2][3][4], [6].The local projection approach project the data in the neighborhood onto the hyper-plane and bring the flow pattern deviated back to the real dynamics system [6].This approach has achieved nice noise reduction effects and has been applied to the speech signals, biomedical signals, mechanical vibration signals, etc.The local models process data in the vicinity of local neighborhoods leading to more detailed and possibly more complex models in comparison with the global projection methods.On the other hand, the local projection noise reduction approach is influenced by the neighborhood selection greatly and usually, neighborhoods merge if all data are contaminated by large amounts of noise, therefore these methods performed well with small or moderate amounts of noise.
The objective of this paper is to investigate the denoising performance of an improved local projection noise reduction approach based on 2D model of neighbors.The neighborhood tensor of 2D neighbors -windows with several consecutive vectors of reconstructed phase-spaceis computed rather than neighborhood matrix of each vector.Tensor approach is compared with its matrix-valued counterpart, which requires stacking the 2D neighbors into one highly structured neighborhood matrix.

Local Projection Noise Reduction Algorithm
A noisy time series ሼ‫ݔ‬ ሽ ୀଵ could be reconstructed to a m-dimensional phase space by selecting the embedding dimension m and time delay ߬, each phase point in the phase space is defined by [10] where ݅ ൌ ͳǡʹǡ ‫ڮ‬ ǡ ‫ܮ‬ െ ሺ݉ െ ͳሻ߬, and ሺȉሻ ் denotes the transpose of a real matrix.At ߬ ൌ ͳ the reconstructed phase space matrix with m rows and ‫ܯ‬ ൌ ‫ܮ‬ െ ݉ ͳ columns (called a trajectory matrix) is defined by The process of finding the neighbors of each point in phase space is the most expensive step in most nonlinear noise reduction algorithms.Usually, the near neighborhood of the reference point is defined as where ߝ -the size of the neighborhood.It becomes a nontrivial problem to identify the correct neighbors if all data are contaminated by large amounts of noise.When the neighborhood radius is too small, the neighborhood is submerged in noise and the fitting direction of this region is nearly random.And while the neighborhood radius is oversized the attractor manifold was also distorted.The high dimension of the embedding space helps partly to identify neighbors also for rather high noise levels [3,7].Furthermore, the neighborhood selection algorithm is proposed as follows: (i) find a given number K of nearest neighbors rather than a neighborhood of fixed radius, (ii) define the near neighborhood of the each no overlapping or overlapping windows ‫܆‬ ‫א‬ ԧ ൈ having P columns ሼx ሽ ୀ ାିଵ of the trajectory matrix rather than the near neighborhood of the each point .The similarity between two matrices and is defined by In consequence the defined neighborhood can be rewritten into a three-mode format where the three modes of tensor are neighbors mode (size of the neighborhood, i. e. number of neighbors), reconstructed phase-space mode (embedding dimension) and window mode (window length, i. e. number of vectors in the window).Finally, by transpose second and third modes, the neighborhood tensor can be defined as In order to perform a decomposition of the neighborhood tensors and split the three mode data into two orthogonal subspaces -the signal and noise subspaces -the higher order singular value decomposition (HOSVD) [11,12] is used.The HOSVD is preferred over other decompositions because it provides an orthonormal basis, allowing an extension of the well-known matrix subspace technique [11].The HOSVD of a three-way array ( 6) is given by where झ ‫א‬ ԧ ൈൈ is the core tensor, which satisfies the all-orthogonality conditions [10] and ‫܃‬ ሺሻ ‫א‬ ԧ ெ ൈெ are the unitary matrices of i-mode singular vectors for i = 1, 2, 3.In ( 7) and ( 8), the notation ×i corresponds to the scalar product along the ith mode.It has been demonstrated in [12] that the estimation of the three singular matrices of a given three-mode array घ ‫א‬ ԧ ൈൈ can be performed by the estimation of left singular matrices of the three possible unfolding matrices, obtained by stacking sub arrays in large matrices This is to say that these unfolding matrices can be decomposed by singular value decomposition (SVD) into The rank of a three-mode array can be defined as ‫ݎ‪ܽ݊݇ሺ‬ݎ‬ ଵ ǡ ‫ݎ‬ ଶ ǡ ‫ݎ‬ ଷ ሻ, composed of the ranks of unfolding matrices [11], i.e. ‫ݎ‬ ൌ ‫ۼ‪ܽ݊݇൫‬ݎ‬ ሺሻ ൯.Similarly to 2D arrays, subspace methods for three-mode arrays are based on a rank approximation of the HOSVD.Consider three-mode array (6) and its decomposition into two three-mode arrays where घ ௦ describes the signal subspace and घ the noise subspace.घ can be expressed in terms of an "economy size" HOSVD in the following way [11,12]: where matrices ‫܃‬ ௦ ሺሻ ‫א‬ ԧ ெ ൈ are obtained by keeping the first ‫ݍ‬ singular vectors (i = 1, 2, 3), associated with the signal subspace.The values ‫ݍ‬ are chosen by finding an abrupt change of slope in the curves of three-mode singular values.Actually, the matrices are the projectors of the three unfolding matrices.The equation ( 12) operates directly on the measurement data and is therefore termed the "direct data approach" [13].In this work, we calculate the signal subspace based on the relation between the HOSVD-based subspace estimate and the SVD-based subspace estimate in the presence of noise [14].SVD is performed on the matrix ‫ۼ‬ ሺଷሻ் ‫א‬ ԧ ൈ , which contains all vectors of neighborhood stacked along the rows where ሾघ ሿ ሺଷሻ denoted a matrix unfolding of the tensor घ along the 3th mode.Applying the 3-mode unfolding to the tensor ट ௦ , that represents the HOSVD-based basis for the signal subspace estimate, and taking its transpose, authors [13] obtain i.e. the matrix-based subspace estimate ‫܃‬ ௦ ሺଷሻ gets premultiplied by a Kronecker product of ‫۾‬ ଵ and ‫۾‬ ଶ , which are the projection matrices ( 13) onto the subspaces spanned be the one-mode and the two-mode vectors of the tensor ट ௦ .The signal subspace घ ௦ defined by matrix unfolding of the tensor ൣघ ௦ ൧ ሺଷሻ along the last dimension can be expressed as Given a tensor as matrix object ൣघ ௦ ൧ ሺଷሻ , we can rearrange its entries back into a tensor घ ௦ , i. e. rebuilt the original data structure and define the projected trajectory matrix for reference matrix of neighborhood.Finally, an enhanced one-dimensional signal is created from the new space, typically by time-aligning and averaging the columns of the trajectory matrix d s X (see [2] for more details).

Denoising performance analyses from simulated data sets
We had applied the reviewed methods to the denoising of an x component of the Rossler system contaminated with additive white, Gaussian and independent from signal noise.The signals of Rossler system may be used for the simulation of some biomedical signals, having a pseudoperiodic character.The Rossler system was simulated having parameters a=0,398, b=2 and c=4.It has been argued that the SVD method can obtain better results for pseudoperiodic signals by over-embedding with time delay Ĳ=1 [3,7].Therefore, the embedding dimension of the reconstructed phase space m=60 (approximately one cycle of the attractor) and time delay Ĳ=1 for all signals were defined.Ten sequences (each 1000 points) are used for evaluating the performance of the denoising at various signal-to noise level (SNR).The neighborhood of each segment n of trajectory matrix is a three-mode array घ ‫א‬ ԧ ଷൈଵൈ , that has neighbors mode (the first 30 nearest neighbors are used for each reference phase window), window mode (number of vectors in the window -10), and embedding mode (embedding dimension -60), which by over-embedding at Ĳ=1 play role of temporally mode.After looking at several combinations of three-mode ranks for the HOSVD subspace method, we have chosen one with rank(2, 2, 2).In the following, we compared multiway tensor data denoising HOSVD method with its matrix-based SVD counterpart, where neighborhood matrix ‫ۼ‬ ‫א‬ ԧ ȉൈ contains all K transposed neighbors matrixes ‫ۼ‬ ‫א‬ ԧ ൈ concatenated along the rows.In kind of criterion evaluating the performance of the denoising of the signals the relative mean square error (MSE) İ between the normalized original signal x and estimated signal ‫ݔ‬ ො is used where ԡ‫ڄ‬ԡ is the Euclidean norm.Data processing and time and frequency analyses were performed using software written in Matlab (The MathWorks, Natick, MA).Fig. 1 shows the difference between the results obtained with the traditional local subspace approach and with the HOSVD local subspace approach.We can see that the HOSVD local subspace approach would lead to a smoother wave form without high-frequency distortions than that of the traditional local subspace approach.

Conclusions
A simulation involving the synthetic and real biomedical signals shows the efficiency, in terms of noise reduction, of the proposed local projection noise reduction approach based on three-mode model of neighborhood compared to the well know local projection approach.Furthermore, for one dimensional signal, signal subspace estimated through the HOSVD of neighborhood tensors घ ‫א‬ ԧ ൈൈ approximately yields the same result as the signal subspace estimated by the applying a SVD on large neighborhood matrixes ‫ۼ‬ ‫א‬ ԧ ȉൈ , containing all K neighbors ‫ۼ‬ ‫א‬ ԧ ȉൈ concatenated along the rows.But tensor-based approach in comparison with the matrix approach is more robust to changes in numbers of singular vectors, where formed the signal subspace.The computational complexity of algorithms based on three-mode model of neighborhood and on highly structured neighborhood matrix is significantly higher than the traditional local subspace approach.Largely the computational time is defined by size of neighbor matrixes.Therefore, the number of consecutive vectors p can't be too large, but approximately not less 10 -otherwise these approaches would degrade to the traditional local subspace approach.

Fig. 1 .
Fig. 1.Denoised signals of Rossler system at SNR=0dB: a) with traditional local subspace method, b) with the HOSVD local subspace method Fig. 2 shows the averaged relative-mean-square errors between the original and denoised signals for HOSVD based technique applied to the neighborhood tensor and for

Fig. 2 .
Fig. 2. The averaged relative-mean-square error versus white Gaussian noise level for denoising the signals of Rossler system with the HOSVD and SVD subspace methods

Fig. 3
Fig.3shows the pulmonary arterial pressure signal contaminated with 5 dB additive white Gaussian noise and denoised signal with the proposed HOSVD based approach.

Fig. 3 .
Fig. 3. Pulmonary Arterial Pressure signal, a) the raw data contaminated with 5 dB additive white Gaussian noise, b) target and denoised signals with the proposed approach