Multi-base compressive sensing procedure with application to ECG signal reconstruction

Standard compressive sensing (CS) scenario assumes a single sparsifying basis used to reconstruct the signals from a small set of incoherent measurements. However, in many cases, the signal cannot be sparsely represented using a single transformation. Particularly, in ECG signal analysis, each signal segment is specific in nature and reflects different physical phenomena. Hence, using the same transformation for all segments may be inappropriate for efficient analysis and reconstruction. Moreover, in the CS scenario, it would be necessary to combine different transforms to achieve compact signal support and to provide successful reconstruction from randomly under-sampled data. This work proposes a hybrid CS reconstruction algorithm that combines different transform basis, based on the concept of orthogonal matching pursuit. The performance of the proposed approach is verified experimentally using the combination of the Fourier and the Hermite transform on the real ECG signals.


Introduction
Compressive sensing (CS) has recently appeared as a new concept that has been used in different applications to collect, store, transmit and reconstruct large amount of data using much fewer measurements [1]. Generally, the CS paradigm has been introduced with the aim to reduce the number of measurements in signal analysis and reduce the memory requirements, as well as the energy consumption and recording time. In certain situations, CS solutions can be also used to recover parts of signals that are discarded due to disturbances and noise [2,3]. The measurements are usually represented as a small subset of original data samples acquired through a random selection process [4][5][6][7][8]. In order to provide a successful and unique signal recovery, CS reconstruction algorithms require a certain basis in which the considered signal is sparse [8]. Namely, the assumption is that a signal s is sparse when represented using linear transform Ψ such that the transform vector x has just a few significant components. However, the problem appears when the signal of interest is not sparse when represented in a single transform domain. Namely, we observe a problem when a signal is composed of different segments and each segment could be sparse in different transform domains, but the signal as a whole is not sparse in any particular domain when observed individually. Particularly, this is the case with the ECG signals composed of several different waves that correspond to electrical and mechanical heart activities. Concerning the biomedical ECG signal representation and compressive sensing, the part of the ECG signal that has been commonly treated so far is the QRS complex [9,10]. Nevertheless, providing a compact representation for the entire ECG signal is still a challenge, and it is the main motivation of this work. To be more specific, there are other parts of ECG signal such as P and T waves, which are very important in medical diagnosis (e.g. long QT syndrome, coronary ischemia), but the optimal transform domain analysis of these segments have not been treated so far. Being able to represent the entire ECG signal with very few transform coefficients would allow efficient compression, CS and classification based on the low number of features.
Exploring the CS theory can be done only if a compact (sparse) support over different ECG signal segments is provided, which cannot be accomplished using a single transformation basis. In that sense, we introduce a combined sparsity model that is suitable for ECG signals represented as a sparse set of combined basis functions (hybrid or multi-base representation). The concept of applying the hybrid approach in CS applications has been used in the literature for speech signals, where the hybrid dictionary is applied for the construction of the basis matrix in CS (linear prediction coding (LPC) and discrete Fourier transform (DFT) basis) [11]. Also, the hybrid transform/prediction sparsity model has been employed for exploiting correlation in different directions for efficient CS and reconstruction of multidimensional signals [12]. It is important to note that most of the work that has been done in this field is related to a particular application or a particular type of signal, due to the fact that the hybrid approaches emerge to deal with challenging signal structures.
An interesting approach that allows to exploit multiple structural features in the reconstruction of biomedical signals has been considered in [13]. This multi-structural signal recovery framework assumes additional a priori information which is included by imposing additional minimization constraints into the optimization model for signal recovery. This approach results in multi-criterion optimization and new convex programming problem. The considered structural properties are signal-dependent such as sparsity, piecewise smoothness and low-rank property, but generally the approach can be applied in its specific form to different types of biomedical data such as MRI, EMG, EEG and ECG signals. Particularly, in the case of ECG signals, the multi-structure optimization problem is defined as a linear combination of TV1, TV2 and L1 norms [13], exploiting piecewise smoothness and sparsity in the wavelet domain. However, as we have mentioned earlier, due to the different nature of ECG signal segments, the wavelet transform could be suitable for transient signal parts while the smooth or quasi-periodic parts could be better suited by other methods. Furthermore, it has been shown in [9] that discrete Hermite transform (DHT) provides much better compression performance and more compact support than the wavelet transform, DFT or DCT, particularly in the case of QRS complexes.
Therefore, in a similar context but without imposing additional optimization constraints, we propose an efficient combined multi-base sparsity model which can be used with the classical optimization methods to provide computationally simple solutions amenable for practical implementations. It is based on the orthogonal matching pursuit (OMP) algorithm applied over different waves of ECG signal using the combination of the two transforms. It is shown that the combined sparsity model produces the most compact representation of ECG signal parts with a highly reduced number of non-zero hybrid coefficients (relevant features). Having an optimal signal representation allows us to apply the CS concept and to recover the signal from a small per cent of randomly acquired samples. The proposed computing method is based on the modified OMP approach [14], which allows low calculation complexity, fast processing time and thus being well-suited for real-time implementations in different mobile devices and server applications [ref]. The proposed approach is also compared with several commonly used classical approaches such as iterative soft-thresholding algorithm based on the Least Absolute Shrinkage and Selection Operator (LASSO) minimization and convex optimization methods [15][16][17], showing that the proposed method represents the most convenient solution when combined with the sparsest multi-base representation. Note that we have observed the single-channel ECG data. In the case of multiple channels that need to be processed simultaneously, the large-scale problem in place requires high processing time, and in these circumstances, some computationally efficient multi-channel modifications of greedy algorithms could be considered, such as Simultaneous Greedy Analysis Pursuit algorithm [18].
The paper is organized as follows. A theory behind CS is given in Section 2. A combined signal sparsity model is introduced in Section 3. The hybrid CS using the combined sparsity model is proposed in Section 4. The application to the ECG signals and the experimental evaluation are provided in Section 5. The concluding remarks are given in Section 6.

