Resonance parameter estimation from spectral data: Cramér–Rao lower bound and stable algorithms with application to liquid sensors

A recently introduced method for robust determination of the parameters of strongly damped resonances is evaluated in terms of achievable accuracy. The method extracts and analyzes the locus of the resonant subsystem of noisy recorded complex spectra, such that the interfering influences of many environmental factors are eliminated. Estimator performance is compared to the absolute lower limit determining the Cramér–Rao lower bound (CRLB) for the variance of the estimated parameters. A generic model that is suitable for representation of a large class of sensors is used and analyzed. It is shown that the proposed robust method converges to the CRLB for low measurement noise.


Introduction
The resonance frequency and the quality factor (Q factor) of resonant sensors are the most important intermediate parameters in the determination of desired physical quantities. For example, quartz crystal resonators (QCRs) can be used as microbalances (QCMs) to measure mass deposits on the resonator surface [1]. The increased resonator mass results in a reduced resonance frequency. More generally, dissipative overlayers can also be detected with such a sensor [2,3]. If the surface of the QCM is in contact with a viscous liquid, the mass drag causes a frequency shift and the viscosity of the fluid reduces the initially high quality factor of the resonator. Examples of resonant sensors are limited to liquid sensors in this work, although the illustrations are just as valid for a variety of other resonant sensors like pressure sensors, force sensors, accelerometers and mass flow sensors [4, p 225]. We consider applications where the resonance can be strongly attenuated and thus determination of the resonance parameters is critical. Since viscosity may vary within a large dynamic range, it is important that the parameter estimation algorithms perform reliable for a wide range of Q factors and resonance frequency shifts, as well. Further important aspects and an overview of the current state of the art with respect to fluid sensors can be found in [5,6].
Unfortunately, there is no strict consensus in the literature about how resonance frequency shifts and Q factors are to be Measurement Science and Technology Resonance parameter estimation from spectral data: Cramér-Rao lower bound and stable algorithms with application to liquid sensors T Voglhuber-Brunnmaier 1,2 , A O Niedermayer 2 , R Beigelbeck 1,3 and B Jakoby 2 measured best. In cases where only the amplitude response is known, the −3 dB method is commonly used. For example, the resonance frequency f r and the resonator bandwidth B are estimated at the maximum amplitude and by the frequency range between the points at which the signal has decreased by a factor of 1 / 2 on both sides of the maximum. The Q factor is calculated from = Q f B^/r , where estimated parameters are marked with a hat symbol. Particularly in the presence of larger spurious components and low Q factors, the −3 dB criterion is not meaningful, as the maximum amplitude point and the resonance point differ. Furthermore, drifts of the spurious effects are difficult to separate from changes in the quantity to be measured. Although alternative approaches like Lorentzian curve fit [7] and resonance curve area [8] are more sophisticated, they also do not perform well under these conditions. When complex spectra are available, the additional phase information should be considered in the signal processing. The approaches used in this work follow [9] and are based on various modified methods summarized in the review paper of Petersan and Anlage [7]. We assume that complex spectra, e.g. from quadrature-demodulated signals are available. The locus plot in the vicinity of an undisturbed resonance closely resembles a circle, which can be modeled by a linear second-order lumped element system, such as an RLC in series or in parallel connection. In agreement with common literature [10][11][12] on quartz crystal resonators, this part is termed a motional arm. The signal of the motional arm is extracted from the measurement data and the resonance parameters are determined. Therefore, the results do not depend on spurious phase shifts and signal components.
In general, estimation of the resonance parameters from complex spectra by direct optimization methods, such as the nonlinear least squares method (nLSQ), shows poor performance. This is mainly caused by ill-conditioned equation systems. The stability of the estimation method can be greatly improved if the geometric properties of resonance curves in the complex plane are considered as proposed in [9]. This, and modified approaches, have been successfully applied to resonances of QCMs [13], Lorentz force actuated double diaphragm sensors [14], pressure wave viscometers [15], cantilevers and U-shaped wire sensors [16] and vibrating platelets [17].
The paper is structured as follows: Starting with acquisition of noisy steady state oscillations, the noise on the spectrum in terms of variances, covariances and the covariance matrix is determined in section 2. In section 3, reasonable signal models for linear second-order resonances are shown. The Cramér-Rao lower bound for this function and present noise is derived in section 4. The realizable and stable estimator as proposed in [9] is outlined and the probabilistic properties of the estimator are derived in section 5.

