Robust and Efficient Parametric Spectral Estimation in Atomic Force Microscopy

TL;DR

This paper introduces a robust two-stage spectral estimation method for atomic force microscopy data, significantly improving accuracy and noise robustness while reducing computational costs.

Contribution

It presents a novel two-stage estimator combining variance stabilization and periodicity testing to enhance spectral parameter estimation in AFM data.

Findings

Two- to ten-fold reduction in mean squared error

Robustness against electronic noise

Efficient parameter estimation with reduced computational cost

Abstract

An atomic force microscope (AFM) is capable of producing ultra-high resolution measurements of nanoscopic objects and forces. It is an indispensable tool for various scientific disciplines such as molecular engineering, solid-state physics, and cell biology. Prior to a given experiment, the AFM must be calibrated by fitting a spectral density model to baseline recordings. However, since AFM experiments typically collect large amounts of data, parameter estimation by maximum likelihood can be prohibitively expensive. Thus, practitioners routinely employ a much faster least-squares estimation method, at the cost of substantially reduced statistical efficiency. Additionally, AFM data is often contaminated by periodic electronic noise, to which parameter estimates are highly sensitive. This article proposes a two-stage estimator to address these issues. Preliminary parameter estimates are first obtained by a variance-stabilizing procedure, by which the simplicity of least-squares combines with the efficiency of maximum likelihood. A test for spectral periodicities then eliminates high-impact outliers, considerably and robustly protecting the second-stage estimator from the effects of electronic noise. Simulation and experimental results indicate that a two- to ten-fold reduction in mean squared error can be expected by applying our methodology.

Tables

Table 1. Table 1 : SHO parameters in baseline environment.

SHO Parameter	Value
Temperature	$T=298K\text{/}$
Stiffness	$k=0.172N\text{/}m$
Resonance Frequency	$f_0=33.533\mathrm{kHz}\text{/}$
Quality Factor	$Q\in \{1,10,100,500\}$

Table 2. Table 2 : Real AFM cantilever parameter estimates and standard errors.

	${\hat{f}}_0$ ($\mathrm{kHz}\text{/}$)	$\hat{Q}$ (unitless)	$\hat{k}$ ($N\text{/}m$)
NLS	312.703 (.0043)	603.01 (14.28)	57.52 (0.91)
LP	312.701 (.0047)	595.40 (12.74)	57.16 (0.81)
MLE	312.700 (.0048)	593.58 (12.66)	57.20 (0.81)

Table 3. Table 3 : Asymptotic relative bias for each parameter resulting from failure to account for the 1 / f 1 𝑓 1/f noise in Figure  10 .

	$f_0$	$Q$	$k$
$Q=1$	1.02	.92	.92
$Q=10$	1.00	.89	.98
$Q=100$	1.00	1.05	1.00
$Q=500$	1.00	1.10	1.00

Table 4. Table 4 : Relative MSE of NLS and LP estimators to MLE in the presence of 1 / f 1 𝑓 1/f noise. For Q = 1 𝑄 1 Q=1 some of the estimators failed to converge, such that ratios are based on 1000 datasets for NLS, 996 for LP, and 768 for MLE.

	Method	$f_0$ MSE Ratio	$Q$ MSE Ratio	$k$ MSE Ratio
Q = 1	NLS	2.23	1.97	2.42
	LP	1.16	1.51	1.78
Q = 10	NLS	2.59	3.74	1.97
	LP	1.01	1.01	1.01
Q = 100	NLS	3.07	5.27	1.48
	LP	1	1	1.01
Q = 500	NLS	3.25	5.94	1.62
	LP	1.01	1	1

Equations

$$\mathcal{S}(f)={\int }_{-\infty }^{\infty }{e}^{-2\pi itf}\gamma (t)\text{operator@font}dt=\mathcal{F}\{\gamma (t)\}.$$

$$m{\ddot{X}}_t=-k{X}_t-\varsigma {\dot{X}}_t+F_t,$$

$$\mathcal{S}(f)=\frac{k_{B}T/(k\cdot\pi f_{0}Q)}{\big{[}(f/f_{0})^{2}-1\big{]}^{2}+\big{[}f/(f_{0}Q)\big{]}^{2}},$$

$${\tilde{X}}_k=\sum _{n=0}^{N-1}{e}^{-2\pi ikn/N}{X}_n.$$

$$\frac{1}{N}{\left|{\tilde{X}}_k\right|}^2\stackrel{\text{operator@font}ind}{\sim }{f}_s\cdot \mathcal{S}(f_k)\times \text{Expo}(1),0\le k<K.$$

$$\ell_{W}({\bm{\theta}}\,|\,{\bm{Y}})=-\sum_{k=1}^{K}\big{(}Y_{k}/\mathcal{S}_{k}({\bm{\theta}})+\log\mathcal{S}_{k}({\bm{\theta}})\big{)},$$

$${\bar{Y}}_m=\frac{1}{B}\sum _{k\in I_m}{Y}_k,I_m=\{k:(m-1)B<k\le mB\}.$$

$${\bar{Y}}_m\stackrel{\text{operator@font}ind}{\sim }{\bar{\mathcal{S}}}_m(\mathbf{\theta })\times \text{Gamma}(B,B),$$

$$\hat{\bm{\theta}}_{\textnormal{NLS}}=\operatorname*{arg\,min}_{\bm{\theta}}\sum_{m=1}^{N_{B}}\big{(}\bar{Y}_{m}-\bar{\mathcal{S}}_{m}({\bm{\theta}})\big{)}^{2},$$

$$\bar{Y}_{m}\stackrel{{\scriptstyle\operator@font{ind}}}{{\sim}}\mathcal{N}\big{(}\bar{\mathcal{S}}_{m}({\bm{\theta}}),\tfrac{1}{B}\bar{\mathcal{S}}_{m}({\bm{\theta}})^{2}\big{)}.$$

