Evaluation of Thoracic Equivalent Multiport Circuits Using an Electrical Impedance Tomography Hardware Simulation Interface

Abstract

Electrical impedance tomography is a low-cost, safe, and high temporal resolution medical imaging modality which finds extensive application in real-time thoracic impedance imaging. Thoracic impedance changes can reveal important information about the physiological condition of patients’ lungs. In this way, electrical impedance tomography can be a valuable tool for monitoring patients. However, this technique is very sensitive to measurement noise or possible minor signal errors, coming from either the hardware, the electrodes, or even particular biological signals. Thus, the design of a good performance electrical impedance tomography hardware setup which properly interacts with the tissue examined is both an essential and a challenging concept. In this paper, we adopt an extensive simulation approach, which combines the system’s analogue and digital hardware, along with equivalent circuits of 3D finite element models that represent thoracic cavities. Each thoracic finite element model is created in MATLAB based on existing CT images, while the tissues’ conductivity and permittivity values for a selected frequency are acquired from a database using Python. The model is transferred to a multiport RLC network, embedded in the system’s hardware which is simulated at LT SPICE. The voltage output data are transferred to MATLAB where the electrical impedance tomography signal sampling and digital processing is also simulated. Finally, image reconstructions are performed in MATLAB, using the EIDORS library tool and considering the signal noise levels and different electrode and signal sampling configurations (ADC bits, sampling frequency, number of taps).

1. Introduction

Electrical impedance tomography (EIT) is a low-cost imaging technique which applies a 1 kHz–1 MHz frequency low amplitude electrical current to a subject under test (SUT) and collects the measured voltages from an electrode cluster [1,2]. Then, an estimation of the SUT internal conductivity distribution is performed through a set of inverse algorithms [3]. Although EIT is characterized by critically low spatial resolution compared to X-ray, CT, and MRI, it has the advantages of the absence of ionizing radiation, patient safety, low hardware cost (20 to 100 times compared to a CT scanner), and high speed. High frame rates lead to significant temporal resolution, which makes EIT suitable for applications that need real-time monitoring, such as thoracic and ventilation imaging [3,4,5].

EIT’s spatial resolution, however, might be further limited by signal noise, since the EIT inverse problem is ill-posed and ill-conditioned [6]. Image artefacts that appear are also related to common signal effects, usually caused from mismatches in the current source and between the electrode channels, combined with limited common mode rejection ratio (CMRR) at the voltage recording circuitry [5,7]. In addition, signal degradation due to the channels’ or IC’s stray capacitances is a factor that reduces the image quality, especially at frequencies above 100 kHz.

In the case of real-time (dynamic) thoracic imaging, these effects are more intense, since higher signal frequencies are required to achieve a sufficient frame rate and thus temporal resolution in order to extract useful clinical information related to the patient’s lung functionality. We note that for perfusion (cardiac-related impedance changes) imaging, even higher frequencies and voltage $SNR$ levels are essential [8,9]. Furthermore, the effect of thoracic boundary change during each breath cycle, as well as the fact that the accurate boundary shape is unknown, leads to the presence of artefacts. Two related studies demonstrated that a mismatch of more than $4\%$ between the actual thoracic boundary shape and the EIT reconstruction domain boundary has major effect in the image quality [10,11].

Considering the many possible sources of imaging errors, the design and implementation of a good performance EIT hardware which is capable of high signal quality excitation is the first essential step for successful EIT imaging. Many different EIT systems have been implemented over time, both generic and application-targeted. Such general purpose EIT systems are the Sheffield MK (v1 1987, v3.5 20021) [12], the ACT 3 and ACT 4 EIT systems [13,14], and the KHU Mark 2.5 system (2014) [15]. In addition, some EIT systems have been developed specifically for dynamic thoracic imaging, such as the Swisstom AG (2012), the ACE1 [16], and the systems presented in [4,17]. These EIT systems make use of active electrodes [18,19] in order to reduce the effect of stray capacitances. Although there has been such improvement in the performance of EIT systems during the last decade, their design and implementation still remains a challenging concept and an open research topic.

Some research has been performed on how to set the proper requirements and design an EIT system. Ref. [20] performs an $SNR$ analysis for the voltage acquisition part of custom EIT systems, deriving a specific model to calculate it. Ref. [21] presents a structured design methodology to achieve a high $SNR$ in EIT systems, focusing on the choice of the instrumentation amplifier at the voltage acquisition part and the analogue-to-digital converter’s (ADC) specifications. Similar analysis has been performed in [22], where a low-power readout front-end is also designed, and in [4,17,23], with integrated circuit designs.

Although research has been done in system design optimization, it mainly focuses on noise compensation. However, the design of an EIT system faces many other challenges, such as common signal effects due to contact impedances and parasitic effects [5]. Hence, a simulation approach which includes these effects could actually assist in examining their possible impact in imaging quality.

Simulation approaches to EIT hardware and impedance modelling have already been presented. In [24], a general 2D SPICE multiresolution impedance method for low frequencies is presented, where a resistor network mesh is created. In [25], a MATLAB interface which transfers a finite element (F.E.) mesh with assigned impedance or conductivity values to each element to a resistive multiport equivalent was presented. This interface is part of the EIDORS library tool, which is commonly used for EIT image reconstructions [26]. A preliminary integration of 2D F.E. mesh circuit equivalents to basic EIT circuitry for SPICE simulation and image reconstruction was presented in [27], where the impact of electrode short-circuit or disconnection was simulated. This approach was extended in [28] where simulations on 3D saline tank equivalents were performed under single ended and mirrored Howland current pumps. The DAC and ADC functions were also simulated in MATLAB, and results were acquired and compared for 2 current signal frequencies and 2 ADC sampling rates. A minor extension was presented in [29], where the thermal, flicker, and quantization signal noise effects, as well as the demodulation part, were added to the previous interface.

