Accurate QTF Sensing Approach by Means of Narrow Band Spectral Estimation

We propose a new approach for the extraction of the equivalent parameters of quartz tuning forks used as sensors by means of noise measurements. Noise is used as the test signal for the determination, by means of spectral analysis, of the frequency response of a circuit including the quartz tuning fork whose parameters need to be determined. A new approach for the analysis of strongly peaked noise spectra was developed in order to allow the correct measurement of the strongly peaked noise spectrum at the output of the system, which is the result of the high-quality factor of any quartz tuning fork-based sensor. With the approach we propose, the best compromise in terms of accuracy and measurement time can be obtained in a single measurement run. The performances of the approach we propose are discussed in comparison with those that can be obtained from a swept spectrum approach in the same operating conditions.


Introduction
AT-cut crystals operating in the MHz range dominate the field of quartz crystal microbalance (QCM) applications, but quartz tuning fork-(QTF-) based sensors have also been demonstrated [1,2]. Notwithstanding that the geometrical shape of QTFs makes their functionalization more challenging with respect to AT-cut crystals, the fact that they operate at much lower frequencies can be advantageous in terms of hardware requirements. While in some cases only a few relevant parameters, such as the resonance frequency f S and the quality factor Q, are taken into consideration for sensing applications, full information on the QFT response may require the extraction of the values of all parameters in the equivalent circuit in Figure 1.
The extraction of all parameters in Figure 1 can be obtained starting from the measurement of the response of circuits such as the one in Figure 2 where assuming the transconductance gain A R to be constant, the voltage gain V O /V I is proportional to the QTF admittance Y Q : However, if we require high accuracy, with the circuit in Figure 2, we face two major challenges. In the first place, deviation of the amplifier response from ideality should be minimized or at least recognized and corrected; in the second place, the frequency response measurement must be carried out ensuring that no artifacts are introduced as a result of the strongly peaked response due to the high-quality factor Q of the QTF.
To address the first issue, we recently proposed an approach, based on a variant of the circuit in Figure 2, that allows to easily detect and correct nonidealities in the frequency response of the amplifier without requiring any adjustment on the circuit [3]. As for the second challenge, there are typically two approaches for the determination of the frequency response of a system such as the one we are concerned with: swept spectrum measurements and inputoutput cross correlation. In the case of swept spectrum measurements, V I in Figure 2 is a sinusoidal signal whose instantaneous frequency is linearly changed with time. Under proper conditions (sufficiently small frequency change rate), the amplitude and phase response of the system can be extracted from the analysis of the sinusoidal signals at the input and at the output. In the case of cross correlationbased measurement, the input signal is a wide band noise and the frequency response is obtained by evaluating the cross spectrum between the input and the output signals. In the case of swept spectrum measurements, errors in the determination of the frequency response may arise from excessive frequency change rates while in the case of cross correlation measurements, window leakage may be responsible for systematic errors [4]. For a common wristwatch QTF exposed to air, the quality factor Q is in the order of 10 4 , resulting in a peaked frequency response, when introduced in the circuit in Figure 2, centered at about 32.7 kHz with a half power bandwidth B of a few Hz. According to the criteria used in swept spectrum measurements [4], systematic errors arise unless the sweep rate is less than 1 Hz/s, leading to long measurement times for exploring a significant bandwidth across the resonance frequency. In the same way, when FFT-based cross correlation measurement approaches are used, the resolution bandwidth Δf should be much smaller than the bandwidth of the peaked response in order to avoid significant spectral leakage [4,5]. Since the length of a time record for conventional spectral estimation using FFT is 1/Δ f , in order to improve accuracy, the measurement time must be increased. In order to better appreciate the type of systematic errors that are the result of an incorrect selection of measurement parameters, we report in Figures 3 and 4 the result of simulated measurements with the circuit in Figure 2 in the case of swept spectrum and cross correlation approaches, respectively. In order to simplify the discussion, we have assumed R P → ∞ and C P = 0. We have also assumed thatA R = R S , so that at the resonance frequency we have, from Equation (1), The results in Figure 3 refer to the case of Q = 104 (virgin QTF exposed to air) with frequency change rates (FCR) for the frequency sweeps ranging from 80 Hz/s down to 1 Hz/s. As it can be clearly verified, the larger the FCR, the larger is the error in the maximum amplitude of the response and in the frequency at which it occurs. If, in actual measurements, we assumed the frequency of the maximum of the response (f MAX ) as an estimate of the QTF resonance frequency f S , we would obtain an estimation error f ERR = f S − f MAX which would depend on FCR (inset in Figure 3). When the QTF is functionalized, Q generally decreases, and the same FCR would result in a smaller error. Conversely, for the same tolerable error, faster sweep rates (larger FCRs) can be used for lower values of Q, thus potentially reducing the time required for the measurement. In actual measurements, however, the problem lies in the fact that, especially in sensing applications, it is not easy to know, beforehand, the measurement conditions that allow to obtain accurate results (within a given maximum error) in the shortest possible time. Of course, reducing the measurement time without affecting the accuracy is particularly important in sensing application as this reflects on the ability of the sensor system to follow rapid changes in the environment. In the case of swept spectrum measurements, comparing the results of repeated measurements obtained with different sweep rates (i.e., different values of FCRs) can provide a way to assess the presence of systematic errors. However, repeated measurements take time and therefore, if possible, they should be avoided. Figure 4 refers to the case of the (simulated) determination of the frequency response in Figure 2 by means of the cross spectrum approach when using a conventional FFT spectrum analyzer. As mentioned before, the relevant parameter, in this case, is the resolution frequency Δf in comparison with the bandwidth B = f S /Q. As it can be noticed in Figure 4, the estimated response is characterized by a significant error unless Δf < <B. For a given sampling frequency f C , the frequency resolution is determined by the number of samples N used for the estimation of the FFT on each single signal record (Δf = f C /N), with 1/Δf being the duration of a record. Clearly, a smaller Δf corresponds to a longer record length, resulting therefore in an overall longer measurement time. As can be deduced from Figure 4, in the case Q = 104, with f C = 200 kHz, we need N in excess of 2 20 (≈10 6 ) in order to ensure minimum error, with a corresponding record duration in excess of 5 s. Moreover, in order to reduce the statistical error, the PSD estimation must be averaged over several records, leading to an overall measurement time in the order of minutes. Similar to the case of the swept spectrum approach, if the value of Q is reduced, Δf can be made larger and the measurement time is proportionally reduced. As before, however, in actual measurements, we do not know beforehand the actual value of Q and therefore it is not easy to set the correct value of Δf for optimizing accuracy and measurement time. In actual measurement sessions, we could proceed by starting with a relatively large Δf (short measurement Figure 1: Equivalent circuit of a quartz tuning fork. Figure 2: Simplified schematic of the circuit that can be used for obtaining a frequency response proportional to the QTF admittance; TIA is an ideal transimpedance amplifier with a constant gain A R . time) and reducing it in successive measurement runs until no significant change is observed in the measured curve, a situation that would indicate that Δf has reached a sufficiently small value for not introducing systematic errors. The important difference with respect to the case of the swept spectrum approach is that while in that case there is no way to perform measurements with different FCRs at the same time, it is possible, in principle, to set up a measurement configuration in which the cross spectrum approach is applied with a number of different values of Δf at the same time, i.e., in a single measurement run. Indeed, one could send the signals V I and V O in Figure 2 to a set of FFT spec-trum analyzers (SAs), all working in parallel but with a different Δf setting from one another. In this way, on the SA with the largest Δf , we would obtain, in a very short time, an estimate of the response of the system; with a delay, due to the longer time record required, the SA with the next smaller Δ f would provide its estimate: if it coincided with the first one, there would be no need to wait any further, otherwise, we wait for the estimate of the SA with the next smaller Δf and so on. Clearly, notwithstanding this potential advantage, resorting to several actual conventional SAs operating in parallel would be unpractical and expensive.
In the case of measurements on QTFs, however, all the information about the device response can be extracted from measurements in a quite limited bandwidth across the resonance frequency. With this observation in mind, we have been able to develop an approach for the correct estimation of the response of a system like the one in Figure 2 (where we expect a strongly peaked response in a limited bandwidth) that is functionally equivalent to the one we have described above (several SAs in parallel) but that can be implemented efficiently and with quite limited hardware resources. The theory underlying such an approach and the measurements demonstrating its effectiveness will be discussed in the following sections.
1.1. Proposed Approach. Conventional spectral estimation can be performed by resorting to the modified periodogram method described by Welch in 1967 [6]. The analog signal whose PSD needs to be estimated is low-pass filtered in order to allow correct sampling at frequency f C (sampling period Δt = 1/f C ). The sequence of samples x i = xðiΔtÞ obtained from the filtered input xðtÞ is divided into records of length N, and each record is multiplied by a proper window function w j ðj = 0, ⋯:,N − 1Þ obtaining a new set of records y h,l : where we have assumed that the first element of each record is x hN (h = 0,1,…) so that the records are adjacent and nonoverlapping. The Discrete Fourier Transform (DFT) of each record is calculated obtaining with jz h,k j 2 representing an estimate, for k < N/2, of the power of the input signal after being filtered by a bandpass filter centred at f k = kf CK /N whose equivalent noise bandwidth (ENB) and detailed frequency response depend on f C , N, and on the window function w. If the window function satisfies the normalization condition and if we assume that the PSD of the input signal can be considered approximately constant within the bandwidth of the  Figure 2 by means of the cross spectrum approach (with C P = 0, R P → ∞). The ideal response is represented by the thick black curve.

