Kramers-Kronig testing
One method for validating immittance spectra involves the use of Kramers-Kronig (KK) transforms. Implementations of the three variants of the linear KK tests (complex, real, and imaginary) described by Boukamp (1995) are included in pyimpspec. These three types of tests have been implemented using least squares fitting. Alternative implementations, which were the default implementation before version 5.0.0, based on matrix inversion are also included. An implementation that uses complex non-linear least squares fitting is also included, but this tends to be significantly slower than any of the other implementations. These tests attempt to fit generally applicable equivalent circuit models (ECM, see the two circuit diagrams below). These ECMs are KK transformable, which means that if they can be fitted to the data with small, random residuals, then the data should also be KK transformable.
The ECM that is used for the impedance representation of the immittance data is shown below and it can contain many parallel RC elements connected in series. The capacitor and inductor connected in series may be necessary for impedance spectra where the imaginary parts do not approach zero at the low- and high-frequency limits, respectively.
For the admittance representation of the immittance data, the following equivalent circuit is used instead. Here, many series RC elements are connected in parallel. Similarly to the circuit shown above, a parallel capacitor and/or a parallel inductor may be needed for some admittance spectra.
A few things to keep in mind about this approach to KK testing:
The fitted circuits have no physical significance and some of the fitted parameters may end up with negative values.
Each parallel/series RC element is replaced with an element where the time constant, \(\tau=RC\), is fixed but either the resistance, \(R\), or the capacitance, \(C\), is still variable for the impedance and admittance representations, respectively.
\(Z(\omega)=\frac{R}{1+j \omega R C}\) becomes \(Z(\omega)=\frac{R}{1+j \omega \tau}\) when operating on the impedance representation.
\(Y(\omega)=\frac{C \omega}{\omega R C - j}\) becomes \(Y(\omega)=\frac{C \omega}{\omega \tau - j}\) when operating on the admittance representation.
Either the complex or the real test should be used. Obtaining good fits with the imaginary test can be challenging even when the immittance spectrum is known to be valid. The error is spread out across the real and imaginary parts when using the complex test, which can make it more difficult to spot issues compared to the real test where the error is concentrated on the imaginary part. However, the real test can be overly sensitive by comparison and one should keep in mind that, e.g., slight increases in the magnitudes of the residuals at the frequency extremes might be resolved by choosing a more appropriate data representation or optimizing the range of time constants.
Many immittance spectra can be validated using the range of time constants within the bounds defined by the inverse of the maximum and the minimum excitation frequencies. However, in some cases it is necessary to extend the range of time constants by some factor \(F_{\rm ext} > 1\) so that \(\tau \in [\frac{1}{F_{\rm ext} \omega_{\rm max}}, \frac{F_{\rm ext}}{\omega_{\rm min}}]\) where \(\omega_{\rm max}\) and \(\omega_{\rm min}\) are the maximum and minimum, respectively, of the measured angular frequencies. The range may also need to be contracted (i.e., \(F_{\rm ext} < 1\)). Pyimpspec includes an implementation for automatically optimizing \(F_{\rm ext}\) and whether or not the suggested \(F_{\rm ext}\) is appropriate can be assessed with the help of a 3D plot of \(\log{\chi^2_{\rm ps}}\) as a function of \(N_\tau\) and \(\log{F_{\rm ext}}\).
In the eaxmple above, the default range of time constants (\(\log{F_{\rm ext}} = 0\)) exhibits a wide range of \(N_\tau\) (\(8 < N_\tau < 45\)) with a gradual decrease of \(\chi^2_{\rm ps}\). An extended range of time constants (\(\log{F_{\rm ext}} = 0.394\), purple markers) is found to be optimal since it achieves a similarly low \(\chi^2_{\rm ps}\) with a lower \(N_\tau\).
An optimum number of parallel/series RC elements (i.e., the number of time constants or \(N_{\tau\rm,opt}\)) should be chosen to avoid over- and underfitting (i.e., fitting to the noise or not fitting to the data, respectively). Pyimpspec implements multiple methods for suggesting \(N_{\tau\rm,opt}\):
Method |
Reference |
|---|---|
1: \(\mu\)-criterion |
|
2: norm of fitted variables |
|
3: norm of curvatures |
|
4: number of sign changes among curvatures |
|
5: mean distance between sign changes among curvatures |
|
6: apex of \(\log{\Sigma_{k=1}^{N_\tau} |\tau_k / R_k|}\) (or \(\log{\Sigma_{k=1}^{N_\tau} |\tau_k / C_k|}\)) versus \(N_\tau\) |
Note
The implementations of methods 1, 3, and 4 include some modifications to make them more robust, but these modifications can be disabled.
The default approach combines the three methods that are based on the curvatures of the immittance spectrum of the fitted ECM in order to:
minimize the number of sign changes of the curvatures (method 4)
minimize the norm of the curvatures (method 3)
maximize the mean distance between sign changes of the curvatures (method 5)
Each method represents a stage that is used to narrow down suitable \(N_\tau\) until one remains. It is also possible to either choose which method(s) to use or to pick a specific number of time constants manually.
Pyimpspec also includes automatic estimation of the lower and upper limits for \(N_{\tau\rm,opt}\) in order to reduce the probability of suggesting an \(N_{\tau\rm,opt}\) that is either too small or too large. The lower limit is estimated using a plot of \(\log{\chi^2_{\rm ps}}\) as a function of \(N_\tau\) while the upper limit is estimated with the help of method 5 (i.e., the mean distances between sign changes of the curvature of the impedance spectra of the fitted ECMs). Either limit, both limits, and/or the difference between the limits can also be specified manually.
How to use
A KK test can be performed by calling the perform_kramers_kronig_test() function, which returns a KramersKronigResult object. This function acts as a wrapper for several other functions that can also be called individually: evaluate_log_F_ext(), suggest_num_RC_limits(), suggest_num_RC(), and suggest_representation().
The evaluate_log_F_ext() function attempts to optimize the range of time constants (i.e., optimize \(F_{\rm ext}\)), but the value of \(F_{\rm ext}\) can also be specified explicitly.
A list of KramersKronigResult can be supplied to the suggest_num_RC_limits() and suggest_num_RC() functions.
The former function will return the estimated lower and upper limits (\(N_{\tau\rm,min}\) and \(N_{\tau\rm,max}\), respectively) of \(N_\tau\) where \(N_{\tau\rm,opt}\) is likely to exist.
The latter function will return a tuple containing the suggested KramersKronigResult instance (i.e., the one that corresponds to \(N_{\tau\rm,opt}\)), a dictionary that maps the numbers of time constants to the scores that were used to suggest \(N_{\tau\rm,opt}\), and the estimated \(N_{\tau\rm,min}\) and \(N_{\tau\rm,max}\).
A list of these tuples, where each tuple corresponds to a KK test that was performed on either the impedance or the admittance representation, can then be provided to the suggest_representation() function.
If perform_kramers_kronig_test() is called with admittance=None, then both the impedance and the admittance representation are tested.
Otherwise, only either the impedance (admittance=False) or the admittance (admittance=True) is tested.
>>> from pyimpspec import (
... DataSet,
... KramersKronigResult,
... generate_mock_data,
... perform_kramers_kronig_test,
... )
>>> from pyimpspec.analysis.kramers_kronig import (
... evaluate_log_F_ext,
... suggest_num_RC,
... suggest_representation,
... )
>>> from typing import Dict, List, Tuple
>>>
>>>
>>> data: DataSet = generate_mock_data("CIRCUIT_1", noise=5e-2, seed=42)[0]
>>>
>>> test: KramersKronigResult # The suggested result
>>> test = perform_kramers_kronig_test(data)
>>> # The line above is equivalent to the lines below
>>> # in terms of the work that is performed
>>>
>>> Z_evaluations: List[Tuple[float, List[KramersKronigResult], float]]
>>> Z_evaluations = evaluate_log_F_ext(data, admittance=False)
>>>
>>> Z_suggested_F_ext: float
>>> Z_tests: List[KramersKronigResult]
>>> Z_minimized_statistic: float
>>> Z_suggested_F_ext, Z_tests, Z_minimized_statistic = Z_evaluations[0]
>>>
>>> Z_suggestion: Tuple[KramersKronigResult, Dict[int, float], int, int]
>>> Z_suggestion = suggest_num_RC(Z_tests)
>>>
>>> Y_evaluations: List[Tuple[float, List[KramersKronigResult], float]]
>>> Y_evaluations = evaluate_log_F_ext(data, admittance=True)
>>>
>>> Y_tests: List[KramersKronigResult] = Y_evaluations[0][1]
>>>
>>> Y_suggestion: Tuple[KramersKronigResult, Dict[int, float], int, int]
>>> Y_suggestion = suggest_num_RC(Y_tests)
>>>
>>> suggestion: Tuple[KramersKronigResult, Dict[int, float], int, int]
>>> suggestion = suggest_representation([Z_suggestion, Y_suggestion])
>>>
>>> scores: Dict[int, float] # Scores for various numbers of RC elements
>>> lower_limit: int
>>> upper_limit: int
>>> test, scores, lower_limit, upper_limit = suggestion
A single KramersKronigResult can be plotted on its own, but it is also possible to plot the suggested KramersKronigResult along with the \(\chi^2_{\rm ps}\) values of all KramersKronigResult instances so that one can see if the suggested KramersKronigResult is indeed the best choice.
From the top-left plot one can see that the estimated lower and upper limits define a range of \(N_\tau\) values (filled circles, the y-axis on the left-hand side) where \(N_{\tau\rm,opt}\) value is estimated to exist. The y-axis on the right-hand side shows the scores assigned based on an approach that makes use of methods 3, 4, and 5. These scores are then used to suggest \(N_{\tau\rm,opt}\) (dashed line).
The perform_kramers_kronig_test() function takes keyword arguments that can be passed on to the suggest_num_RC() function.
This can be used to, e.g., select which method(s) to use or to adjust any method-specific settings such as the \(\mu\)-criterion of method 1.
>>> from pyimpspec import (
... DataSet,
... KramersKronigResult,
... generate_mock_data,
... perform_kramers_kronig_test,
... )
>>> from pyimpspec.analysis.kramers_kronig import (
... evaluate_log_F_ext,
... suggest_num_RC,
... )
>>> from typing import List, Tuple
>>>
>>> mu_criterion: float = 0.85
>>>
>>> data: DataSet = generate_mock_data("CIRCUIT_1", noise=5e-2, seed=42)[0]
>>>
>>> test: KramersKronigResult
>>> test = perform_kramers_kronig_test(data, mu_criterion=mu_criterion)
>>> # The above is equivalent to the following lines
>>> # in terms of the work that is performed
>>>
>>> evaluations: List[Tuple[float, List[KramersKronigResult], float]]
>>> evaluations = evaluate_log_F_ext(data)
>>>
>>> optimum_log_Fext: Tuple[float, List[KramersKronigResult], float]
>>> optimum_log_Fext = evaluations[0]
>>>
>>> tests: List[KramersKronigResult] = optimum_log_Fext[1]
>>> suggestion: Tuple[KramersKronigResult, Dict[int, float], int, int] = suggest_num_RC(
... tests,
... mu_criterion=mu_criterion,
... )
>>>
>>> scores: Dict[int, float]
>>> lower_limit: int
>>> upper_limit: int
>>> test, scores, lower_limit, upper_limit = suggestion
The plot of relative residuals is typically used to interpret the validity of the immittance spectrum that was tested. Alternatively, statistical tests performed on the residuals can also be used.
>>> from pyimpspec import (
... DataSet,
... KramersKronigResult,
... generate_mock_data,
... perform_kramers_kronig_test,
... )
>>> from pyimpspec.analysis.kramers_kronig import (
... evaluate_log_F_ext,
... suggest_num_RC,
... suggest_representation,
... )
>>> from typing import List, Tuple
>>>
>>>
>>> data: DataSet = generate_mock_data("CIRCUIT_1", noise=5e-2, seed=42)[0]
>>>
>>> test: KramersKronigResult # The suggested result
>>> test = perform_kramers_kronig_test(data)
>>> statistics: str = test.to_statistics_dataframe().to_markdown(index=False)
The contents of statistics would look something like:
| Label | Value |
|:----------------------------------------------------|--------------:|
| Log pseudo chi-squared | -4.46966 |
| Number of RC elements | 13 |
| Log Fext (extension factor for time constant range) | -0.100386 |
| Series resistance (ohm) | 103.525 |
| Series capacitance (F) | 0.0120676 |
| Series inductance (H) | -2.24707e-06 |
| Mean of residuals, real (% of |Z|) | -5.25249e-06 |
| Mean of residuals, imag. (% of |Z|) | 0.00220288 |
| SD of residuals, real (% of |Z|) | 0.0523757 |
| SD of residuals, imag. (% of |Z|) | 0.075694 |
| Residuals within 1 SD, real (%) | 68.2927 |
| Residuals within 1 SD, imag. (%) | 58.5366 |
| Residuals within 2 SD, real (%) | 97.561 |
| Residuals within 2 SD, imag. (%) | 95.122 |
| Residuals within 3 SD, real (%) | 100 |
| Residuals within 3 SD, imag. (%) | 100 |
| Lilliefors test p-value, real | 0.252269 |
| Lilliefors test p-value, imag. | 0.513698 |
| Shapiro-Wilk test p-value, real | 0.591578 |
| Shapiro-Wilk test p-value, imag. | 0.168292 |
| Estimated SD of Gaussian noise (% of |Z|) | 0.0643079 |
| One-sample Kolmogorov-Smirnov test p-value, real | 0.871214 |
| One-sample Kolmogorov-Smirnov test p-value, imag. | 0.60763 |
All three statistical tests (Lilliefors, Shapiro-Wilk, and Kolmogorov-Smirnov) return \(p\)-values greater than 0.05 (our chosen threshold) for the residuals of both the real and the imaginary parts. The means of the residuals are close to zero as well. All of this indicates that the tested immittance spectrum is likely to be valid. This is also in agreement with the interpretation based on inspecting the plot of the relative residuals.
Some immittance spectra might not be possible to validate based on testing the impedance representation. For example, the Nyquist plot below shows a synthetic impedance spectrum that includes a negative differential resistance (the larger, outer loop that goes from the right-hand side to the left-hand side as the frequency is decreased). Similar impedance spectra have been reported when measuring, e.g., in the passive region of a system with a tantalum working electrode in hydrofluoric acid (Fig. 3b in the reference).
Attempting to perform KK tests on this impedance data as shown in the previous example incorrectly indicates that the spectrum is not linear, causal, and stable.
However, there are two approaches that can be used to successfully validate this impedance spectrum.
The first approach is to perform the KK tests on the admittance data either explicitly (i.e., by specifying admittance=True when calling the perform_kramers_kronig_test() and evaluate_log_F_ext() functions) or by calling perform_kramers_kronig_test() with admittance=None (default value). The latter should then test both the impedance and the admittance representation before ultimately suggesting the result for the most appropriate represention, which in this case is the admittance representation.
The second approach is to add a parallel resistance of a suitable magnitude to the impedance data and to perform the KK tests on the resulting impedance data.
The resistance, \(R_{\rm par}\), is known a priori to be KK transformable. Adding the resistance in parallel to the experimental data, which is represented in this circuit diagram as \(Z_{\rm data}\), does not negatively affect the compliance of the resulting circuit. Thus, the KK compliance of the resulting circuit is dependent on whether or not \(Z_{\rm data}\) is KK compliant.
Note
The magnitude of the resistance to choose depends on the original impedance data. In this example, the real part of the impedance at the lowest frequency in the original data is approximately \(-100\) \(\Omega\). A value of 50 \(\Omega\) was chosen for the parallel resistance after testing a few different values.
As can be seen from the results below, the new, and thus also the original, impedance data has been validated successfully.
References: