Very Low Resource Digital Implementation of Bioimpedance Analysis ^†

Abstract

Bioimpedance spectroscopy consists of measuring the complex impedance of biological tissues over a large frequency domain. This method is particularly convenient for physiological studies or health monitoring systems. For a wide range of applications, devices need to be portable, wearable or even implantable. Next generation of bioimpedance sensing systems thus require to be implemented with power and resource savings in mind. Impedance measurement methods are divided into two main categories. Some are based on “single-tone” signals while the others use “multi-tone” signals. The firsts benefit from a very simple analysis that may consist of synchronous demodulation. However, due to necessary frequency sweep, the total measurement may take a long time. On the other hand, generating a multi-frequency signal allows the seconds to cover the whole frequency range simultaneously. This is at the cost of a more complex analysis algorithm. This makes both approaches hardly suitable for embedded applications. In this paper, we propose an intermediate approach that combines the speed of multi-tone systems with a low-resource analysis algorithm. This results in a minimal implementation using only adders and synchronous adc. For optimal performances, this small footprint digital processing can be synthesized and embedded on a mixed-mode integrated circuit together with the analog front-end. Moreover, the proposed implementation is easily scalable to fit an arbitrary frequency range. We also show that the resulting impact on noise sensitivity can be mitigated.

1. Introduction

Bioimpedance spectroscopy consists of measuring the complex impedance of biological tissues over a large frequency domain [1]. This method is convenient in particular for studying body composition [2], blood characterization [3] and even cancer detection [4]. This wide range of applications makes it suitable as a part of health monitoring systems. Today’s self-monitoring devices tend to be portable, wearable or even implantable. Next, generation bioimpedance sensing systems thus require to be implemented with power and resource savings in mind.

Impedance measurement methods are divided into two main categories. Some are based on “single-tone” signals while the others use “multi-tone” signals. The firsts use a pure frequency sine wave to make the measurement [5,6,7,8]. They benefit from a very simple analysis that can consist of synchronous demodulation or sampling. However, the operation must be repeated for each frequency over the domain of interest. Due to this necessary frequency sweep, the total measurement may take a long time. On the other hand, generating a multi-frequency signal allows the analysis to cover the whole frequency range simultaneously [9,10,11]. This is at the cost of a more complex analysis algorithm (discrete cosine transform—dct, typically). Unfortunately, both methods result in excess power consumption: a long time of measurement for single-tone frequency sweep, hardware and computational resources for multi-tone. This makes both approaches hardly suitable for embedded applications. In 2008, Ronk and Toomessoo proposed an implementation of bioimpedance measurement using multi-frequencies on a field-programmable gate array (fpga) [12]. Their method is similar to dct, but use square-waves instead of sines. This avoids the need of multiplication in the processing.

Some intermediate methods of frequency analysis have been developed with resource savings in mind. They combine the speed of multitone measurement with a much simpler analysis algorithm than dct or fast Fourier transform (fft). For instance, Goertzel filters can be used to calculate frequency components [13]. These methods are particularly useful for applications like the built-in self test of integrated circuits [14].

Last year, we presented a similar approach with further simplifications thanks to the power-of-two frequency distribution [15] at the International Conference on Sensing Technology (icst 2018, Limerick, Ireland). The present paper proposes a slightly improved implementation of the concept resulting in a reduced (divided by two) clock frequency in most of the digital circuitry. As a consequence, the circuit itself avoids a few frequency dividers, all resulting in further power savings. Compared to the conference paper, the present article also adds a new frequency domain analysis of the concept, studies the noise impact on impedance values and gives an extended bibliography.

The rest of the paper is organized as follows. First, we consider some a priori facts about bioimpedance. Those allow us to simplify the resolution needs for a bioimpedance-specific sensor compared to a generic purpose impedance analyzer. Then, we compare classical approaches to impedance analysis for sparse logarithmic frequency distribution and we propose a new methodology. In the fourth section, we analyze the proposed method based on synchronous sampling demodulation and present the processing steps. The fifth section gives the proposal of a digital implementation of the algorithm using only adders and frequency dividers. In the sixth section, we study the impact of noise on impedance estimation. Finally, the last section opens some discussions about the limits of the proposed approach.