Journal of Sensors
filter (white noise approximation), the quantity jz h,k j 2 /ENB can be regarded as an estimate of the PSD of the input signal at f k as obtained from the elaboration of the h th record. The estimate from a single record, however, is quite a crude one, as the standard deviation of the error in the estimate can be as large as the PSD to be estimated. It is for this reason that the results from M records need to be averaged as this reduces the statistical error by a factor √M (assuming that the estimates obtained from nonoverlapping records are uncorrelated). It must be noted that while the ENB changes depending on the shape of the window, for the most common window types it ranges from f C /N up to 1:5 f C /N [7] so that, in a first approximation, it essentially coincides with the frequency resolution Δf = f C /N. The success of the approach outlined above mostly depends on the fact that if N is chosen as a power of 2, the DFT in Equation (3) can be calculated by resorting to the very efficient fast Fourier transform (FFT) algorithm [8]. As we have noted before, in the case of the measurements on QTF, we should be working with Δf < <f s /Q that may correspond to record lengths in the order of 10 6 for the largest expected values of Q. On the other hand, if we are interested in the extraction of the QTF parameters from the measurement of its admittance, we do not need to explore the entire frequency range from 0 to f C /2, but we can limit our focus to a few bandwidths B (B = f S /Q) across the resonance frequency. The chirp Z-transform (CZT) algorithm initially developed by Rabiner et al. [9] can be set up for evaluating the DFT of the input sequence of N samples at N equally spaced frequency values with an arbitrary frequency start and arbitrary frequency step at the cost of 3 FFT computations on sequences of length 2N (only 2 FFT calculations are required in the case of repeated estimations). Note, however, that the ENB still depends on N and therefore, the chirp-z transform cannot, in itself, provide any advantage for the problem at hand. To work around this issue, we have devised an approach for obtaining an arbitrarily small ENB by a proper elaboration of the DFTs calculated on adjacent input records regardless of the record length N. When this approach is employed together with the CZT, we can efficiently elaborate records of N samples at a time obtaining the desired frequency resolution close to f S and combine the results of the elaboration on the DFTs on adjacent records in order to obtain any (and more than one at the same time) ENB we may require. In order to better understand the approach we propose, it can be useful to refer to Figure 5.
The block diagram in Figure 5 can be regarded as the idealized continuous time counterpart of the process that in the modified periodogram method leads to the estimation of the power of the process when filtered by a bandpass filter centered at f k . The output of the filter is sampled at intervals NΔt, and the squared moduli of the samples are averaged for obtaining an estimation of the PSD of the process at f k . In the conventional modified periodogram approach, the filter is a finite response filter (window) with an impulse response lasting exactly NΔt. This property ensures, in the numerical domain, that only the last N signal samples are required for each new estimation of the quantities z h,k . At the same time, as noted before, the length of the impulse response sets a limit to the selectivity of the filter. Suppose now that we maintain the configuration in Figure 5, still sampling the output of the filter at intervals NΔt, but no longer enforcing a limit for the impulse response duration of the filter and, hence, for the selectivity that can be obtained. While in principle we can use any impulse response, for the computation algorithm to remain efficient, we must maintain the condition, if possible, that each new power estimate (sampled value at the output of the filter) can be obtained by maintaining in memory only a finite (an possibly small) number of sampled input values. A filter with the impulse response wðtÞ sketched in Figure 6 does indeed satisfy this constraint since, as we shall presently demonstrate, the power estimate at the end of each interval of duration NΔt can be obtained as a function of the power estimate in the preceding interval and of the last N sampled input values. In the discrete time domain, we replace the impulse response wðtÞ with its sampled counterpart w i = αw ðiΔtÞ. Note that the parameter α has the physical dimension of a time, and it is required in order to obtain a dimensionless w i . Note that the values of α and of the parameter A defining the amplitude of wðtÞ in Figure 6 must be chosen so that their product allows to satisfy the normalization condition with (N → ∞) in Equation (4). In order to simplify the discussion, assume thatx i = 0fori < 0. In the discrete time domain, the output z h , k of the filter in Figure 6 at i = ðh + 1ÞNΔt is

