Time-frequency design of a multi-sine excitation with random phase and controllable amplitude for (bio)impedance measurements

Impedance spectroscopy has become a standard electroanalytical technique to study (bio)electrochemical and physiological systems. From an instrumentation point of view, the measurement of impedance can be carried out either in the frequency domain using the classical frequency sweep method or in the time domain using a variety of broadband signals. While time-domain techniques can be implemented with relatively simple hardware and can achieve faster acquisition time, they are still not that popular because of their lower accuracy and modularity. In this work we present a method and an algorithm for computing the impedance spectrum of samples using an arbitrary time-domain signal excitation initially engineered from a multi-sine summation. The signal is designed in such a way that its spectral phase is derived from a discrete uniform distribution (obtained from a physically-true random number generator), but its time-domain maximum amplitude remains controllable. As such, samples can be excited with adjustable amplitude signals (as low as 50mV to 2V in the current setup) in order to avoid exciting the nonlinearities in the samples under test, which is critical for the very definition of the impedance. Furthermore, the method is proved to combine both speed of measurement, given that it permits to excite at once a large number of frequency components of the sample, as well as accuracy given that the signal’s power spectrum is relatively flat. We verified our method and algorithm for monitoring yogurt quality using (bio)impedance measurements over a period of 48 hours.


I. INTRODUCTION
T HE spectral impedance Z(ω) (or admittance Y (ω) = 1/Z(ω)) is a complex-valued characteristic property of solids, liquids, semi-liquids as well as organic and inorganic materials, from which many mechanistic and physicochemical information can be deducted [1]- [3]. Its measurement is carried out via small-amplitude perturbations around an equilibrium steady-state point (current-voltage profiles in most (bio)electrochemical systems are in general not linear but can be considered linear for small enough perturbations), and is capable of probing relaxation processes over several orders of magnitude of frequencies (e.g. 10 −3 to 10 6 Hz) [1]. In principle, if the impedance function is sampled over an infinite bandwidth, it should be able to provide all necessary information about the sample from an electrical point of view [1], something not possible with transient techniques such as potentiodynamic polarization or dc excitations [4]. In addition, the measured data can be readily checked and validated against the Kramers-Kronig transforms for subsequent analysis. These features have made impedance spectroscopy a well-established technique that is frequently used for the study of many (bio)electrochemical materials, devices and systems [5] . For instance, the capacitiveresistive behavior in electrochemical energy devices [6]- [8], kinetic data of semiconductive materials and solar cells [9]- [11], microstructural properties of mixed phase systems [?], interfacial corrosion reactions and quality control of coatings [12], as well as some biochemical mechanisms in organic systems [13], have all been extensively studied using impedance spectroscopy. In particular, the agricultural industry and (bio)impedance sensing in general have seen a steady rise over the past few years in the use of this technique to assess for instance the ripeness of produces [14]- [18], the overall quality of dairy products [19] and vegetable oils [20], as well as for detecting fungal and bacterial infections [21]. However, in terms of instrumentation and measurement, further research is still needed in order to speed-up the data acquisition time without picking up external noises, which is especially important when monitoring the very low-frequency behavior of fast-evolving systems (< 10 −2 Hz) [22]. Fast measurements are also desired in portable impedance analyzer systems to save power consumption and to allow long-lifetime during continuous monitoring of samples in the field. High-frequency measurements (> 10 6 Hz) present also their own challenges such as stray capacitance and transmission line effects in the leads and the sample [3].
Measurement methodologies of impedance are most conveniently classified based on the type of the applied excitation signal and its independent variable. For frequencydomain techniques, a small-amplitude single-frequency sine voltage (potentiostatic control) is applied on the sample in the form of: v(t) = V 0 sin(ωt) (1) and the current response is measured. The resulting current should be a single-frequency sine function of the form: provided that the sample under test is linear, stable and obeys the causality principle. However, in practical situations the response may show additional harmonics and noise contents that can be represented by [3]: Impedance measurement systems based on single-frequency sine excitations and measurement at the same frequency are able to minimize these effect especially when longer acquisition times are used. Similarly, under the so-called galvanostatic control, the sample is excited with a smallamplitude sinusoidal current input and the voltage response is measured. For most cases, the two types of excitations should provide the same results. The impedance function of the sample is then computed from: for each frequency ω, where V (jω) and I(jω) are the frequency transforms of time-dependent voltage excitation v(t) and current response i(t), respectively, and |Z(ω)| and φ(ω) are the magnitude and phase of impedance, respectively. They are given by: |Z(ω)| = Z (ω) 2 + Z (ω) 2 and φ(ω) = tan −1 (Z (ω)/Z (ω)) where the real and imaginary components are computed from the integrals given by: and respectively. The measurement process is repeated at several frequencies in order to draw the characteristic impedance spectrum of the sample across the frequency range of interest.
On the other hand, the applied signal and the response of the sample can be recorded in the time domain with time being the explicit independent variable in this case. The applied voltage can be any arbitrary function containing multiple frequency components, which can be represented by a Fourier integral: The resulting current passing through the sample would be in this case: This is a two-step approach different from direct frequencyby-frequency measurements of impedance in the sense that at first the input voltage and response current are sampled and recorded in a given time window. Then, the computation of the impedance is carried out using hardware or software timeto-frequency conversion techniques such as Fourier transform and similar signal processing operations. Examples of such excitations are multi-sine, sinc, chirp, and pseudorandom signals which are all broad bandwidth signals having the advantage of shortening the measurement time compared to frequency-domain approaches [3]. Most of these signals are designed in the time domain, therefore, there is minimum controllability on the resulting frequency-domain behaviour. A number of authors have previously introduced multi-sine impedance measurement techniques with a random phase in the time domain (see [23] and the references therein), but without a clear method on controlling the maximum amplitude in the generated time-domain signal.
In a recent paper [22], the authors of this work introduced and verified a fast system for impedance measurement relying on the use of a voltage excitation with constant spectral power but random phase. The signal was constructed from a Gaussian noise signal in the time domain, which in the frequency domain presents random spectral magnitude and phase. Then, calibration and conditioning loops were applied in order to ensure a relatively flat magnitude while maintaining the phase random as-is. This technique, while having the advantage of being faster than classical frequency sweep methods with comparable accuracy, does not allow controlling the voltage level of the signal to be applied on the sample; an issue which also exists in other random phase, multi-sine based techniques [23]. This is a very critical aspect in impedance measurement since nonlinearities in some electrochemical and biological samples may be involved in their response when high enough amplitude signals are applied, which in turn invalidates the measured data. On the hand, some hard, less conductive materials may require a sufficiently larger voltage in order to induce a measurable current for the computation of their impedance spectrum.
The aim of the work is to address this limitation of amplitude controllability by presenting a new algorithm and method for impedance measurement while keeping the main original features of time-domain signals with random phase, i.e. fast measurement time. The initial signal is constructed from a multi-sine summation in the time domain, but its spectral phase (computed with discrete-time Fourier transform (DTFT) algorithm) is modified with a discrete series of uniformly-distributed numbers processed from a stationary, continuous random fluctuations. This is done iteratively in two steps: (i) appending the phase spectrum by random numbers extracted from a True Random Number Generator (TRNG), and (ii) re-constructing the time domain signal and comparing its maximum amplitude with the user-defined value, as explained in detail later below. The impedance system is made compact and cost-effective, and consists simply of an off-the-shelf data acquisition board interfaced with a PC. No additional hardware is required. (Bio)impedance measurements carried out on fruits and vegetable samples are compared to those obtained with a state-of-the-art benchtop impedance spectroscopy station and show that our method is equally precise and accurate.