Method
The aim of this study is to prove that the combined sparsity model can provide a compact transform domain representation of ECG signal parts which is used for an efficient compressive sensing scenario with a large amount of missing data. Namely, the ECG signals are not sparse when observed in a single transform domain and therefore it is important to provide a comprehensive approach for dealing with this type of data. Moreover, the achieved compact support representation can be important for classification, feature description and compression of ECG signals which is highly important in biomedical applications. The paper presents the general combined theoretical model for sparse signal representation, which in the case of ECG signals combines the Hermite transform and discrete Fourier transform, as well as the adapted, simple and computationally efficient compressive sensing scenario based on the combined sparsity model.
The theoretical developments are supported by the experimental evaluation on the real ECG signals from the MIT-BIH ECG Compression Test Database, while the simulations have been carried out with Matlab R2015b. twice the maximal signal frequency. In practice, it means that we are faced with a large amount of data to be sensed, stored and transmitted. Alternatively, the CS allows acquiring fewer samples while still being able to reconstruct the entire signal afterwards [1]. Thereby, the assumption is that the signal samples (called observations or measurements) are chosen randomly, and there exists a certain domain in which the signal of interest exhibits sparse representation. It is important to mention that biomedical imaging such as MRI was among the first and widely used applications of CS concept based on the advantages of using 2D DFT domain. Still, the applications on other types of biomedical signals are much less frequent in the literature, due to their specific nature [19][20][21].
Let us observe a signal s(n), with a full dataset of length N, that can be represented in basis Ψ using basis vectors ψ i : or in the matrix form s = Ψx.If most of the transform coefficients in vector x are zero (or close to zero), then we can say that x is a sparse representation of s [8]. Some of the frequently used transformation domains for CS applications are DFT, Hermite transform (HT) [22,23], discrete cosine transform (DCT), etc. The process of acquiring random signal measurements is often modelled by a certain observation matrix Φ of order M×N , where M < < N is a small number of random measurements selected out of N using Φ: while y is the vector of measurements of length M. Following (1) and (2), we have: where A = ΦΨ. Under the assumption that the sparsity, i.e., the number of non-zero elements in x is K (where K < M < < N), the linear system of equations y=Ax can be solved using different mathematical algorithms leading to the reconstructed signal. The problem is then reduced to finding the sparsest transform vector corresponding to the measurements y: where ‖ * ‖ 0 is the ℓ 0 -norm. The ℓ 0 -norm is non-convex and which, in terms of computation, makes (4) difficult to solve. The problem of finding the sparsest solution is commonly solved by using the greedy approaches such as the OMP. Firstly, a convex relaxation of the problem is applied using the ℓ 1 -norm instead of the ℓ 0 -norm. Hence, under certain conditions, one can use the corresponding convex ℓ 1 -norm optimization problem as a suitable convex approximation of (4), i.e., The problem defined in (5) can be recast as a linear program [24], and it is known as the basis pursuit problem. Another commonly used approach that was introduced in the statistics literature is based on the LASSO using the ℓ 1 -penalty to promote sparsity [25]: where t is a nonnegative real parameter, while (6) represents a quadratic program. If x is not sparse enough meaning that K>M or even K close to N, then the observed system of equations is underdetermined. Therefore, it is of high importance to identify the sparsity domain for the signal to be reconstructed. Obviously, for certain signals such as the ECG signal, this condition may be a severe limitation, since it is not appropriately sparse in any of the mentioned domains.