Journal of Sensors
With the index change m = l − hN in the last sum and exploiting the following property for w i , we have Note that the last sum over m, except for a normalization factor, is nothing but the DFT of the sequence of the last N samples when a uniform window is used. Equation (7), therefore, offers a quite efficient approach for evaluating z h,k recursively, starting from the evaluation of the DFT over the last N samples. What makes the function w i especially useful is the fact that the values of p (0 ≤ p < 1) is set the bandwidth of the filter and hence the ENB for spectral estimation.
Using quite standard mathematical techniques, it can be easily demonstrated that, in the continuous time domain, the frequency response Wð f Þ corresponding to the impulse response w(t) is with an ENB: From Equation (7), we notice that the same DFT computed over the last N samples can be used for updating the values of as many z h,k as desired (each one corresponding to a different p and, hence, to a different ENB) with a very small computational overhead. In short, thanks to the combination of the chirp z-transform and the approach we have just described, we are now in the condition to implement the multiple spectrum analyser approach discussed at the end of the previous section with acceptable cost in terms of memory allocation and computational complexity. We can proceed as follows: (1) We start by setting the number of frequency points that are required (N) (2) We set up the parameters of the chirp z-transform algorithm in such a way as to obtain the spectra estimation at N closely spaced frequencies across the (expected) series resonance frequency (3) We set up a number of values for p leading to geometrically decreasing values of ENB (typically with a common ratio of 2) (4) Each time a new input record (length N) is acquired, we perform a chirp-z transform on such a record in order to obtain the rightmost sum in Equation (7). Note that this transform needs to be calculated only once, regardless of the number of selected values for p About the 6 th step, it must be remarked the fact that while the z h,k corresponding to different p are all updated at time intervals equal to the duration of the input record, subsequent values of z h,k for a given p are strongly correlated over time distances shorter that the inverse of the corresponding ENB and therefore, in power spectra estimation or in cross spectra estimation, we obtain a reduction of the statistical uncertainty (according to the inverse of the number of averages law) only if the corresponding z h,k are sampled at time intervals at least equal to the inverse of the ENB.
In order to have an estimate of the resources required to perform the elaboration we propose, let us assume N = 8192. With f CK = 100 kHz, the ENB corresponding to p = 0 is about 12 Hz. AssumingQ MAX = 10000, we should be able to obtain a minimum ENB significantly smaller thanf s /Q MAX ≈ 3 Hz. Let us assume that we require to reach a minimum ENB of 0.1 Hz. It follows that 7 different values of p (including p = 0) are sufficient to cover the entire range of ENB values form 12 Hz down to less than 0.1 Hz in geometric progression with a common ratio of 2. Since we operate with two channels and we want to estimate the power spectra of each channel and the cross spectra for all values of p, it follows that, assuming 8 bytes for the storage of a double precision value (one complex value is equivalent, in terms of memory occupation, to two double precision values), the total memory occupation is less than 10 Mbytes, that is almost negligible for modern standards. As far as the computational cost is concerned, the most time-consuming operation is by far the calculation of the two chirp z-transforms on both channels, with all the other computations (including the z h,k update for all ps and spectral averaging) having a negligible impact. The computational cost, therefore, reduces to that of the calculation of 4 FFT (2 for each channel) of length 2N for each record of length N. Such computational cost is completely manageable even by low-end CPUs by today's standards.
In order to favour the experimentation of the proposed approach, we extended the public domain QLSA library [10] that now includes all the required data structures and computation algorithms for the implementation of narrow band spectral estimation approach discussed above. While developed with specific reference to the application in this 5 Journal of Sensors paper, the implementation in QLSA allows to extend the approach we propose to the estimation of narrow band spectra and cross spectra involving any number of channels [11]. A link to the web page where the library QLSA is described, with indications on how it can be obtained, is reported in Data Availability.