II. EXCITATION SIGNAL A. SIGNAL DESIGN AND ALGORITHM
The proposed time-domain signal is derived from the realvalued discrete sum of a group of sine waves, initially all of constant amplitude A 0 and zero phase φ 0 , and with frequencies ranging from f i to f f such that: The N frequency points can be chosen equidistant from f i to . ., f i + (N − 1)∆f , or log-spaced and t = 0, 1, 2, . . . , M − 1 are the consecutive discrete time instants which depend on the sampling frequency f s .
Here f s is taken as 8 to 10 times the maximum frequency of the sample interval. The flowchart summarizing the steps followed for the measurement method is depicted in Fig.1.
The user starts by setting the four parameters, f i to f f , N , and the maximum voltage that can be allowed on the sample, V max . The number of segments that each sub-interval of frequency (e.g one decade if log-spacing is selected) will split into requires another constant integer N s that we will introduce later on. Discrete-time Fourier transform (DTFT) using fast Fourier transform (FFT) algorithm is then applied on the initial signal {x i (t)} (i = 0) such that: from which the frequency-domain magnitude and phase, |X i (jω k )| and arg(X i (jω k )), are computed. The phase vector at this point is discarded and replaced by randomlygenerated and independent phase values that we denoted {φ m } in the flowchart.
The random numbers to be appended onto the phase of the signal of this first iteration is produced using a physical random number generator (Random Bit Generator of RBG). The RBG system is based on the highly-fluctuating current time series observed between the tip of two electrodes when atmospheric pressure air microplasma is ignited [24]. The data can be either collected as a digital current using a current probe clipped around one of the electrodes [24], or in the form of binaries directly from the ADC register of the AD2 board. In this work, we collected the series of bit 3 of each stored ADC sample, which were then packed into 32-bit binary numbers. The binary numbers are then converted to decimals and scaled to the phase range between 0 and 360 deg. Fig. 1(b) and 1(c) present a box plot and a polar histogram of sample of 3.072 × 10 6 datapoints processed from the RBG system which show that the values are uniformlydistributed on the interval [0, 360].
The signal of magnitude |X i (jω k )| and arg(X i (jω k )) appended with the random phase values is then reconstructed (back to the time domain) by applying the inverse DTFT. To achieve voltage controllability which is an important contribution of this method and algorithm, the maximum voltage in the newly-constructed signal is checked and compared with the set value of V max . If the maximum voltage of the reconstructed signal is larger than V max , the same process is repeated with another sequence of random numbers. If after trying enough batches (3000 trials were done in this work) without reaching an acceptable maximum voltage, the initial amplitude A is reduced (a step of 0.01 is used in this work) and then the whole process is repeated again. Once V max is reached the loop quits and the target signal {x c (t)} ( Fig. 1(a)) is constructed. This procedure guarantees reaching a random-phased signal with a maximum pre-set value of V max , with a typical execution time in the order of the second (using MATLAB v2020 running on Apple Mac Book Pro 2.2 GHz Quad-Core Intel Core i7).
Again, we emphasize here that the designed signal, which is performed in the frequency domain with time-domain constraints, has the following advantages: 1) The random phase sequence is appended onto the frequency domain spectrum of the signal and not in the time domain as commonly performed with traditional multi-sine techniques [23]. 2) Therefore, the time-domain maximum amplitude of the signal can be user-controlled, which is very important VOLUME 4, 2016 User-defined parameters when characterizing samples such as soft tissues and avoiding deviation from linearity.