2. Conceptualization

Let’s put into evidence some general characteristics of bioimpedance. All biological tissues present similar frequency behaviors. Typically, the impedance is globally decreasing, with one or several relaxation domains [16] as we can see in Figure 1. Variations are smooth, with noticeably no resonance and only require sparse frequency resolution.

Typical bioimpedance can be modeled with resistors, capacitors and so-called constant phase elements (cpe) of impedance

$$\begin{array}{ccc}\hfill Z_{\mathrm{CPE}}(f)& =\frac{1}{q_0{(2\mathrm{i}\pi f)}^{\alpha }},\hfill & \hfill \mathrm{with}{q}_0\mathrm{arbitrary}\mathrm{parameter},\mathrm{and}0<\alpha <1.\end{array}$$

(1)

Physiological parameters (fat content, body water, tissue characteristics) are related to objective markers such as: relaxation frequencies, amplitude, and phase of cpe. Estimating such values is better achieved with logarithmically growing frequencies, as illustrated Figure 2. Designing a bioimpedance-specific sensor, we can adopt a frequency distribution of the form

$$\begin{array}{ccc}\hfill f_i& =\frac{f_0}{2^i},\hfill & \hfill 0\le i<n_{\mathrm{f}},\end{array}$$

(2)

for the multitone excitation signal. This gives $n_{\mathrm{f}}$ frequencies with a maximum of $f_0$ and a resolution of ${log}_2(10)\approx 3.32$ points per decade.

3. Methodology

The basic idea behind the proposed method of impedance estimation is the use of coherent sampling together with a multitone excitation signal. Classical approaches either use coherent sampling of a single-tone (pure frequency sine) or Fourier analysis (sine values multiplications) of a multi-tone signal.

The first case is summarized in Figure 3 where analysis is performed by down-sampling, then low-pass filtering the voltage signal. In the frequency domain, it consists of convoluting the acquired signal by a Dirac comb, then estimating the dc (zero-frequency) content. The signal must be acquired for at least one period of the lowest frequency signal, i.e., $\frac{2^{n_{\mathrm{f}}-1}}{f_0}$. Because the operation must be repeated for each frequency, the total acquisition time reaches $n_{\mathrm{f}}\times \frac{2^{n_{\mathrm{f}}-1}}{f_0}$.

The second case is summarized in Figure 4 where analysis is performed by multiplying the multi-frequency voltage signal by pure sine values, then low-pass filtering. In the frequency domain, it consists of convoluting the acquired signal by a unique Dirac, then estimating the dc content. In this case, the signal still must be acquired for at least the inverse of the lowest frequency, i.e., $\frac{2^{n_{\mathrm{f}}-1}}{f_0}$), but only one time.

Coherent sampling only requires one addition per period, and no multiplication. Thus, the total number of additions is ${\sum }_{i=0}^{n_{\mathrm{f}}-1}{2}^i=2^{n_{\mathrm{f}}}-1$. This number must be multiplied by two to get real and imaginary parts of the impedance.

As said before, bioimpedance spectroscopy can be performed using a somehow sparse number of frequency points. Therefore, dct reveals itself to be more appropriate than the fft algorithm. The processing requires as additions and multiplications as the number of samples, $2^{n_{\mathrm{f}}}$ (with the adc at the Shannon rate, $2f_0$), repeated for each of the $n_{\mathrm{f}}$ frequencies of interest. Again, this number must be multiplied by two for real and imaginary parts, i.e., $n_{\mathrm{f}}\times 2^{n_{\mathrm{f}}+1}$.

Mixing coherent sampling with multi-tone excitation results in short acquisition time while avoiding the need for multipliers (see

Table 1. Expected characteristics of impedance estimation methods.

			(Proposed Method)
	Single Freq.	Multi Freq.	Multi Freq.
	& Coherent Sampling	& DCT	& Coherent Sampling
Acquisition time (mini.)	$n_{\mathrm{f}}\times \frac{2^{n_{\mathrm{f}}-1}}{f_0}$	$\frac{2^{n_{\mathrm{f}}-1}}{f_0}$	$\frac{2^{n_{\mathrm{f}}-1}}{f_0}$
Additions	$2\times (2^{n_{\mathrm{f}}}-1)$	$n_{\mathrm{f}}\times 2^{n_{\mathrm{f}}+1}$	$2\times (2^{n_{\mathrm{f}}}-1)$
Multiplications	0	$n_{\mathrm{f}}\times 2^{n_{\mathrm{f}}+1}$	0

). This method should be implemented using the lowest power consumption. Unfortunately, Figure 5 shows that aliasing occurs in this case at null frequency, making it harder to process. However, we show in the following section that it is still possible to discriminate valuable information from the measured dc component.