In this paper, the SPICE-MATLAB EIT simulation interface described in [28,29] is applied to three-dimensional thoracic structures. In specific, two fine F.E. thoracic structures are extracted, based on CT images and representing the deflated (full-exhalation) and inflated (full-inhalation) states respectively. Each tissue’s relative conductivity and permittivity are loaded from a web database using Python, according to the current signal selected frequency. Both structures are transferred to 16-port RLC networks and added in the SPICE EIT circuitry for transient simulations, resulting in two corresponding measurement frames. Signal noise is included, and the sampling and digital processing are simulated in MATLAB, in order to acquire the expected image reconstructions. Simulations are carried out assuming 2 electrode sizes and various ADC specifications and number of taps per voltage channel measurement. It is noted that instead of the previously presented interface which uses only resistive equivalent networks, in this work, both SUT conductivity and permittivity are considered, leading to a frequency-dependent RLC multiport equivalent network. Thus, a standard in-phase and quadrature (IQ) demodulation technique is also simulated in MATLAB. Finally, qualitative and quantitative evaluation of the reconstructed images is performed in order to compare the corresponding image results of each simulated configuration.

The remainder of this paper is written as follows. In Section 2, the EIT measurement principle is briefly described. In Section 3, the thoracic structures evaluated in this study as well as the tissues admittance data are presented. Section 4 synopsizes the SPICE-MATLAB interface used for the simulations. Furthermore, in Section 5, the reconstruction domain used, the algorithm with its parameters, and the image evaluation method are described. In Section 6, the hardware and electrode configurations are presented, while presentation and comparison of the results are performed. Finally, Section 7 concludes this paper.

2. EIT Principle

In EIT, a low or medium frequency AC current (usually up to 1 MHz, depending on the application) of a small amplitude (up to $9m{A}_{p-p}$) is injected to the SUT through 2 selected electrodes from an electrode array attached on its surface. At the same time, differential voltages are acquired through other electrode pairs, rendering a set of tetrapolar measurements. The process is repeated for a number of current electrode pairs. The final set of tetrapolar measurements acquired consists of a measurement frame $\mathit{V}$. In time-difference EIT (td-EIT) imaging, which is usually used for real-time thoracic imaging, two measurement frames are needed to reconstruct an image. Their difference $\delta \mathit{V}$ is used as data in the reconstruction algorithm.

Although many measurement patterns have been used in EIT, the most commonly utilized in lung imaging is the adjacent pattern [30,31]. It is characterized by the current excitation of two neighbouring electrodes and voltage measurement between two other neighbouring electrodes at each time [32]. Although the adjacent pattern presents low field sensitivity to the current injected near the centre of the SUT, it is often preferred due to the fact that it offers more independent measurements compared to other patterns [32,33]. In terms of this work, we make use of the adjacent current pattern. However, the simulation interface can be easily modified in order to include other measurement patterns (as in [29], where the skip-3 current pattern is also activated [34]).

The basic EIT circuitry consists of a voltage controlled current source (VCCS) of sufficient transconductance and large output impedance in the largest possible frequency range [5,35]. The signal waveform (usually sinusoidal, but sometimes pulse signal as superposition of multifrequency signals is used) can be digitally produced from a digital-to-analogue converter (DAC), controlled by a direct digital synthesizer (DDS). The voltages are acquired by an analogue front-end (AFE) which usually includes one or more instrumentation amplifiers (IA), filters, programmable gain amplifiers (PGA), and an ADC driver [1,5,21]. Both the current source and the voltage recording parts are connected through switching multiplexers to an electrode array which is attached on the SUT surface. The analogue voltages are sampled from an ADC (usually successive-approximation, SAR), and the samples are sent to a processing unit (MCU, DSP, or FPGA), which performs the synchronization with the input signal, the demodulation, and the switch control. Finally, the processed measurements are sent to a PC for the image reconstruction. A brief schematic of the description above is presented in Figure 1. It is noted that this is a basic EIT configuration principle and not necessarily the unique or optimal one. For example, some systems use active electrodes [4,17], placing the multiplexers before the VCCS (current source) and after the first IA. Other systems drive multiple electrode pairs with more than one current frequency [36] (these systems apply frequency-difference EIT usually used in brain imaging applications). Moreover, some systems perform simultaneous measurements at the voltage recording front-end [16], while many systems enact analogue demodulation before the ADC part [37,38].

3. Thoracic Structures

The thoracic SUT was simulated by implementing two fine (dense) 3D F.E. models, representing the full-exhalation (deflated) and the full-inhalation (inflated) lung states respectively (Figure 2). Both models were implemented in MATLAB using the NETGEN tool [39]. Each model includes the following tissues: left lung, right lung, heart, vertebra, skin, and the (muscle-plasma) background. Each tissue’s conductivity $\sigma$ and permittivity $ϵ$ were loaded from a database demonstrated in [40,41,42], according to the current signal frequency selection (i.e., the frequency of the current signal produced from the VCCS). The corresponding conductivities and permittivities per frequency between 1 kHz and 1 MHz are demonstrated in Figure 3. In this particular work, we have examined measurement (current signal) frequencies of 15 kHz (low frequency case) and 100 kHz (medium frequency case). For each one of these two frequencies, we assigned the corresponding $\sigma$ and $ϵ$ values to each 3D F.E. model’s element, according to which one of the mentioned tissues corresponds. To take into consideration possible inhomogeneity in each particular tissue, a standard deviation ($std$) of $2\%$ for the muscle-plasma background and the heart, $3\%$ for the lungs, and $1\%$ for the skin and bones were assumed. The nominal values assigned to each tissue along with the $std$ are shown in

Table 1. Assigned conductivity and permittivity values to the thoracic models’ tissues for $f=15$ kHz and $f=100$ kHz. The admittance is estimated as $\gamma =\sigma +j\omega ϵ{ϵ}_o$. For the skin and fat case, the average values of wet skin and breast fat admittances (see Figure 3) were used.

Tissue	$\mathit{\sigma }$ at 15 kHz (S/m)	$\mathit{\omega }\mathit{ϵ}{\mathit{ϵ}}_{\mathit{o}}$ at 15 kHz (F·Hz/m)	$\mathit{\sigma }$ at 100 kHz (S/m)	$\mathit{\omega }\mathit{ϵ}{\mathit{ϵ}}_{\mathit{o}}$ at 100 kHz (F·Hz/m)
Heart	$0.164\pm 0.003$	$0.041\pm 0.001$	$0.215\pm 0.004$	$0.0548\pm 0.001$
Inflated Lung	$0.095\pm 0.003$	$0.010\pm 0.000$	$0.107\pm 0.002$	$0.014\pm 0.000$
Deflated Lung	$0.247\pm 0.007$	$0.020\pm 0.001$	$0.272\pm 0.003$	$0.029\pm 0.001$
Bones	$0.021\pm 0.000$	$0.000\pm 0.000$	$0.021\pm 0.000$	$0.001\pm 0.000$
Skin and Fat	$0.015\pm 0.000$	$0.012\pm 0.000$	$0.045\pm 0.000$	$0.043\pm 0.000$
Muscle and Plasma	$0.350\pm 0.007$	$0.017\pm 0.001$	$0.380\pm 0.008$	$0.024\pm 0.001$

