Fast Impedance Measurement of Li-Ion Battery Using Discrete Random Binary Excitation and Wavelet Transform

Electrochemical impedance spectroscopy (EIS) is a widely used means for characterization of the dynamics of batteries and electrochemical energy conversion systems in general. EIS is useful for on-line condition monitoring since it contains valuable information on the internal condition of the batteries. The conventional approach to the EIS of batteries relies on their successive excitation with mono-component sinusoidal signals at the pre-defined frequencies. When inferring about the battery’s state-of-health, the low-frequency part of the impedance characteristic is of particular interest. Excitation in the low-frequency region can take an excessively long time. Moreover, maintaining stable experimental conditions over long time intervals i na way that the external disturbances will not affect the estimated impedance, might be demanding, especially in the in-field applications. To alleviate the said limitations, and to minimally intrude the process operation, in this paper we apply broadband electrochemical impedance spectroscopy based on a discrete random binary sequence for perturbation of the battery input. The impedance is evaluated by processing voltage and current signals with continuous Morse wavelet transform. The main contributions of the paper refer to (i) the accurate evaluation of the impedance spectra from $\mu $ Hz to kHz range with high-frequency resolution (more than 200 points per decade) and (ii) the evaluation of the uncertainty region of the impedance characteristic. The entire characterisation takes only a fraction of the time required by the classical sine-based electrochemical impedance spectroscopy. The algorithm is successfully demonstrated on a commercial Li-ion battery, which, together with the all datasets are available for download at https://repo.ijs.si/gnusev/supplementary_material.git.


I. INTRODUCTION
Thanks to their high energy and power density, Li-ion batteries have become indispensable in a broad range of stationary and automotive applications as well as in portable devices [1]. Reliable indication of the change in the internal state of health of all the main battery components (cathode, anode, electrolyte) and remaining useful life could be gained by EIS since change in the patterns of the impedance spectra could be related with the change in the internal condition. In the recent The associate editor coordinating the review of this manuscript and approving it for publication was Bernardo Tellini . paper from Zhang et al. [2], it was demonstrated on a huge set of long-term experiment data that information contained in EIS is very useful in diagnosing the degradation modes as well as assessing the remaining useful life.
EIS is performed by applying small amplitude perturbation signal to the battery input and then measuring its response on the output. Typically EIS measurements use mono-component sinusoidal excitation [3]- [9]. Hence the impedance is calculated only at a limited number of frequency points (usually ten frequencies per decade).
It has been shown that in the sub-millihertz region the Nyquist curve contains valuable information for estimating the actual capacity of the battery, i.e. state of health (SOH), and its state of charge (SOC) [3], [9]. Therefore, performing accurate measurement in the low frequency region is of great interest.
EIS using sinusoidal excitation usually requires 3 to 10 periods [10]- [12]. Consequently, for low frequencies the time required to perform excitation might become notoriously long. For example, the time required to measure EIS at 1 mHz, by following the recommended 3-10 periods, would take 3000 to 10000 seconds. During long excitation period it is non-trivial to avoid drifts and environmental disturbances that might significantly impact the accuracy of the final result. Moreover, in the ultra-low frequency region (<10 mHz), a reasonable resolution would require excessive overall measurement times. On top of that, the information about the uncertainty of the measured impedance spectra is mainly ignored.
To evaluate EIS with high resolution over a dense set of frequency points, while, at the same time, reducing the excitation times to the shortest possible duration, we propose the application of broadband excitation signals [13]- [19]. From the reported cases, the lowest frequency at which EIS was measured was 100 mHz requiring approximately 90 seconds of excitation time [13], [15].
The contribution of this paper to the EIS evaluation for Li-ion batteries is twofold. The first is adoption of discrete random binary sequence (DRBS) for battery excitation. DRBS could be viewed as composited of an infinite number of sinusoids over a continuum of frequencies. Second, the approach enables straightforward estimation of the confidence region of the evaluated EIS.
Small-amplitude DRBS guarantees that linearity assumption, the experiment is based on, is valid. However, instead of exciting the system with consecutive mono-component sinusoids, simultaneous excitation over a range of frequencies is performed. The linearity assumption guarantees that the output is sum of responses obtained at the particular sinusoidal inputs. DRBS was successfully applied to the EIS evaluation of fuel cells [20], [21]. The same approach was used by Li et al. [22] for estimating Warburg-like impedance spectra. All of them used Morlet wavelet as mother wavelet. However, the main disadvantage of the Morlet wavelet is that it looses its analytical properties for low values of the central frequency parameter [23], which makes it ineffective for calculating impedance at low frequencies (below 10 mHz).
In the proposed approach, particular attention is paid to the accurate estimation of impedance at low frequencies.
The original approach from [20], [21] is further improved so that the loss of analytical properties of Morlet wavelet is overcome by using Morse wavelet for impedance calculation.
Validation of the proposed EIS approach was done on a commercial 3.0 Ah 811 NMC-Graphite 18650 cylindrical Li-ion battery type LG18650HG2 (LG Chem). The same battery was then characterised using the DRBS excitation and evaluation of the Nyquist characteristic by processing voltage and current signals with the Morse wavelet transform.
The main goal is to show that by using time-domain signals we get results comparable with those obtained with the established high-end laboratory equipment. We go beyond the state of the art as our approach is capable to estimate the confidence intervals of the measured impedance. The signal processing code is available in the supplementary material.
The organization of the paper is as follows. The main properties of the excitation waveform and the signal processing approach are presented in section III. Validation of the approach by using a simulated equivalent circuit model is presented in section IV. Description of the experiment is given in section V. Finally, section VI presents the experimental validation of the approach.