B. TIME-FREQUENCY ANALYSIS
The signal {x c (t)} generated to compute the spectral impedance of the samples is non-stationary when considering its entirety from f i to f f (an example for a short sample is shown in Fig. 1(a)) which also means that from the frequency-domain perspective its spectrum keeps changing along with the duration of the measurement. In Fig. 1 we show typical time, frequency, and time-frequency domains representations of the signal that results in one frequency decade from 1 to 10 Hz with N = 40. Similar results were observed and verified in all other decades from f i = 1 Hz to f f = 1 M Hz. The signal has been segmented into four portions of N s = 10 pts each as shown in Fig. 1(b) in order to allocate a relatively higher power for each targeted frequency to run the measurements when compared to having a single segment covering the whole decade from 1 to 10 Hz. This has been shown to be of practical importance when conducting the measurements (see below). In the timedomain plot in Fig. 1(b), one can visualize the four bursts of excitation which are spread over the frequency range from 1 Hz to at least 10 Hz according to its frequencydomain representation in Fig. 1(c). Higher power level in the low frequencies are due to the longer duration of its corresponding time-domain signal. Because the frequencydomain plot does not provide any type of information on the time, and the time-domain plot does not provide any type of information on the frequency, the two representations can be combined in a joint time-frequency spectrogram that shows the concentration of energy at certain time instant or certain frequency band or more generally, in some particular time and frequency region [25]. The spectrum is estimated over sliding windows (within which stationarity is assumed) using short-time Fourier transform (STFT), and the results are depicted in Fig. 1(d). Four frequency ridges in bright yellow color can be visualized indicating energy concentration areas in this time-frequency representation, which are in line with the four burst of energy seen in Fig. 1(b). The ridges are, however, relatively spread over several adjacent frequency bins (see frequency power levels in the color map) due to the leakage of the windowing method used in both time and frequency analysis. The frequency resolution was set to 571 mHz and the time resolution to 4.5 s in this case. It is understood that improving the time and frequency resolutions simultaneously is not possible given that the two variables are inversely proportional to each other.