.

The deflated (i.e., state where the lungs are empty) 3D F.E. model’s geometry and tissues’ boundaries have been taken from a CT-scan of a healthy adult male at the 4th intercostal muscle height. The CT boundary and the corresponding 3D model are available online in the EIDORS library tool [26]. The F.E. model uses this geometry as a cross-section to create a 3D fine structure, with the x-axis normalized between $-1$ and 1, and the height set to 1.

For the inflated 3D F.E. thoracic model (i.e., state where the lungs are air-filled), we assumed a chest boundary movement of $10\%$ and a lung cross-section area increase of $15\%$ (see Figure 2a)). These values fall into the ranges described in [43]. The 3D F.E. model was created accordingly, as shown in Figure 2c). In both models, we have assumed two different electrode placements. In the first, $N_{el}=16$ electrodes were accurately placed on the 4th intercostal height (i.e., in the middle of each F.E. model’s total height, $z=1/2$). In the second one, an electrode position error has been added, assuming that the electrodes present a $5\%$ height $std$ and a $3\%$ angle std around their nominal values used in the first configuration. Furthermore, for both placements, circular electrodes of radius 0.05 (arbitrary unit—AU) and 0.03 AU are considered. These configurations were added in order to examine the effect of electrodes’ placement error and their size in the signal and the reconstructed images’ quality. The numbers of nodes and elements per model and electrode configuration are shown in

Table 2. Number of tetrahedral elements and nodes of each 3D model case.

Model	No of Elements (${\mathit{L}}_{\mathit{e}}$)	No of Nodes (${\mathit{n}}_{\mathit{e}}$)
Inflated, Uniform Electrodes, $R_{el}=0.05$	145,900	29,507
Inflated, Uniform Electrodes, $R_{el}=0.03$	133,756	26,861
Inflated, Non-Uniform Electrodes, $R_{el}=0.05$	146,000	29,542
Inflated, Non-Uniform Electrodes, $R_{el}=0.03$	135,330	27,120
Deflated, Uniform Electrodes, $R_{el}=0.05$	134,200	27,460
Deflated, Uniform Electrodes, $R_{el}=0.03$	133,756	24,849
Deflated, Non-Uniform Electrodes, $R_{el}=0.05$	133,529	27,328
Deflated, Non-Uniform Electrodes, $R_{el}=0.03$	119,654	23,965

.

4. Simulation Interface

In this section, we briefly describe the simulation approach used for this work. The interface used includes a MATLAB part which creates (with the assistance of NETGEN) the 3D models and assigns the tissue admittance values, loaded via Python from the database [41] to each element. MATLAB also uses the EIDORS library tool [26] with minor modifications to transfer the 3D thoracic models to RLC 16-port equivalent networks. These networks are integrated in LT SPICE with the whole EIT analogue circuitry, where transient simulations are executed. The transient measurements are imported in MATLAB, where the expected white noise is added and sampling and digital processing are simulated. From the final measurements acquired, the image reconstruction process is performed with the EIDORS tool. The overall process is described in [28,29]; however, for clarification, we synopsize the main parts. Furthermore, emphasis is given in the active electrode configuration, which was not simulated in previous works.

4.1. F.E.M. to RLC Equivalent Circuit Transformation

Considering the 3D thoracic F.E. model, we obtain the following system equation:

$$\left[\begin{array}{cc}{\mathit{A}}_{\mathit{m}}+{\mathit{A}}_{\mathit{z}}& {\mathit{A}}_{\mathit{v}}\\ {{\mathit{A}}_{\mathit{v}}}^{*}& {\mathit{A}}_{\mathit{D}}\end{array}\right]\left[\begin{array}{c}\mathit{U}\\ {\mathit{V}}_{\mathit{l}}\end{array}\right]=\left[\begin{array}{c}\mathbf{0}\\ {\mathit{I}}^{\mathit{d}}\end{array}\right]$$

(1)

where $\mathit{U}={\left[U_i\right]}_{i=1}^{n_e}$ are the nodal potentials for the corresponding model, ${\mathit{A}}_{\mathit{m}}\in R^{n_e\times n_e}$ is the element-assembled matrix, ${\mathit{A}}_{\mathit{z}}\in R^{n_e\times n_e}$ and ${\mathit{A}}_{\mathit{v}}\in R^{n_e\times N_{el}}$ express the complete electrode model (CEM) boundary conditions [44], and ${\mathit{A}}_{\mathit{d}}\in R^{N_{el}\times N_{el}}$ refers to the nodes where the measurements are taken [45]. Furthermore, ${\mathit{V}}_{\mathit{l}}\mathbf{\in }{\mathit{N}}_{\mathit{el}}$ is the voltage measurement vector, and ${\mathit{I}}^{\mathit{d}}$ is the current injected per electrode. By writing ${\mathit{A}}_{\mathit{k}}={\mathit{A}}_{\mathit{m}}+{\mathit{A}}_{\mathit{z}}$, it can be proved that [25,28]

$$({\mathit{A}}_{\mathit{d}}-{{\mathit{A}}_{\mathit{v}}}^{*}{{\mathit{A}}_{\mathit{k}}}^{-1}{\mathit{A}}_{\mathit{v}}){\mathit{V}}_{\mathit{l}}={\mathit{I}}_{\mathit{d}}$$

(2)

By setting $\mathit{F}\left(\mathit{\sigma }\right)={\mathit{A}}_{\mathit{d}}-{{\mathit{A}}_{\mathit{v}}}^{*}{{\mathit{A}}_{\mathit{k}}}^{-1}{\mathit{A}}_{\mathit{v}}$, we obtain a sparse, almost symmetric conductivity (or admittance) matrix. Therefore, the impedance between two nodes i and j (where voltage measurements are taken) is computed as [25,28,46]

$$Z_{ij}=-1/F{\left(\sigma \right)}_{ij}$$

(3)