$$\mathrm{var}(\log {\bar{Y}}_m)\approx {\left({\frac{\mathrm{d}\log y}{\mathrm{d}y}|}_{y=E[{\bar{Y}}_m]}\right)}^2\times \mathrm{var}({\bar{Y}}_m)=\frac{1}{{\bar{\mathcal{S}}}_m(\mathbf{\theta }{)}^2}\times \frac{{\bar{\mathcal{S}}}_m(\mathbf{\theta }{)}^2}{B}=1/B,$$

$$Z_{m}=\log(\bar{Y}_{m})\stackrel{{\scriptstyle\operator@font{ind}}}{{\sim}}\mathcal{N}\big{(}\log\bar{\mathcal{S}}_{m}({\bm{\theta}}),B^{-1}\big{)}.$$

$$\hat{\bm{\theta}}_{{\textnormal{LP}}}=\operatorname*{arg\,min}_{\bm{\theta}}\sum_{m=1}^{N_{B}}\big{(}Z_{m}-\log\bar{\mathcal{S}}_{m}({\bm{\theta}})\big{)}^{2}.$$

$$H_0:W_k=\frac{Y_k}{{\mathcal{S}}_k({\mathbf{\theta }}_0)}\stackrel{\text{operator@font}iid}{\sim }\text{Expo}(1),$$

$$\mathcal{M}=\underset{1\le k\le K}{\max }\frac{W_k}{{\sum }_{j=1}^K{W}_j}.$$

$$\mathrm{Pr}(\mathcal{M}>a|H_0)=\sum _{k=1}^K(-1{)}^{k+1}\left(\genfrac{}{}{0ex}{}{K}{k}\right)(1-k\cdot a{)}_{+}^{K-1},x_{+}:=\max (x,0).$$

$$\mathcal{S}(f\,|\,{\bm{\theta}})=A_{\textrm{w}}+\frac{k_{B}T/(k\cdot\pi f_{0}Q)}{\big{[}(f/f_{0})^{2}-1\big{]}^{2}+\big{[}f/(f_{0}Q)\big{]}^{2}},$$

$${\mathcal{R}}_j(\phi )=\frac{{\sum }_{i=1}^M({\hat{\phi }}_j^{(i)}-{\phi }_0{)}^2}{{\sum }_{i=1}^M({\hat{\phi }}_W^{(i)}-{\phi }_0{)}^2},$$

$${\mathcal{R}}_j(\phi )=\frac{{\sum }_{i=1}^M({\hat{\phi }}_{j,\text{noise}}^{(i)}-{\phi }_0{)}^2}{{\sum }_{i=1}^M({\hat{\phi }}_{j,\text{base}}^{(i)}-{\phi }_0{)}^2},$$

$${\mathcal{R}}_j(\phi )=\frac{{\sum }_{i=1}^M({\hat{\phi }}_{j,\text{corr}}^{(i)}-{\phi }_0{)}^2}{{\sum }_{i=1}^M({\hat{\phi }}_{W,\text{base}}^{(i)}-{\phi }_0{)}^2},$$

$$\mathcal{S}(f\,|\,k,f_{0},Q,A_{\textrm{w}},A_{\textrm{f}},\alpha)=A_{\textrm{w}}+\frac{A_{\textrm{f}}}{f^{\alpha}}+\frac{k_{B}T/(k\cdot\pi f_{0}Q)}{\big{[}(f/f_{0})^{2}-1\big{]}^{2}+\big{[}f/(f_{0}Q)\big{]}^{2}}.$$

$$\mathcal{S}(f|\mathbf{\theta })$$

$$\displaystyle=\tau\times\left\{R_{\textrm{w}}+\frac{1}{\big{[}(f/f_{0})^{2}-1\big{]}^{2}+(f/\gamma)^{2}}\right\}=\tau\cdot\mathcal{G}(f\,|\,{\bm{\eta}}),$$

$$\mathcal{Q}_{\textnormal{NLS}}(\tau,{\bm{\eta}})=\sum_{m=1}^{N_{B}}\big{(}\bar{Y}_{m}-\tau\cdot\bar{\mathcal{G}}_{m}({\bm{\eta}})\big{)}^{2},$$

$$\hat{\tau }(\mathbf{\eta })=\underset{\tau }{\mathrm{argmin}}{\mathcal{Q}}_{\text{NLS}}(\tau ,\mathbf{\eta })=\frac{{\sum }_{m=1}^{N_B}{\bar{\mathcal{G}}}_m(\mathbf{\eta })\cdot {\bar{Y}}_m}{{\sum }_{m=1}^{N_B}{\bar{\mathcal{G}}}_m(\mathbf{\eta }{)}^2}.$$

$$\mathcal{Q}_{\textnormal{LP}}({\bm{\theta}},{\bm{\eta}})=\sum_{m=1}^{N_{B}}\big{(}Z_{m}-\log(\tau)-\log\bar{\mathcal{G}}_{m}({\bm{\eta}})\big{)}^{2},$$

$$\hat{\tau}({\bm{\eta}})=\operatorname*{arg\,min}_{\tau}\mathcal{Q}_{\textnormal{LP}}({\bm{\theta}},{\bm{\eta}})=\exp\left\{\frac{1}{N_{B}}\sum_{m=1}^{N_{B}}\big{(}Z_{m}-\log\bar{\mathcal{G}}_{m}({\bm{\eta}})\big{)}\right\}.$$

$${\mathcal{Q}}_{\text{MLE}}(\mathbf{\theta },\mathbf{\eta })=\sum _{k=1}^K\left(\frac{Y_k}{\tau \cdot {\mathcal{G}}_k(\mathbf{\eta })}+\log (\tau )+\log {\mathcal{G}}_k(\mathbf{\eta })\right),$$

$$\hat{\tau }(\mathbf{\eta })=\underset{\tau }{\mathrm{argmin}}{\mathcal{Q}}_{\text{MLE}}(\mathbf{\theta },\mathbf{\eta })=\frac{1}{K}\sum _{k=1}^K\frac{Y_k}{{\mathcal{G}}_k(\mathbf{\eta })}.$$

$$\mathcal{S}(f\,|\,k,f_{0},Q,A_{\textrm{w}},A_{\textrm{f}},\alpha)=A_{\textrm{w}}+\frac{A_{\textrm{f}}}{f^{\alpha}}+\frac{k_{B}T/(k\cdot\pi f_{0}Q)}{\big{[}(f/f_{0})^{2}-1\big{]}^{2}+\big{[}f/(f_{0}Q)\big{]}^{2}},$$

Full text

Key Words: Cantilever Calibration; Periodogram; Whittle Likelihood; Variance-Stabilizing Transformation; Fisher’s $g$-Statistic.

1 Introduction

An atomic force microscope (AFM) is a scientific instrument producing high-frequency ultra-precise displacement readings of a minuscule ($\sim 100\text{\,}\mathrm{\SIUnitSymbolMicro m}\text{/}$-long) pliable beam referred to as a cantilever. The cantilever bends in response to various forces exerted by its surrounding environment, the recording of which has been immensely useful for the study of e.g., the composition of polymers and other chemical compounds (Sugimoto et al., 2007; García et al., 2007), interatomic and intramolecular forces (Radmacher, 1997; Hoffmann et al., 2001), pathogen-drug interactions (Alsteens et al., 2008), cell adhesion (Evans and Calderwood, 2007), and the dynamics of protein folding (Yu et al., 2017).

In a typical AFM experiment, the cantilever’s bending response is measured in opposition to its spring-like restoring force, which requires proper calibration of the cantilever stiffness in order to convert measured displacement readings into force (Cleveland et al., 1993; Burnham et al., 2002; Clarke et al., 2006; Sader et al., 2011). This calibration is accomplished by fitting various parametric models to a baseline spectral density recording, i.e., to a cantilever driven by thermal noise alone. A representative baseline spectrum calculated from experimental data is displayed in Figure 1.

This data serves to illustrate two outstanding challenges in AFM parametric spectral density estimation. First, the experiments produce massive amounts of data, for which maximum likelihood estimation can be prohibitively expensive. A much faster least-squares method is routinely employed in practice (Nørrelykke and Flyvbjerg, 2010) – at the cost of substantially reduced statistical efficiency. Second, parametric estimates by either method are severely affected by electronic noise, due to periodic fluctuations in the AFM’s circuitry. Such noise is evidenced by the presence of sharp peaks (i.e., vertical lines) in the baseline spectrum of Figure 1.

In this article, we propose a two-stage estimator addressing both of these issues. A preliminary estimator first applies a variance-stabilizing transformation which renders the least-squares estimator virtually as efficient as the MLE. After the preliminary fit, an automated denoising procedure, based on a well-known statistical test for hidden periodicities, robustly protects the second-stage estimator from most of the effects of electronic noise. Extensive simulations and experimental results indicate that a two- to ten-fold reduction in mean squared error can be expected by applying our methodology.

The remainder of the paper is organized as follows. Section 2 provides an overview of parametric spectral density estimation in the AFM context. Section 3 describes our proposed two-step estimator. Section 4 presents a detailed simulation study comparing our proposal to existing methods. Section 5 applies the methodology to calibration of a real AFM cantilever and we close in Section 6 with a discussion of future work.

2 Parametric Spectral Density Estimation in AFM

Let $X_{t}=X(t)$ denote a continuous stationary Gaussian stochastic process with mean $E[X_{t}]=0$ and autocorrelation $\gamma(t)=\operatorname{cov}(X_{s},X_{s+t})$. The power spectral density (PSD) of $X_{t}$ is then defined as the Fourier transform of its autocorrelation,

[TABLE]

Spectral densities can be used to express the solutions of various differential equations and are thus commonly employed in many areas of physics. A particularly important example for AFM applications is that of a simple harmonic oscillator (SHO). This model for the thermally-driven tip position $X_{t}$ of the AFM cantilever is

[TABLE]

where $\dot{X}_{t}$ and $\ddot{X}_{t}$ are velocity and acceleration, $m$ is the tip mass, $k$ is the cantilever stiffness, $\varsigma$ is the viscous damping from the surrounding medium (e.g., air, water), and $F_{t}$ is the thermal force which drives the cantilever motion. It is a stationary white noise process with $E[F_{t}]=0$ and $\operatorname{cov}(F_{s},F_{s+t})=2k_{B}T\varsigma\cdot\delta(t)$, where $T$ is temperature and $k_{B}$ is Boltzmann’s constant. While the autocorrelation of $X_{t}$ has no simple form, a straightforward calculation in the Fourier domain obtains the spectral density

[TABLE]

where f_{0}=\mathchoice{{\hbox{\displaystyle\sqrt{k/m,}}\lower 0.4pt\hbox{\vrule height=7.5pt,depth=-6.00003pt}}}{{\hbox{\textstyle\sqrt{k/m,}}\lower 0.4pt\hbox{\vrule height=7.5pt,depth=-6.00003pt}}}{{\hbox{\scriptstyle\sqrt{k/m,}}\lower 0.4pt\hbox{\vrule height=5.25pt,depth=-4.20003pt}}}{{\hbox{\scriptscriptstyle\sqrt{k/m,}}\lower 0.4pt\hbox{\vrule height=3.75pt,depth=-3.00002pt}}}/(2\pi) is the cantilever’s resonance frequency and Q=\mathchoice{{\hbox{\displaystyle\sqrt{km,}}\lower 0.4pt\hbox{\vrule height=6.94444pt,depth=-5.55559pt}}}{{\hbox{\textstyle\sqrt{km,}}\lower 0.4pt\hbox{\vrule height=6.94444pt,depth=-5.55559pt}}}{{\hbox{\scriptstyle\sqrt{km,}}\lower 0.4pt\hbox{\vrule height=4.8611pt,depth=-3.8889pt}}}{{\hbox{\scriptscriptstyle\sqrt{km,}}\lower 0.4pt\hbox{\vrule height=3.47221pt,depth=-2.77779pt}}}/\varsigma is its “quality factor”, which measures the width of the PSD amplitude peak around $f_{0}$ (see Figure 2).

Remark 1**.**

The presentation above glosses over several technical details, e.g., non-integrability of many well-defined autocorrelations in (1) (such as that of the SHO), and limited interpretability of the white noise process $F_{t}$ in (2) as a function of $t$. For a rigorous treatment of these issues see Itō (1954).

2.1 Parametric Inference

In a parametric setting, the PSD is expressed as $\mathcal{S}(f,{\bm{\theta}})$ and the goal is to estimate the unknown parameters ${\bm{\theta}}$ from discrete observations ${\bm{X}}=(X_{0},\ldots,X_{N-1})$ recorded at sampling frequency $f_{s}$, such that $X_{n}=X(n\cdot{\Delta t})$ and $f_{s}=1/{\Delta t}$. Ideally one would work directly with the loglikelihood in the time domain, $\ell({\bm{\theta}}\,|\,{\bm{X}})$. However, this approach is inviable in practice since it (i) requires Fourier inversion of the PSD to obtain the variance of ${\bm{X}}$, and (ii) scales quadratically in the number of observations (e.g., Brockwell and Davis, 2013, Proposition 8.2.1). Instead, parametric inference can be considerably simplified by making use of the following result.

Proposition 1**.**

Let $N=2K+1$ and define the finite Fourier transform $\tilde{{\bm{X}}}=(\tilde{X}_{-K},\ldots,\tilde{X}_{0},\ldots,\tilde{X}_{K})$ of ${\bm{X}}$ as

[TABLE]

For each $\tilde{X}_{k}$, let $f_{k}=\tfrac{k}{N}f_{s}$ denote the corresponding frequency. Then if $\mathcal{S}(f)$ is the PSD of $X_{t}$, under suitable conditions on $\mathcal{S}(f)$, and as $N\to\infty$ and ${\Delta t}\to 0$, we have

[TABLE]

Proposition 1 leads to the so-called Whittle loglikelihood function (Whittle, 1957)

[TABLE]

where $Y_{k}=\tfrac{1}{N}\left|\tilde{X}_{k}\right|^{2}$ and $\mathcal{S}_{k}({\bm{\theta}})=f_{s}\cdot\mathcal{S}(f_{k},{\bm{\theta}})$. Since the periodogram ${\bm{Y}}=(Y_{1},\ldots,Y_{K})$ can be computed in $\mathcal{O}(N\log N)$ time using the Fast Fourier Transform, maximization of the Whittle loglikelihood (4) is considerably easier than of the original likelihood $\ell({\bm{\theta}}\,|\,{\bm{X}})$. Conditions for the convergence of the Whittle MLE $\hat{\bm{\theta}}_{W}=\operatorname*{arg\,max}_{\bm{\theta}}\ell_{W}({\bm{\theta}}\,|\,{\bm{Y}})$ to the true MLE $\hat{\bm{\theta}}=\operatorname*{arg\,max}_{\bm{\theta}}\ell({\bm{\theta}}\,|\,{\bm{X}})$ have been established by Fox and Taqqu (1986); Dahlhaus (1989). Since the true MLE is typically unavailable, we shall refer to Whittle’s $\hat{\bm{\theta}}_{W}$ simply as the MLE in the developments to follow.

2.2 Periodogram Binning

Despite its computational advantages relative to exact maximum likelihood, obtaining $\hat{\bm{\theta}}_{W}$ often remains a practical challenge, due to the enormous size of typical AFM datasets and the difficult numerical optimization of $\ell_{W}({\bm{\theta}}\,|\,{\bm{Y}})$. A common technique to overcome these issues is to group the periodogram frequencies into consecutive bins (e.g., Daniell, 1946; Brockwell and Davis, 2013, Section 10.4). That is, assume that $K=B\cdot N_{B}$ is a multiple of the bin size $B$, and consider the average periodogram value in bin $m$,

[TABLE]

It then follows from Proposition 1 that if $\mathcal{S}_{i}(f)$ is relatively constant within bins, the distribution of $\bar{\bm{Y}}=(\bar{Y}_{1},\ldots,\bar{Y}_{N_{B}})$ can be well approximated by

[TABLE]

where $\bar{\mathcal{S}}_{m}({\bm{\theta}})={f_{s}}\cdot\mathcal{S}(\bar{f}_{m},{\bm{\theta}})$, $\bar{f}_{m}=\tfrac{1}{B}\sum_{k\in I_{m}}f_{k}$, and $\textnormal{Gamma}(B,B)$ is a Gamma distribution with mean 1 and variance $1/B$. This leads to the non-linear least-squares (NLS) estimator

[TABLE]

which is a consistent estimator of ${\bm{\theta}}$ (Nørrelykke and Flyvbjerg, 2010). The sum-of-squares criterion (7) can be minimized using specialized algorithms such as Levenberg-Marquardt (Levenberg, 1944; Marquardt, 1963), rendering the calculation of $\hat{\bm{\theta}}_{\textnormal{NLS}}$ considerably simpler than that of $\hat{\bm{\theta}}_{W}$. However, this gain often incurs a significant loss in statistical precision.

3 Robust and Efficient Parametric PSD Inference

The choice between NLS and MLE estimators imposes a trade-off between computational and statistical efficiency. In addition, both estimators are highly sensitive to periodic noise which commonly plagues AFM spectral data (Section 3.2). Here we describe a two-stage parametric spectral estimator designed to overcome these issues.

3.1 Variance Stabilizing Transformation

To see why the NLS estimator is sub-optimally efficient, note that the approximate Gamma distribution of the binned periodogram (6) can itself be approximated by a Normal with matching mean and variance:

[TABLE]

Substituting any constant for the parameter-dependent variance in (8) then gives rise to $\hat{\bm{\theta}}_{\textnormal{NLS}}$ in (7). However, by a straightforward application of the statistical delta method (also known as propagation of errors, e.g., Bevington and Robinson, 2003), we note that taking the logarithm of the binned periodogram is a variance-stabilizing transformation:

[TABLE]

such that

[TABLE]

Maximizing the likelihood resulting from (9) leads to the log-periodogram (LP) estimator

[TABLE]

This simple sum-of-squares can be effectively minimized by the methods of Section 2.2, yet with $\hat{\bm{\theta}}_{\textnormal{LP}}$ achieving nearly the same precision as $\hat{\bm{\theta}}_{W}$. The LP estimator is commonly used in statistics to estimate long-range dependence (Geweke and Porter-Hudak, 1983; Robinson, 1995). Its asymptotic properties have been derived by Fay et al. (2002) and compared favorably therein to the efficient estimators of Taniguchi (1987).

Remark 2**.**

A different variance-stabilizing transformation of the periodogram ${\bm{Y}}$ is to take logarithms before binning, i.e., let $\tilde{Z}_{m}=\frac{1}{B}\sum_{k\in I_{m}}\log(Y_{k})$. Since the log-Exponential distribution is more Normal than the Exponential itself, a smaller $B$ is required for $\tilde{Z}_{m}$ than $Z_{m}$ for within-bin normality to hold. However, assuming that approximately $Y_{k}\stackrel{{\scriptstyle\operator@font{iid}}}{{\sim}}\bar{\mathcal{S}}_{m}({\bm{\theta}})\times\textnormal{Expo}(1)$ for $k\in I_{m}$, it can be shown that $\operatorname{var}(\tilde{Z}_{m})=\pi^{2}/6\times\frac{1}{B}>\operatorname{var}(Z_{m})$. Therefore, the analogous estimator to (10) with $\tilde{Z}_{m}$ in place of $Z_{m}$ is expected to be less efficient, and indeed this was the case in our numerical experiments.

3.2 Periodic Noise Removal

The PSD as defined in (1) tacitly assumes that the data are “purely stochastic” (in the sense of their Wold decomposition, e.g., Lindquist and Picci, 2015, Theorem 5.1.1). However, the periodogram in Figure 1 has several vertical lines in the ${10}^{5}{10}^{6}\text{\,}\mathrm{Hz}\text{/}$ range, suggesting the presence of periodic terms which cannot be explained by a PSD alone. Indeed, the AFM is a complex instrument operated by extensive electronics, which inevitably leads to periodic noise from various electrical components and power sources. Careful engineering can significantly reduce the effects of electronic noise on the final cantilever displacement readings. However, the residual periodic components shown in Figure 1 can gravely impact PSD parameter estimates as will be demonstrated shortly. Fortunately, the more severe electronic noise can be easily and automatically removed from the PSD by the following method due to Fisher (1929).

Suppose that the periodogram data ${\bm{Y}}=(Y_{1},\ldots,Y_{K})$ contain no periodic components. Under this null hypothesis, we have

[TABLE]

where ${\bm{\theta}}_{0}$ is the true parameter value. Now consider the maximum jump of the normalized cumulative periodogram density, also known as Fisher’s $g$-statistic:

[TABLE]

Under $H_{0}$, ${\mathcal{M}}$ is distributed as the maximum distance between the order statistics of $K$ iid uniform random variables, of which the distribution is given by (e.g., Brockwell and Davis, 2013, Corollary 10.2.2)

[TABLE]

3.3 Proposed Estimator

The developments above motivate a two-stage parametric spectral density estimator consisting of the following steps:

Preliminary Estimation. Calculate a preliminary estimate $\hat{\bm{\theta}}_{{\textnormal{LP}}}^{(1)}$ using the log-periodogram likelihood function (10). 2. 2.

Periodicity Removal. Calculate ${\mathcal{M}}$ upon substituting $\hat{\bm{\theta}}_{{\textnormal{LP}}}^{(1)}$ for the unknown value of ${\bm{\theta}}_{0}$, and the $p$-value against large ${\mathcal{M}}$ using (11). If the $p$-value is small – say less than 1% – replace the corresponding periodogram ordinate $Y_{k}$ by a random draw from $\mathcal{S}_{k}(\hat{\bm{\theta}}_{{\textnormal{LP}}}^{(1)})\times\textnormal{Expo}(1)$. Repeat this procedure until Fisher’s $g$-test does not reject $H_{0}$. 3. 3.

Final Estimation. Calculate $\hat{\bm{\theta}}_{{\textnormal{LP}}}$ on the periodogram obtained from Step 2, from which the unwanted periodicities have been removed.

Remark 3**.**

We have opted in Step 2 to replace the periodic outliers by random draws, instead of simply deleting them and repeating Fisher’s $g$-test with $K-1$ variables. This is because the largest of these $K-1$ variables is in fact the second largest of the original $K$, for which (11) does not give the right distribution under $H_{0}$.

4 Simulation Study

In order to evaluate the parametric spectral estimator proposed in Section 3.3, we consider the following simulation study reflecting a broad range of AFM calibration scenarios. Each simulation run consisted of a $5\text{\,}\mathrm{s}\text{/}$ time series sampled at $10\text{\,}\mathrm{MHz}\text{/}$ ($N=$5\text{\times}{10}^{6}$$ data points) from the SHO model (3) with added white noise,

[TABLE]

where ${\bm{\theta}}=(k,f_{0},Q,A_{\textrm{w}})$. Data was generated using a standard FFT-based algorithm (e.g., Labuda et al., 2012b). For all simulations, the baseline parameters are displayed in Table 1.

All parameters being fixed except $Q\in\{1,10,100,500\}$, the corresponding SHO spectra are displayed in Figure 2.

For each of the four baseline settings, $M=1000$ datasets were generated, and for each dataset we calculated the three estimators $\hat{\bm{\theta}}_{W}$, $\hat{\bm{\theta}}_{\textnormal{NLS}}$, and $\hat{\bm{\theta}}_{\textnormal{LP}}$. This was done using only periodogram frequencies in the range $f_{0}\pm f_{0}/\sqrt{2}$, a range typically provided by the cantilever manufacturer, and outside of which the remaining frequencies provide little additional information about ${\bm{\theta}}$. For the NLS and LP estimators, the bin size was set to $B=100$. For all estimators, the optimization was reduced from four to three parameters by the method of profile likelihood described in Appendix A.

4.1 Baseline Environment

Figure 3 displays boxplots for each estimator of each parameter estimate relative to its true value. The numbers on top of each boxplot correspond to the mean squared error (MSE) ratios between each estimator and the MLE $\hat{\bm{\theta}}_{W}$. That is, for each of the SHO parameters $\varphi\in(k,f_{0},Q)$ and estimator $j\in\{\textnormal{NLS},{\textnormal{LP}},\textnormal{MLE}\}$, the corresponding MSE ratio in Figure 3 is calculated as

[TABLE]

where $\varphi_{0}$ is the true parameter value and $\hat{\varphi}_{j}^{(i)}$ is its estimate by method $j$ for dataset $i$.

For low quality factor $Q=1$, the NLS method has roughly 1.5-2 times higher MSE than the MLE. For higher values of $Q$, the MSE of NLS increases to roughly 3-5 times that of MLE. In contrast, the LP estimator achieves virtually the same MSE as the MLE at a small fraction of the computational cost.

4.2 Electronic Noise Contamination

In order to assess the impact of electronic noise, a random sine wave of the form $D\cdot\sin(2\pi\zeta t+\phi)$ was added to each of the baseline datasets from the simulations above. The parameters of each sine wave were chosen to mimic the electronic noise in the real AFM data of Figure 1, a particularly difficult scenario for SHO parameter estimation due to the proximity of the electronic noise to the resonance frequency $f_{0}$. Specifically, the frequency $\zeta$ of each sine wave was generated from a Normal with mean $f_{0}=$33.533\text{,}\mathrm{kHz}\text{/}$$ and standard deviation $10\text{\,}\mathrm{Hz}\text{/}$, the phase $\phi$ was drawn uniformly between 0 and $2\pi$, and the amplitude $D$ was set to achieve an approximately ten-fold increase from the maximum value of the baseline PSD near $f_{0}$. The small jitter in the sine wave parameters was added both to mimic the small variations measured in real AFM data, and to investigate the impact of spectral leakage.

Figure 4 displays a simulated dataset with electronic noise contamination. Also displayed is the 1% threshold for periodic noise detection by Fisher’s $g$-test (Section 3.2). This is calculated by solving numerically for $\Pr({\mathcal{M}}>a_{\textnormal{cut}}\,|\,H_{0})=.01$ using (11), then setting the threshold for frequency $f_{k}$ to $a_{\textnormal{cut}}\times\mathcal{S}_{k}(\hat{\bm{\theta}}_{\textnormal{LP}}^{(1)})\cdot\sum_{j=1}^{K}(Y_{j}/\mathcal{S}_{j}(\hat{\bm{\theta}}_{\textnormal{LP}}^{(1)}))$. The threshold in Figure 4 indicates that most electronic noise detectable to the naked eye can be easily removed by the denoising procedure.

Figure 5 displays boxplots of each parameter estimate relative to its true value in the presence of electronic noise. To assess the impact of the noise corruption, these estimates do not include the denoising step of Section 3.3. The numbers in the plot correspond to MSE ratios between the estimator with noise corruption, relative to its own performance in the baseline dataset. The ratios are thus calculated as

[TABLE]

where $\hat{\varphi}_{j,{\textnormal{base}}}^{(i)}$ and $\hat{\varphi}_{j,{\textnormal{noise}}}^{(i)}$ are parameter estimates with method $j$ for dataset $i$ under baseline and noise-contaminated settings, respectively.

At low $Q$, the MSE ratios are close to one, indicating that the estimators are relatively insensitive to the electronic noise. However, for high $Q$ the effect of the noise is considerably more detrimental, particularly for NLS. In all cases, the performance of the LP estimator is affected the least, indicating it is naturally more robust than NLS and MLE to periodic noise contamination, even before the denoising technique is applied.

Figure 6 displays boxplots for the second-stage parameter estimates, after electronic noise removal. Each estimator (NLS, LP, and MLE) used its own preliminary fit to determine the noise cutoff value. Here, the MSE ratios are calculated relative to an “ideal” estimator: the MLE with perfect denoising. That is, the MSE ratios are

[TABLE]

where $\hat{\varphi}_{j,{\textnormal{corr}}}^{(i)}$ is the noise-corrected estimate of method $j$ for dataset $i$.

In general, the denoising procedure is extremely effective for both LP and MLE, but somewhat less so for NLS (for example, at $Q=100$ the MSE relative to $\hat{Q}_{W,{\textnormal{base}}}$ for $\hat{Q}_{\textnormal{NLS},{\textnormal{base}}}$ is 5.30, whereas for $\hat{Q}_{\textnormal{NLS},{\textnormal{corr}}}$ it is 8.90). However, for very high $Q=500$, the denoising procedure for LP and MLE fails to completely remove the upward bias in $\hat{Q}$ and the downward bias in $\hat{k}$. Upon further investigation, Figure 7 reveals that this is due to spectral leakage. Indeed, a close look at the 50 frequencies on either side of $f_{0}$ (Figure 7) shows that several periodogram variables adjacent to the electronic noise at $33.5490\text{\,}\mathrm{kHz}\text{/}$ have been pushed upward by its presence. The denoising procedure is able to remove the noise at $33.5490\text{\,}\mathrm{kHz}\text{/}$, but not in the neighboring frequencies. The net effect after noise correction is a slight upward bias in the binned periodogram (Figure 7), which, due to the high curvature of the SHO at $Q=500$, causes an upward bias in $\hat{Q}_{{\textnormal{corr}}}$. However, the overall amplitude of the SHO remains unaffected, and since by (12) this amplitude is proportional to $k_{B}T/(k\cdot\pi f_{0}Q)$, the upward bias in $\hat{Q}_{{\textnormal{corr}}}$ is accompanied by a downward bias in $\hat{k}_{{\textnormal{corr}}}$.

5 Application to Experimental AFM Data

We now turn to the problem of calibrating the AFM cantilever for which the periodogram is displayed in Figure 1. The data consist of $5\text{\,}\mathrm{s}\text{/}$ of an AC160 Olympus cantilever recorded at $5\text{\,}\mathrm{MHz}\text{/}$ ($N=$25\text{\times}{10}^{6}$$ observations). The objective is to determine the parameters of the best-fitting SHO model to the first cantilever eigenmode (Figure 1).

Calibration of a real AFM cantilever is subject to at least two complications not addressed in the simulations of Section 4:

While the PSDs used in simulation are dominated at low frequencies by white noise, those measured in the real data of Figure 1 exhibit power-law behavior, $\mathcal{S}(f)\sim 1/f^{\alpha}$ as $f\to 0$. This is referred to as “$1/f$ noise”; it features prominently in AFM experiments (e.g., Harkey and Kenny, 2000; Giessibl, 2003; Heerema et al., 2015), and is due in this case to slow fluctuations of the measurement sensor. Depending on the exponent, $1/f$ noise induces long-range dependence in the cantilever displacement ($0<\alpha<1$), or even lack of stationarity ($\alpha\geq 1$). Failing to account for it can significantly bias SHO parameter estimates. Fortunately, $1/f$ noise can be dealt with readily by adding a correction term to the SHO model, which becomes

[TABLE]

While important at low frequencies, the $1/f$ noise around the first eigenmode (Figure 1) is nearly imperceptible. Consequently, we estimated the first eigenmode’s SHO parameters using the simpler model (12). We have constructed a simulation in which $1/f$ noise severely affects SHO parameter estimation in Appendix B. Relative performance of LP to NLS and MLE estimators was similar to Section 4. 2. 2.

In addition to the first eigenmode at roughly $313\text{\,}\mathrm{kHz}\text{/}$, the data contain higher eigenmodes corresponding to flexural oscillations of the clamped cantilever beam (Sader, 1998). The first of these higher eigenmodes is displayed in Figure 1. Calibration of higher eigenmodes is of essential importance for popular bimodal and multifrequency AFM imaging techniques (e.g., Martinez et al., 2008; García and Proksch, 2013; Herruzo et al., 2014; Labuda et al., 2016b), on which we elaborate in the Discussion (Section 6).

Figure 8 displays the periodogram of the AFM data from Figure 1 over the frequency range used for parameter estimation. The electronic noise at $312.5\text{\,}\mathrm{kHz}\text{/}$ and $300.0\text{\,}\mathrm{kHz}\text{/}$ was easily removed with Fisher’s $g$-statistic. Table 2 displays parameter estimates and standard errors for NLS, LP, and MLE methods, the first two being calculated with bin size $B=100$. For LP and MLE, standard errors are calculated by inverting the observed Fisher information matrices corresponding to (10) and (4). For NLS, standard errors are obtained by the sandwich method (e.g., Freedman, 2006).

For this particular dataset, the NLS, LP, and MLE estimators are fairly similar, all being within one standard error of each other. This is because the difference between the estimators is largely driven by the relative amplitude of the SHO peak to its base. Here this ratio is about 10, which is similar to the $Q=10$ scenario examined in Section 4. Indeed, repeating the simulations of Section 4 with true parameters values taken as the MLE estimates in Table 2 produced similar results to the aformentioned scenario, i.e., indistinguishable LP and MLE estimators having three times smaller MSE than NLS.

5.1 Bin Size

While for this dataset there is little difference between the various estimators, NLS and LP can be substantially faster than MLE due to periodogram frequency binning. In practice, the choice of bin size affects both computational efficiency and approximation accuracy. Large bin sizes can group periodogram variables with very different amplitudes, thus invalidating the Gamma approximation to $\bar{Y}_{m}$ in (6). On the other hand, small bin sizes can strain the Normal approximations to $\bar{Y}_{m}$ and $\log(\bar{Y}_{m})$ in (8) and (9).

To investigate the impact of bin size, Figure 9 plots NLS and LP estimators for the values of B=\textrm{50250}. The behavior of NLS is considerably more erratic, presumably due to small changes in the bin end points having larger impact on $\bar{\mathcal{S}}_{m}({\bm{\theta}})$ than $\log\bar{\mathcal{S}}_{m}({\bm{\theta}})$. Note that the downward trend in $\hat{Q}$ is caused by increased flattening of the periodogram curvature as bin size increases.

6 Discussion

Parametric spectral density estimation plays a key role in AFM cantilever calibration. We have proposed a two-stage parameteric spectral estimator having statistical efficiency comparable to MLE at a small fraction of the computational cost, robust to most adverse effects of periodic noise contamination (except perhaps for very sharply peaked spectra). As spectral leakage due to binning affects the choice of bin size, a possible direction for future work is the construction of variable bin sizes, to be determined after the preliminary fit. Another line of future investigation is the calibration of higher eigenmodes. In principle, this can be done by fitting separate SHO models to each successive eigenmode. However, as the peak amplitude of these higher modes gets closer and closer to the noise floor, the accuracy of separate SHO estimators rapidly deteriorates. Instead one might wish to combine SHO models on the basis of hydrodynamic principles (Van Eysden and Sader, 2006; Clark et al., 2010) and other scaling laws (Labuda et al., 2016a).

Appendix A Profile Likelihood for SHO Fitting

We begin by reparametrizing the SHO model (12) as

[TABLE]

where $\tau=k_{B}T/(k\cdot\pi f_{0}Q)$, $R_{\textrm{w}}=A_{\textrm{w}}/\tau$, $\gamma=f_{0}Q$, and ${\bm{\eta}}=(f_{0},\gamma,R_{\textrm{w}})$. The objective function for the NLS estimator then becomes

[TABLE]

where $\bar{\mathcal{G}}_{m}({\bm{\eta}})={f_{s}}\cdot\mathcal{G}(\bar{f}_{m},{\bm{\eta}})$. For any fixed value of ${\bm{\eta}}$ the value of $\tau$ which minimizes $\mathcal{Q}_{\textnormal{NLS}}(\tau,{\bm{\eta}})$ is

[TABLE]

It follows that by setting $\hat{\bm{\eta}}_{\textnormal{NLS}}=\operatorname*{arg\,min}_{{\bm{\eta}}}\mathcal{Q}_{\textnormal{NLS}}(\hat{\tau}({\bm{\eta}}),{\bm{\eta}})$ and $\hat{\tau}_{\textnormal{NLS}}=\hat{\tau}(\hat{\bm{\eta}}_{\textnormal{NLS}})$, we have $(\hat{\tau}_{\textnormal{NLS}},\hat{\bm{\eta}}_{\textnormal{NLS}})=\operatorname*{arg\,min}_{(\tau,{\bm{\eta}})}\mathcal{Q}_{\textnormal{NLS}}(\tau,{\bm{\eta}})$. We can then recover the corresponding estimator $\hat{\bm{\theta}}_{\textnormal{NLS}}=\operatorname*{arg\,min}_{\bm{\theta}}\mathcal{Q}_{\textnormal{NLS}}({\bm{\theta}})$ by applying the inverse transformation $Q=\gamma/f_{0}$, $k=k_{B}T/(\tau\cdot\pi\gamma)$, and $A_{\textrm{w}}=R_{\textrm{w}}\cdot\tau$. Thus, we have obtained $\hat{\bm{\theta}}_{\textnormal{NLS}}$ at the cost of the three parameter optimization of $\mathcal{Q}_{\textnormal{NLS}}(\hat{\tau}({\bm{\eta}}),{\bm{\eta}})$, rather than the four parameter direct optimization of $\mathcal{Q}_{\textnormal{NLS}}({\bm{\theta}})$.

An analogous “profiling” procedure can be applied to the LP and MLE estimators, in order to reduce the optimization problem from four parameters to three. For LP, the objective function is

[TABLE]

for which

[TABLE]

Similarly, the objective function for the MLE estimator is

[TABLE]

where $\mathcal{G}_{k}({\bm{\eta}})=\mathcal{G}(f_{k}\,|\,{\bm{\eta}})$, for which

[TABLE]

Appendix B $\bm{1\,/f}$ Noise

The presence of $1/f$ noise is a common feature of AFM power spectra. This type of noise typically arises from slow fluctuations of the laser and photodiode sensor (Labuda et al., 2012a) and other long-term cantilever instabilities (Paolino and Bellon, 2009). It is manifested by a power law behavior at low frequencies, $\mathcal{S}(f)\sim 1/f^{\alpha}$ as $f\to 0$. Thus, a PSD model for the SHO with both white noise and $1/f$ noise contamination is

[TABLE]

where $\alpha$ and $A_{\textrm{f}}$ are the $1/f$ noise exponent and amplitude parameters, respectively. While the SHO estimates for the real AFM data in Figure 1 were not impacted by the $1/f$ noise, here we construct a simulation study in which they are. Namely, we use the baseline parameters described in Section 4, to which we add $1/f$ noise with parameters $\alpha=0.55$ and $A_{\textrm{f}}=$1.0\text{\times}{10}^{7}\text{,}{\mathrm{fm}}^{2}\text{/}\mathrm{Hz}$$. Baseline and noise contaminated power spectra are displayed in Figure 10.

To quantify the severity of the $1/f$ noise, Table 3 displays the asymptotic relative bias (i.e., as $N\to\infty$) due to fitting the SHO + white noise model (12) without accounting for the $1/f$ noise in Figure 10. This was calculated by a direct curve-fitting procedure.

While the bias on $f_{0}$ and $k$ is relatively small, for $Q$ it is on the order of $510$%.

To evaluate the different estimators, $M=1000$ datasets are generated under each setting as in Section 4, and NLS, LP, and MLE parameter estimates are calculated for each dataset. For the NLS and LP estimators the bin size was $B=100$. Table 4 displays the parameter-wise MSE ratio for NLS and LP estimators relative to the MLE. For moderate $Q\geq 10$, the performance of the LP estimator is virtually the same as the MLE, and $1.55$ times superior than that of NLS. For very low $Q=1$, the $1/f$ noise in Figure 10 is almost undetectable, leading to parameter identifiability issues in the fitting algorithms. In such a setting we recommend to first estimate the $1/f$ parameters separately from the low frequency periodogram values, then estimate the SHO and white noise parameters with $\hat{\alpha}$ and $\hat{A}_{\textrm{f}}$ fixed.

Supplementary Materials

Software: All code for the various PSD fitting algorithms is available at

https://github.com/mlysy/realSHO.