Experimental Results.
In order to test the proposed approach, we resorted to the circuit configuration in Figure 7. The complete list of the components is reported in Table 1. The first block in Figure 7 is a white noise source (WNS) obtained starting from two nominally identical 10V Zener diodes (D Z1 , D Z2 ) biased in breakdown with nominally identical currents (V DD = −V EE = 12 V; R Z1 = R Z2 = 200 kΩ). I A1 is an instrumentation amplifier (INA128) that rejects the (nominally) identical DC voltages at the cathodes of D Z1 and D Z2 and amplifies the uncorrelated noise signals across R Z1 and R Z2 . Note that notwithstanding that the Zener diodes and the resistances R Z are mounted in direct thermal contact with each other, a DC voltage difference in the order of few mV is present at the input of I A1 , so that a high-pass filter (C A1 , R A1 ) is used in order to reject the DC before the connection to the QTF and to the first input (V I ) of the data acquisition system. The transimpedance amplifier (TIA) stage in Figure 5 is the same as the one used in [3]. The resistance R R , required to provide DC coupling to the inverting input of O A1 , is such that at frequencies close to the resonance of the QTF, it can be completely neglected with respect to the effect of the feedback capacitor C R . By including the voltage gain stage (VGS), the overall transimpedance gain, from the inverting input of O A1 to the output V O , close to the QTF resonance frequency, is As discussed in [3], this solution has the advantage of allowing a quite straightforward recognition and correction of the deviation of the gain of the system from the approximated expression in Equation (10) obtained in the assump-tion of virtual ground at the input of O A1 . In particular, in the assumption of the ideal response in Equation (10) for the amplifying chain, the overall frequency response where [3] ,   Journal of Sensors with c PR = C P C R ; Note that H IF , as defined in Equation (12), has the physical dimensions of a frequency. Thanks to the peculiar circuit configuration used in Figure 7, H R and H IF have the following symmetry properties: From a graphical point of view, on a log scale for the frequency, Equation (14) corresponds to an even symmetry of H IF with respect to f S , while Equation (15) corresponds to an odd symmetry with respect to the point (f S , −A V c PR ) as shown in Figure 8 (solid lines). As mentioned before, Equation (10) is an approximation of the actual response obtained in the assumption of virtual ground. When high accuracy is required, however, the full response of the amplifying chain must be taken into account, which requires accounting for the effects of poles at higher frequencies [12]. As an example, suppose that a pole is present in the actual response at a frequency of, say, 20 times the resonance frequency of the QTF. This pole would introduce, at frequencies close to f S , a gain error of less than 0.2% in the response amplitude and a phase error in the order of -3 o . While the amplitude error can be usually neglected, the phase error may have a significant impact in the extraction of the QTF parameters [13]. The presence of a phase error is especially apparent with reference to H IF , as shown in Figure 8 (red curve) and can be corrected by introducing (numerically) the phase shift required for restoring symmetry in the measured response H IF [4]. The way in which the QTF parameters can be obtained from the measured H R and H IF , once symmetry is restored, is discussed in [3]. Here, we are interested in demonstrating how the new method we propose for spectra estimation can help in reaching, in a short time and in a reliable way, an estimate of the response H V ð f Þ by means of noise measurements. In order to do so, we will mostly focus on the properties of H IF and on the values of the resonance frequency f S and of the quality factor Q, which are the most important parameters that are taken into consideration in sensing applications.
When employing cross correlation, the frequency response H V is obtained as where S II is the PSD of the noise at the input V I and S IO is the cross spectrum between the noise at the input V I and the noise at the output V O in Figure 7. In our experiments, signal acquisition was performed using a National Instruments PCI-4462 four-channel dynamic signal acquisition board while a dedicated software was developed for PSD estimation according to the approach we propose. As we mentioned above, the software we have developed is based on an upgraded version of the public domain library QLSA [10]. Acquisition frequency was set to 204.8 kHz, and the record length was 8192 (2 13 ). We employed a wristwatch QTF exposed to air as a DUT to work with the highest possible Q (functionalization for sensing applications is expected to reduce the quality factor). A set of 10 different values of p was used in order to obtain geometrically decreasing (with a factor of 2) ENB from 12.5 Hz down to about 24 mHz in 7 Journal of Sensors a frequency range of 100 Hz across the resonance frequency (expected to be about 32760 Hz). A phase correction is performed in order to compensate for the distortion introduced by the nonideal behaviour of the amplifier. The correction, as discussed before, is performed numerically, and the correction angle is chosen as the one that ensures the best symmetry (Equations (14) and (15)). Figure 9 represents the situation for H IF after 60 s (top graph) and 240 s (bottom graph) since the measurement was started. As we have discussed before, curves corresponding to different values of ENB are all obtained at the same time.
The fact that the curves corresponding to ENB down to about 400 mHz do not superimpose to one another is a clear indication that these curves cannot represent the actual H IF . On the other hand, curves obtained with ENB below 200 mHz appear to be superimposed to one another, a clear indication that the corresponding ENB is sufficiently small not to induce systematic error in the estimate of H IF . As it was expected, the curves corresponding to larger ENBs appear smoother because a large number of uncorrelated averages were accumulated in the available time. As we move toward smaller and smaller ENBs, fewer uncorrelated averages were possible resulting in larger residual statistical errors. Note that, for the same ENB, as the measurement time increases, more uncorrelated averages are accumulated and the statistical fluctuations in the estimates are reduced, as it is apparent by comparing the situation after 60 s and after 240 s in Figure 9. When looking at a representation like the one in Figure 9 in real time (during experiments), the measurement can be stopped as soon as one detects the situation in which a curve from a given ENB superimposes to the one corresponding to the next smaller ENB. To evidence the effect of an incorrect (i.e., too large) ENB on the extraction of the QTF parameters, we fitted all the curves in Figure 9 against the expression of H IF in Equation (12). Because of the symmetry ofH IF (once the phase error is compensated), distortion due to large ENBs has a small or no effect on the frequency at whichH IF reaches its maximum (the resonance frequency). However, as it is apparent from Figure 10, a consistent (and therefore correct) estimate of the quality factor can only be obtained with ENBs below 200 mHz. As observed before, the threshold of 200 mHz for obtaining, in our test experiment, the correct estimate of Q is clearly detectable by observing the plots in Figure 9 in real time during measurements: all plots obtained with ENBs below 0.39 Hz essentially superimpose.
It is worth noticing that while curves obtained with low ENBs require longer measurement time for assuming a stable and highly regular shape, the error due to the low number of averages is statistical in nature. Therefore, if curve fitting over a frequency range is used for the estimation of the parameters we are interested in, such statistical errors tend to cancel out, even if they are still quite evident as in the plots in Figure 9. This is clearly demonstrated by the results reported in Figure 11, where the estimated value of f S , obtained from fitting the measured H IF for ENB = 195 mHz against the expression in Equation (12) Figure 10: Values of the quality factor Q extracted from the curves in Figure 9 as a function of the ENB. The values of Q (together with all other parameters) are obtained by fitting each curve in Figure 9 against the function H IF in Equation (12) in a frequency interval of 100 Hz centered across 32760 Hz.      We have already observed that the distortion due to spectral leakage preserves the symmetry of the estimated H IF . Indeed, when we plot the results of the extraction of f S for curves with different ENBs, as in Figure 12, we obtain that the estimated values of f s are almost the same (within a few ppm) regardless of the systematic error introduced by large values of ENBs. This is a clear advantage over the swept spectrum approach. If we indeed plot the results of the estimation of H IF as can be obtained with the swept spectrum approach on the very same QTF used for obtaining the results in Figure 9, we obtain the plots in Figure 13, where it is clear that even with a FCR of 1 Hz/s we still cannot conclude, with an uncertainty comparable to what is obtained in Figure 12, that the correct value of f s is obtained at the frequency for which the curve has a maximum.
When dealing with functionalized QTF in sensing applications, a significant reduction of Q is generally expected. A lower Q, as observed before, translates into a large bandwidth B. Since, as a consequence of a larger B, the ENB below which the curves in Figure 9 will superimpose increase, this means that the required measurement time for accurate results will decrease as well, so that a complete characterization of the QTF parameters can be typically obtained in a matter of a few seconds or, at most, a few tens of seconds.

Discussion and Conclusions
We have developed an effective approach for the accurate determination of QTF parameters based on noise measurements. In order to reach this goal, we had to address and solve the problem of the accurate estimation of the spectra of noise signals whose PSD changes significantly over a narrow bandwidth. The approach we have developed relies on a proper elaboration of the sampled noise signals that allows spectral and cross spectral estimation to be performed using several values of ENB at the same time. This allows the user to easily select the maximum ENB that do not cause systematic errors and that can therefore lead to the correct estimate of the QTF parameters in the shortest time possible. Since the results corresponding to each ENBs are all available at the same time, the need for repeated measurements is eliminated. A software library for the easy integration of the approach we propose in virtual instrumentation applications has been developed, and it is part of the QLSA public domain package [10]. While we used a high-end signal acquisition board for the measurements reported in this paper, lowcost sound boards for PC are widely available that may allow the implementation of the approach we propose at a very low cost [13]. We believe that the ability to accurately estimate the PSD over a narrow bandwidth across the resonance frequency of a QTF-based sensor may facilitate the exploration of fluctuation enhanced sensing (FES) with QTF-based sensors. Evidence of the potentialities of FES approaches in the case of quartz crystal microbalances (QCM) based on AT-cut quartz crystals can be found in the literature. In [14], for instance, it is stipulated that part of the noise observed in a 10 MHz oscillator driven by a functionalized crystal is due to the interaction of the sensing material with the chemicals in the environment. In this context, the availability of the library we have developed may contribute to the development of sensitive experiments for the exploration of FES approaches in QTF-based sensors.

Data Availability
The data used to support the findings of this study are in part included within the article and available from the corresponding author upon request. Details on the operation and usage of QLSA can be found online at http://iee1 .unime.it/QLSA.

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.  Figure 12: Values of the resonance frequency f S extracted from the curves in Figure 9 as a function of the ENB. The values of f S (together with all other parameters) are obtained by fitting each curve in Figure 9 against the function H IF in Equation (12) Figure 13: Estimation of H IF obtained with the swept spectrum approach with different sweep speeds (FCR). For these measurements, the WNS block in Figure 7 was replaced with a swept sinusoidal source. 9 Journal of Sensors