In the case where admittances are considered, $\mathit{Z}\in C^{N_{el}\times N_{el}}$. Thus, each impedance $Z_{ij}$ can be implemented by a parallel combination of a resistor and either a capacitor or an inductor. In this way, we compose an RLC equivalent circuit, where each terminal corresponds to an EIT electrode.

4.2. EIT SPICE Circuitry

The SPICE EIT circuitry used for the simulations is mainly based on the basic structure presented in Section 2. It is implemented according to the passive electrode configuration (i.e., the analogue circuitry is connected to the electrode channels via cables and switches with parasitic capacitors) and the partially active electrode configuration (i.e., the first voltage readout stage—usually a buffer—is placed on each particular electrode, significantly reducing the parasitic effects).

The input signal is created in MATLAB in discretized (vector) form, demonstrating a $L_{DAC}=16$-bit DAC look-up-table (LUT) [28]. LUTs for 15 kHz and 100 kHz sinusoidal waveforms are considered, while the DAC sampling rate is simulated at 16 times higher than the corresponding input signal’s frequency. Each input waveform is stored in a PWL file, which is transferred to LT SPICE. In addition, the multiplexer digital inputs, according to the measurement pattern adopted, are also stored in PWL files, read by LT SPICE.

The analogue SPICE circuitry current injecting front-end consists of a second-order multiple feedback wideband band pass filter (BPF), a fully differential voltage output driver (THS413-Texas Instruments), and an enhanced mirrored Howland current pump (HPC) which acts as a VCCS. The voltage recording front-end is implemented as a partially active electrode, exclusively in the 100 kHz input signal case, and as a typical passive electrode for both the 15 kHz and 100 kHz input signal cases. In the active electrode configuration, each electrode is directly connected to a buffer (which is placed very close to the electrode). Then, each buffer’s output is connected to the IA through analogue (bi-directional) multiplexer switches (ADG426 analogue devices) and two first-order high-pass $RC$ filters (one per IA input channel). In the passive electrode configuration, each electrode is directly connected to two ADG426 analogue multiplexers, one for each differential measurement channel. Each ADG426 output is connected to a buffer and each buffer’s output is driven to the corresponding input of the AD8421 IA through a first-order high-pass $RC$ filter. The mirrored HPC VCCS is connected to the electrodes through two ADG426 multiplexers (one for the source and one for the sink electrode) in both the active and passive electrode configurations. The multiplexer switches on resistors are included in the ADG426 SPICE model, while parasitic capacitances (randomly chosen in the range of 100–200 pF) between the electrode channels have been manually added in SPICE, in order to simulate parasitic effects. Furthermore, to include channel imbalance effects at the output (i.e., common signal effects), unequal resistors (20–60 $\mathsf{\Omega }$) have been introduced to each measurement channel. Finally, the electrodes are modelled as in [28], according to the descriptions in [47,48,49] (assuming AgCl electrodes). A brief schematic of both configurations, including the parasitic effects and channel impedances, is presented in Figure 4.

4.3. Sampling and Digital Signal Processing

A transient simulation is performed in LT SPICE for each image frame. Its total duration depends on the input signal’s frequency as well as the selected time between two particular differential voltage measurements. The simulation stops when all the measurements defined by the selected pattern (adjacent) have been acquired. The transient results are written on a text file and sent to MATLAB for processing. The input reference signal is also measured from an IA and sent to MATLAB, in order to sample it and demodulate it with the output signal. Finally, a noise simulation at the SPICE circuit produces a power spectral density diagram (PSD). The total amount of noise, resulting from the integration of PSD, is added to the IA signal output in MATLAB [29].

A MATLAB software program simulates the ADC function in the following way. First, the number of bits $L_{ADC}$ and the sampling frequency $f_s$ are selected. Assuming negligible conversion times, the time points where the signal must be sampled are defined for each measurement window (time space between two switch changes) [28]. Then, for each defined time point, the nearest corresponding SPICE time value is found and the corresponding transient voltage value is sampled (note that the longest timestep in LT SPICE simulations is set to 100 ns, many times smaller than the ADC sampling period which is 625 ns at maximum). The sampling process is explained in detail in [28].

The sampled values from the input and output signals are stored in binary form (dropping the two LSBs which are assumed as noisy), and IQ demodulation is performed (see Figure 1). The multiplication results are driven to an adder, which actually acts as a digital low pass filter (LPF). From the in-phase and quadrature results, the signal amplitude of each measurement window is estimated. The final set of amplitudes is then used as data for the difference image reconstruction.

The whole simulation interface is shown in Figure 5. It includes the MATLAB/NETGEN modelling interface, where the tissues’ electrical properties are loaded from the mentioned database to the two 3D simulation structures (deflated and inflated). These structures are transformed to an RLC equivalent circuit, which is simulated at LT SPICE, along with the whole analogue EIT circuitry, which is shown in a “white box”. Any noise, parasitic effects, channel mismatches, and components’ non-idealities are added in the SPICE circuitry. The SPICE transient measurements are transferred in MATLAB to be processed as described before, and the reconstruction is performed with the EIDORS tool.

It is important to note that both the simulated analogue circuitry and the digital part are not optimal and are not necessarily recommended for novel EIT system designs, since this is beyond the scope of this work. Instead, the simulations presented compare the expected performance of state-of-the-art EIT hardware designs at different configurations in lung imaging. However, the simulation interface proposed can be easily modified to simulate the performance of new hardware designs.

5. Reconstruction and Evaluation Method

In this section, a description of the image reconstruction approach is written, presenting the selected algorithm and its parameters. Furthermore, a method for quantitative evaluation of the produced images is presented, in order to perform a more considered comparison.

5.1. Image Reconstruction

Assume that $\mathit{\delta }\mathit{V}$ is the difference voltage data, acquired from two continuous measurement frames (e.g., from a thoracic deflated and inflated state, respectively). We also assume a 2D EIT coarse-element reconstruction domain $\mathsf{\Omega }$, consisting of $L_c$ triangular elements. This domain’s elliptic boundary $\partial \mathsf{\Omega }$ differs from the 3D thoracic F.E. (both deflated and inflated states) boundaries since (A) the SUT boundary is not accurately known in real EIT lung imaging applications, and (B) to obtain a robust reconstruction result, we need to avoid the well-known inverse crime [6]. Then, the conventional difference-EIT reconstruction is considered as an optimization problem, where we have to find the optimal $\mathit{\delta }\mathit{\gamma }$ (admittance change) distribution that minimizes the following quantity [50]:

$$\mathit{\delta }{\mathit{\gamma }}_{\mathbf{*}}=\underset{\gamma \in C^{L_c}}{\mathrm{argmin}}\left\{F\left(\mathit{\delta }\mathit{\gamma }\right)\right\}$$

(4)

where

$$F\left(\mathit{\gamma }\right)={∥{\mathit{U}}_{\mathbf{2}}({\mathit{\gamma }}_{\mathit{o}}+\mathit{\delta }\mathit{\gamma })-{\mathit{U}}_{\mathbf{1}}\left({\mathit{\gamma }}_{\mathit{o}}\right)-\mathit{\delta }\mathit{V}∥}_{\mathit{W}}^2+{\lambda }^{2}P\left(\mathit{\delta }\mathit{\gamma }\right).$$

(5)

In the expression above, ${\mathit{U}}_{\mathbf{1}}\left({\mathit{\gamma }}_{\mathit{o}}\right)$ and ${\mathit{U}}_{\mathbf{2}}({\mathit{\gamma }}_{\mathit{o}}+\mathit{\delta }\mathit{\gamma })$ are the simulated electrode voltage vectors at the two states (assuming admittance distributions ${\mathit{\gamma }}_{\mathit{o}}$ and ${\mathit{\gamma }}_{\mathit{o}}+\mathit{\delta }\mathit{\gamma }$ in $\mathsf{\Omega }$, respectively), $\mathit{W}$ is a weighting matrix, $\lambda$ is the regularization hyperparameter, and $P\left(\mathit{\delta }\mathit{\gamma }\right)$ is the regularization term [1,3,50]. We note that regularization is essential, since the EIT reconstruction problem is both ill-posed and ill-conditioned [1,6]. It is also a heavily non-linear problem.

For this particular work, we use the linearized difference-EIT reconstruction with a smooth regularization term. This is the most common approach in dynamic thoracic imaging, since it enacts only one reconstruction step, allowing real-time imaging. However, it often suffers from accuracy, since the admittance changes during the breath cycle are usually significant [43]. The linearized difference-EIT reconstruction assumes small admittance changes near ${\mathit{\gamma }}_{\mathit{o}}$; hence, $\mathit{J}\mathit{\delta }\mathit{\gamma }\simeq {\mathit{U}}_{\mathbf{2}}({\mathit{\gamma }}_{\mathit{o}}+\mathit{\delta }\mathit{\gamma })-{\mathit{U}}_{\mathbf{1}}\left({\mathit{\gamma }}_{\mathit{o}}\right)$, where

$$\mathit{J}=\frac{\partial \left({\mathit{U}}_{\mathbf{2}}({\mathit{\gamma }}_{\mathit{o}}+\mathit{\delta }\mathit{\gamma })-{\mathit{U}}_{\mathbf{1}}\left({\mathit{\gamma }}_{\mathit{o}}\right)\right)}{\partial \mathit{\delta }\mathit{\gamma }}$$

(6)

is the Jacobian matrix. The optimization problem (4)–(5) then becomes

$$F\left(\mathit{\gamma }\right)={∥\mathit{J}\mathit{\delta }\mathit{\gamma }-\mathit{\delta }\mathit{V}∥}_{\mathit{W}}^2+{\lambda }^2{∥\mathit{\delta }\mathit{\gamma }∥}_{\mathit{Q}}^2.$$

(7)

where $\mathit{Q}$ is the prior regularization matrix [50]. We enact generalized Tikhonov regularization, obtaining the following single-step solution [3,50]:

$$\mathit{\delta }\mathit{\gamma }={\left({\mathit{J}}^{\mathit{T}}\mathit{W}\mathit{J}+{\lambda }^2\mathit{Q}\right)}^{-1}\left({\mathit{J}}^{\mathit{T}}\mathit{W}\mathit{\delta }\mathit{V}\right)$$

(8)

In terms of this work, the Laplace smoothness prior is used for the reconstruction [50]. However, other priors, such as the NOSER, standard Tikhnonov, or Gaussian could be also used to obtain the EIT images without loss of generality. The hyperparameter has been heuristically selected to $\lambda =0.1$ for all the cases. We note that the selection of single-step generalized Tikhonov regularization is indicative, and more recent approaches, such as learning approaches [51,52,53], could also be applied and can be considered for future work.

The reconstruction domain $\mathsf{\Omega }$ is a $L_c=1024$ triangular element 2D thoracic shape mesh, shown in Figure 6a. The gap electrode model has been applied in this domain to model the electrodes [54,55]. Shape mismatch between the true (cross-section 3D) boundary and the reconstruction domain’s boundary ($\partial \mathsf{\Omega }$, Figure 6b) introduces a minor modelling error. However, the exact thoracic shape is never accurately known, and it also changes during the breath cycle [4,11].

5.2. Image Evaluation Method

To evaluate the reconstructed images, we have to set a reference image for the comparisons on the same domain $\mathsf{\Omega }$. $\partial \mathsf{\Omega }$, however, differs from both the 3D models’ boundary shapes. As a result, it is necessary to primarily set a “true boundary” 2D reference, based on the cross-section of the 3D models. At this point, we take into consideration that (A) the 3D model’s shape is not constant, and (B) we perform difference-EIT imaging; thus, we reconstruct a single EIT image which represents the conductivity (and permittivity) differences between the two states. To this end, we assume two “true boundary” cross-section 2D references, one for the deflated and one for the inflated case. Then, the $\mathsf{\Omega }$ domain reference image is created in the following way. First, we scale the $\mathsf{\Omega }$ and each “true boundary” domain in the y-axis, matching them as shown in Figure 6b. Their horizontal limits are also normalized between $-1$ and 1. Secondly, we assume $A_i$ to be the $\mathsf{\Omega }$ domain’s ith-element area ($i=\{1,2,\dots ,L_e\}$). For each element i, we find the following weight vector:

$${\mathit{w}}_{\mathit{d}}^{\mathit{i}}={\left[w_{j,d}^i\right]}_{j=1}^6\in R^{6\times 1}$$

(9)

which represents the percentage of $A_i$ which is included in the following corresponding six tissue curves of the deflated state “true” domain: left lung, right lung, vertebra, heart, skin, and muscle-plasma. The corresponding vector which represents the tissues’ (mean) admittances at the deflated state is defined by