C. EXPERIMENTAL SETUP
The setup used of (bio)impedance measurements consists of a cost-effective and compact in size Digilent Analog Discovery 2 board equipped with a 14-bit analog-to-digital conversion module (ADC), and capable of a maximum sampling rate of 100 MS/s. The AD2 power supply has a voltage range of ±2.5 V (more than enough for (bio)impedance measurements) and a maximum current driving capability of 50 mA. One channel of the AD2 board was used for both measurement and data acquisition via two gold-plated electrodes.

D. IMPEDANCE MEASUREMENT
The computation of impedance Z w (jω) of a sample is conducted the same way as in our previous work [22].

E. COMPARISON WITH A COMMERICAL IMPEDANCE ANALYZER
Prior to impedance measurements on real samples, data calibration was performed in order to eliminate trend errors (e.g. resulting from cabling). The calibration was carried out on a standard RC circuit: R 0 + (R 1 C 1 ) + (R 2 C 2 ) with R 0 = 1 kΩ, R 1 = 1 kΩ, R 2 = 3.57 kΩ, C 1 = 10 nF , and C 2 = 2.2 µF [22]. Its complex impedance function is given by: The impedance results obtained using the signal {v c (t)} = {x c (t)} are compared to those obtained from frequencysweep method on a Hioki IM-3536 LCR meter (measurement range: 100 mΩ to 100 MΩ with accuracy of ±0.05% on |Z H | and ±0.03 • on φ(Z H ); measurement frequency capability: 4 Hz up to 8 MHz). Fig. 2 summarizes the data calibration results conducted with V max set to 0.5 V over the frequency range 1 Hz to 1 MHz. The error trends in Figs. 2(a) and 2(b) (black scattered dots) are calculated by dividing the as-is measured impedance and phase values of the calibration box by those of  The average of the magnitude relative error is 1.24% with standard deviation of 1.02% compared to the impedance analyzer and 1.00% with standard deviation of 1.02% compared to the ideal magnitude response (Fig. 2(e)), whereas for the phase these values are 3.89% with standard variation of 4.3% and 2.74% with standard deviation of 3.61%, respectively ( Fig. 2(f)). In general, in any impedance measurement system the distortion of the data can be attributed to several factors including leakage error and broadening of the Fourier spectrum, aliasing error because of the discrete nature of the VOLUME 4, 2016 data, and to the discrete-time approximation of the infinite Fourier transformation using a finite number of datapoints [3]. Here, we tried to minimize these effects by applying the DTFT to several segments (N s = 4 for these results) within each single decade of frequency from 1 Hz to 1 MHz, and for which the sampling frequency was chosen to be greater enough than the highest frequency in each segment (8 to 10 times). However, a compromise between speed and quality of measurements should be made by the end-user via the value of N s (see Fig. 1(a)). The larger is the value of N s the better is the quality of the data, but that comes at the expense of longer measurement time, and vice-versa. We also verified the calibration procedure at different voltage levels from 50 mV to 1.8 V using the same polynomial calibration function, and the magnitude and phase results are shown in Fig. 1(g) and Fig. 1(h), respectively. It is clear that the fact all spectra coincide with each other that the desired voltage controllability for this method and algorithm is successfully achieved.

F. APPLICATION
The calibrated setup described in Sect. II-C has been used to characterize by impedance measurement a sample of natural, low fat, unsweetened yogurt purchased fresh from the local market. This is solely shown for the purpose of demonstrating the proposed method in a useful application. The measurements were carried out at room temperature once every three hours for a total duration of 48 hours, shown in Fig. 3(a). The experiment was repeated four times on four different samples of the same yogurt type and brand to check for consistency. Studying the chemical changes in parallel with the impedance measurement is outside the scope of this work. The sample is assumed to remain biochemically stable during the period of each measurement, which is here just 79.9 s for the frequency range from 1 Hz to 1 MHz. Otherwise, as mentioned above, if it was nonlinear or drifts with time (something frequently encountered with (bio)electrochemical and physiological samples), then the Fourier transform will incorrectly ascribe the harmonics due to system nonlinearities to input signal components which may or may not be there. Consequently, the measured impedance may be corrupt and invalid, which can be screened using the Kramers-Kronig integral transforms. Fig. 3(b) shows the Nyquist plot of the conducted experiments over 48 hrs covering frequency range from 1 Hz, to 1 MHz. In Figs., 3(c), and 3(d), we show the (bio)impedance results in terms of magnitude and phase spectra vs. frequency. The shown amplitude and phase are captured within the frequency range between 2 Hz and 100 KHz, at which the amplitude shows an interesting behavior. Note that below 2 Hz the readings started to exhibit unstable and noisy values, on the other hand, there was no significant increase in the amplitude beyond 100 KHz.
The figures indicate that there are changes in the electrical properties of the yogurt over time. In particular, as shown in Figs. 3(e), 3(f), 3(g), and 3(g) in which the impedance  magnitude at 2 Hz, 100 Hz, 1 KHz and 100 KHz, respectively, is plotted vs. time (in log scale), one can notice a linear trend. Linear fitting of the data are shown in all figures in red dashed lines. As shown, the impedance increases at a higher rate in the lower frequencies, the slope is equal to 53 at 2 Hz, and a lower rate at the higher frequencies, the slope is equal to 18 at 100K Hz. The phase behavior is not shown at individual frequencies because it did not show substantial change compared to the magnitude as can be deduced from Fig. 3(d). The four samples used in the experiment manifested similar behaviour.
This experiment indicates that the magnitude response can be used as a potential monitoring tool for yogurt quality. Note that we started observing the local growth of fungus culture around the positive electrode around 30 hours onward for the four samples. Further detailed investigation (including samples from other manufacturers, effect of temperature and other environmental conditions) is required in order to capitalize on the advantages of this (bio)impedance sensing technique and is beyond the scope of this work [26].

III. CONCLUSION
In this work, we presented an engineered excitation voltage produced from a multi-sine signal with random phase data (obtained from a true random number generator) appended to this signal in the frequency domain. The advantage of using this type of excitation signal in impedance measurements, compared to the classical frequency-sweep technique, is mainly the significantly reduced measurement time at low frequencies arising from its wide-band noise-like power spectrum. More importantly, amplitude controllability, which is an important feature in practical applications, can be maintained. The proposed signal was used in an impedance measurement setup and compared with a commercial highend impedance analyzer over the frequency range 4 Hz to 1 MHz on a standard calibration circuit box showing less than 1.24% magnitude relative error and 3.89% phase relative error. Furthermore, we measured the impedance of a sample of natural yogurt over 48 hours using the same setup as a demonstrative application. This is particularly done to show the portability and low-cost of the technique. Further future improvements will focus on optimizing the frequency segmentation by finding the optimal number of frequency points within a segment, while keeping the power level high enough at these points to reduce the error. Also, alignment of the target frequency points with the resolution of the FFT can be further optimized.