II. EXCITATION SIGNAL
Conventional EIS in batteries is done by using a sequence of mono-or multi-component sine waves, thus obtaining the EIS values at a discrete set of frequencies. From the system identification perspective, a broad-band noise-like signal can excite the system dynamics over a broad frequency range. Such a signal must fulfil several properties [24], [25]: • it must be stationary, • its bandwidth must include the highest frequency of interest, • the power spectral density must be large enough to guarantee an appropriate signal-to-noise ratio.
A signal that satisfies the above properties and whose generation is fairly simple is the random binary sequence. Its value switches between −a and +a at random time instances. The number of changes over a period of time is Poisson-distributed, where the intensity parameter of the Poisson distribution determines the signal bandwidth.
For the case of discrete time signals, changes can not occur at arbitrary time instances, but only at the discrete time points kλ (k ∈ N 0 ), where λ is minimal time between the two switchings [24, pp. 161-162]. Such a signal is referred to as discrete random binary sequence (DRBS). A time realisation of the DRBS and its frequency spectrum are shown in FIGURE 1.
The main idea of using DRBS for system excitation is to simultaneously excite the system with almost all frequencies of interest. In an ideal scenario this would mean using band-pass limited white noise. DRBS is a signal with close statistical resemblance with white noise. The power spectral density d X (ω) of the DRBS is shown in FIGURE 1b and reads: (1) where λ represents the minimal time between the two switchings, ω represents the angular frequency and a represents the amplitude of the generated signal. Power spectrum has zeros exactly at integer multiples of 1/λ. Useful, near-flat part of the frequency band f B is approximately determined by VOLUME 9, 2021 FIGURE 1. (a) Discrete random binary sequence waveform generated with λ = 3.3 seconds with the effective band f B = 0.1Hz. (b) Power spectral density of the discrete random binary sequence. Numerical plot is calculated using a 300 seconds long realisation of the DRBS signal and Welch method for power spectral density estimation.
the −3 dB crossover frequency f B = 1 3λ [26]. The complete derivation of (1) is given in Appendix A.

III. CONTINUOUS WAVELET TRANSFORM
Continuous wavelet transform (CWT) belongs to the group of time-frequency analysis methods. The simplest example of the time-frequency analysis is the short-time Fourier transform. It segments the measured signal f m (t) by using sliding window g(t) in the time domain. The short-time Fourier transform is defined as: where g(t) is window function. As such, (2) can be regarded as a time localised Fourier transform. Fourier transform is performed at each position of the sliding window g(t), which results in a sequence of spectra, hence the name time-frequency analysis.
A challenge with time-frequency analysis concerns the limits of resolution expressed by the Heisenberg-Gabor inequality [27]. We assume f (t) = f m (t) · g(t) and its Fourier transform F(ω) are well localized functions. The center in time of the function f is defined as The concentration of the function around the center is expressed with the spread t In the same spirit the center µ ω and the spread ω are defined Then the Heisenberg-Gabor inequality reads [28] and γ stands for the area of the Heisenberg box described with the rectangle [ω, Message contained in (5) is illustrated in FIGURE 2 by indicating the bounds of the selected time-frequency resolution. For more detailed definition of the joint time-frequency resolution one may refer to [29,Chapter 2].
The left-hand plot represents the case of the short-time Fourier transform (2) where the time-frequency resolution is kept constant. The shape of the rectangle is determined by the selected window function g(t). The right-hand graph depicts the so-called adaptive time-frequency resolution. In this case, at high frequencies, higher precision in time localisation is achieved at the cost of lower frequency resolution. The opposite is true for the low-frequency region, where high frequency resolution is traded for the low time resolution. Regardless of the partitioning, the area of the ''rectangles'' remains constant as guided by (5).
Fixed length time window in short-time Fourier transform (2) results in an inflexible time-frequency partitioning (plot on the left in FIGURE 2). Better results can be achieved by using more sophisticated functions that define orthogonal space capable of describing the transient signals in a more efficient manner [30]. One such transform is continuous wavelet transform, in which any finite energy signal f (t) can be represented with a set of orthogonal functions ψ(t) referred to as wavelets: where u represents time shift parameter and s is a scaling parameter. The scaling parameter contracts and expands the mother wavelet thus achieving adaptive time-frequency resolution, cf. the right-hand part of FIGURE 2.
There are numerous families of mother wavelets each with their own strengths and weaknesses. Since the goal is to analyse impedance data, only complex mother wavelets can be considered for this task. The most notable mother wavelets for performing CWT are Morlet [31], Bump wavelet, the Morse wavelet [32], [33] and the Lognormal wavelet [34].

A. SELECTION OF MOTHER WAVELET
The most common criterion for selecting the mother wavelet is its time-frequency resolution [29], [30], [34]. The time and frequency resolution of a wavelet function defines the minimal time and frequency difference for which two mono-component sinusoidal signals can be still reliably distinguished [29]. In the sequel we will analyse the time and frequency resolution separately as well as joint time-frequency resolution.
The Morlet wavelet and its Fourier transform are given respectively by where ω 0 is the central frequency.
The Morlet wavelet has optimal joint time-frequency concentration i.e. the Heisenberg box area reaches its lower bound t ω = 2π √ 2 . However, in spite of strong properties, the Morlet wavelet suffers from some weaknesses. In particular, it depends only on one parameter ω 0 . Also, we are also limited with the parameter choices it cannot be considered analytic for ω 0 < 5 s −1 .
The Morse wavelet is designed as a two-parameter family of wavelets aimed to achieve optimal localization in the sense that the eigenvalues of a joint time-frequency localization operator are maximized. The general form of the Morse wavelet in frequency domain is defined as follows [35]: where U (ω) is Heaviside unit step function and K α,β is normalizing factor. However, for more computationally efficient calculation the following form of the Morse wavelet is used in the analysis that follows [29]: where parameter q is related to the central frequency ω 0 as ω 0 = q a (1/a) and e represents the Euler's number.
Due to its complexity, the wavelet transform is performed by a simple multiplication in the frequency domain by using FFT algorithm. The computational complexity is n log(n), where n is the number of samples in the signal. For that reason the time-domain waveform of the Morse wavelet is not required.
In can be noticed that in the cases of Morlet as well as Morse wavelet the parameter ω 0 = 2πf 0 directly affects the wavelet's frequency resolution ω. A comparison of the features of the Morlet and Morse wavelet is presented in FIGURE 3. Three aspects are analysed: a) the frequency resolution depending on the central frequency parameter ω 0 = 2πf 0 , b) time resolution (minimal time τ ) in which two consecutive impulses can be distinguished depending on the frequency resolution, and c) the variation of the joint time-frequency resolution γ with respect to the frequency resolution. The value of γ is related to the minimal resolvable area in which two signals can be reliably distinguished both in time and in frequency. Regarding the first aspect, addressed in FIGURE 3(a), the Morlet wavelet reaches a plateau for the frequencies below ω 0 = 3 s −1 , which is not the case for the Morse wavelet. This effect is mainly due to the fact that for the low values of central frequency ω 0 the Morlet wavelet looses the analytical properties [23]. The time resolution of both wavelets is similar, as shown in FIGURE 3(b). The joint time-frequency resolution γ , shown in FIGURE 3(c), is slightly lower in the case of Morlet wavelet.
For the case of battery EIS, the goal is to achieve sufficiently accurate results at the low frequencies. Since the Morse wavelet preserves the analytical properties even at the ultra-low frequencies [23], all of the subsequent analysis is performed using Morse wavelet as the mother wavelet.

B. IMPEDANCE EVALUATION USING CWT
The impedance values can be obtained by CWT in a pretty much similar way as with the Fourier transform. CWT analysis of the voltage u(t) and current i(t) will result in a set of complex wavelet coefficients: (10) VOLUME 9, 2021 FIGURE 3. Comparison of time-frequency characteristics between the Morlet and the Morse wavelet. The notation is adopted from [34]. (a) the minimal resolvable frequency depending on the central frequency parameter ω 0 , (b) the minimal time τ required for distinguishing two impulses and (c) the joint time-frequency resolution, which is related to the minimal resolvable area in which two signals can be reliably distinguished both in time and in frequency.
where Wx(t, f ) denotes wavelet of a signal x(t) and is function of time t and frequency f . The impedance is the ratio of the wavelet coefficients: It should be noted that the impedance calculated in (11) is defined by time and frequency and provides instantaneous amplitude and phase at every time-frequency point. This means that one can track the evolution of the phase and amplitude of the impedance for each frequency over time.
Conversely, the impedance calculated through FFT provides only time-averaged values of the amplitude and phase for each frequency. Hence, having the time evolution of the complex impedance, it becomes possible to perform more detailed analysis on the impedance values, a task that renders impossible when the impedance is calculated through FFT.

IV. A SIMULATED EXAMPLE
Method described in Section III is validated on a simulated Randles circuit, shown in FIGURE 4. Such a circuit is commonly used to simulate impedance of lithium-ion batteries [36]- [38]. The theoretical impedance of the Randles circuit reads: where R s is the series resistance, R 1 , Q 1 and α 1 are the parameters of the pole, where The goal is to estimate the impedance characteristic of Randles circuit (12) in the frequency band between 600 µHz and 2 kHz using time domain excitation and response signals. Since the DRBS has almost flat power spectrum approximately up to 1 3λ (see FIGURE 1b), one might assume that only one waveform with large enough λ will suffice. However, an almost flat power spectrum can be achieved only with an infinitely long signal. We suggest to split the required bandwidth into decades and for each decade use a DRBS excitation with different λ. In such a way, each frequency band will be excited with sufficient energy allowing for accurate spectral estimation. So, for the required numerical validation, the system was excited with 6 DRBS signals whose bandwidths and duration are listed in Table 1. The simulation of the model (12) was performed in potentiostatic mode. The DRBS excitation was applied to the voltage u(t) and the response was recorded through the output current i(t). Since the simulation was done in time domain, the time response of (12) was simulated using closed-form solutions of linear fractional-order differential equations [39]. The sampling frequency was different for each simulated signal and is given in Table 1. The waveforms of the voltage u(t) and current i(t) are shown in FIGURE 5.
With 6 excitation signals, the overall impedance was obtained by merging results of the CWT from each DRBS excitation and the corresponding response i(t), see FIGURE 6. The estimated characteristic almost ideally matches the theoretical impedance, shown with grey markers.
From the Table 1 it can be seen that approximately 97 minutes of excitation is required to evaluate the whole impedance. Also, it should be noted that by using our approach the impedance spectrum is evaluated with a frequency resolution of 122 frequency points per decade. The length of the lowest DRBS bandwidth dictates the lowest frequency, at which the 46156 VOLUME 9, 2021  Comparison between the theoretical impedance and calculated impedance by applying the method proposed in Section III. In order to calculate the impedance 6 different DRBS signals (Table 1) were used.
impedance can be calculated. In the particular simulation the length of the DRBS with the bandwidth 0.01 Hz corresponds to the 3 periods of the lowest frequency for which the impedance is calculated (600 µHz). If the same number of points had to be evaluated by using the classical single sine-wave excitation signal with only 2 periods, that would overall take around 54 hours of excitation. For comparison, by using the standard 10 points per decade the time needed for impedance evaluation by 10 mono-component sinusoidal signals in the same frequency bandwidth from 600 µHz to 2 kHz, will take around 4.7 hours of excitation time. This example clearly indicates that the method described in Section III accurately determines the low frequency spectra in spite of a shorter excitation time.

V. EXPERIMENTAL VALIDATION
The method for fast EIS measurement was experimentally validated using 811 NMC-Graphite 18650 cylindrical Li-ion battery (LG18650HG2, LG Chem). During measurements the battery was held fixed under constant applied force as shown in the lower part of the FIGURE 7.

A. EXPERIMENTAL SET-UP
Since the application of the method described in Section III requires access to the time domain data, the excitation and measurement equipment include high-precision data acquisition system, electronic load with built-in arbitrary waveform generator and control unit. We used the Keysight N6705C DC power analyser together with the N6781A Source/Measure Unit module, which allows access to the time-domain signals. This device is a multi-functional power supply system that combines the functions of a DC voltage/current source together with built-in data acquisition system. The output of the device also includes arbitrary waveform generator. In order not to surpass the 5 mV voltage change of the battery, the experiment was performed in potentiostatic mode. The requirement of the restricted magnitude of imposed voltage perturbation when determining impedance response of a Li-ion battery is explained and clarified in the Supplementary Appendix D A block scheme of the set-up is shown in FIGURE 7.

B. EXPERIMENTAL PROTOCOL
Prior to impedance measurements the battery cell was preconditioned. At the side of positive terminal and at the lower part of the battery steel housing (negative pole) we have carefully spot-welded two nickel tab strips that served for connecting the sense (U+, U−) ports (FIGURE 7). Battery cell was placed in a specially designed measurement rack that exerts constant mechanical force of 150 N at the positive terminal and bottom of the housing via two massive brass electrodes that are connected with the current (I+, I−) ports. This means we can perform true 4-electrode measurements where all the resistances of the measuring wires (cables) and contact resistances at the two battery terminals are effectively eliminated. We placed the battery together with the rack in a temperature chamber and allowed for 2 hours for the temperature of the battery to equilibrate at 23 • C. We performed five initial galvanostatic ±3 A charge/discharge cycles in voltage range from 2.5 V up to 4.2 V, followed by a half-charge up to 3.705 V and 3 hours of voltage hold VOLUME 9, 2021 whereby the current dropped below 1.5 mA. By applying this pre-conditioning procedure the battery was driven to be in thermal and electrochemical equilibrium at State of Charge 0.5. The impedance of the battery was measured in the frequency region 600 µHz to 1.1 kHz. Experiments were performed in potentiostatic mode by using DRBS excitation signals with amplitude V AC = 5 mV. The amplitude was chosen in a way to guarantee linearity, but still big enough to maintain sufficient signal-to-noise ratio. The EIS was measured using 6 DRBS excitation signals as listed in TABLE 1. The pause between each DRBS sequences was 30 seconds.
As stated in the previous section, the lowest observable frequency in the spectra is determined by the duration of the excitation signal (f min = 3 T ), where T represents the length of the acquired signal. In our case the lower observable frequency is f min = 600 µHz. On the other hand, the highest observable frequency is determined by the sampling frequency of the acquired signal.

VI. RESULTS AND DISCUSSION
The same battery sample was measured using three different approaches: 1. Initial measurements Initially, the battery was measured using standard laboratory equipment BioLogic SP-200 by performing mono-component sinusoidal signal excitation in the frequency range between 1 mHz and 10 kHz with frequency resolution of 10 points per decade. These measurements are only used as a reference.

Single sine measurements
For validation purposes using the equipment described in Section V, the EIS curve was measured using single sine excitation in the frequency interval between 1 mHz and 1 kHz. This might seem as a repetition of the first step. However, the goal was to reconstruct the original results using standard method so that we can check validity of the measurement equipment.

DRBS based EIS
Using the same equipment as in the second step, the battery's EIS was measured in the same frequency band between 600µ Hz and 1.1 kHz using the proposed DRBS and CWT based approach.
The implementation of the approaches is available as supplementary material the details of which are presented in Appendix E.

A. INITIAL REFERENCE MEASUREMENTS
The EIS curve of the initial measurements is shown with black × in FIGURE 9. The EIS was measured after 5 initial charge/discharge cycles and comprises 68 frequency points (10 points per decade) in the frequency region between 1 mHz and 10 kHz. Excitation at each frequency point took 1 period. The whole characterisation process took 96 minutes together with pauses between the frequency points. The impedance characteristics, as expected, contain one semi-circle in the frequency region between 3 Hz and 5 kHz and constant-phase capacitive properties at frequencies below 3 Hz.

B. IMPEDANCE CALCULATION USING MONO-COMPONENT SINUSOIDAL EXCITATION
In order to validate the proposed approach, additional measurements were performed using the classical monocomponent sinusoidal excitation in the same frequency band 1 mHz to 2 kHz with the resolution of 10 points per decade. The length of the sinusoidal excitation signal was 3 periods for the frequencies below 1 Hz and 10 periods for the frequencies above 1 Hz. The overall excitation time, without pauses between the frequency points, was 157.53 minutes. The obtained EIS curve is shown in FIGURE 9, along with the results of the DRBS based approach.

C. IMPEDANCE EVALUATION USING DRBS EXCITATION
With the DRBS excitation the EIS curve was evaluated in the frequency range 600 µHz and 1.1 kHz by using 6 different DRBS signals, as described in TABLE 1. The voltage excitation signal and the battery current response are shown in FIGURE 8. The characterisation was performed at V DC = 3.705 V with amplitude V AC = ± 5 mV. The main rationale for using potentiostatic mode was to fulfill the requirement not to exceed the 5 mV voltage change, when performing characterization on battery. The amplitude of the noise in the voltage signal was estimated V noise = 114 µV and the signal-to-noise ratio SNR = 90 dB. When calculating the impedance values, the parameters of the Morse wavelet in equation (9) were a = 3 and q = 1.224. The obtained EIS curve is shown in FIGURE 9. By using the proposed method, the impedance is calculated at 862 frequency points, which amounts to approximately 100 points per decade. It should be noted that this can be easily increased by up to 3000 points per decade without changing the duration of the excitation. The total length of the excitation signals is 97 minutes. In order to accurately estimate impedance for the frequencies lower than 10 mHz, the length of the DRBS signal has to be at least 3 periods of the lowest frequency. Hence, in our example the length of the longest DRBS signal with 10 mHz bandwidth is 5000 seconds, which corresponds to 3 periods of the lowest observable frequency f = 600 µHz. FIGURE 9 shows the comparison between initial impedance measurement performed with BioLogic SP-200, the impedance measured using DRBS signals and the impedance measured using classical single-sine excitation signals. It can be seen that the impedances obtained by the two approaches are identical.
It should be noted that the inductive part of the impedance exhibits a higher variance. This is due to a hardware's attenuation at higher frequencies. More details about the measured impedance spectra at frequencies higher than 1 kHz are given in Appendix C.
The key result is that by using the method described in section III, it becomes possible to estimate the probability distribution of the impedance values at the  calculated frequency points. This provides information regarding the confidence interval of the measurements. The uncertainty region of the estimated impedance is represented with brighter colours (yellow) in FIGURE 9. The details of the method for estimating impedance distribution are presented in Appendix B.

D. DISCUSSION
With classical sine-based EIS measurements the impedance is calculated at a discrete set of frequencies. Broadband excitation, combined with a sophisticated signal processing, allows the evaluation of the impedance at a high frequency resolution. The theoretical limit of the number of points N , in which the impedance can be calculated using the proposed approach, is T = 1 f s × N, where T is duration of excitation and f s is the sampling frequency. For example, in our case the sampling frequency is f s = 24.45 kHz and the duration of measurement is 5834 seconds. That means with only 5834 seconds of measurement the impedance can be evaluated for the frequency range 600 µHz to 1.5 kHz at almost 200 million points using only single excitation signal.
As shown in FIGURE 9, the uncertainty of the impedance estimation at lower frequencies increases i.e. the spread of the confidence is getting larger for lower frequencies. Instead of analysing only time-averaged impedance values when using sine-based excitation, with the proposed approach it becomes possible to estimate the uncertainty of the measured impedance as explained in Appendix B. It should be pointed out that this uncertainty includes both unavoidable measurement errors, as well as the inherent randomness of the processes occurring during battery's operation. The uncertainty region is a valuable information reflecting the quality of the measurement, which is in field of EIS measurements usually neglected.

VII. CONCLUSION
The paper demonstrates the advantages of using DRBS excitation with Morse wavelet for broadband EIS compared to the conventional sine-based approach. The EIS is estimated at almost the continuum of frequency points without the need of additional excitation cycles. This greatly reduces the excitation time while in the same time improves the quality of the results.
Furthermore, having broad-band excitation as input, it becomes possible to assess the uncertainty of the measured impedance at each frequency point. This is a feature that is unavailable in classical sine-based approaches and has not been exploited yet.
Having a tool capable to get fast and accurate impedance values at ultra-low frequencies it will become possible to address the issues of SOH estimation in more details. The complete approach has computationally efficient implementation. The proposed CWT based impedance can be easily calculated as series of Fourier transforms. All these properties show the practical merit of proposed CWT based EIS.

APPENDIX A AUTOCORRELATION AND POWER SPECTRAL DENSITY OF THE DRBS
The imposed DRBS voltage variation is denoted as u(t) and its realisation reads [40] where U n is the amplitude of the n th pulse, which can be either −a or a, D is the initial displacement (delay) of the initial pulse, λ is the minimal width of the pulse and w(t) is a rectangular pulse defined as: The autocorrelation function of (13) can be written as: FIGURE 10. Schematic of a general measurement system form performing EIS with n u (t ) is the input noise and system uncertainties and n i (t ) is the measurement noise.
Since the amplitudes of different pulses l = n are independent, and are also independent of the displacement D then the above relation becomes: For equiprobable values of −a and a the second term is zero. The first term is an autocorrelation of two rectangular signals. Consequently, it can be easily shown that the autocorrelation of the DRBS with equiprobable states −a and a reads: where τ = t 2 −t 1 . The power spectral density (1) can be easily derived as a Fourier transform of the autocorrelation (17).

APPENDIX B ESTIMATING THE PROBABILITY DISTRIBUTION OF THE MEASURED EIS VALUES
The schematic of the measurement setup is shown in FIGURE 10. The battery under test can be regarded as an ideal system. When performing characterisation we have two sources of disturbances. The first one are the disturbances imposed on the excitation signal denoted as n u (t) ∼ N (0, σ 2 u ). These disturbances include any noise from the excitation circuit as well as any influence of inherent random processes within the battery. The second source is the measurement noise n i (t) ∼ N (0, σ 2 i ). Consequently, the time evolution of the wavelet coefficients (10) at each frequency follow the statistical properties of a zero-mean Gaussian circular complex random variable [41]. Since the wavelet coefficients of the voltage Wu(t, f ) and the current Wi(t, f ) are linked through a linear system, 1 they are correlated. The resulting bivariate complex Gaussian distribution has the following covariance matrix: where ρ = ρ r + jρ i is the complex correlation coefficient such that |ρ| ≤ 1. The resulting distribution can be derived in a closed form as [20], [42] where z r and z i are real and imaginary components of the random variable Z . The location of the mode of f Z (z) depends on the correlation coefficient ρ.

APPENDIX C FREQUENCY DOMAIN ANALYSIS OF THE OBTAINED RESULTS
The Nyquist plot of the measured impedance spectra of the battery shown in FIGURE 9 describes the joint behaviour of the real and imaginary parts. To get a better overview of the accuracy of the proposed approach, the variance of the real and imaginary impedance components over frequency are shown in FIGURE 11 and FIGURE 11, respectively. FIGURE 12, respectively. As already shown in the analysis of the Nyquist plot, the uncertainty of the measured impedance is higher at lower frequencies. There is an additional effect of increased uncertainty at higher frequencies (1kHz). This is a direct result of the frequency characteristics of the data acquisition system. FIGURE 13 shows the comparison between the bandwidths of the theoretical and calculated power spectral density (PSD) of each excitation signal DRBS used to perform EIS. It is clear that the measured PSD for all DRBS excitation signals are identical to the theoretical PSD, except for DRBS 1000 Hz (FIGURE 13f). Above 1 kHz, the effective amplitude of the measured signal is significantly lower. This leads to a reduced SNR and thus to an increase of the measurement uncertainty. This effect can be clearly seen in FIGURE 14, which is a zoomed part of the FIGURE 9 at high frequencies.

APPENDIX D RESTRICTION OF THE MAGNITUDE OF THE INPUT EXCITATION SIGNAL
During the impedance response measurements of Li-ion batteries, we are subject to the constraint of the maximum magnitude of the voltage (and current) perturbation. The limitation is of electrochemical in origin. First, it is important to distinguish between two fundamentally different aspects of impedance measurements. The first is the accuracy of the imposed (excitation) signal and the obtained response signal, which determines the ''quality'' of the obtained impedance data and is directly related to the measurement parameters such as SNR, etc. The second aspect, which is not related to instrumentation and signal processing, is the aspect of electrochemical stability of the tested system (battery) during the measurement. In general, tests on Li-ion battery systems (e.g. the LG18650HG2 used here) have to be considered with great care when performing good quality impedance measurements over a wide frequency range (down to very low frequencies).
Moreover, it should be noted that Li-ion batteries can exhibit significant (long-lasting) relaxation times when approaching equilibrium (for selected SOC). In our particular case of the LG18650HG2 battery, which consists of an 811 Ni rich NMC cathode and a graphite anode, corresponding relaxation times of the two active electrode materials have been observed in the range of (at least) several hours. For example, in a recent study of the 811 NMC cathode at Galvanostatic Intermittent Titration Technique (GITT), it was observed that the typical relaxation times of this material ranged from 2 hours to more than 4 hours (depending on SOC) [43]. Therefore, a reasonable relaxation time must be allowed prior to the impedance measurement. There are two ways to implement this: either by introducing a classical open-circuit period or by applying a suitable selected constant voltage to the battery cell. In this work, we used the second approach because we wanted the battery to operate at the selected SOC = 0.5, which in our particular case of the 811 NMC graphite battery corresponds to 3.705 V. The decision that FIGURE 13. Comparison between the bandwidths of the theoretical and measured power spectral density (PSD) of each excitation DRBS signal used to perform EIS. In the figure a) represents DRBS 0.01 Hz, b) represents DRBS 0.1 Hz, c) represents DRBS 1 Hz, d) represents DRBS 10 Hz, e) represents DRBS 100 Hz and f) represents DRBS 1000 Hz. Above 1 kHz, the effective amplitude of the measured signal is significantly lower leading to a reduced SNR which corresponds to an increased measurement uncertainty.  the SOC should be close to 0.5 was made based on the fact that at this SOC the battery (with all its internal components) is electrochemically and chemically most stable. In other words, for the case of LG18650HG2, the choice of a SOC = 0.5 ensures that the system will not change over long periods of time, allowing a systematic and reliable study of the corresponding impedance. In contrast, Ni-rich NMC active materials are subject to strong side electrochemical processes when approaching the upper limit voltage (e.g., 4.2 V), leading to degradation of the cathode during long-term charge/discharge cycles [44].
With the above considerations in mind, we can move on to the important point of whether to use a (small) applied voltage perturbation or a current-controlled perturbation. Li-ion insertion materials possess exactly defined (relaxed) open-circuit voltage (V OC ) as a function of SOC. The value of V OC is determined by the difference in the chemical potential of lithium in the two host materials (cathode-anode). In other words, when a battery is allowed to relax for a sufficiently long period of time at open-circuit voltage (or, alternatively, for a sufficiently long period at constant voltage so that the corresponding current drops to a very low level), it is considered to be very close to the true thermodynamic equilibrium state. In this case, the corresponding voltage V OC is very close to V eq (SOC, T ) (this notation shows that the equilibrium voltage is a function of both SOC and temperature T). Thus, for accurate electrochemical measurements, both SOC and temperature have to be accurately controlled (temperature control is realized by placing the tested battery in an appropriate temperature chamber or dipping it in a thermostatic bath). FIGURE 15 shows measured voltage profiles of the galvanostatic charge/discharge cycle for the 811 NMC-Li cell (C/20 obtained by GITT sequence) and the LG18650HG2 battery cell, obtained at a current density of C/50 (the C rate corresponds to reaching a full nominal capacity in 50 hours of charging or discharging, i.e., nominally 60 mA). The measured voltage of the battery LG is about 100 mV lower compared to the 811 NMC-Li cell, which is due to the fact that the main operating voltage of the graphite anode is about +100 mV compared to the metallic lithium anode. We can observe that for the selected SOC = 0.5, the voltage hysteresis of the NMC-Li cell is ±9 mV relative to the V eq (3.805 V) and the voltage hysteresis of the LG18650HG2 battery is ±12 mV relative to the V eq (3.705 V). The intrinsic voltage hysteresis of the two compared cells is indeed even lower and could be determined, for example, by ''quasistatic'' measurements at extremely low rates [45]. Let us define the maximum voltage deviation from the equilibrium voltage as: V max = V −V eq , where V is the applied voltage. The measured voltage hysteresis at finite current density sets the upper bound on the limited (maximum) voltage deviation we can afford in impedance measurements. Namely, if we increase or decrease the applied voltage V by more than ±δV max relative to V eq , we push the system out of voltage hysteresis. So in this particular case, if the applied voltage perturbation is greater than about ±9 mV, the system (battery) slips out of linearity. Thus, one of the basic conditions for impedance determination will no longer be satisfied. In voltage-controlled experiments, we must be careful not to exceed the range of voltage hysteresis. In the present work, we decided to limit the maximum voltage deviation to ±5 mV. In this way we have ensured that the linearity of the system is maintained at all tested frequencies.
An exception to the above rule is possible for cases where experiments are of very short duration (mainly high frequency). In particular, if the time/frequency is short/high enough that the storage of energy (charge) in the interior of solid-state particles of the active material has not yet started, then higher voltage perturbations can be used. In the present case of NMC graphite battery, this cut-off frequency is about 0.3 Hz [43].
On the other hand, when performing measurements by using current-driven (so-called galvanostatic) perturbations, we have to be very careful not to exceed the minimum frequency (f min ) at which we would push the system under test out of voltage hysteresis. As described and discussed above for the currently under tested LG18650HG2 battery, the largest input perturbation at SOC = 0.5 should definitely be less than V max = ±12 mV. Also, to ensure that we stay within the linear response of the system, we choose the maximum perturbation smaller than V max of the 811 NMC cathode material (±9 mV) by taking the narrower bound of ±5 mV. We directly obtain the maximum current amplitude (I a,max ) that can be applied to the battery under test by the relation: where |Z (f )| is the magnitude of the impedance at a chosen frequency f. Specifically, in the present work, for the case of LG18650HG2 battery at SOC = 0.5, we found the following values for |Z  In summary, for impedance measurements performed down to very low frequencies, we have defined a constraint on the maximum magnitude of the current perturbation that we can apply without violating the linearity of the system. For the battery studied in this work (LG18650HG2), the current perturbation has to be limited to below 90 mA for measurements below 1 mHz. Since we were also interested in impedance behavior at low (or even ultra-low) frequencies in this work, we decided to use voltage-controlled perturbations to safely avoid problems with the limited current. Also, as mentioned earlier, most measurements were made using the Keysight N6705C DC power analyser. The accuracy of the system output voltage is 0.025% + 600 µV, while the accuracy of the current is 0.04% + 300 µA.

