Compressive sensing in electrical impedance tomography for breathing monitoring

Objective: Electrical impedance tomography (EIT) is a functional imaging technique in which cross-sectional images of structures are reconstructed based on boundary trans-impedance measurements. Continuous functional thorax monitoring using EIT has been extensively researched. Increasing the number of electrodes, number of planes and frame rate may improve clinical decision making. Thus, a limiting factor in high temporal resolution, 3D and fast EIT is the handling of the volume of raw impedance data produced for transmission and its subsequent storage. Owing to the periodicity (i.e. sparsity in frequency domain) of breathing and other physiological variations that may be reflected in EIT boundary measurements, data dimensionality may be reduced efficiently at the time of sampling using compressed sensing techniques. This way, a fewer number of samples may be taken. Approach: Measurements using a 32-electrode, 48-frames-per-second EIT system from 30 neonates were post-processed to simulate random demodulation acquisition method on 2000 frames (each consisting of 544 measurements) for compression ratios (CRs) ranging from 2 to 100. Sparse reconstruction was performed by solving the basis pursuit problem using SPGL1 package. The global impedance data (i.e. sum of all 544 measurements in each frame) was used in the subsequent studies. The signal to noise ratio (SNR) for the entire frequency band (0 Hz–24 Hz) and three local frequency bands were analysed. A breath detection algorithm was applied to traces and the subsequent error-rates were calculated while considering the outcome of the algorithm applied to a down-sampled and linearly interpolated version of the traces as the baseline. Main results: SNR degradation was generally proportional with CR. The mean degradation for 0 Hz–8 Hz (of interest for the target physiological variations) was below ~15 dB for all CRs. The error-rates in the outcome of the breath detection algorithm in the case of decompressed traces were lower than those associated with the corresponding down-sampled traces for CR  ⩾  25, corresponding to sub-Nyquist rate for breathing frequency. For instance, the mean error-rate associated with CR  =  50 was ~60% lower than that of the corresponding down-sampled traces. Significance: To the best of our knowledge, no other study has evaluated the applicability of compressive sensing techniques on raw boundary impedance data in EIT. While further research should be directed at optimising the acquisition and decompression techniques for this application, this contribution serves as the baseline for future efforts.


Introduction
The electrical conductivity (σ) and the permittivity (ε) of materials determine their admittivity (γ = σ + j ωε), where j = √ −1 and ω is the angular frequency (Cao and Xu 2013). Assuming γ is constant, for a given electric field (E), the current density (J) in a material may be found as J = γE. For a given geometry, γ is inversely proportional to the impedance (Z) of the sample which is identified as the ratio of the voltage across and current passing through the sample. An indication of this Z may be measured by applying electrical current through a pair of electrodes and recording the corresponding induced electrical voltage across another two in the vicinity to measure a transfer impedance (Grimnes and Martinsen 2006). The relative positions of current carrying and recording electrodes determine the segment of the sample whose impedance has the highest impact on the measurement. In conventional electrical impedance tomography (EIT), alternating current is sequentially applied via different electrode pairs and voltage is recorded via other pairs in each case. The real and imaginary components of boundary trans-impedances can then be calculated. These boundary impedances can then be used to estimate the internal impedance changes of the structure and reconstruct an image. Solving for impedance from boundary trans-impedance measurements is an ill-posed, non-linear inverse problem.
Proposed biomedical applications of EIT include gastric (Podczeck et al 2007), brain (Li et al 2017) and thoracic (Frerichs et al 2016) monitoring as well as breast (Polydorides 2018) and prostate (Murphy et al 2018) imaging. Amongst all applications, due to a significant contrast between the conductivity of lung tissue and that of other structures in human thorax and considerable variations of this conductivity during a breathing cycle, one of the most researched applications of EIT has been that of breathing monitoring (Hsu et al 2017). The breathing cycle is composed of two phases of inspiration, and expiration. These phases and their key physiological attributes may be monitored using EIT.
There is great interest in high temporal resolution for 2D and 3D EIT which will increase the number of measurements for each frame. This will produce a large volume of data which will be a limiting factor for data transfer and storage. In a typical 32 channel EIT system implementing adjacent electrode current injection and corresponding adjacent voltage measurements, as many as 1024 measured impedances for each frame are produced, leading to in excess of 1.5 GB data at a sampling rate of ~50 frames per second depending on the data format. This is particularly a limiting factor when raw measurements should be transmitted wirelessly to a base unit. This presents a strong motivation for exploring compression strategies for raw impedance data in EIT. Compressive sensing (CS) is a technique through which it may be possible to sense the desired signal directly into a compressed form and recover the original signal if certain conditions are met. This efficient compression strategy has been used in electroencephalogram (Moy et al 2017) and magnetic resonance imaging (Lustig et al 2008) acquisition amongst other applications (Rani et al 2018).
In this paper, we investigate the feasibility of employing conventional CS techniques in compressively sensing raw boundary impedance in EIT for breathing monitoring. Since breathing is a relatively low-frequency event (⩽2 Hz) and existing EIT systems operate at frequencies much higher than the Nyquist rate for breathing, it is possible to assess CS, in terms of information loss, over a wide range of compression ratios (CRs). CR is defined as the ratio of the number of samples in the original signal over the number of samples in the compressed signal. Deductions may be extended to predict performance for higher rate physiological events such as heart rate. Random demodulation acquisition technique was applied on appropriate channels in a number of frames for 30 subjects at different CRs. We then applied a reconstruction algorithm to the compressed traces and compared the decompressed traces with the non-compressed signal in terms of their signal to noise ratio (SNR). We further applied a breath detection algorithm to these traces and compared the number of breaths detected while monitoring the correlation coefficient between decompressed traces and the original signal. We also tested the outcome on a downsampled version of traces. While the outcome may not be quantitatively pertinent at this stage, this study qualitatively indicates that there is a strong motivation for further research in the use of CS in EIT and future research should focus on optimising the CS processes.
In section 2, the theory of CS in general and as related to EIT is summarised. In section 3 methods of data acquisition and post processing including acquisition methods, decompression algorithm and breath detection algorithm are described. In section 4, the results in terms of changes in the SNR as a result of compressive sensing and relative errors in the outcome of breath detection are presented. The ensuing discussion and concluding remarks are presented in sections 5 and 6, respectively.