Acquisition of complex spectra
In this section, the transfer of measurement noise to the frequency spectra is described. The spectra are calculated by a Fourier transform of noisy sampled time signals. The basic finding is that the real and imaginary parts of the spectra are approximately Gaussian, linearly independent and of equal variance.
The resonance parameters are determined from complex spectra obtained with an impedance analyzer (e.g. Agilent 4294A) or the YAPIA circuit [18]. Both systems generate harmonic oscillations at discrete and known frequencies f j and analyze the oscillation after sufficient settling time. The commercial impedance analyzer uses analogous demodulation, whereas the YAPIA design samples the electrical current signal at the instants t i and calculates the admittance y m [t i ; f j ] 4 The signal y m [t i ; f j ] is inevitably affected by noise δ[t j ; f j ]. In order to determine the spectral properties A j and φ j from y m [t i ; f j ], the Goertzel algorithm 5 [19, p 633] is applied to each discrete time signal vector y m [t 0 , ..., t L − 1 ; f j ], yielding the admittance spectrum 6 The spectral transform can be represented by A scaling factor of 1/L is used in (2d) to obtain spectral amplitudes independent of the sample number L. The readout electronics records the oscillations of different frequencies f j with integer periods, such that inner products vanish and consequently the amplitude and phase in (1) are related to Y by Additive noise on the data yields additive complex noise on the spectrum, which is not correlated with the parameters (see (2a)-(2d)). The variances and the covariance are (4c) 4 Sampled data are indicated by square brackets. 5 The Goertzel algorithm is a computationally efficient implementation of the discrete Fourier transform of a single spectral line. 6 The underlines denote complex valued quantities.
For a large number of sampling points, the restrictions on the measurement noise δ can be relaxed to independent, identically distributed random variables with finite variance, as a result of the central limit theorem [20, p 44]. The definition of variance and covariance can be given in terms of expected values E{x} (i.e. E{x} is an average of x) by 7 A unified representation of variances and covariances for a vector of complex random variables x is given by the covari- where the superscript H indicates Hermitian transposition. On the main diagonal of C are the variances (real valued) and on the off-diagonal elements are the covariances in complex conjugate pairs.