$${\mathit{\gamma }}_{\mathit{t},\mathit{d}}={\left[{\gamma }_{j,d}\right]}_{j=1}^6\in C^{1\times 6}$$

(10)

Then, we compute the corresponding reference admittance as (see Figure 7a)

$${\gamma }_{r,d,i}={\mathit{\gamma }}_{\mathit{t},\mathit{d}}·{\mathit{w}}_{\mathit{d}}^{\mathit{i}}\in C$$

(11)

In the same way, we compute each i element’s reference admittance for the inflated state:

$${\gamma }_{r,f,i}={\mathit{\gamma }}_{\mathit{t},\mathit{f}}·{\mathit{w}}_{\mathit{f}}^{\mathit{i}}\in C$$

(12)

where

$${\mathit{w}}_{\mathit{f}}^{\mathit{i}}={\left[w_{j,f}^i\right]}_{j=1}^6\in C^{6\times 1}$$

(13)

is the ith element’s weight vector at the inflated state and

$${\mathit{\gamma }}_{\mathit{t},\mathit{f}}={\left[{\gamma }_{j,f}\right]}_{j=1}^6\in C^{1\times 6}$$

(14)

the vector which represents the tissues’ (mean) admittances at the inflated state. We set the deflated and inflated state reference admittances for the domain $\mathsf{\Omega }$:

$$\begin{array}{cc}\hfill {\mathit{\gamma }}_{\mathit{r},\mathit{d}}& ={\left[{\gamma }_{r,d,i}\right]}_{i=1}^{L_e}\in C^{L_c\times 1}\hfill \\ \hfill {\mathit{\gamma }}_{\mathit{r},\mathit{f}}& ={\left[{\gamma }_{r,f,i}\right]}_{i=1}^{L_e}\in C^{L_c\times 1}\hfill \end{array}$$

(15)

The final $\mathsf{\Omega }$ domain reference has the following admittance difference values:

$${\mathit{\gamma }}_{\mathit{r}}={\mathit{\gamma }}_{\mathit{r},\mathit{f}}-{\mathit{\gamma }}_{\mathit{r},\mathit{d}}$$

(16)

where ${\mathit{\gamma }}_{\mathit{r}}\in C^{L_e\times 1}$. The reference image obtained is shown in Figure 7b.

At both the reference and the reconstructed images, the $\mathit{\delta }\mathit{s}={\left[∥\delta {\gamma }_i∥\right]}_{i=1}^{L_c}\in R^{L_c\times 1}$ (amplitude) term is demonstrated. For the quantitative evaluation and comparisons, the correlation coefficient $CC$ metric has been considered, defined as follows:

$$CC=\frac{\mathrm{Cov}(\mathit{\delta }{\mathit{s}}_{\mathit{r}},\widehat{\mathit{\delta }\mathit{s}})}{\mathrm{Std}\left(\mathit{\delta }{\mathit{s}}_{\mathit{r}}\right)\mathrm{Std}\left(\widehat{\mathit{\delta }\mathit{s}}\right)}$$

(17)

where $\mathit{\delta }{\mathit{s}}_{\mathit{r}}$ is the ground truth admittance amplitude, $\widehat{\mathit{\delta }\mathit{s}}$ is the admittance amplitude distribution estimated from the image reconstruction, and Cov is the covariance between reference and estimated values. The value of $CC$ ranges between 0 and 1, and higher $CC$ indicates more similarity between the reference and the reconstructed image.

In addition, to enhance the quantitative evaluation, the relative reconstruction error $RRE$ metric has also been used, defined as follows [52,53]:

$$RRE=\frac{{∥\mathit{\delta }{\mathit{s}}_{\mathit{r}}-\widehat{\mathit{\delta }\mathit{s}}∥}_2}{{∥\mathit{\delta }{\mathit{s}}_{\mathit{r}}∥}_2}$$

(18)

We note that for the calculation of $RRE$, both $\mathit{\delta }{\mathit{s}}_{\mathit{r}}$ and $\widehat{\mathit{\delta }\mathit{s}}$ were normalized between $-1$ and 1. Higher $RRE$ indicates a larger error between the reference and the reconstructed image.

6. Results and Discussion

In this section, the simulation cases are firstly defined, based on the electrodes’ properties, dimensions, and deviation and the digital acquisition part’s configurations. The theoretically expected $SNR$s are also computed for each ADC configuration case. In addition, the obtained images per case are presented, as well as the resulting $CC$ metrics. Furthermore, both qualitative and quantitative results are discussed.

6.1. Simulation Cases

For the simulations, we consider two input sinusoidal frequencies: $f=15$ kHz and $f=100$ kHz (see Section 4.2). For the $f=15$ kHz case, SPICE simulations are carried out using the passive electrode configuration, while for the $f=100$ kHz, both passive and active electrode configurations are considered. This is because parasitic effects are much more significant at $f=100$ kHz than at $f=15$ kHz. Furthermore, for each of the previous three cases, we examine (A) accurate electrode placement on the 4th intercostal height, (B) electrode deviation ($5\%$ height and $3\%$ angle), (C) circular electrodes of $0.03$ (A.U.) radius, and (D) circular electrodes of $0.05$ (A.U.) radius (for further details, see Section 3). Moreover, we simulate the following ADC cases: (A) resolution of $L_{ADC}=12$ bits and sampling rate of $f_s=4f$, (B) $L_{ADC}=16$ bits and sampling rate of $f_s=4f$, (C) $L_{ADC}=12$ bits and sampling rate of $f_s=16f$, and (D) $L_{ADC}=16$ bits and sampling rate of $f_s=16f$. Finally, we also consider $N_T=2$ and $N_T=4$ sine periods to be sampled per voltage channel measurement (see [29] and Figure 8). Hence, we obtain a total of 96 sub-cases for imaging, evaluation, and comparison.

The expected voltage amplitude $SNR$ can be computed according to the following formula [20]:

$$SNR=10lo{g}_{10}\left(\frac{A^2(K/2)}{\frac{V_{ref}^2}{12·2^{2L_{ADC}}}+{\sigma }_n^2}\right)$$

(19)

where A is the mean acquired voltage signal amplitude, K the number of matched digital LPF filter taps (samples) at the IQ demodulation part, $V_{ref}$ the ADC voltage reference ($3.3$ V in the simulated cases), $L_{ADC}$ the ADC resolution (bits), and ${\sigma }_n$ the circuit noise std. The ratio $V_{ref}/2^{L_{ADC}}$ denotes the LSB voltage range. Furthermore, the number of taps per measurement is defined by $K=N_T·f_s/f$. As voltage signal amplitude A, we take the mean amplitude between all the channel measurements (208 for the adjacent pattern when 16 electrodes are attached [1,56]). For the $SNR$ estimation, we also assume that it is not considerably affected by the electrode configurations or the patient’s state. However, if we consider different topologies than in Figure 4, possible changes in noise levels between active and passive electrode configuration should be considered (see [4]). Moreover, an increase in signal frequency has a substantial role in $SNR$ reduction, firstly due to the need for higher sampling frequencies $f_s$ and secondly due to the lower voltages acquired, resulting from the increased tissue admittances (see Figure 3 and Table 1). Therefore, we examine 16 different $SNR$ cases, as demonstrated in

Table 3. Calculated mean SNR values (dB) for each measuring parameter case, according to (19).

f (kHz)	${\mathit{L}}_{\mathit{ADC}}$ (bits)	${\mathit{f}}_{\mathit{s}}$	${\mathit{N}}_{\mathit{T}}$	$\mathit{SNR}$ (dB)
15	12	$4f$	2	$27.5$
15	12	$4f$	4	$30.5$
15	12	$16f$	2	$33.5$
15	12	$16f$	4	$36.5$
15	16	$4f$	2	$27.5$
15	16	$4f$	4	$30.5$
15	16	$16f$	2	$33.5$
15	16	$16f$	4	$36.6$
100	12	$4f$	2	$20.8$
100	12	$4f$	4	$23.8$
100	12	$16f$	2	$26.8$
100	12	$16f$	4	$29.8$
100	16	$4f$	2	$20.8$
100	16	$4f$	4	$23.8$
100	16	$16f$	2	$26.8$
100	16	$16f$	4	$29.9$

. A ${\sigma }_n$ value of $4.2m{V}_{p-p}$ was extracted from LT SPICE noise simulations. This value strongly depends on the input current amplitude and the instrumentation amplifier’s gain (set as $2m{A}_{p-p}$ and 200 V/V, respectively, for these particular simulations). Table 3 shows the calculated voltage amplitude $SNR$ per case.

As demonstrated, the expected voltage $SNR$ is reduced as the input signal frequency increases from 15 kHz to 100 kHz. At the same time, increasing the sampling frequency from 4 to 16 times the signal frequency and the sampled periods from 2 to 4 (i.e., the number of taps/samples per measurement) leads to improved $SNR$ values. In addition, an increase in the ADC resolution did not have a noticeable effect on $SNR$ in the cases tested. This can be explained by observing (19) denominator term. It is demonstrated that if noise level ${\sigma }_n^2$ from the analogue circuitry is significantly higher than the first term, then the ADC resolution does not have a substantial role in the SNR. Nonetheless, this noise can be filtered from the digital LPF by increasing the samples per measurement (K).

We note that these $SNR$s are the mean ones per each particular measurement frame. In fact, there is a large variation between the voltage channels’ measurement $SNR$s, due to the wide voltage amplitude range. For example, voltage electrode pairs close to the current electrode pair drive significantly higher voltages than the voltage electrode pairs on the other side of the setup, resulting in higher $SNR$ levels. It is also noted that the selected hardware circuitry, as well as the ratio between the voltages acquired and $V_{ref}$, is not optimized in this particular work. However, the simulation interface presented can be an effective tool for this direction.

6.2. Simulation Results

The images obtained for $f=15$ kHz and passive electrode configuration are shown in Figure 9. The corresponding $CC$ and $RRE$ values are shown for each subcase. In almost all the images, the admittance reduction due to the change between the end-inspiration (deflated) and end-expiration (inflated) states is visible. However, some artefacts are present, including a ringing effect between the lungs (which is often misinterpreted as “heart” [57]) or near the boundary. Furthermore, in some cases (especially when sampling with lower rates), the lung-related admittance change regions appear very close to the reconstruction domain’s centre.

The presence of artefacts is a result of hardware noise and issues (electrode’s contact impedances, common signals between the channels, etc.), boundary and tissues’ motion, and the ill-posed and ill-conditioned nature of the inverse problem (low noise levels can lead to strong artefacts) [5,6,58]. When lower sampling rates $f_s$ or when a lower number of sampling periods ($N_T$) is set, the intense change of the lungs’ admittance near the chest is often not properly detected. This is verified quantitatively, where for $f_s=60$ ksps and $N_T=2$, we obtain a $CC$ between 0.58 and 0.71, while for $f_s=240$ ksps and $N_T=4$, $CC$ takes values between 0.58 and 0.77. The electrode’s surface area seems to be an important factor, since smaller surfaces lead to higher contact impedances. The contact impedance effect is more intense at lower frequencies [59,60]. Mismatches between the channels’ contact impedances may lead to common signals, especially at the voltage recording circuitry part. This can explain the improved performance ($CC$ of 0.68–0.77 and $RRE$ of 0.50–0.64) observed for an electrode radius of $0.05$ compared to this for the $0.03$ radius electrode configuration ($CC$ of 0.58–0.70 and $RRE$ of 0.61–0.78). Furthermore, a small electrode position deviation overall leads to minor additional image artefacts, assuming the same gap electrode model for each reconstruction. This is anticipated, since we enact difference-EIT imaging, where minor modelling errors are compensated [61]. Finally, as expected, an increase in the ADC resolution $L_{ADC}$ for 12 to 16 bits does not lead to significant reconstruction improvement with the presence of relatively high white noise levels (see Section 6.1 and Table 3).

Overall, for the 15 kHz passive electrode case and high white noise levels, we conclude that the electrode radius and the total number of taps per measurement (proportional to $f_s$ and $N_T$) are two key factors for the image reconstruction performance. In addition, according to (19), for lower analogue circuitry noise levels, the ADC resolution is also expected to boost the performance. However, such a case has not been simulated in this particular work.