4. Formal Analysis

The chosen principle of bioimpedance sensing is illustrated in Figure 6. With the above considerations, generated current $i(t)$ can be expressed as the superposition of $n_{\mathrm{f}}$ sine waves. Current being real, hermitian symmetry applies to complex amplitudes ($I_{-i}=I_i^{*}$) such that

$$i(t)=\sum _{i=0}^{n_{\mathrm{f}}-1}{I}_i{\mathrm{e}}^{2\mathrm{i}\pi \frac{f_0}{2^i}t}+I_i^{*}{\mathrm{e}}^{-2\mathrm{i}\pi \frac{f_0}{2^i}t}.$$

(3)

For the sake of simplicity, we can let all $I_i$ be equal to $\frac{1}{2}$ (normalized amplitudes and null phase for all frequency components). If not so, we can get back to this case by resistor calibration:

$$i(t)=\frac{1}{2}\sum _{i=0}^{n_{\mathrm{f}}-1}{\mathrm{e}}^{2\mathrm{i}\pi \frac{f_0}{2^i}t}+{\mathrm{e}}^{-2\mathrm{i}\pi \frac{f_0}{2^i}t}.$$

(4)

Let us express the Fourier transform $I(f)=\mathcal{F}\{i(t)\}$ in the frequency domain

$$I(f)=\frac{1}{2}\sum _{i=0}^{n_{\mathrm{f}}-1}\delta (f-\frac{f_0}{2^i})+\delta (f+\frac{f_0}{2^i}),$$

(5)

with $\delta (t)$ denoting the Dirac distribution. We can deduce from the last equation the voltage expression $V(f)=\mathcal{F}\{v(t)\}$ in the sinusoidal steady state:

$$\begin{array}{cc}\hfill V(f)& =Z(f)I(f)=\frac{1}{2}\sum _{i=0}^{n_{\mathrm{f}}-1}{Z}_i\delta (f-\frac{f_0}{2^i})+Z_i^{*}\delta (f+\frac{f_0}{2^i}),\hfill \end{array}$$

(6)

with $Z_i=Z(f_i)$ and $Z_i^{*}=Z(f_{-i})$.

The objective of the following processing is to extract the values $Z_i=Z(f_i)$ from the acquired $v(t)$ signal. The main idea is to use synchronous sampling demodulation. Indeed, $Z(f_0)$ can be obtained by sampling $v(t)$ at the higher frequency $f_0$ and by low-pass filtering (averaging) the resulting samples. However, sampling the signal at lower frequencies results in aliasing. Let us explicitly express the result of sampling $v(t)$ at $f_j=\frac{f_0}{2^j}$. In the frequency domain, sampling is equivalent to a convolution product by a Dirac comb:

$$V(f)\ast \sum _{k\in \mathbb{Z}}\delta (f-k\frac{f_0}{2^j})=\frac{1}{2}\sum _{k\in \mathbb{Z}}\sum _{i=0}^{n_{\mathrm{f}}-1}{Z}_i\delta \left(f-f_0(\frac{k}{2^j}+\frac{1}{2^i})\right)+Z_i^{*}\delta \left(f-f_0(\frac{k}{2^j}-\frac{1}{2^i})\right).$$

(7)

We isolate from the last equation the dc component $S_j$, i.e., the coefficient of $\delta (f)$. For all i between 0 and j, there are two integers $k=\pm 2^{j-i}$ that satisfy

$$\frac{1}{2^i}\pm \frac{k}{2^j}=0.$$

(8)

Summing all contributions, we get

$$S_j=\frac{1}{2}\sum _{i=0}^j(Z_i+Z_i^{*})=\sum _{i=0}^j\Re (Z_i).$$

(9)