APPENDIX E DESCRIPTION OF THE SUPPLEMENTARY MATERIAL
The supplementary material contains all the necessary data and Python notebook required for the reconstruction of the presented results. The entry point is the notebook Supplementary_material.ipynb. This notebook contains all the necessary guidelines and referenced libraries, which are readily available as open source versions. The corresponding git references are given in the notebook. All scripts and datasets are available for download at https://repo.ijs.si/gnusev/ supplementary_material.git.
Voltage was used as an excitation signal to test the battery. The amplitude of the DRBS signal was U AC = 5 mV The dataset used for calculation of impedance spectra consist of 6 different files, named: − DRBS_0.01hz.mat − DRBS_0.1hz.mat − DRBS_1hz.mat − DRBS_10hz.mat − DRBS_100hz.mat − DRBS_1000hz.mat Each file consists of the time-domain voltage (u) and current signal (i), their sampling frequency (f s ) and minimum (f min ) and maximum (f max ) frequency in which the CWT is calculated. Additionally the measured impedance spectra from BioLogic SP-200 is also provided in order to compare measured impedance spectra. Since 1981, he has been affiliated with the Jožef Stefan Institute as the Deputy Head of the Department of Systems and Control. Since 2013, he has also been a Full Professor with the University of Nova Gorica, Nova Gorica, Slovenia. He is the (co)author of more than 200 journal and conference papers with peer review and has delivered ten invited talks. He has been the coordinator or person responsible for nearly 40 research projects and projects for industry. His research interests include condition monitoring, with particular focus on fault detection, diagnosis, and prognosis; nonlinear dynamic system identification; statistical signal processing; and optimal control. He is also a member of IFAC Technical Committee SAFEPROCESS and the General Assembly and Council of the European Control Association. He received several awards, including DAAD Fellowship in 2001 and the ISA Transactions Best Paper Award for 2009.