Combined sparsity model
Sparse signal reconstruction in CS assumes that a signal has a compact (or compressible) representation in certain transform domain. In other words, the under-sampled signal can be reconstructed from its incoherent measurements only if it can be represented by using a small number of significant coefficients in an appropriate transform basis. When exploring the sparsity property, signals are usually observed in a single transformation domain. Here, we will observe the case when a signal is not sparse in any known transform basis, but can be represented as a sparse set of combined basis functions.
Let us firstly observe the additive signal model which is common in practice: where s 1 (n), s 2 (n),…, s P (n) are of the same length N, and s 1 (n) is sparse in Ψ1, s 2 (n) is sparse in Ψ2, etc. The signal is represented using P different transform basisΨpdefined by the basis vectors: fψ 1 p; ψ 2 p; …; ψ N pg, p=1,…, P. The signal of length N contains K= K 1 +K 2 +…+K P components (K p <<N and K<N), such that: -K 1 components belong to the basis Ψ1, -K 2 components belong to the basis Ψ2, -while K P components belong to the basis ΨP.
In that sense, we might observe that a signal s(n) is not sparse in any of the transform basis Ψp when observed individually, but some parts of the signal could be observed as sparse in appropriate basis Ψp . This signal model corresponds to the concatenation concept introduced in [26], with an example of a linear combination of spikes and sines that will be sparse when concatenating coordinate and Fourier bases. However, the standard CS framework cannot be applied if there is a certain correlation between columns of concatenated basis [26].
Nevertheless, in the case of the ECG signals, we are dealing with a slightly different signal model obtained as a combination of different signal segments (combined signal model): where s 1 ðn 1 Þ ¼ i¼i 1 x i 1ψ i 1ðn 1 Þ;…, s P ðn P Þ ¼ P l K P l¼l 1 x l Pψ l Pðn P Þ; and each segment has its own duration: s 1 (n 1 ) is N 1 samples long, while s P (n P ) is N P samples long. Namely, each cardiac cycle within the ECG signal is composed of consecutive waves, where the most prominent part is the QRS complex, while the other parts are known as P, T and U waves.
In analogy with the additive model, we might assume that in the combined model defined by (8), the signal is composed of K 1 components belonging to the basis Ψ1, K 2 components belonging to the basis Ψ2, etc.