In the same manner, we can get the quadrature components by delaying the sample time by a quarter of period $\frac{2^j}{4f_0}$:

$$V(f)\ast \sum _{k\in \mathbb{Z}}\delta (f-k\frac{f_0}{2^j}){\mathrm{e}}^{-2\mathrm{i}\pi \frac{2^j}{4f_0}f}=\frac{1}{2}\sum _{k\in \mathbb{Z}}\sum _{i=0}^{n_{\mathrm{f}}-1}{Z}_i\delta \left(f-f_0(\frac{k}{2^j}+\frac{1}{2^i})\right){\mathrm{e}}^{-\mathrm{i}\pi \frac{k}{2}}+Z_i^{*}\delta \left(f-f_0(\frac{k}{2^j}-\frac{1}{2^i})\right){\mathrm{e}}^{-\mathrm{i}\pi \frac{k}{2}}.$$

(10)

Again, for $k=\pm 2^{j-i}$, we get the dc component

$$Q_j=\frac{1}{2}\sum _{i=0}^j{Z}_i{\mathrm{e}}^{-\mathrm{i}\pi \frac{-2^{j-i}}{2}}+Z_i^{*}{\mathrm{e}}^{-\mathrm{i}\pi \frac{2^{j-i}}{2}}.$$

(11)

If $j\ge 2$, we can isolate the terms corresponding to $i=j$ and $i=j-1$:

$$Q_j=\frac{1}{2}\sum _{i=0}^{j-2}{Z}_i{\mathrm{e}}^{\mathrm{i}\pi 2^{(j-i-1)}}+Z_i^{*}{\mathrm{e}}^{-\mathrm{i}\pi 2^{(j-i-1)}}+\frac{1}{2}(Z_{j-1}{\mathrm{e}}^{\mathrm{i}\pi }+Z_{j-1}^{*}{\mathrm{e}}^{-\mathrm{i}\pi })+\frac{1}{2}(\mathrm{i}{Z}_j-\mathrm{i}{Z}_j^{*}).$$

(12)

Thus,

$$Q_j=\sum _{i=0}^{j-2}\Re (Z_i)-\Re (Z_{j-1})-\Im (Z_j).$$

(13)

From Equations (9), (11) and (13), we can recursively calculate real and imaginary parts of $Z_j$.

$$\begin{array}{cccc}\hfill \Re (Z_0)& =S_0,\hfill & \hfill \Im (Z_0)& =-Q_0,\hfill \end{array}$$

(14)

$$\begin{array}{cccc}\hfill \Re (Z_1)& =S_1-S_0,\hfill & \hfill \Im (Z_1)& =-Q_1-S_0,\hfill \end{array}$$

(15)

$$\begin{array}{cccc}\hfill \Re (Z_j)& =S_j-S_{j-1},\hfill & \hfill \Im (Z_j)& =2S_{j-2}-S_{j-1}-Q_j.\hfill \end{array}$$

(16)

5. Numerical Simulation

In order to validate the formal analysis described in Section 4, we have performed some numerical simulations using the open-source computational software Scilab. The code provided in Listing 1 tests the algorithm for eight frequencies and a sample rate of $1\mathrm{MHz}$. To make results more realistic, we use measurements from tuna muscle bioimpedance (Figure 1) as the reference model. As expected, estimated values by decimation–accumulation and post-processing are in perfect agreement with the reference (

Table 2. Results of the numerical simulation performed with Scilab compared with reference values.

Frequencies	Reference	Estimated	Reference	Estimated
(Hz)	Modulus (Ohm)	Modulus (Ohm)	Phase (Rad)	Phase (Rad)
250,000	235.55872	235.55872	−0.0740422	−0.0740422
125,000	242.33682	242.33682	−0.0852266	−0.0852266
62,500	251.06561	251.06561	−0.1040394	−0.1040394
31,250	262.43457	262.43457	−0.1302557	−0.1302557
15,625	277.09086	277.09086	−0.1660522	−0.1660522
7812.5	295.96917	295.96917	−0.2183331	−0.2183331
3906.25	321.75482	321.75482	−0.3006170	−0.3006170
1953.125	362.5306	362.5306	−0.4326722	−0.4326722

).

6. Implementation