Sensor models of fluid sensors
Statements about the principal achievable accuracies of estimation algorithms are based on a specific signal model. The models used are reduced representations of the actual systems (partial differential equations in general) and in the form of electrical lumped element circuits. Although, the fundamental sensor typically exhibits many different vibrational modes, these equivalent circuits are valid representations. For each mode, the potential and kinetic energy can be calculated, from which in turn equivalent mass-spring-damper systems for a generalized coordinate can be deduced (see e.g. [21, p 22]). Equivalent electrical circuits consisting of passive elements can be determined for every system, such that the proposed lumped element circuits are well founded. The frequency response of, e.g. a liquid loaded QCM or a piezoelectric tuning fork sensor can be modeled by the lumped element Butterworth-Van Dyke equivalent circuit [22,23], as shown in figure 1 (left). The parameters f r and Q denote the undamped resonance frequency and the quality factor of the motional arm (L m , C m , R m ), respectively. Changes of the mass-density and viscosity of the liquid under test change these parameters. The shunt capacitance C 0 = ϵA/h is basically determined by the electrode area A, permittivity ϵ and resonator thickness h. This parameter can vary strongly if there are changes of the electrical parameters (permittivity and conductivity) of the liquid under test as shown e.g. in [24]. On the right side of figure 1, an equivalent circuit, valid for many Lorentz force actuated and inductively read out sensors (e.g. [16,[25][26][27]), is shown. In many setups, the inductive crosstalk between excitation and readout path is significant and is modeled by the inductor L 0 . The ideal transformer couples the electrical and mechanical systems. The coupling strength is represented by the transformer ratio n and depends on the external magnetic field required for excitation and readout.
Due to wiring, sensor interface and calibration errors a phase shift ϕ can be present. The comparison of the admittance and (trans)impedance expressions given in figure 1 shows that they are structurally similar. This means that for the piezoelectric sensors the admittance spectrum and for inductive sensors the impedance spectrum are processed using the same signal model. Since both models are equivalent, we restrict the analysis to piezoelectric sensors with additive noise ε f ( ). The equation is rearranged into 7 The superscript asterisks denote complex conjugates.
where the spurious components are modeled by a com- I and a phase shift ϕ. One fact that will be exploited is that the locus plot of the motional arm alone resembles a circle with diameter ∼ Y centered at + ∼ Y / 2 0i. The additional part Y l causes a significant distortion of the circular shape at Q factors lower than approximately 1000. In fluid sensing especially, the Q factors are much lower and the impact of the shunt capacitance, as shown in figure 2, or the inductive crosstalk is relevant.

Cramér-Rao lower bound
In this section, the limits of the achievable accuracies of the estimated resonance parameters for the model from the previous section are shown. It is demonstrated that the spurious elements lead to larger errors compared to the spurious-free case, and that the choice of the frequency sampling points of the spectrum also has a certain influence.
Since resonant sensors may be represented by different models, and the parameters of each model react differently to measurement noise, the case is considered where the sensor can be described by a series RLC arm with additional phase shift (e.g. caused by wiring or calibration errors) and broadband spurious elements, which may be due to the sensor as for the QCR and owing to signal conditioning (e.g. signal  . Right: application of the −3 dB rule to Y (red) cannot be used in this case, because the signal does not decay on both sides to −3 dB. In general, this method yields largely wrong results for significant spurious signal components, such as for many liquid sensors. transformers [28]). The Cramér-Rao lower bound (CRLB) defines a lower bound for the variance of the estimated parameters [29, p 27]. This limit can always be found for the considered class of sensors, even though the existence of an estimator achieving this bound is not assured. For the signal model in (7a), additive Gaussian noise ε f [ ] j and multiple real parameters θ ϕ = ∼ f Q Y k k l l [ , , , , , , , ] r R I R I T , the CRLB (C θ ) can be specified by (see [29, p 525 where I(θ) represents the so-called Fisher information and C Y the covariance matrix of the spectral data. This form can only be used for independent and equal variances of the real and imaginary parts of ε, which was verified in section 2. Equation (8b) can be simplified in the present case where the covariance matrix of the data is a diagonal matrix with identical entries independent of the parameters (C Y = σ 2 E, with identity matrix E). The vector and matrix multiplications in (8b) are replaced by summation, yielding The achievable variance depends on the distribution of the samples around the locus plot. The lowest variances are obtained for frequency distributions with uniform arc lengths in (10a) and shown in figure 3.
((4 cot ( / 2)) cot ( / 2)) , In [9], an equidistant frequency distribution centered around the resonance f r is proposed. Using this frequency vector does not yield closed form expressions for the CRLB and generally larger values. The differences between equidistant arc lengths and equidistant frequency distributions are examined in appendix A. It can be argued that for this frequency vector the resonance parameters Q and f r must be known a priori. In order to obtain usable measurement results, it is evident that the resonance frequency and the Q factor must be coarsely known. The results in appendix A show no drastic changes for deviations from the sample frequencies in (10a).
The entries of the Fisher information are determined by (9) and are given by   The matrix is symmetric (i.e. I ij = I ji ) and only the non-zero values of the upper triangular part are shown. The thick printed expressions are approximated with relative errors below 2% at Q = 2 and below 0.02% at Q = 10. To the best of our knowledge, no suitable approximation for I 57 and I 68 is known.
The Fisher-information matrix in (11) can also be used to determine the CRLB for various reduced models. If some param- If there are dependencies between parameters (i.e. correlation between the columns of H ), then the variances on the remaining estimated parameters are decreased. In the following, the CRLB for a single RLC series arm, (i.e. free from spurious elements) is derived and compared to the full signal model in (7a).

CRLB for spurious-free second-order resonance
In this case, only the upper left 3 × 3 submatrix of I(θ) is relevant The expressions in (12) are exact, but only valid for the given frequency distribution in (10a). 8 The inverse of I(θ) yields the covariance matrix for this problem with the single variances The approximation in (13a) is valid for 12Q 2 ≫ 1. These results are assigned to σ σ , f Q 2 2 and σ Y 2 and define the minimum variances achievable in the best spurious-free case. Because the center of the locus plot is at ∼ Y / 2 and not at 0 (see figure 2), a covariance between Q and ∼ Ŷ results. For better interpretation of the results, the standard deviations of the parameters relative to their absolute values are shown These results are intuitive, because a sharper resonance (higher Q) allows a more accurate estimation of the frequency. Furthermore the noise amplitude with respect to the circle diameter ∼ Y is relevant for these parameters. It is therefore convenient to define a characteristic signal-tonoise ratio by σ ∼ Y / .

CRLB of second-order resonance with spurious elements
In a first step the linear frequency parameters l R and l I are neglected and thus the upper left 6 × 6 submatrix of I(θ) is considered. Compact expressions can be obtained in this case The variance on the frequency estimate is approximately six times higher and the variance on the Q factor is twice as large as in the spurious-free case. The estimation of the offset introduces also covariances with all parameters, which are not further discussed because the main interest is in the resonance frequency and quality factor, as they are in close relationship with the relevant fluid parameters. Although the spurious components are of minor interest for sensing applications, they have to be taken into account when present and unknown as they increase the variances of the relevant parameters.
Adding lines/rows 7 and 8 of the matrix influences all parameters and the expressions become very bulky. However, the resulting variances of the resonance parameters, as shown in figure 4, are close to the previous results in (15a)-(15e), such that these expressions provide good approximations.
The CRLB for the parameter estimation problem in (7a) is determined. The calculation of the CRLB in this manner is continued in Appendix B for various submodels and additional spurious components. In the next section, the variances of the stable estimation algorithm are calculated and compared with the lower bound.

Performance of the stable estimator
In this section, a stable estimation procedure, as outlined in [9], is analyzed. It will be shown that the variances of this estimator converge to the CRLB for small noise amplitudes. The robust method is based on the investigation of the locus plots of the complex spectra. In a first step, additive, linearly frequency dependent components are removed using so called off-resonance measurements (ORMs) [9]. The ORMs are acquired at both sides and at a significant distance from the resonance and are used for linear regression. The accuracy of the regression depends on the separation and the number of ORMs. (The effect of higherorder background functions are discussed in appendix C in more detail.) Subsequently, the center of the resonance circle is estimated using the hyperaccurate circle fit [30,31], which is free from essential bias. We used the implementation of the algorithm in MATLAB available from the web page of Nikolai Chernov. 9 The phase of the residual function has a smooth profile and the resonance parameters are derived using a nLSQ method. For convergence, a reasonable initial guess of the parameters is still mandatory, but their range is increased such that simple heuristic methods are sufficient.
The model function in (7a) is rearranged such that the function in brackets is centered at 0 + 0i For angles equally distributed (see (10a)) around the circle center, the minimum variances of estimating the circle center and the diameter are (18a) which are reached by using the hyperaccurate circle fit. The estimate of the circle diameter ∼ Y is therefore optimal (compare with (15c)). After subtracting the circle center estimate, a random bias The observation matrix H φ contains first derivatives and S φ second derivatives. It is observed that S φ can be neglected for the present case, which results in the Gauss-Newton method. Furthermore, the convergence behavior of the iteration is not enhanced notably by S φ . For a poor initial guess, it is appropriate to replace S φ by a scaled identity matrix μE where the scalar value μ is reduced towards zero during the iteration. By choosing a suitable positive μ, the search direction of the iteration is continuously directed to a local minimum of the residuum (i.e. θ φ φ − (^)) n m , resulting in a stable iteration. This approach is also known as the Levenberg-Marquardt modification [32, p 145]. For small second derivatives, the essential bias can be assumed small compared to the parameter noise, and the estimation yields values close to the correct parameter. Therefore, the deviation from the correct parameter can be estimated for small errors (Δ θ ≪ θ) by Phase deviations are caused by errors of the center estimation and by the measurement noise The bias and phase distortion function resemble scaled and phase shifted sinusoidal functions for Consequently, the covariance matrix for E{Δθ} ≈ 0 reduces to The result for f r in (23c) seems to fall below the CRLB in (15a) for Q < 1.8. At such low Q, the approximations used for the derivation of the CRLB are inaccurate. By using accurate numerical calculations, the CRLB is shown to be always lower for any constellation. The interesting variances are approximately and thus close to the CRLB in (15a) and (15b). Although this holds only for small noise amplitudes, the convergence to the CRLB is remarkable, because this implies that the method uses the information provided by the data efficiently. The critical operation, which could reduce the information content of the data, is the calculation of the phase angle of (19) yielding (20a). In general, also the magnitude spectrum must be considered to prevent loss of information. However, in the present case, the circle center was subtracted such that the magnitude spectrum is constant and contains no information but noise. All of the information is therefore passed to the phase and it is possible to achieve an optimum result with the benefit of a largely increased numerical stability compared to direct numerical approaches. Figure 5 (left) shows the standard deviation of the parameters relative to the optimum standard deviation. The empirical standard deviations are calculated based on 100 000 simulated measurements with M = 100 frequency points per locus plot. Furthermore, locus plots at different noise levels are shown. The right graph shows the bias relative to the standard deviation. Although the results show a slight bias at higher noise levels, this bias is practically irrelevant as long as it is much smaller than the standard deviation of the parameters.

Conclusions
We demonstrated that the approach for circle estimation in the locus plot and subsequent data processing is not only numerically stable but also an efficient approach to determine the resonance parameters of the motional arm of resonators. The method was successfully applied to piezoelectric transducers and Lorentz force actuated and inductively read out sensors in [14,16,17]. The approach in [9] uses more advanced model functions with arbitrary and polynomial frequency dependencies. The same approach, as shown in this paper, can be used to assess the performance of these methods. When the discussed method is employed for fluid sensing, an interesting coincidence is observed in that the real and imaginary parts of the acoustic impedance of the fluid layer can be estimated with equal accuracy.