Multi-base compressive sensing using combined sparsity model
In the CS context, instead of a full-length signal s(n), we are dealing with a random set of M measurements, where M < < N. The measurement process can be defined as: where y denotes measurement vector, while s 1 , s 2,…, s P are different parts of a signal vector. Having in mind that s p is sparse when represented inΨp, for p=1,…, P, then we can write s p ¼ Ψpxp. Consequently, (9) becomes: y¼Φ Ψ1x1 Ψ2x2 … ΨPxP with Ψp being N P ×N P transform matrix and xpis an N P × 1 vector of transform coefficients and Ap ¼ ΦΨp. Alternatively, we may write it in the matrix form as follows: i.e., y = A cs X. (11) The combined transform domain representation is formed as: Note that in the case of additive signal model the concatenated CS matrix would be A CS ¼ ½A1 A2…AP , where each of the sub-matrices is of size N×N. The reconstruction problem can now be defined as follows: where the symbol ∧ is used to denote conjunction, i.e., logical and operation. This problem can be solved using the OMP algorithm as follows: The core of the algorithm is a standard OMP algorithm which is used in the context of multi-base approach and reconstruction of different signal parts x1 , x2 ,…, xP (extracted at line 11). Note that the problem defined by (13) can be solved by applying other existing reconstruction algorithms, such as convex optimization algorithms or thresholding methods. However, it will be shown that in the application with the ECG signal reconstruction, the OMP provides more accurate results comparing to other classical approaches. Furthermore, the OMP is more convenient due to its realization simplicity and the execution time when used in engineering applications [27].
In the application with ECG signals, the combination of two basis sets will be considered in the sequel, namely the set of discrete Hermite functions combined with the discrete Fourier basis. It has been proven in [9,10] that the discrete Hermite functions provide the most compact representation of the most prominent ECG segments, i.e., QRS complexes having transient signal characteristics, because of the high similarity between QRS complexes and Hermite functions. Moreover, for this purpose, the discrete Hermite transform outperforms the DFT, DCT, but also the DWT [9]. The remaining parts of ECG signal are much smoother and locally quasi-periodic in nature as can be seen in Fig. 1, and consequently, the DFT basis appears as the most suitable choice for these signal parts.
Finally, we would like to remark that the proposed multi-base approach combined with the multi-structural signal recovery [13], could be considered as an interesting future topic and extension of this work. Namely, the reconstruction problem can be further extended as a linear combination of L1 and TV-norm minimization [13], to explore the advantages of additional structures (at the cost of higher algorithm complexity). It could be also of particular importance in the sense of generalization for different types of biomedical signals. However, the scalars in the linear combination of different constraints need to be determined optimally.

Application to ECG signals
In order to demonstrate the importance of using the proposed combined sparsity domain for ECG signal analysis, we observed the real-world signals obtained from the database given in [28]. The signals are sampled with frequency 1/250 Hz. The amplitudes are scaled by 1/400 in order to have vertical axes expressed in [mV], as specified in [28]. Selected part of ECG signal is composed of different segments appearing consecutively within the signal: P, Q, R, S, T. Note that it is sufficient to observe only one combination of segments to test the proposed method, since this combination repeats over the signal (Fig. 1a). Therefore, we have (1) the segment of ECG signal consisting of T and P wave (Fig. 1b), and (2) QRS complex/segment itself (Fig. 1c). The QRS segments can be identified using some of the common algorithms [29]. In the simplified form, after detecting R peaks as the local maxima, one can use an empirically determined number of samples on the left and right side of R peaks to identify the QRS complex (e.g. 25 samples on each side for the set of considered signals). The remaining parts consist of T and P waves. Next, we apply the combined sparsity model to highlight the importance of the proposed approach. For the simplicity, in this example, we used the a priori assumption that the sparse representation of QRS complex can be achieved using the HT [9,10], while for the remaining parts the testing was performed using other standard basis functions (DFT, DCT, etc.). Consequently, it is shown that the first segment has a compact support in the DFT domain (Fig. 1b), while the second segment is sparse in the HT domain (Fig. 1c).
This means that the analysed signal can be represented by a small set of combined DFT and HT coefficients, while the remaining coefficients are zeros. In that sense, we have concentrated the signal description into few coefficients that can be used as relevant signal features in other applications. For the comparison of achieved results, we can observe the selected part of the signal in separated domains (in DFT or HT, respectively). Namely, the individual representations (Fig. 2a, b) contain large number of non-zero coefficients and cannot be considered sparse as it is the case with the proposed hybrid (combined) representation (Fig. 2c). Consequently, in order to provide compact support representation (for compression or CS purpose), the hybrid combined DFT-HT should be employed, as shown in Fig. 2c. One can thus conclude that the hybrid combined DFT-HT transform provides the most suitable representation for ECG segments ensuring the lowest possible number of non-zero coefficients. Moreover, in order to explore the concept of CS, instead of full dataset we observe the signal with missing samples. According to (11), the measurement vector is given by: where A1 denotes a random partial DFT matrix and A2denotes a random partial HT matrix. Since x1 and x2 are approximately of the same length N/2 (N is the total segment length), it follows that A1 and A2 are of the same size M/2×N/2. The measurement vector is supplemented by zeros at the positions of missing values and it is shown in Fig. 3a.
The initial hybrid transform is shown in Fig. 3b, where the sparsity is ruined as a consequence of dealing with the missing samples (60% of samples are missing).
As a reference signal (desired signal) we use the original signal with full set of samples and the corresponding reference hybrid transform (Fig. 4, upper row). The results of applying the proposed signal reconstruction algorithm are shown in Fig. 4, bottom row (reconstructed time domain signal and the corresponding hybrid transform). It can be observed from the figure that the reconstructed data highly match the original data, confirming the high reconstruction accuracy.

Results and discussion
The influence of the number of available samples to the quality of reconstruction is expressed in terms of the mean squared error (MSE) calculated for different realizations, i.e., different random combinations of available samples. The MSE between the original (desired) and the reconstructed signals are shown in Fig. 5 for different percentage of available samples. In general, when observing the averaged results from different realizations, we may conclude that the number of available samples should be at least 40% in order to assure the acceptable quality of reconstructed signal. However, in CS scenarios it is usually possible to apply the most efficient random combination of available samples that would allow the high-quality reconstruction even for the lower percentage of available measurements.

Comparison results
In order to assess the performance of the proposed approach, in this part, we examine the reconstruction using a single transform, namely, the HT (Fig. 6a) or the DFT (Fig.  6b). The signal with missing samples and the corresponding transform domain representations are shown in the upper row in Fig. 6. Also, the original signal with full set of samples and the corresponding transforms are shown in the middle row in Fig. 6, while the reconstructed signal and the corresponding transforms (HT and the DFT) are shown in the bottom row in Fig. 6. The results indicate that the reconstructed signal differs significantly from the original one in both observed cases. Next, we provide a comparison between the proposed OMP based method (with hybrid transform) and other commonly used algorithms, such as (1) the primal-dual interior point method implemented within the package l 1 -magic [16], and (2) the iterative soft-thresholding algorithm (ISTA) [30] (which is based on the LASSO minimization). The reconstruction results using the two considered algorithms are  shown in Fig. 7. The original (desired) signal and its transform are given in Fig. 7a. The signal reconstructed using l 1 -magic algorithm is shown in Fig. 7b, while the reconstructed signal obtained using LASSO-ISTA algorithm is shown in Fig. 7c. Based on the reconstructed signals, one can observe that the results obtained using these two algorithms are significantly worse that the result obtained by employing the proposed method (Fig. 4). Moreover, it is important to highlight that, since the ECG signals are used in diagnosis of different heart disorders, the two algorithms (l 1 -magic and LASSO-ISTA) cannot be used because they produce errors higher than the acceptable threshold in biomedical applications. To support this conclusion, the numerical comparison is also provided in terms of the MSE and it is given in Table 1. Note that the MSE for the proposed algorithm is very low, meaning that the reconstructed signal is approximately the same to the original (desired) signal, while the same cannot be said for l 1 -magic nor LASSO-ISTA. The CPU times for each of considered algorithms is given in Table 1 as well (the evaluation is performed per ECG signal segments of 100 samples, corresponding approximately to one cardiac cycle).

Conclusion
In this work, the CS algorithm for data reconstruction was introduced based on the combined (hybrid or multi-base) sparse representation. The main objective here was to enhance the sparsity of ECG signals representation and to open up new perspectives for CS not only as an emerging technology, but also for standard signal compression and classification using a lower number of relevant features. The general concept of the combined sparsity model was introduced, while the particular methodology was developed to identify signal components belonging to different sources. In this way, we were able to provide an optimal representation and CS reconstruction of ECG signals. The main contributions of the proposed solution are twofold. Firstly, it provides the sparsest (compact) representation of ECG signals that cannot be sparsely represented in any of the commonly used transformation domains when observed individually. Namely, the state of the art methods, observe only QRS complexes when considering signal compression, while the other parts of the signal are not treated. This is because these parts cannot be appropriately represented in the same transform basis as the QRS complexes.
In that sense, we used the combined multi-base transform, providing improved representation compared to the commonly used domains. In particular, the combined transform is concentrated in just a few coefficients that can represent the signal accurately, while the rest of the coefficients are zeros. Secondly, the proposed solution allows an efficient CS scenario for ECG signals, where it was shown that the exact signal