The main advantage of the proposed method is the simplicity of its digital implementation. Like classical synchronous demodulation, it only requires an adc sharing the same clock as the signal generation (Figure 6). Samples are fed by the adc at the rate of $f_{\mathrm{clk}}=4f_0$. Then, they are re-sampled at each frequency $f_i$ composing the multitone signal. The extraction of the dc component can be easily done by averaging the samples over a period of the whole multitone signal. Then, sample decimation is performed between each stage by dividing the frequency by two (Figure 7).

Practically, these functions are grouped in a block called Demod for demodulator that can be easily instantiated from a generic description (Figure 8). In more detail, at each stage, we generate two clock signals in quadrature (Figure 9). They control two identical accumulators (Figure 10). Samples are accumulated at the rate corresponding to a particular frequency $f_i$. The adder size is adapted to avoid overflow. Then, dividing by $\frac{1}{2^{n_{\mathrm{f}}-j-1}}$ is done by keeping only the $n_b$ most significant bits (msb). Thus, all $S_j$ and $Q_j$ accumulation are eventually encoded with the same precision.

Values of $S_j$ (in-phase) and $Q_j$ (in-quadrature) are saved in the output register at the end of the global period ($\frac{2^{n_{\mathrm{f}}}}{f_0}$). This is achieved using a simple frequency divider on the clock output of the last demodulator. We choose not to embed the computation of actual real and imaginary parts of the impedance as described by Equations (14) to (16). It can be done without transmission penalty (same size of data) as a post-process. The objective is to keep the digital circuit as simple as possible in a low-power constraint.

It must be noticed that the proposed algorithm is mathematically equivalent to a bank of finite impulse response (fir) filters. This is in particular made explicit in the implementation proposed by [17] that uses a similar multiplication-less approach. However, the straightforward implementation of such filters increases a lot the complexity of embedded computations and lacks genericity due to differing coefficients of each filter.

In our design, the parallel architecture can smoothly scale with the number $n_{\mathrm{f}}$ of frequency components. Moreover, each stage only uses two adders together with a few flip-flops to divide the frequency. For optimal performances, it is totally realistic to embed this digital processing directly with the analog front-end (current source and sensing amplifier) in a mixed-mode application-specific integrated circuit (asic).

7. Noise Sensitivity

The proposed architecture raises the question of the effect of noise on measurements. In particular, low frequencies suffer from averaging low numbers of samples. The worst occurring for $f_{n_{\mathrm{f}}-1}=\frac{f_0}{2^{n_{\mathrm{f}}-1}}$ where only one sample is used to calculate each of $S_{n_{\mathrm{f}}-1}$ and $Q_{n_{\mathrm{f}}-1}$. In this section, we show that, in the case of bioimpedance spectroscopy, this effect can be mitigated.

We now model perturbations by additive white noise n. Let ${\sigma }^2$ be the noise power. Voltage samples from the adc can be written $v[i]+n[i]$. At the end of the period, in-phase and quadrature accumulators of rank j have summed $2^{n_{\mathrm{f}}-j-1}$ samples:

$$\begin{array}{cc}\hfill \widehat{S_j}& =\frac{1}{2^{n_{\mathrm{f}}-j-1}}\sum _{i=0}^{2^{n_{\mathrm{f}}-j-1}-1}(v[2^{j+2}i]+n[2^{j+2}i]),\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill \widehat{Q_j}& =\frac{1}{2^{n_{\mathrm{f}}-j-1}}\sum _{i=0}^{2^{n_{\mathrm{f}}-j-1}-1}(v[2^j(4i+1)]+n[2^j(4i+1)]).\hfill \end{array}$$

(18)

From Equation (15), the estimated value of the impedance real part is

$$\Re ({\widehat{Z}}_j)={\widehat{S}}_j-{\widehat{S}}_{j-1}.$$

(19)

This combination contains $2^{n_{\mathrm{f}}-j-1}$ noise samples, thus the signal power

$$\mathrm{P}\{\Re ({\widehat{Z}}_j)\}=\Re {(Z_j)}^2+2^{n_{\mathrm{f}}-j-1}\frac{{\sigma }^2}{{(2^{n_{\mathrm{f}}-j-1})}^2}=\Re {(Z_j)}^2+\frac{{\sigma }^2}{(2^{n_{\mathrm{f}}-j-1})}.$$