Compressive sensing
Assuming x N×1 is an arbitrary discrete signal in R N , it may be represented using a basis matrix . Thus, α N×1 represents the signal in Ψ N×N (e.g. frequency for Fourier basis). The signal is referred to as K-sparse in a target domain when only K coefficients are needed to represent the signal. If K N ( N may reflect the number of samples collected at Nyquist frequency), the signal becomes compressible (Baraniuk 2007, Candès andWakin 2008). This is utilised in conventional post-sampling compression by transforming the signal to the domain in which it is compressible, identifying the relatively large (depending on application) coefficients and only storing the indices and values of those.
It is more efficient to sense the signal directly in a compressed form, from which the original signal may be recovered. A linear sensing system may be thought of as a process that computes M < N inner products between x N×1 and a vector set Baraniuk 2007). It is noted that in a conventional sampling scenario {φ j N×1 } M j=1 is a set of vectors that have only one element as 1 and the rest as 0 and the position of 1 s are shifted. Arranging the M measurement vectors, In CS, the primary tasks are to find a nonadaptive Φ M×N to sample the compressible signal while reducing its dimensionality without loss of necessary information and to design a decompression algorithm to reconstruct the original signal. The sufficient condition in designing the measurement matrix (Φ M×N ) operating on a signal sparse in Ψ N×N is that its rows (Candès and Wakin 2008) or they should not sparsely represent each other (Baraniuk 2007). This condition is referred to as incoherence which is a special case of restricted isometry property (Candès and Wakin 2008). By selecting Φ M×N as a random matrix, the incoherence condition may be met with high probability.
Various acquisition strategies in hardware have been proposed to have this randomisation in the dataset and appropriate reduction of the dimensionality. Namely, random demodulation (Kirolos et al 2006), modulated wideband conversion (Sun et al 2013), random modulation pre-integration (Yoo et al 2012), random filtering (Tropp et al 2006), random convolution (Romberg 2009), compressive multiplexer (Slavinsky et al 2011), random equivalent sampling (Zhao et al 2012), random triggering based modulated wideband compressive sampling (Zhao et al 2017) and quadrature analogue to information conversion (Haque et al 2015) are some of the acquisition strategies proposed in the literature. A comprehensive review of these techniques may be found in Rani et al (2018).
The second problem is to reconstruct the original signal. Since M < N , y M×1 = Θ M×N α N×1 will have many solutions. p -norm of a vector is defined as p » N i=1 |υ i | p and 0 -norm is an operator that counts the number of non-zero elements. Inverse problems are classically solved by minimising 2 -norm (distance). However, this will not give the sparse solution. The problem may be solved by minimising 0 -norm which counts the number of non-zero entries in α N×1 as the first intuitive method. However, this method is numerically unstable (Baraniuk 2007) and NP-hard (Rani et al 2018). Considering 1 -norm reconstruction, the convex optimisation problem may be solved using linear basis pursuit shown in (1): where α is the estimate of α. Inputs to algorithms used to solve this problem are y M×1 and Θ M×N , and the output is an estimation of sparse representation of x N×1 (i.e. α) which can be converted to x by applying inverse transform from the said domain to the original domain. The algorithms used for reconstruction apart from convex optimisation (Donoho 2006) are based on greedy (Eldar and Kutyniok 2012), thresholding (Blumensath and Davies 2009), combinatorial (Cormode and Muthukrishnan 2006), non-convex (Chartrand and Staneva 2008) and Bayesian (Ji et al 2008) approaches. A comprehensive review of the acquisition methods and reconstruction algorithms may be found in Rani et al (2018).

