# Accurate QTF Sensing Approach by Means of Narrow Band Spectral Estimation

Authors: Graziella Scandurra; Gino Giusi; Carmine Ciofi
Journal: Journal of Sensors
 (2020)

---

## Abstract

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.


---

## Body

## 1. 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 fS 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.Figure 1
Equivalent circuit of a quartz tuning fork.The extraction of all parameters in Figure1 can be obtained starting from the measurement of the response of circuits such as the one in Figure 2 where assuming the transconductance gain AR to be constant, the voltage gain VO/VI is proportional to the QTF admittance YQ:
(1)HVf=VOVI=ARYQFigure 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 gainAR.However, if we require high accuracy, with the circuit in Figure2, 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 Figure2, 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 input-output cross correlation. In the case of swept spectrum measurements, VI 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 correlation-based 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 104, 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 RP→∞ and CP=0. We have also assumed thatAR=RS, so that at the resonance frequency we have, from Equation (1), HVfS=YQfSRS=1.Figure 3
Simulated measurement of the frequency responseHVf with the swept spectrum approach assuming an ideal TIA with AR=RS, fS=32.7kHz, Q=104, and CP=0; RP→∞. FCR is the frequency change rate, ranging from 1 Hz/s up to 80 Hz/s. The error fERR in assuming the frequency of the maximum of the response as the resonance frequency is reported in the inset, normalized with respect to B=fS/Q, as a function of the FCR.Figure 4
Simulated results of the effect of the choice of the resolution bandwidthΔf on the measurement of the response HVf of the circuit in Figure 2 by means of the cross spectrum approach (with CP=0, RP→∞). The ideal response is represented by the thick black curve.The results in Figure3 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 (fMAX) as an estimate of the QTF resonance frequency fS, we would obtain an estimation error fERR=fS−fMAX 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.Figure4 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=fS/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 fC, the frequency resolution is determined by the number of samples N used for the estimation of the FFT on each single signal record (Δf=fC/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 fC=200kHz, we need N in excess of 220 (≈106) 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 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 VI and VO in Figure 2 to a set of FFT spectrum 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 Figure2 (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 fC (sampling period Δt=1/fC). The sequence of samples xi=xiΔt obtained from the filtered input xt is divided into records of length N, and each record is multiplied by a proper window function wjj=0,⋯.,N−1 obtaining a new set of records yh,l:
(2)yh,l=xhN+lwl,l=0,1,⋯,N−1,where we have assumed that the first element of each record is xhN (h = 0,1,…) so that the records are adjacent and non-overlapping. The Discrete Fourier Transform (DFT) of each record is calculated obtaining
(3)zh,k=∑l=0N−1yh,le−j2πkNl;k=0,1,⋯,N−1with zh,k2 representing an estimate, for k<N/2, of the power of the input signal after being filtered by a bandpass filter centred at fk=kfCK/N whose equivalent noise bandwidth (ENB) and detailed frequency response depend on fC, N, and on the window function w. If the window function satisfies the normalization condition
(4)∑l=0N−1wl2=1,and if we assume that the PSD of the input signal can be considered approximately constant within the bandwidth of the filter (white noise approximation), the quantity zh,k2/ENB can be regarded as an estimate of the PSD of the input signal at fk as obtained from the elaboration of the hth 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 fC/N up to 1.5fC/N [7] so that, in a first approximation, it essentially coincides with the frequency resolution Δf=fC/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<<fs/Q that may correspond to record lengths in the order of 106 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 fC/2, but we can limit our focus to a few bandwidths B (B=fS/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 fS 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.Figure 5
Continuous time equivalent of the approach used for the estimation of the PSD in DFT spectrum analyzers.The block diagram in Figure5 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 fk. 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 fk. 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 zh,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 wt 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 wt with its sampled counterpart wi=αwiΔt. Note that the parameter α has the physical dimension of a time, and it is required in order to obtain a dimensionless wi. Note that the values of α and of the parameter A defining the amplitude of wt 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 thatxi=0fori<0. In the discrete time domain, the output zh,k of the filter in Figure 6 at i=h+1NΔt is
(5)zh+1,k=∑l=0h+1N−1xlwh+1N−1−le−j2πkNl==∑l=0hN−1xlwh+1N−1−le−j2πkNl++∑l=hNh+1N−1xlwh+1N−1−le−j2πkNl.Figure 6
Impulse response of the continuous time domain filter from which the actual window function used in the proposed approach is derived. The parameterp ranges from 0 to 1.With the index changem=l−hN in the last sum and exploiting the following property for wi,
(6)wi+N=pwi,we have
(7)zh+1,k=p∑l=0hN−1xlwh+1N−1−le−j2πkNl++e−j2πkh∑m=0N−1xhN+mwN−1−me−j2πkNm==pzh,k+αAe−j2πkh∑m=0N−1xhN+me−j2πkNm.Note that the last sum overm, 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 zh,k recursively, starting from the evaluation of the DFT over the last N samples. What makes the function wi 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 responseWf corresponding to the impulse response w(t) is
(8)Wf=ATsinπfTπfTe−jπfT1−pe−j2πfT;T=N∆t=NfC,with an ENB:
(9)ENB=1−p1+pfCN.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 zh,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 atN closely spaced frequencies across the (expected) series resonance frequency(3)
We set up a number of values forp leading to geometrically decreasing values of ENB (typically with a common ratio of 2)(4)
Each time a new input record (lengthN) 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(5)
For each value ofp, we update the corresponding vectors of values zh,k (one vector of N complex values for each value of p)(6)
The values ofzh,k corresponding to each p are used for spectral and cross spectral estimationAbout the 6th step, it must be remarked the fact that while the zh,k corresponding to different p are all updated at time intervals equal to the duration of the input record, subsequent values of zh,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 zh,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 assumeN=8192. With fCK=100kHz, the ENB corresponding to p=0 is about 12 Hz. AssumingQMAX=10000, we should be able to obtain a minimum ENB significantly smaller thanfs/QMAX≈3Hz. 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 zh,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 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.


### 1.2. Experimental Results

In order to test the proposed approach, we resorted to the circuit configuration in Figure7. 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 (DZ1, DZ2) biased in breakdown with nominally identical currents (VDD=−VEE=12V; RZ1=RZ2=200kΩ). IA1 is an instrumentation amplifier (INA128) that rejects the (nominally) identical DC voltages at the cathodes of DZ1 and DZ2 and amplifies the uncorrelated noise signals across RZ1 and RZ2. Note that notwithstanding that the Zener diodes and the resistances RZ 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 IA1, so that a high-pass filter (CA1, RA1) is used in order to reject the DC before the connection to the QTF and to the first input (VI) of the data acquisition system. The transimpedance amplifier (TIA) stage in Figure 5 is the same as the one used in [3]. The resistance RR, required to provide DC coupling to the inverting input of OA1, 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 CR. By including the voltage gain stage (VGS), the overall transimpedance gain, from the inverting input of OA1 to the output VO, close to the QTF resonance frequency, is
(10)AR=−AVj2πfCR,AV=11.Figure 7
Detailed schematic of the circuit used for testing the approach we propose. The main blocks are a white noise source (WNS), a transimpedance amplifier (TIA) with capacitive feedback, and a voltage gain stage (VGS).Table 1
Component list for the circuit in Figure5.
ComponentTypeValueIA1Instrumentation amplifierIF3601OA1, OA2Operational amplifiersOP27YQQTF under test—RZ1, RZ20.1% metallic film resistor200 kΩRGMulti turn trimmer20 kΩRA11% resistor10 kΩCA1Polyester, 66 V4.7μFRA21% resistor10 kΩCA2Polyester, 66 V1μFRR1% resistor10 MΩCR1%100 pFR10.1% metallic film resistor1 kΩR20.1% metallic film resistor10 kΩ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 assumption of virtual ground at the input of OA1. In particular, in the assumption of the ideal response in Equation (10) for the amplifying chain, the overall frequency response HV from VI to VO would be
(11)HVf=VOVI=YQAR=HR+jHI,where [3]
(12)HRf=−AVcPR+cSR×Q21−f/fS2Q21−f/fS22+f/fS2,HIFf=fHIf=AVQfScSRrSP+cSR×f/fS2Q21−f/fS22+f/fS2,with
(13)cPR=CPCR;cSR=CSCR;rSP=RSRP.Note thatHIF, as defined in Equation (12), has the physical dimensions of a frequency. Thanks to the peculiar circuit configuration used in Figure 7, HR and HIF have the following symmetry properties:
(14)HIFffS=HIFfSfS,(15)HRf/fS+AVcPR=−HRfSf−AVcPR.From a graphical point of view, on a log scale for the frequency, Equation (14) corresponds to an even symmetry of HIF with respect to fS, while Equation (15) corresponds to an odd symmetry with respect to the point (fS, −AVcPR) 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 fS, a gain error of less than 0.2% in the response amplitude and a phase error in the order of -3o. 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 HIF, 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 HIF [4]. The way in which the QTF parameters can be obtained from the measured HR and HIF, 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 HVf by means of noise measurements. In order to do so, we will mostly focus on the properties of HIF and on the values of the resonance frequency fS and of the quality factor Q, which are the most important parameters that are taken into consideration in sensing applications.Figure 8
Simulated behaviour ofHR and HIF for a QTF operating at 32.8 kHz. The red curve is obtained in the presence of a phase error of -3 degrees in the frequency response.When employing cross correlation, the frequency responseHV is obtained as
(16)HVf=SIOfSIIf,where SII is the PSD of the noise at the input VI and SIO is the cross spectrum between the noise at the input VI and the noise at the output VO 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 (213). 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 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 HIF 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.Figure 9
Result of the estimation ofHIF at TM=60s and at TM=240s since the measurement was started. All curves are obtained simultaneously in real 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 actualHIF. 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 HIF. 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 HIF in Equation (12). Because of the symmetry ofHIF(once the phase error is compensated), distortion due to large ENBs has a small or no effect on the frequency at whichHIFreaches 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.Figure 10
Values of the quality factorQ 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 HIF in Equation (12) in a frequency interval of 100 Hz centered across 32760 Hz.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 Figure9. This is clearly demonstrated by the results reported in Figure 11, where the estimated value of fS, obtained from fitting the measured HIF for ENB=195mHz against the expression in Equation (12), is reported vs the measurement time. As it can be observed, the estimate changes by less than 1 ppm from 30 s up to 480 s.Figure 11
Values of resonance frequencyfS extracted from the estimate of HIF with ENB=195mHz vs measurement time. The values of fS (together with all other parameters) are obtained by fitting the estimated HIF against the function HIF 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 estimatedHIF. Indeed, when we plot the results of the extraction of fS for curves with different ENBs, as in Figure 12, we obtain that the estimated values of fs 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 HIF 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 fs is obtained at the frequency for which the curve has a maximum.Figure 12
Values of the resonance frequencyfS extracted from the curves in Figure 9 as a function of the ENB. The values of fS (together with all other parameters) are obtained by fitting each curve in Figure 9 against the function HIF in Equation (12) in a frequency interval of 100 Hz centered across 32760 Hz.Figure 13
Estimation ofHIF 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.When dealing with functionalized QTF in sensing applications, a significant reduction ofQ 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.


## 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 fC (sampling period Δt=1/fC). The sequence of samples xi=xiΔt obtained from the filtered input xt is divided into records of length N, and each record is multiplied by a proper window function wjj=0,⋯.,N−1 obtaining a new set of records yh,l:
(2)yh,l=xhN+lwl,l=0,1,⋯,N−1,where we have assumed that the first element of each record is xhN (h = 0,1,…) so that the records are adjacent and non-overlapping. The Discrete Fourier Transform (DFT) of each record is calculated obtaining
(3)zh,k=∑l=0N−1yh,le−j2πkNl;k=0,1,⋯,N−1with zh,k2 representing an estimate, for k<N/2, of the power of the input signal after being filtered by a bandpass filter centred at fk=kfCK/N whose equivalent noise bandwidth (ENB) and detailed frequency response depend on fC, N, and on the window function w. If the window function satisfies the normalization condition
(4)∑l=0N−1wl2=1,and if we assume that the PSD of the input signal can be considered approximately constant within the bandwidth of the filter (white noise approximation), the quantity zh,k2/ENB can be regarded as an estimate of the PSD of the input signal at fk as obtained from the elaboration of the hth 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 fC/N up to 1.5fC/N [7] so that, in a first approximation, it essentially coincides with the frequency resolution Δf=fC/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<<fs/Q that may correspond to record lengths in the order of 106 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 fC/2, but we can limit our focus to a few bandwidths B (B=fS/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 fS 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.Figure 5
Continuous time equivalent of the approach used for the estimation of the PSD in DFT spectrum analyzers.The block diagram in Figure5 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 fk. 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 fk. 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 zh,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 wt 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 wt with its sampled counterpart wi=αwiΔt. Note that the parameter α has the physical dimension of a time, and it is required in order to obtain a dimensionless wi. Note that the values of α and of the parameter A defining the amplitude of wt 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 thatxi=0fori<0. In the discrete time domain, the output zh,k of the filter in Figure 6 at i=h+1NΔt is
(5)zh+1,k=∑l=0h+1N−1xlwh+1N−1−le−j2πkNl==∑l=0hN−1xlwh+1N−1−le−j2πkNl++∑l=hNh+1N−1xlwh+1N−1−le−j2πkNl.Figure 6
Impulse response of the continuous time domain filter from which the actual window function used in the proposed approach is derived. The parameterp ranges from 0 to 1.With the index changem=l−hN in the last sum and exploiting the following property for wi,
(6)wi+N=pwi,we have
(7)zh+1,k=p∑l=0hN−1xlwh+1N−1−le−j2πkNl++e−j2πkh∑m=0N−1xhN+mwN−1−me−j2πkNm==pzh,k+αAe−j2πkh∑m=0N−1xhN+me−j2πkNm.Note that the last sum overm, 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 zh,k recursively, starting from the evaluation of the DFT over the last N samples. What makes the function wi 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 responseWf corresponding to the impulse response w(t) is
(8)Wf=ATsinπfTπfTe−jπfT1−pe−j2πfT;T=N∆t=NfC,with an ENB:
(9)ENB=1−p1+pfCN.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 zh,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 atN closely spaced frequencies across the (expected) series resonance frequency(3)
We set up a number of values forp leading to geometrically decreasing values of ENB (typically with a common ratio of 2)(4)
Each time a new input record (lengthN) 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(5)
For each value ofp, we update the corresponding vectors of values zh,k (one vector of N complex values for each value of p)(6)
The values ofzh,k corresponding to each p are used for spectral and cross spectral estimationAbout the 6th step, it must be remarked the fact that while the zh,k corresponding to different p are all updated at time intervals equal to the duration of the input record, subsequent values of zh,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 zh,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 assumeN=8192. With fCK=100kHz, the ENB corresponding to p=0 is about 12 Hz. AssumingQMAX=10000, we should be able to obtain a minimum ENB significantly smaller thanfs/QMAX≈3Hz. 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 zh,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 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.


## 1.2. Experimental Results

In order to test the proposed approach, we resorted to the circuit configuration in Figure7. 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 (DZ1, DZ2) biased in breakdown with nominally identical currents (VDD=−VEE=12V; RZ1=RZ2=200kΩ). IA1 is an instrumentation amplifier (INA128) that rejects the (nominally) identical DC voltages at the cathodes of DZ1 and DZ2 and amplifies the uncorrelated noise signals across RZ1 and RZ2. Note that notwithstanding that the Zener diodes and the resistances RZ 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 IA1, so that a high-pass filter (CA1, RA1) is used in order to reject the DC before the connection to the QTF and to the first input (VI) of the data acquisition system. The transimpedance amplifier (TIA) stage in Figure 5 is the same as the one used in [3]. The resistance RR, required to provide DC coupling to the inverting input of OA1, 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 CR. By including the voltage gain stage (VGS), the overall transimpedance gain, from the inverting input of OA1 to the output VO, close to the QTF resonance frequency, is
(10)AR=−AVj2πfCR,AV=11.Figure 7
Detailed schematic of the circuit used for testing the approach we propose. The main blocks are a white noise source (WNS), a transimpedance amplifier (TIA) with capacitive feedback, and a voltage gain stage (VGS).Table 1
Component list for the circuit in Figure5.
ComponentTypeValueIA1Instrumentation amplifierIF3601OA1, OA2Operational amplifiersOP27YQQTF under test—RZ1, RZ20.1% metallic film resistor200 kΩRGMulti turn trimmer20 kΩRA11% resistor10 kΩCA1Polyester, 66 V4.7μFRA21% resistor10 kΩCA2Polyester, 66 V1μFRR1% resistor10 MΩCR1%100 pFR10.1% metallic film resistor1 kΩR20.1% metallic film resistor10 kΩ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 assumption of virtual ground at the input of OA1. In particular, in the assumption of the ideal response in Equation (10) for the amplifying chain, the overall frequency response HV from VI to VO would be
(11)HVf=VOVI=YQAR=HR+jHI,where [3]
(12)HRf=−AVcPR+cSR×Q21−f/fS2Q21−f/fS22+f/fS2,HIFf=fHIf=AVQfScSRrSP+cSR×f/fS2Q21−f/fS22+f/fS2,with
(13)cPR=CPCR;cSR=CSCR;rSP=RSRP.Note thatHIF, as defined in Equation (12), has the physical dimensions of a frequency. Thanks to the peculiar circuit configuration used in Figure 7, HR and HIF have the following symmetry properties:
(14)HIFffS=HIFfSfS,(15)HRf/fS+AVcPR=−HRfSf−AVcPR.From a graphical point of view, on a log scale for the frequency, Equation (14) corresponds to an even symmetry of HIF with respect to fS, while Equation (15) corresponds to an odd symmetry with respect to the point (fS, −AVcPR) 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 fS, a gain error of less than 0.2% in the response amplitude and a phase error in the order of -3o. 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 HIF, 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 HIF [4]. The way in which the QTF parameters can be obtained from the measured HR and HIF, 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 HVf by means of noise measurements. In order to do so, we will mostly focus on the properties of HIF and on the values of the resonance frequency fS and of the quality factor Q, which are the most important parameters that are taken into consideration in sensing applications.Figure 8
Simulated behaviour ofHR and HIF for a QTF operating at 32.8 kHz. The red curve is obtained in the presence of a phase error of -3 degrees in the frequency response.When employing cross correlation, the frequency responseHV is obtained as
(16)HVf=SIOfSIIf,where SII is the PSD of the noise at the input VI and SIO is the cross spectrum between the noise at the input VI and the noise at the output VO 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 (213). 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 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 HIF 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.Figure 9
Result of the estimation ofHIF at TM=60s and at TM=240s since the measurement was started. All curves are obtained simultaneously in real 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 actualHIF. 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 HIF. 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 HIF in Equation (12). Because of the symmetry ofHIF(once the phase error is compensated), distortion due to large ENBs has a small or no effect on the frequency at whichHIFreaches 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.Figure 10
Values of the quality factorQ 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 HIF in Equation (12) in a frequency interval of 100 Hz centered across 32760 Hz.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 Figure9. This is clearly demonstrated by the results reported in Figure 11, where the estimated value of fS, obtained from fitting the measured HIF for ENB=195mHz against the expression in Equation (12), is reported vs the measurement time. As it can be observed, the estimate changes by less than 1 ppm from 30 s up to 480 s.Figure 11
Values of resonance frequencyfS extracted from the estimate of HIF with ENB=195mHz vs measurement time. The values of fS (together with all other parameters) are obtained by fitting the estimated HIF against the function HIF 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 estimatedHIF. Indeed, when we plot the results of the extraction of fS for curves with different ENBs, as in Figure 12, we obtain that the estimated values of fs 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 HIF 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 fs is obtained at the frequency for which the curve has a maximum.Figure 12
Values of the resonance frequencyfS extracted from the curves in Figure 9 as a function of the ENB. The values of fS (together with all other parameters) are obtained by fitting each curve in Figure 9 against the function HIF in Equation (12) in a frequency interval of 100 Hz centered across 32760 Hz.Figure 13
Estimation ofHIF 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.When dealing with functionalized QTF in sensing applications, a significant reduction ofQ 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.


## 2. 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, low-cost 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.



---

*Source: 8951340-2020-06-08.xml*