(20)

In the same manner,

$$\Im ({\widehat{Z}}_j)=2{\widehat{S}}_{j-2}-{\widehat{S}}_{j-1}-{\widehat{Q}}_j.$$

(21)

Consider $2^{n_{\mathrm{f}}-j+1}$ noise samples, thus the power

$$\begin{array}{c}\hfill \mathrm{P}\{\Im ({\widehat{Z}}_j)\}=\Im {(Z_j)}^2+2^{n_{\mathrm{f}}-j+1}\frac{{\sigma }^2}{{(2^{n_{\mathrm{f}}-j+1})}^2}=\Im {(Z_j)}^2+\frac{{\sigma }^2}{(2^{n_{\mathrm{f}}-j+1})}.\end{array}$$

(22)

Now, the signal-to-noise ratio (snr) can be expressed as

$$\mathrm{SNR}=2^{n_{\mathrm{f}}-j-1}\frac{\Re {(Z_j)}^2}{{\sigma }^2}.$$

(23)

Because we study bioimpedance, $Z(f)$ can be assimilated to cpe. Using Equation (1), we assume that $\alpha \approx 1/2$:

$$\Re {(Z_j)}^2=\frac{1}{2}{\left|Z\left(\frac{f_0}{2^j}\right)\right|}^2=\frac{2^j}{2q_0\pi f_0}.$$

(24)

Leading to the signal-to-noise ratio

$$\mathrm{SNR}=\frac{2^{n_{\mathrm{f}}-2}}{{\sigma }^2{q}_0\pi f_0}.$$

(25)

We can notice that the expression does not depend on rank j. Obviously, the same conclusion applies for the imaginary part. Therefore, with the previous assumptions, the snr is almost constant along the whole frequency range. This is due to typical bioimpendances decreasing at the same rate as the squared root variance of the proposed estimator (Figure 1). A numerical simulation has been performed to illustrate this behavior (Figure 11).

Finally, in the case where the impedance does not decrease as $1/\sqrt{f}$, a solution may be to adapt the generated current signal by increasing the amplitudes for frequencies where the snr worsens.

8. Discussion

In the presented bioimpedance spectroscopy system, all frequency components are processed in parallel. Thus, mismatch between channels might impact output accuracy. However, with the proposed architecture, all signal paths (real, imaginary, for each frequency) are fed from the same ADC. Considering a particular channel, data are just sub-samples of a unique sample flow. All further processing being digital, the only source of mismatch appears to be the ADC clock jitter. Under the reasonable assumption that jitter is statistically independent from the sub-sampling operation, we hope that this effect will be moderated by the accumulation operation. From this point of view, jitter may be treated as an additional noise (or error) on sample values.

Another concern may be the very-low frequency noise due to amplifier offset, electronics $1/f$ noise, electrode interface potential, etc. It would induce dc component that would impact all $S_j$ and $Q_j$ values. However, it appears relatively easy to get rid of this effect by averaging all samples at $f_{\mathrm{clk}}$, i.e., before decimation and subtracting this dc value during post-processing.

The last issue could be the fixed resolution due to frequency distribution as powers of $\frac{1}{2}$. If more resolution is needed, a non-power-of-two frequency division can be introduced, e.g., $\frac{f_{\mathrm{clk}}}{3}$ and then used instead of $f_{\mathrm{clk}}$ in a duplicated structure. Of course, expressions of $S_j$ and $Q_j$ will be affected and need to be re-evaluated. However, only post-processing computation will be modified.

9. Conclusions

We have proposed a very efficient implementation of multitone analysis for bioimpedance sensing. Thanks to general characteristics of bioimpedance, we have chosen a logarithmic frequency distribution. With this hypothesis, we can explicit the aliasing resulting from successive synchronous demodulations and decimations. This expression can be used to post-process the measurements and estimate real and imaginary parts of the impedance for each frequency. The digital implementation consists only in accumulators and frequency dividers and does not make any use of multipliers. Moreover, the structure can be easily scaled to wider frequency ranges. Despite expected degradation of snr for low frequencies, we show that this effect is mitigated by the increase of impedance in this range.