The reconstructed images for $f=100$ kHz and passive electrode configuration are demonstrated in Figure 10. First of all, it is observed that the overall spatial resolution and image quality are reduced, compared to the corresponding 15 kHz case. Although the electrodes appear to have reduced contact impedances at this frequency, the parasitic effects taking place between the channel cables become more critical. These effects, combined with the reduced VCCS output impedance and instrumentation amplifier’s CMRR, lead to image quality degradation [5,7]. However, due to the low contact impedances, the electrodes’ surface does not significantly affect the performance ($CC$ between $0.54$ and $0.74$, $RRE$ between $0.54$ and $0.74$ for the $0.03$ radius cases and $CC$ between $0.51$ and $0.76$, $RRE$ between $0.54$ and $0.71$ for the $0.05$ radius cases). As in the 15 kHz cases, the electrodes’ position deviation has a small impact on the images’ quality. The best results are obtained by increasing the number of taps per measurement (maximum $CC$ is observed for $f_s=1.6$ Msps and $N_T=4$ sampling sine periods. A small improvement is observed when increasing the ADC resolution from 12 to 16 bits.

The EIT images for $f=100$ kHz and the active (on the voltage recording part) electrode configuration are finally shown in Figure 11. A slight improvement of the image quality is detected, mainly due to the reduction of cable parasitic effects at the voltage measurement channels. The maximum $CC$ of $0.77$ (corresponding minimum $RRE$ of $0.52$) is achieved for non-deviated, $0.05$ radius electrodes, when a 16 bit resolution ADC is utilized, sampling at $1.6$ Msps for $N_T=4$ sine periods. As in the passive electrode cases, the performance is less affected by the electrode deviation and the ADC resolution and mostly affected by the total number of taps per measurement channel. In addition, the electrode surface has minor effect on the performance, due to the lower contact impedances presented at $f=100$ kHz.

6.3. Discussion

From the results presented above, we can observe the effects of contact impedances, parasitic effects, noise, and shape mismatch under different hardware configurations and measurement frequencies. Specifically, the contact impedance effect becomes more intense at lower frequencies. However, the absence of capacitive effects and attenuation allows a moderate to good performance with a passive electrode configuration, which is achieved at 15 kHz (see Figure 9).

At higher input sinusoidal signal frequencies, the use of an active electrode configuration is essential for achieving a sufficient image quality by reducing the capacitive effects. Such systems have been already implemented in [4,16,18]. Fully active electrode EIT systems that implement the current source near each corresponding electrode pair have shown even better performance, since apart from the voltage channels’ parasitic effects, they compensate for the current channels’ ones. In [4], a comparison between passive and fully active electrode configurations on a lung-imaging EIT system has shown superior performance for the second case over a high frequency span. Hence, higher signal frequency spans can be achieved with active electrodes that allow higher temporal analysis which is essential for real-time lung imaging applications.

Furthermore, high noise levels from the circuitry can be taken care of by increasing the ADC sampling frequency $f_s$ and the number of sine periods $N_T$ that are sampled per measurement window. This leads to the increase of the number of taps K per measurement and benefits the $SNR$ values, according to (19). In this way, the matched LPFs’ performance is improved (after IQ demodulation) by actually performing signal averaging. Some systems, such as [16], use a high K per measurement window and achieve good voltage $SNR$ levels. However, increasing $N_T$ also increases the total time needed to take a measurement frame. As a result, lower image frame rates are obtained, significantly reducing the imaging temporal resolution. Therefore, lung imaging has to be performed at frequency rates above 100 kHz and relatively high sampling rates $f_s$. In addition, these frequency ranges usually include valuable information about the thoracic tissues’ electrical behaviour [4]. Nonetheless, there is always a trade-off between signal noise and frame speed. Refs. [4,17], for example, achieve an imaging speed of more than 100 frames per second (fps) at a 1 MHz signal frequency; however, the $SNR$ rates fairly drop between 125 kHz and 1 MHz. On the other hand, ref. [16], which, as mentioned, uses a large number of filter taps, keeps a high $SNR$ by dropping the speed to lower than 20 fps. Moreover, another way to increase the $SNR$ is by enacting higher amplitude current signals. However, apart from the patients’ safety consideration, this solution demands higher power supplies and thus power consumption, which would be impractical in the case of active electrodes. The usage of a programmable gain amplifier (PGA), which amplifies the voltage signals collected from electrodes that are far from the current source pair, could also improve signal resistance to noise [7,23]. Finally, the effect of enacting different current patterns instead of the adjacent one to $SNR$ has been partially studied in [29]. Skip current patterns might increase the signal amplitude; however, sometimes they have a negative effect at the imaging reconstruction process, since they lead to a more ill-posed problem than the adjacent pattern. This happens due to the reduced number of electrode boundary measurements and the reduced number of independent measurements [34].

Finally, the shape mismatch between the chest (surface attached by electrode array) boundary and the reconstruction domain, as well as the boundary motion, leads to further image artefacts. Boundary tracking methods, such as the use of accelerometers, have recently been proposed to overcome this effect [4,11]. Although these methods can improve the image quality, they are often more expensive, due to the extra hardware needed. Furthermore, since the boundary motion due to the breath (or even the subject’s movement) is continuous, artefacts will be still present due to the fact that the reconstruction domain does not change through the imaging frames. Until now, such change in real-time thoracic EIT imaging applications is not applicable due to the complexity and time the mesh refinement needs. Therefore, compensation of motion-related imaging errors in EIT imaging applications is still a field under research.

A prior code implementation of the proposed interface can be found online on the following website: https://github.com/chdim100/SPICE-MATLAB-Interface-for-Electrical-Impedance-Tomography-Simulation-SPICEIT- (accessed on 4 August 2021).

7. Conclusions

An EIT simulation interface which makes use of Python, MATLAB, and LT SPICE software is applied in order to evaluate the impact of modelling, parasitic, and noise effects in thoracic imaging. The SUT is modelled as a fine F.E. structure, which is transformed to a multiport RLC circuitry block and merged in the analogue SPICE EIT circuitry. Transient and noise simulations are carried out in LT SPICE, while MATLAB creates the input data (input sinusoidal waveform and digital signals for the multiplexers’ control) and handles the output SPICE transient data to simulate the digital EIT hardware part. Configurations of the electrodes and the signal sampling are considered, and comparison of the results is performed. Simulations showed the significant effect of the channels’ impedances at low frequencies and passive-electrode configurations in reconstructions as well as the effect of the parasitic channels’ capacitors at higher frequencies. The active electrode configuration can compensate for these effects, while increasing the demodulation filter taps can increase the acquired signal $SNR$ and, as a result, the reconstructed image quality. This interface can be used as a valuable tool in EIT system design for prior estimation of custom EIT system performance under different analogue and digital part characteristics.