Compressive sensing in EIT
Previous implementations of CS techniques in EIT in the literature , Tehrani and McEwan 2010, Gehre et al 2012 were confined to the use of similar concepts (e.g. 1 -norm minimisation method) in solving the ill-posed EIT problem. However, no other study has applied CS techniques in acquiring boundary impedance data. It was noted that a single frame of EIT is reconstructed based on many boundary impedance measurements. The partial periodicity of physiological variations is reflected in the frequency spectra of these measurements at a bin (or a range of bins) associated with the frequency of these variations. Since a measurement reflects the induced voltage across a pair of electrodes due to current injected across another pair, the extent to which each of these variations are portrayed in the measurements is varied across different measurements. Figure 1 shows the frequency spectra of selected channels in 2000 frames (~42 s). There are sharp spikes in the frequency spectra of all measurements although it should be noted that obviously traces are not completely sparse in frequency domain and this may affect the results. Each channel, which is a measurement due to a specific current injection, may have a different spectrum with a different set of dominant spikes. Thus, these may indicate that conventional CS techniques may be employed to compressively sense raw impedance data and recover original traces after transmission prior to image reconstruction.

Data collection
Data collection was performed within an ongoing prospective clinical study (CRADL) 11 . The study is registered in a clinical trials registry (ClinicalTrials.gov, NCT02962505). It was approved by the ethics committees at three clinical study sites: the Emma Children's Hospital, Amsterdam, Netherlands (Ethics number: METC 2016/184), the Archbishop Makarios III Hospital, Nicosia, Cyprus (Ethics number: EEBK/EP/2016/32) and the Oulu University Hospital, Oulu, Finland (Ethics number: EETTMK 35/2017). The data associated with a group of 30 subjects, simply based on the order of recruitment, were used in this study.
Infants with a body weight less than 600 g, those with postmenstrual age less than 25 weeks at inclusion, infants with electrically active implants, or those suffering from thoracic skin lesions were excluded from the study. Informed written consent was obtained from the parents of the neonatal study participants. The raw EIT data were acquired by the CRADL study EIT device (SenTec AG, Therwil, Switzerland, www.sentec.com/) specially designed for small thoracic diameters with 32 textile electrodes at a frame rate of 48 frames per second (Waldmann et al 2017, Sophocleous et al 2018. Current injections with amplitude of 3 mA (root mean square) at a frequency of 200 kHz were applied using a skip-4 pattern. The resulting voltage differences were measured by the remaining electrode pairs after each current injection and the impedance data were derived. Data associated with current carrying electrodes and two neighbouring electrodes were discarded due to the saturation of the recording system leading to 544 measurements for each frame. 2000 frames (i.e. ~42 s) from each of the 30 subjects were used. Only the real part of measurements was used.

Dimensionality reduction acquisition matrix
All the subsequent data processing was performed by post-processing in Matlab R2017a (Mathworks, Natick, USA) using a computer with Intel Core i7-6700 CPU @ 3.4 GHz and 32 GB RAM.
Each of the 544 measurements was baseline corrected by subtracting the output of a 250 sample wide moving average function from measurements to suppress low frequency drift and DC levels. Random demodulation (Tropp et al 2010) acquisition method which may be readily implemented in hardware was implemented for dimensionality reduction. The 2000 frames for each of the 544 channels were compressed at seven different CRs defined as the ratio of the size of the original samples and the compressed vector (i.e. CRs of 2, 4, 8, 16, 25, 50, 100 corresponding to 1000, 500, 250, 125, 80, 40 and 20 frame samples in the associated compressed vectors, respectively).
Random demodulation can be applied by multiplying the signal by a random sequence of ±1 followed by low pass filtering which has been shown to appropriately acquire signal in a reduced dimension for a signal sparse in frequency domain (Tropp et al 2010). In terms of matrix implementation, random demodulation may be thought of as a diagonal matrix ( D N×N ) with the random +1 or −1 as its diagonal elements and an accumulateand-dump matrix ( H M×N ). Assuming that the dimensionality is to be reduced from N to M and N/M has no remainder, H M×N has N/M 1s in each row. The first row of H M×N sums the first N/M samples, the second row will sum the second N/M samples and so forth. The overall measurement matrix can be formulated as Φ M×N = H M×N D N×N . The first operation spreads the frequency spectrum by randomising the trace and the second captures the low frequency portion of the spectrum (Rani et al 2018). An example matrix is shown in (2) for 6 to 2 reduction: A unique random acquisition matrix was generated for each subject to diversify the outcome across all subjects but the same matrix was used for the 544 traces associated with each subject to be consistent in each case.

Decompression algorithm
SPGL1 (van den Berg and Friedlander 2007) package was used to solve the basis pursuit problem in Matlab. Our preliminary tests did not indicate substantial change in results when implementing basis pursuit denoise (BPDN) or least absolute shrinkage and selection operator (LASSO) in SPGL1 to decompress and this package was relatively faster compared to 1 -Magic (Candes and Romberg 2005). The arguments passed on to the algorithm were the measured dimensionality reduced vector y M×1 and Θ M×N . Ψ N×N in each case was formed by calculating the complex conjugate of discrete Fourier transform matrix (∆ N×N ) divided by N. ∆ N×N was found by performing a fast Fourier transform on an identity matrix. ∆ −1 N×N was applied to the output of the algorithm to find the temporal trace.

SNR comparison
For the subsequent processing of data, the total impedance signal was formed by summing all the measurements associated with each frame (i.e. n = 544). The power spectral density of the total impedance signal was calculated. The local peaks of the power spectral density where they were at least 10 −4 times the maximum peak in terms of the height and prominence with at least 6-bin spacing between the peaks were identified. To avoid excessive number of peaks while including all the prominent peaks, the power spectral density was passed through a moving average filter 4 bin wide twice. The bins associated with these peaks and 6 bins before and after them were considered to be those containing signal while the rest of the spectrum was noise. The SNR was estimated by forming the ratio of the integral of the power spectral density over the bins associated with signal and its integral over those associated with noise. This was done for the entire frequency band and individual 8 Hz bands.

Breath detection
A breath detection algorithm (Khodadad et al 2018) based on zero crossing algorithm with amplitude threshold was used. The fact that the algorithm is dependent on the excursions of traces is important in assessing the way CS affects signal fidelity. The total impedance signal in each case was passed through a moving average filter 25 samples wide in non-compressed and decompressed cases. The number of breaths detected in each case was noted while the error with respect to the number of breaths detected in non-compressed signal was calculated as |NB C − NB NC | /NB NC , where NB C is the number of breaths detected for decompressed trace and NB NC is the number of breaths detected for the associated original trace. We further downsampled the uncompressed traces to a level corresponding to the CRs considered, then linearly interpolated the resulting traces to decompress them. This linear piecewise representation would facilitate the subsequent processing while providing a baseline for the minimum performance. These traces were also passed to the breath detection algorithm and the corresponding errors were noted. The median error and its 25th and 75th percentiles were calculated for each CR. Since the excursions of the temporal traces are of interest, we also monitored the correlation coefficient between decompressed traces and the corresponding non-compressed traces. Figure 2 shows a 15 s long example trace of total impedance for decompressed signals corresponding to different CRs and the associated frequency spectra. The traces have been normalised in each case for a clearer presentation. From the example traces it is clear that the signal fidelity is clearly inversely proportional with the CR. Heart rate which appears as a peak in the frequency spectrum just below 3 Hz is only distinctively recovered for CR = 2 and 4 and for higher CRs it may be buried in spurious components introduced after reconstruction. Higher CRs appear to lead to more sparsely reconstructed traces although more spurious components are introduced.

Results
A clearer representation of the outcome across different frequency bands can be observed in figure 3 where ∆SNR defined as the difference between the SNRs of the uncompressed and decompressed traces is introduced. Understandably, at higher frequency bands SNR is degraded more severely. Another observation is that at higher frequency bands and at higher CRs the SNR degradation has a higher standard deviation across subjects. The SNR degradation in the middle band is less dependent on the CR. It is noted that these SNR values are only presented for comparison between different CRs and different bands and do not necessarily reflect the effects on clinical decision making. However, SNR may be considered to be proportional with the accuracy of clinical decision making as it may be considered as a measure of fidelity. The extent to which this is important should be assessed in a case by case basis depending on the features of a given target marker as it is shown for breath detection in this paper and shown in figure 4. Figure 4 shows the median and respective percentiles of the error in the number of breaths detected for different CRs for compressively sensed traces and downsampled-interpolated ones across all subjects. The latter may  be used as the baseline. While the median error for decompressed traces is comparable and slightly higher than that of downsampled-interpolated traces, for CR ⩾ 25 the median error associated with decompressed traces becomes clearly lower than the median error of the number of breaths detected for downsampled-interpolated traces due to the deterioration of the outcome in the latter. It is noted that CR ⩾ 16 correspond to sub-Nyquist sampling.
As an estimate of phase error, the mean and standard deviation of the correlation coefficients between downsampled-interpolated and decompressed traces, each with the corresponding original signals, across all subjects for various CRs are shown in figure 5. This shows that only for CR ⩾ 25 is the mean coefficient associated with decompressed traces in a similar range to those associated with downsampled-interpolated traces. The exact value of the correlation coefficient may not be quantitatively relevant in terms of the way clinical decision making is impacted; however, this can give an estimate of how temporal excursions of the signal are preserved. Nonetheless, certain excursions may or may not affect a specific clinical deduction. For instance, in the case of the breath detection algorithm used, zero crossing is of value. Thus, not all forms of dissimilarity may lead to an error in that case.

Discussion
CS has been implemented in a variety of biomedical applications where either Nyquist rate is too high or data volume is excessive. In the case of EIT while the Nyquist rate with respect to target physiological events may not be comparatively high but the volume of data produced prior to image reconstruction is excessive. This is a problem for transmission and storage if the solution is to be deployed clinically. The problem will be more severe if multi-plane, multi-frequency and faster EIT systems are needed for an application. In this study, it was shown that based on the acquisition and reconstruction methods used and comparing the outcome with that of downsampled-interpolated traces, the ensuing errors are moderate and the performance is better for breath detection at sub-Nyquist sampling rates (relative to downsampling). The reason for the non-ideal outcome is that impedance signals may not be completely sparse which may be caused by aperiodic organ movement. It is noted that any periodic movement of organs still leads to sparsity in impedance signals.
A key point here is to note the changes required in hardware to ensure compressive sensing is possible. EIT instrumentation often implements synchronous detection in which recorded signal is multiplied by in-phase and quadrature versions of the applied signal to derive real and imaginary parts of impedance. Random demodulator may be incorporated in the multiplier that already exists in the hardware. Other existing acquisition methods may require more elaborate changes to the existing hardware.
We implemented a simple sensing matrix readily implementable in hardware (i.e. random demodulation) and used an algorithm that solved the basic basis pursuit problem for reconstruction. While this study sets the baseline for future efforts, a more elaborate sensing matrix, taking into account the multichannel nature of the signal, and more suitable reconstruction algorithms may lead to significantly lower error-rates. Also, it should be noted that the decompression process may slow down the real time performance. Thus, efficient computational techniques should be implemented.

Conclusions
As discussed, the large volume of data generated for every frame of EIT may impede its clinical deployment. This paper demonstrated a proof of concept study in using CS techniques to reduce data volume in EIT prior to image reconstruction by post-processing impedance data taken from 30 neonates and assessing the outcome of a breath detection algorithm on decompressed traces. While perfect reconstruction is not possible at this stage since the breathing signal in the raw impedance data is not completely sparse in frequency, further improvements to acquisition and decompression techniques as well as the signal processing paradigm may lead to a significantly better outcome. In this proof of concept study, breath detection and the way it may be affected by CS incorporation were assessed. Since the algorithm used in this study for breathing monitoring uses the total impedance data, this feature served well as a single indicator of possible future prospects. However, future studies should identify how other clinical information such as regional ventilation can be affected.
The outcome of this study may be extended to deduce the possible outcome in detecting events of higher rates such as heart rate. Future studies may aim to use more optimised acquisition and decompression techniques while systematically analysing the effects of extra complexity added to the system due to the implementation of CS techniques. In terms of the clinical decision making, also, the CS embedded EIT system should be optimised while taking into account the more nuanced features of data as mentioned above.