# Data-driven Computing in Elasticity
via Chebyshev Approximation

Rahul-Vigneswaran K*† *Website : rahulvigneswaran.github.io†Work done while
interning at CEN Department of Mechanical Engineering
Amrita Vishwa Vidyapeetham
Amritapuri, India
<EMAIL_ADDRESS>Neethu Mohan, Soman KP  Center for Computational
Engineering and Networking (CEN)
Amrita School of Engineering, Coimbatore
Amrita Vishwa Vidyapeetham, India
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

This paper proposes a data-driven approach for computing elasticity by means
of a non-parametric regression approach rather than an optimization approach.
The Chebyshev approximation is utilized for tackling the material data-sets
non-linearity of the elasticity. Also, additional efforts have been taken to
compare the results with several other state-of-the-art methodologies.

###### Index Terms:

Data-driven computational mechanics, Model-free method, Nonparametric method,
Chebyshev polynomials, elasticity, Chebyshev approximation, chebfun

##  I Introduction

In the era where machine learning is catapulting the conventional
methodologies into the future, computational mechanics is a field where there
was not much of development until recently. Data-driven solvers have become an
emerging field lately due to its obvious advantages over conventional
heuristic solvers. Data-driven solvers attempt to utilize the material data-
set, without empirically modelling a constitutive law based on the stress-
strain relation of the material, upon which the conventional methodologies of
computational mechanics were heavily dependent on. Conventional approaches get
the job done but are not good enough because of the unpredictable mapping
between the equations and the real world behaviours of the systems. One cannot
convert every possible real-time scenario into an appropriate equation. This
calls for the need of Data-driven solvers that function based on real-time
data rather than a set of heuristically created equations.

Elastic systems like shown in Fig. 1 (a), are easy to solve due to its linear
nature by using Hooke’s law [1] while a system like shown in Fig. 1 (b) are
inherently difficult to solve due to its non-linearity of stress-strain
relationship even at infinitesimally small strains inputs. There are several
models used for solving nonlinear elasticity namely Ramberg-Osgood [2],
Hyperbolic law [3], Hardin-Drnevich [4], Duncan-Chang [5], Duncan-Selig [6]
and the one that is commonly used is called Power law or also called as
Ludwik’s law [7]. Check [19] for an extensive review. Even though the above-
mentioned model capture the essence of nonlinear elasticity of the system to a
certain degree, they are still a borderline approximation of the actual non-
linearity. This is where date driven solvers like [8, 9, 10, 11, 12, 13, 14,
15, 16, 17, 18] become superior over the conventional methods mentioned
earlier.

Section I provides a basic introduction about the research work, Section II
explores and gives a detailed analysis of already existing literature and
Section III gives a background knowledge of topics such as nonlinear
elasticity and Chebyshev polynomials for a better understanding of the work.
Section IV provides an in-depth knowledge of the experimental setup and the
methodologies utilized. Finally Section V contains all the results obtained
from the experiments and Section VI gives a strong closing statement by
pressing on the importance of this work, the significance of the results
obtained and future scope

##  II Related Works

Several methods like [8, 9] rely on minimizing the distance between a provided
data set of the material and the subset of the appropriate stress-strain
fields in equilibrium conditions. [8] considers the Mahalanobis distance [28]
in order to incorporate the statistical uncertainty (2nd order) of data in
computations. It is the distance between two groups’ multiple means
(centroids) used in the analysis of discriminant. In [18], the method of
minimizing the distance is extrapolated for nonlinearity in elasticity and
Lagrange multipliers enforce the physical constraints. Data-driven solvers
like [10], irrespective of the ideology of following a pure approach guided
the given data in its entirety, heuristically adds a few constraints in order
to develop a solver that is thermodynamically consistent during its solving
phase. This directly contributes to fulfilling the second and first principles
of thermodynamics which is independent of the experimental result’s quality.
While most solvers lack generality due to its nature of the mathematical
formulation, they fail at adapting themselves to new experimental data. This
poses a real issue in real-time embedded applications as the solvers have to
be remodelled periodically. [11] overcomes this downfall by directly
connecting the data to computers for performing simulations. The so-performed
numerical simulations will employ rules based on universal laws while
minimizing the need for specifically engineered models by making use of
manifold learning [29] methodologies.

##  III Background

###  III-A Nonlinear Elasticity

In order to understand nonlinear elasticity in detail, it becomes vital to
have a basic understanding of elasticity in general. The key pillars for a
material to be elastic in nature are – when the application of the load is
stopped, the deformation vanishes, the deformation is independent of its
history and the potential-energy exist as a function of deformation. When a
material satisfies the above criteria for elasticity and its stress-strain
relation is nonlinear like in Fig. 1 (b), then they are called as nonlinear
elastic materials.

Figure 1: Stress-strain relation graphs of (a) Linearly Elastic and (b) Non-
linearly elastic systems.

This nonlinearity between the stress and strain is exactly the hurdle that
prevents conventional methodologies from functioning at their fullest and
makes it hard to map to a set of predefined equations.

###  III-B Chebyshev Polynomials

Chebyshev polynomials play a vital role in approximation theory [20] due to
its roots of the Chebyshev polynomials of the first kind and are also known as
Chebyshev nodes [21]. During the polynomial interpolation, these nodes are
used. As a result, the resulting polynomial of interpolation reduces the
problem of Runge’s phenomenon. It also delivers an approximation that is
extremely near to the optimum approximation polynomial to a function that is
continuous under the maximum norm. This approximation inherently leads to the
Clenshaw–Curtis quadrature method [22].

Figure 2: Illustration of The Chebyshev polynomials of the first kind
${T_{n}}\left(x\right)$ for $x\in\left[{0,1}\right]$ and $n=1,2,\ldots,5$.

The Chebyshev polynomials of the first kind ${T_{n}}\left(x\right)$ are
illustrated in Fig. 2 for $x\in\left[{0,1}\right]$ and $n=1,2,\ldots,5$. Using
a Chebyshev Polynomial of the First Kind $T$, we can define

$\displaystyle{{c}_{j}}$
$\displaystyle\equiv\frac{2}{N}\sum\limits_{k=1}^{N}{f\left({{x}_{k}}\right)}{{T}_{j}}\left({{x}_{k}}\right)$
(1)
$\displaystyle=\frac{2}{N}\sum\limits_{k=1}^{N}{f\left[\cos\left\\{\frac{\pi\left(k-\frac{1}{2}\right)}{N}\right\\}\right]}\cos\left\\{\frac{\pi
j\left(k-\frac{1}{2}\right)}{N}\right\\}$ (1)

where, $k=1,2,\ldots,n$. Then,

$f\left(x\right)\approx\sum\limits_{k=0}^{N-1}{{{c}_{k}}{{T}_{k}}\left(x\right)}-\frac{1}{2}{{c}_{0}}$

For the $N$ zeros of ${T_{N}}\left(x\right)$ it is exact and this way of
approximation is essential because, in case of truncation, the error is
distributed evenly and smoothly over $\left[{-1,1}\right]$. The Chebyshev
approximation [23, 30] is very near to the Minimax Polynomial.

##  IV Experiments

###  IV-A Experimental Setup

We use MATLAB as the test platform and take advantage of chebfun [27] package
for carrying out calculations of Chebyshev approximation. For exponentially
increasing the agility of processing of data, a computer with Intel Core
i7-6700HQ Skylake CPU running at 2.60Hz and a DDR4 16 GB ram memory. The truss
configuration used for the results shown in this paper is illustrated in Fig.
6 but the methodology can be extended to any truss configuration. GitHub
Repository for all the codes used in this paper:
https://github.com/rahulvigneswaran/Data-Driven-Computing-in-Elasticity-via-
Chebyshev-Approximation

(a)

(b)

Figure 3: Obtained paths of equilibrium by varying load multiplier: (a) Method
proposed in this paper; (b) Method proposed in [26].

(a)

(b)

Figure 4: Solutions obtained for the load multiplier $\lambda=10$: (a) Method
proposed in this paper; (b) Method proposed in [26].

(a)

(b)

Figure 5: Obtained stress-strain pairs of member (A). (a) The proposed method
and (b) Method from [24] applied to the formulation in [14]. $\bigcirc\to$
material data point, $\Delta\to$ obtained strain-stress pair.

Figure 6: 10-bar truss configuration [26]

###  IV-B Methodology

Even though this methodology can be extended to any nonlinear elasticity
problem unlike [24], for the purpose of simplicity, we have limited our
explanation to the scenario of truss undergoing small deformations. This
problem when formulated like in [15], inherently becomes a nonconvex
optimization [25] problem. In this paper, we have utilized the method used in
[26], which proposes a novel method by converting the problem into a set of
nonlinear equations that can be solved easily.

The compatibility conditions and the force balance relations are as follows,

$\begin{gathered}{\varepsilon_{i}}=b_{i}^{T}u\\\
\sum\limits_{i=1}^{m}{{\upsilon_{i}}{\sigma_{i}}{b_{i}}}=p\\\ \end{gathered}$

where, ${\varepsilon_{i}}$ is axial strain, ${{\sigma_{i}}}$ is axial stress
and $\upsilon_{i}$ is the volume of the member $i(i=1,2,\ldots,m)$
respectively. In equation (IV-B), $b_{1},b_{2},\ldots,b_{m}$ are constant
vectors. Lets suppose the material dataset has been provided and is
represented as $D=\left\\{{\left({{{\overset{\lower
5.0pt\hbox{$\smash{\scriptscriptstyle\smile}$}}{\varepsilon}}_{1}},{{\overset{\lower
5.0pt\hbox{$\smash{\scriptscriptstyle\smile}$}}{\sigma}}_{1}}}\right),\ldots,\left({{{\overset{\lower
5.0pt\hbox{$\smash{\scriptscriptstyle\smile}$}}{\varepsilon}}_{d}},{{\overset{\lower
5.0pt\hbox{$\smash{\scriptscriptstyle\smile}$}}{\sigma}}_{d}}}\right)}\right\\}$
where, the observed uniaxial strains and stresses are denoted by
${{{\overset{\lower
5.0pt\hbox{$\smash{\scriptscriptstyle\smile}$}}{\varepsilon}}_{j}}}$ and
${{{\overset{\lower
5.0pt\hbox{$\smash{\scriptscriptstyle\smile}$}}{\sigma}}_{j}}}$ respectively.
For each member $i$, the stress estimated $(S_{i})$, for a given
${{\varepsilon_{i}}}$ is calculated by using the kernel regression [31] in
[26] using the Gaussian kernel of the following form,

${s_{i}}=\frac{{\sum\limits_{j=1}^{d}{\exp\left[{-\alpha{{\left({{\varepsilon_{i}}-{{\overset{\lower
5.0pt\hbox{$\smash{\scriptscriptstyle\smile}$}}{\varepsilon}}_{j}}}\right)}^{2}}}\right]}{{\overset{\lower
5.0pt\hbox{$\smash{\scriptscriptstyle\smile}$}}{\sigma}}_{j}}}}{{\sum\limits_{j=1}^{d}{\exp\left[{-\alpha{{\left({{\varepsilon_{i}}-{{\overset{\lower
5.0pt\hbox{$\smash{\scriptscriptstyle\smile}$}}{\varepsilon}}_{j}}}\right)}^{2}}}\right]}}}$

where, $\alpha>0$ is acquired by cross validation of D. A data-driven solving
approach can be formulated as,

$minimize\left\|{\sigma-s}\right\|$

with constraints (IV-B) and (IV-B). Here, equation (IV-B) in conjunction with
the given material data, is basically used as a lookup table but the drawback
is that the equation (IV-B) has to be calculated every time the material data
has to be referred by the algorithm.

In this paper, we have experimented by replacing the equation (IV-B), that is,
the Kernel regression with Chebyshev approximation (Equation (III-B)). The
advantage is that, by doing so, we create a one-time smooth data curve from
the given coarse experimental data which can be accessed easily like simple
function $f(x)$ where $x$ is the input strain and $f(x)$ will provide the
corresponding stress output.

##  V Results

An extensive comparison has been carried out for better assessment of the
results. Fig. 7 is the stress-strain relation of the material that is
considered for the assessment. It is evident from it that the dataset is
coarse in nature. The obtained equilibrium paths of the method proposed in
[26] and this paper are depicted in Fig. 3. The solutions obtained for the
load multiplier $\lambda=10$ are shown in Fig. 4. The obtained stress-strain
pairs of member (A) in the truss (Fig. 6) are shown in Fig. 5.

Table I: Computational Costs Formulation | Solving Process | Time (s)
---|---|---
Proposed Method | Chebyshev Approx. | 0.3
| Nonlinear Eqn. | 3.2
| Total | 3.5
Kernal Regression [26] | Cross Validation | 2.3
| Nonlinear Eqn. | 2.64
| Total | 7.18
Formulation in [14] | MIP [24] | 309.4
Formulation in [14] | Heuristic [14] | 0.6
Formulation in [15] | NLP | 13.8

Table I compares the computation time of the proposed method with methods
proposed in [14, 15, 24, 16]. It is evident from the result that the proposed
method has significantly outperformed the previously proposed state-of-the-art
model by $51.3\%$. Even though the model proposed in [14] has a very low
computation time compared to other models, it is heuristically designed. It
does not necessarily yield an optimal solution. For example, for a load
multiplier of 10, the value of the solution is 1.316 J, while the actual
solution is 0.251 J. Thus, having 81% deviation from the actual solution,
proving that the heuristic converged an inexact solution, irrespective of that
fact that it is computationally cheap. As a result, the method proposed in
this paper has a significant advantage over the existing models, making it
suitable for the application of nonlinear elasticity. A much more extensive
analysis with varying load multipliers have been done and can be viewed at
https://rahulvigneswaran.github.io/Data-Driven-Computing-in-Elasticity-via-
Chebyshev-Approximation/.

Figure 7: Material dataset used for the experiments.

##  VI Conclusion

The same experiment was carried out several other algorithms including neural
network based curve-fitting which showed significantly faster results compared
to Chebyshev approximation. Yet, they were not considered as suitable
candidates due to its inherent nature that requires it to fit a curve, which
might become time consuming when the distribution of the data is sparse and
the curve that is fit becomes increasingly erroneous and might not necessarily
mimic the behaviour of the actual system while Chebyshev approximation can
handle it with ease.

From the results, it becomes evident that the proposed methodology is
significantly faster than the existing cutting edge models while preserving
the generality approach in elasticity. As we can define the linear elasticity
as a special subset of nonlinear elasticity, the proposed method holds good
for both linear and nonlinear elastic cases, making it a universal solver for
elastic bodies. Deep learning [34] has shown its effectiveness in many of the
emerging fields ([32, 33]) which were traditionally completely based on
heuristic methods. Therefore, using Deep learning in conjunction with the
methodology proposed in this paper can be considered as a future endeavour.

## Acknowledgment

The authors would like to show their acknowledgement to Dr Soman KP, head of
CEN lab for exchanging his wisdom and providing with continuous support during
the course of this research. Also, an immense gratitude is conveyed to the
anonymous reviewers for their valuable insights.

## References

  * [1] Rychlewski, J., 1984. On Hooke’s law. Journal of Applied Mathematics and Mechanics, 48(3), pp.303-314.
  * [2] Ramberg, W. and Osgood, W.R., 1943. Description of stress-strain curves by three parameters.
  * [3] Keyfitz, B.L. and Kranzer, H.C., 1980. A system of non-strictly hyperbolic conservation laws arising in elasticity theory. Archive for Rational Mechanics and Analysis, 72(3), pp.219-241.
  * [4] Fahey, M. and Carter, J.P., 1993. A finite element study of the pressuremeter test in sand using a nonlinear elastic plastic model. Canadian Geotechnical Journal, 30(2), pp.348-362.
  * [5] Huang, B., Bathurst, R.J. and Hatami, K., 2009. Numerical study of reinforced soil segmental walls using three different constitutive soil models. Journal of Geotechnical and Geoenvironmental engineering, 135(10), pp.1486-1498.
  * [6] Katona, M.G., 2015. Influence of soil models on performance of buried culverts (No. 15-0560).
  * [7] Isselin, J., Iost, A., Golek, J., Najjar, D. and Bigerelle, M., 2006. Assessment of the constitutive law by inverse methodology: Small punch test and hardness. Journal of Nuclear Materials, 352(1-3), pp.97-106.
  * [8] Ayensa-Jiménez, J., Doweidar, M.H., Sanz-Herrera, J.A. and Doblaré, M., 2018. A new reliability-based data-driven approach for noisy experimental data with physical constraints. Computer Methods in Applied Mechanics and Engineering, 328, pp.752-774.
  * [9] Conti, S., Müller, S. and Ortiz, M., 2018. Data-driven problems in elasticity. Archive for Rational Mechanics and Analysis, pp.1-45.
  * [10] González, D., Chinesta, F. and Cueto, E., 2019. Thermodynamically consistent data-driven computational mechanics. Continuum Mechanics and Thermodynamics, 31(1), pp.239-253.
  * [11] Ibanez, R., Abisset-Chavanne, E., Aguado, J.V., Gonzalez, D., Cueto, E. and Chinesta, F., 2018. A manifold learning approach to data-driven computational elasticity and inelasticity. Archives of Computational Methods in Engineering, 25(1), pp.47-57.
  * [12] Ibañez, R., Borzacchiello, D., Aguado, J.V., Abisset-Chavanne, E., Cueto, E., Ladevèze, P. and Chinesta, F., 2017. Data-driven non-linear elasticity: constitutive manifold construction and problem discretization. Computational Mechanics, 60(5), pp.813-826.
  * [13] Kanno, Y., 2018. Simple heuristic for data-driven computational elasticity with material data involving noise and outliers: a local robust regression approach. Japan Journal of Industrial and Applied Mathematics, 35(3), pp.1085-1101.
  * [14] Kirchdoerfer, T. and Ortiz, M., 2016. Data-driven computational mechanics. Computer Methods in Applied Mechanics and Engineering, 304, pp.81-101.
  * [15] Kirchdoerfer, T. and Ortiz, M., 2017. Data driven computing with noisy material data sets. Computer Methods in Applied Mechanics and Engineering, 326, pp.622-641.
  * [16] Latorre, M. and Montáns, F.J., 2018. Experimental data reduction for hyperelasticity. Computers & Structures.
  * [17] Leygue, A., Coret, M., Réthoré, J., Stainier, L. and Verron, E., 2018. Data-based derivation of material response. Computer Methods in Applied Mechanics and Engineering, 331, pp.184-196.
  * [18] Nguyen, L.T.K. and Keip, M.A., 2018. A data-driven approach to nonlinear elasticity. Computers & Structures, 194, pp.97-115.
  * [19] Zeidler, E., 1988. Basic Equations of Nonlinear Elasticity Theory. In Nonlinear Functional Analysis and its Applications (pp. 158-232). Springer, New York, NY.
  * [20] Butzer, P.L. and Nessel, R.J., 1971. Fourier analysis and approximation, Vol. 1. Reviews in Group Representation Theory, Part A (Pure and Applied Mathematics Series, Vol. 7).
  * [21] Adibi, H. and Assari, P., 2010. Chebyshev wavelet method for numerical solution of Fredholm integral equations of the first kind. Mathematical problems in Engineering, 2010.
  * [22] Gentleman, W.M., 1972. Implementing Clenshaw-Curtis quadrature, II computing the cosine transformation. Communications of the ACM, 15(5), pp.343-346.
  * [23] Ames, W.F. and Brezinski, C., 1993. Numerical recipes in Fortran (The art of scientific computing): WH Press, SA Teukolsky, WT Vetterling and BP Flannery, Cambridge Univ. Press, Cambridge, 1992. 963 pp., US $49.95, ISBN 0-521-43064-X.
  * [24] Kanno, Y., 2018. Mixed-integer programming formulation of a data-driven solver in computational elasticity. Optimization Letters, pp.1-10.
  * [25] Jain, P. and Kar, P., 2017. Non-convex optimization for machine learning. Foundations and Trends® in Machine Learning, 10(3-4), pp.142-336.
  * [26] Kanno, Y., 2018. Data-driven computing in elasticity via kernel regression. Theoretical and Applied Mechanics Letters, 8(6), pp.361-365.
  * [27] Driscoll, T.A., Hale, N. and Trefethen, L.N., 2014. Chebfun guide.
  * [28] De Maesschalck, R., Jouan-Rimbaud, D. and Massart, D.L., 2000. The mahalanobis distance. Chemometrics and intelligent laboratory systems, 50(1), pp.1-18.
  * [29] Wang, J., Zhang, Z. and Zha, H., 2005. Adaptive manifold learning. In Advances in neural information processing systems (pp. 1473-1480).
  * [30] Mohan, N. and Soman, K.P., 2018, February. Power system frequency and amplitude estimation using variational mode decomposition and chebfun approximation system. In 2018 Twenty Fourth National Conference on Communications (NCC)(pp. 1-6). IEEE.
  * [31] Bishop, C.M., 2006. Pattern recognition and machine learning. springer.
  * [32] Rahul-Vigneswaran K, Vinayakumar R, Soman KP and Prabaharan Poornachandran. Evaluating Shallow and Deep Neural Networks for Network Intrusion Detection Systems in Cyber Security. In 2018 9th International Conference on Computing, Communication and Networking Technologies (ICCCNT) (pp. 1-6). IEEE.
  * [33] Rahul-Vigneswaran K, Prabhaharan Poornachandran and Soman KP. A Compendium on Network and Host Based Intrusion Detection Systems. In 2019 International Conference on Data Science, Machine learning & Applications (ICDSMLA) (pp. 1-8). Springer
  * [34] LeCun, Y., Bengio, Y. and Hinton, G., 2015. Deep learning. nature, 521(7553), p.436.


# Precise sinusoidal signal extraction from noisy waveform in vibration
calibration

Tomofumi Shimoda, Wataru Kokuyama, Hideaki Nozato  National Metrology
Institute of Japan, National Institute of Advanced Industrial Science and
Technology, 1-1-1 Umezono, Tsukuba, Ibaraki 305-8563, Japan
<EMAIL_ADDRESS>
###### Abstract

Precise extraction of sinusoidal vibration parameters is essential for the
dynamic calibration of vibration sensors, such as accelerometers. However,
several standard methods have not yet been optimized for large background
noise. In this work, signal processing methods to extract small vibration
signals from noisy data in the case of accelerometer calibration is discussed.
The results show that spectral leakage degrades calibration accuracy. Three
methods based on the use of a filter, window function, and numerical
differentiation are investigated with theoretical calculations, simulations,
and experiments. These methods can effectively reduce the contribution of the
calibration system noise. The uncertainty of micro vibration calibration in
the National Metrology Institute of Japan is reduced by two orders of
magnitudes using the proposed methods. The theoretical analyses in this work
can lay the foundation for the optimization of signal processing in vibration
calibration, and can be applied to other dynamic calibration fields.

##  1 Introduction

Micro vibration measurement is required in various fields, such as
infrastructure health monitoring [1] or satellite performance analysis [2].
Many types of sensors have been developed and used, such as high-sensitivity
accelerometers, broadband seismometers, or low-noise microelectromechanical
system (MEMS) accelerometers (e.g., [3, 4, 5]). The calibration of the sensor
frequency response is essential for the reliability of micro vibration
measurements. As the applications of micro vibration measurement increase,
measuring responses to small input vibrations is becoming increasingly
important.

Accurate extraction of amplitude and phase from a sinusoidal waveform is
required for the calibration of accelerometer sensitivity. The target
accelerometer is sinusoidally vibrated by a vibration exciter, and the
amplitudes and phases of the sensor voltage signal and reference displacement
signal are estimated and compared to calibrate the sensitivity and phase shift
of the accelerometer. The estiation accuracy is essential for the calibration
uncertainty. The requirement for the primary calibration is about 0.1 % and
0.1∘ for the amplitude and phase, respectively [6, 7]. In these studies, large
vibration that was nearly $10$ m/s2 was applied. To obtain the response to a
micro vibration with an amplitude down to about $10^{-3}$ m/s2, the applied
vibration amplitude should also be small (on the same order of magnitude)
because the linearity of the response is not ensured in general. The
extraction of such a small vibration signal usually suffers from the
background noise of the calibration system, which originates from the
background vibration or electrical noise.

For accelerometer calibration, amplitude and phase extraction are performed by
the sine approximation method (SAM) in ISO16063-11 [8]. However, the SAM is
not optimized for real cases with large background noise. In addition to the
background noise reduction of the calibration system, the optimization of
signal processing is important. Some of the calibration institutes empirically
apply a digital filter to deal with the problem, although an unified approach
has not yet been established. Another signal processing method using
correlation has also been proposed [9]. In this paper, we discuss the
limitations of the SAM and its optimization through theoretical
investigations, simulations, and experiments. These investigations are
important to calibrate accelerometers under a large noise or with a small
vibration amplitude. The results can be applied to not only accelerometer
calibration but also other signal processing settings that require accurate
sinusoidal parameter estimation.

The remainder of this paper is organized as follows: Section 2 summarizes the
mathematical background of the SAM and its verification with simulations.
Section 3 proposes optimization methods, including filtering, changing the
window function, and numerical differentiation. Section 4 applies the proposed
method to the vibration calibration in the National Metrology Institute of
Japan (NMIJ).

##  2 Contribution of system noise to sinusoidal signal extraction

###  2.1 Mathematical framework of the conventional SAM

Figure 1: Primary calibration system of an accelerometer.

The overview of the accelerometer calibration system is depicted in Fig. 1.
The accelerometer is sinusoidally vibrated by the vibration exciter, and its
output signal $V_{\mathrm{s}}$ is recorded along with the reference
displacement signal $x_{\mathrm{r}}$ measured by the laser interferometer.
Throughout this paper, signal processing is assumed to be digital signal
processing using discretely sampled data. The SAM specified in ISO16063-11 [8]
determines the sensitivity modulus and phase shift from two waveforms. The
amplitude and phase of each waveform at vibration frequency $f_{\mathrm{v}}$
are extracted in the following process. Here, a waveform $x(t)$ is sampled
within $0<t<T$ at $t_{n}=nT/N$ ($n=$0, 1, …, $N-1$) ($N$: the number of data).
The recorded data $x_{n}$ is modeled as

$x_{n}=b_{0}+b_{1}\cos(2\pi f_{\mathrm{v}}t_{n})+b_{2}\sin(2\pi
f_{\mathrm{v}}t_{n})+\epsilon_{n},$ (1)

where $b_{0}$, $b_{1}$, and $b_{2}$ are the fitting parameters, and
$\epsilon_{n}$ is the residual from the model. In this work, we assume that
the sampling is sufficiently faster than the vibration frequency, $N/T\gg
f_{\mathrm{v}}$, and the data length is an integer multiple of the vibration
period, $T=N_{\mathrm{c}}/f_{\mathrm{v}}$ ($N_{\mathrm{c}}$: integer). These
conditions are easily achieved using a commercially available digitizer and by
setting the proper record length. Under these conditions, the optimal
parameters that minimize the residual $\sum\epsilon_{n}^{2}$ are given by

$b_{1}=\frac{2}{N}\sum_{n=0}^{N-1}x_{n}\cos(2\pi f_{\mathrm{v}}t_{n}),\ \
b_{2}=\frac{2}{N}\sum_{n=0}^{N-1}x_{n}\sin(2\pi f_{\mathrm{v}}t_{n}),$ (2)

which are identical to the Fourier transform of the data. Here, the estimated
complex amplitude for the variable $x$ is defined as
$\hat{x}_{\mathrm{est}}\equiv b_{1}+ib_{2}$ for convenience. The amplitudes
are computed for both the sensor signal $V_{\mathrm{s}}$ and the reference
signal $x_{\mathrm{r}}$ to calculate the sensitivity modulus and phase delay
as

$S_{\mathrm{cal}}=\frac{|\hat{V}_{\mathrm{s,est}}|}{(2\pi
f_{v})^{2}|\hat{x}_{\mathrm{r,est}}|},\hskip
20.0pt\Delta\phi_{\mathrm{cal}}=\arg(\hat{V}_{\mathrm{s,est}})-\arg(\hat{x}_{\mathrm{r,est}})-\pi,$
(3)

respectively.

Using Eq. (2), the complex amplitude $\hat{x}_{\mathrm{est}}$ can be modified
as

$\displaystyle\hat{x}_{\mathrm{est}}(t_{0})$ $\displaystyle=$
$\displaystyle\frac{2}{N}\sum_{n=0}^{N-1}x_{n}e^{2\pi if_{\mathrm{v}}t_{n}}$
(4) $\displaystyle\simeq$ $\displaystyle\frac{2}{T}\int_{0}^{T}x(t)e^{2\pi
if_{\mathrm{v}}t}dt$ $\displaystyle=$
$\displaystyle\frac{2}{T}\int_{-\infty}^{\infty}w(t+t_{0})x(t)e^{2\pi
if_{\mathrm{v}}t}dt$ $\displaystyle=$
$\displaystyle\frac{2}{T}\int_{-\infty}^{\infty}\tilde{W}(f)\tilde{X}(f_{\mathrm{v}}-f)e^{-2\pi
ift_{0}}df$

Here, $t_{0}$ is the start time of the measurement, and

$w(t)=w_{\mathrm{r}}(t)=\left\\{\begin{array}[]{@{\kern 2.5pt}lL}1\ (0<t<T)\\\
0\ (\rm{otherwise})\end{array}\right.$ (5)

is a rectangular window function. $\tilde{W}(f)$ and $\tilde{X}(f)$ are the
Fourier transform of $w(t)$ and $x(t)$, respectively. In summary, the
conventional SAM is the Fourier transformation with the rectangular window
function.

As Eq. (4) indicates, the parameter estimation at $f_{\mathrm{v}}$ is affected
by the frequency component $X(f)\ (f\neq f_{\mathrm{v}})$, which is known as
the spectral leakage. For the rectangular window, the leakage is determined by

$\tilde{W}_{\mathrm{r}}(f)=T\mathrm{sinc}(fT)e^{\pi ifT},$ (6)

where sinc function is defined as $\mathrm{sinc}(z)=\sin(\pi z)/\pi z$.
Although the amplitude of the purely sinusoidal vibration
$x_{0}(t)=\hat{x}_{0}\sin(2\pi f_{\mathrm{v}}t)$ is calculated to be
$\hat{x}_{0\mathrm{,est}}(f_{\mathrm{v}})=i\hat{x}_{0}$, the background noise
at $f_{\mathrm{v}}$ or the leakage from $f\neq f_{\mathrm{v}}$ can be a
problem in real-world cases.

###  2.2 SAM under noise

Figure 2: Model of the signal flow in accelerometer calibration.

Fig. 2 shows the signal flow model for accelerometer calibration shown in Fig.
1. A similar discussion is applicable to other experiments conducted to
measure the amplitude ratio between two vibration timeseries. The error
sources considered in Fig. 2 are divided into four types

  * •

independent random background noise of each device ($n_{\mathrm{s}}$,
$n_{\mathrm{r}}$)

  * •

common random noise to both sensor and reference signals ($n_{x}$)

  * •

independent line noise or harmonics of each device ($l_{\mathrm{s}}$,
$l_{\mathrm{r}}$)

  * •

common line noise or harmonics to both sensor and reference signals ($l_{x}$)

In the calibration process, sinusoidal vibration $x_{0}=\hat{x}_{0}\sin(2\pi
f_{\mathrm{v}}t)$ is applied. The random vibration noise $n_{x}$ and line
noise $l_{x}$ are added due to the background vibration, electrical noise of
the system, or the distortion of the waveform. Here, $l_{x}$ includes the line
noise, such as the power supply noise, which appears at constant frequencies,
and the harmonics of the input vibration, which appear at $2f_{\mathrm{v}}$,
$3f_{\mathrm{v}}$, and so on. Then, the reference signal measures the
displacement of the waveform, while the sensor to be calibrated responds to
the second derivative of it. Each output signal contains the independent
random and line noise; $n_{\mathrm{r}}$, $l_{\mathrm{r}}$, $n_{\mathrm{s}}$,
and $l_{\mathrm{s}}$. The sensor outputs the signal with the sensitivity
modulus of $S$. The recorded signals are the sum of these contributions, as
shown in Fig. 2. The calculated amplitudes of the recorded signals
$\hat{V}_{\mathrm{s,est}}$ and $\hat{x}_{\mathrm{r,est}}$ are affected by the
noise components as

$\displaystyle\hat{V}_{\mathrm{s,est}}=S(\hat{\ddot{x}}_{\mathrm{0,est}}+\hat{\ddot{n}}_{x\mathrm{,est}}+\hat{\ddot{l}}_{x\mathrm{,est}}+\hat{n}_{\mathrm{s,est}}+\hat{l}_{\mathrm{s,est}})$
(7)
$\displaystyle\hat{x}_{\mathrm{r,est}}=\hat{x}_{0\mathrm{,est}}+\hat{n}_{x\mathrm{,est}}+\hat{l}_{x\mathrm{,est}}+\hat{n}_{\mathrm{r,est}}+\hat{l}_{\mathrm{r,est}},$
(8)

which results in $S_{\mathrm{cal}}\neq S$ due to the second or later terms.
Note that the harmonics generated through a nonlinear process, such as the
nonlinearity of the sensor, are not independent of the input signal; for
example, $l_{\mathrm{s}}$ can be correlated to $l_{x}$ in general. However, we
assumed that the error sources in Fig. 2 were independent of each other for
simplicity. The random noise is discussed in Section 2.3, and the line noise
and harmonics are discussed in Section 2.4.

###  2.3 Effect of random noise

We consider the general case for random noise $n(t)$, which is characterized
with (one-sided) power spectral density (PSD) $G(f)$. The standard deviation
of the estimated amplitude for the noise,
$\sqrt{\langle|\hat{n}_{\mathrm{est}}|^{2}\rangle}$, is the standard
uncertainty of the amplitude of $x(t)=x_{0}(t)+n(t)$. From Eq. (4), the real
and imaginary parts of $\hat{n}_{\mathrm{est}}(t_{0})$ can be written as

$\displaystyle\mathrm{Re}[\hat{n}_{\mathrm{est}}(t_{0})]$ $\displaystyle=$
$\displaystyle\frac{1}{T}\int_{-\infty}^{\infty}\tilde{W}(f)\left\\{\tilde{N}(f_{\mathrm{v}}-f)+\tilde{N}^{\ast}(f_{\mathrm{v}}+f)\right\\}e^{-2\pi
ift_{0}}df,$ (9) $\displaystyle\mathrm{Im}[\hat{n}_{\mathrm{est}}(t_{0})]$
$\displaystyle=$
$\displaystyle\frac{1}{iT}\int_{-\infty}^{\infty}\tilde{W}(f)\left\\{\tilde{N}(f_{\mathrm{v}}-f)-\tilde{N}^{\ast}(f_{\mathrm{v}}+f)\right\\}e^{-2\pi
ift_{0}}df.$ (10)

$\tilde{N}(f)$ is the Fourier spectrum of $n(t)$. The PSD corresponding to the
Fourier spectrum
$\tilde{N}(f_{\mathrm{v}}-f)\pm\tilde{N}^{\ast}(f_{\mathrm{v}}+f)$ is given by
$G(|f_{\mathrm{v}}-f|)+G(|f_{\mathrm{v}}+f|)$. Therefore, the standard
deviations of Eq. (9) and (10) are both given by the integral of the PSD
$|\tilde{W}(f)|^{2}(G(|f_{\mathrm{v}}-f|)+G(|f_{\mathrm{v}}+f|))/T^{2}$. Since
$\mathrm{Re}[\hat{n}_{\mathrm{est}}(t_{0})]$ and
$\mathrm{Im}[\hat{n}_{\mathrm{est}}(t_{0})]$ are independent from each other
because of the randomness of $n(t)$, the amplitude estimation uncertainty
under the random noise is derived as

$\displaystyle\sqrt{\langle|\hat{n}_{\mathrm{est}}|^{2}\rangle}$
$\displaystyle=$
$\displaystyle\sqrt{\langle\mathrm{Re}[\hat{n}_{\mathrm{est}}]^{2}\rangle}=\sqrt{\langle\mathrm{Im}[\hat{n}_{\mathrm{est}}]^{2}\rangle}$
(11) $\displaystyle=$
$\displaystyle\sqrt{\int_{0}^{\infty}\frac{|\tilde{W}(f)|^{2}}{T^{2}}\left\\{G(|f_{\mathrm{v}}-f|)+G(|f_{\mathrm{v}}+f|)\right\\}df}$
(12) $\displaystyle=$
$\displaystyle\sqrt{\int_{-\infty}^{\infty}\frac{|\tilde{W}(f_{\mathrm{v}}-f)|^{2}}{T^{2}}G(|f|)df}$
(13)

This equation represents the uncertainty of vibration amplitude estimation
under the random noise with the PSD of $G(f)$. The uncertainty follows a
Gaussian distribution.

Note that the phase estimation uncertainty (in radian) under the random noise
is identical to the relative uncertainty of the amplitude:

$u(\arg[\hat{x}_{\mathrm{est}}])=u\left(\arctan\left(\frac{\mathrm{Im}[\hat{x}_{0}+\hat{n}_{\mathrm{est}}]}{\mathrm{Re}[\hat{x}_{0}+\hat{n}_{\mathrm{est}}]}\right)\right)=\frac{\sqrt{\langle|\hat{n}_{\mathrm{est}}|^{2}\rangle}}{\hat{x}_{0}}=\frac{u(|\hat{x}_{\mathrm{est}}|)}{\hat{x}_{0}}.$
(14)

This is because the distribution of $\hat{n}_{\mathrm{est}}(t_{0})$ is
isotropic in the complex plane. The component with the same argument as
$\hat{x}_{0}$ is the amplitude error, and the orthogonal component is the
phase error. Therefore, we investigated only the amplitude uncertainty using
Eq. (13) for the random noise in the following subsections and Section 3 for
simplicity as the same results are applicable to the phase. Their equivalence
is confirmed in Section 4.

####  2.3.1 Independent random noise for the sensor and reference

The sensor and reference signals contain independent random background noise
$n_{\mathrm{s}}$ and $n_{\mathrm{r}}$. They contribute to the measurement
result via $\hat{n}_{\mathrm{s,est}}$ and $\hat{n}_{\mathrm{r,est}}$, which
have random values for different measurements. Their standard deviations are
the standard uncertainty of the amplitude $V_{\mathrm{s}}$ and
$x_{\mathrm{r}}$. Their expressions
$\langle|\hat{n}_{\mathrm{s,est}}|^{2}\rangle$ and
$\langle|\hat{n}_{\mathrm{r,est}}|^{2}\rangle$ are given by Eq. (13) using the
PSDs $G_{\mathrm{s}}(f)$ and $G_{\mathrm{r}}(f)$. Consequently, the relative
standard uncertainty of the calibration sensitivity in the absence of the
other noise sources is

$\displaystyle\frac{u(S_{\mathrm{cal}})}{S}$ $\displaystyle=$
$\displaystyle\sqrt{\left(\frac{u(|\hat{V}_{\mathrm{s,est}}|)}{S|\hat{\ddot{x}}_{\mathrm{0,est}}|}\right)^{2}+\left(\frac{u(|\hat{x}_{\mathrm{r,est}}|)}{|\hat{x}_{0\mathrm{,est}}|}\right)^{2}}$
(15) $\displaystyle=$ $\displaystyle\frac{1}{(2\pi
f_{\mathrm{v}})^{2}\hat{x}_{0}}\sqrt{\int_{-\infty}^{\infty}\frac{|\tilde{W}(f_{\mathrm{v}}-f)|^{2}}{T^{2}}\left\\{G_{\mathrm{s}}(|f|)+(2\pi
f_{\mathrm{v}})^{4}G_{\mathrm{r}}(|f|)\right\\}df}.$ (16)

In the limit of long measurement time $T$, the window function is
asymptotically identical to the Dirac delta function as
$\lim_{T\to\infty}|\tilde{W}(f)|^{2}=T\delta(f)$; hence, the amplitude
estimation uncertainty is

$\sqrt{\langle|\hat{n}_{\mathrm{s,est}}|^{2}\rangle}\rightarrow\sqrt{\frac{G_{\mathrm{s}}(f_{\mathrm{v}})}{T}}.$
(17)

This relation indicates that the uncertainty is determined by the noise
spectrum at the vibration frequency and is inversely proportional to the
square-root of the measurement time. Eq. (17) gives the theoretical limit of
measurements, which cannot be avoided unless the background noise of the
system is reduced.

In reality, spectral leakage happens as shown in Eq. (13). This means that
even if the background noise is small at $f_{\mathrm{v}}$, the overall signal-
to-noise ratio (S/N) can be degraded by the noise in the other frequency band.
The amount of spectral leakage is determined by the shape of the window
function. Therefore, the leakage can be reduced by changing the window
function $w(t)$ or filtering the signal $V_{\mathrm{s}}$ and $x_{\mathrm{r}}$
to suppress the noise at $f\neq f_{\mathrm{v}}$. The details of these
modifications are discussed in Section 3.

Figure 3: Simulation results of the standard deviation of amplitude estimation
for two different noise spectrum cases; flat noise (left column) and
frequency-dependent noise (right column). The simulated waveforms (top row),
their PSDs (center row), and histograms of the estimated amplitudes (bottom
row). The standard deviations expected from Eq. (13) are shown in the legend
of the PSD figures, and the simulated values are shown in the histogram
figures.

To validate the calculations above, a simulation of amplitude estimation was
performed. The excited waveform was $x_{0}=\hat{x}_{0}\sin(2\pi
f_{\mathrm{v}}t)$ with $(2\pi f_{\mathrm{v}})^{2}\hat{x}_{0}=1$ m/s2 and
$f_{\mathrm{v}}=1$ Hz. The data length was set to $T=100$ s. The noise
$n_{\mathrm{r}}$ was randomly generated 300 times, and the vibration amplitude
was estimated from each $x_{0}+n_{\mathrm{r}}$ using Eq. (2). The standard
deviation of the estimated amplitudes $\hat{x}_{\mathrm{est}}$ relative to the
true amplitude $\hat{x}_{0}$ was calculated and compared to the theoretical
expectation from Eq. (13). For the PSD of the noise, $G_{\mathrm{r}}(f)$, two
cases were considered as examples: the flat spectrum $(2\pi
f)^{2}\sqrt{G_{\mathrm{r}}(f)}\simeq 10^{-2}$ (m/s2)/Hz1/2 and the frequency-
dependent spectrum. The latter had the same noise levels as the former at the
vibration frequency, $(2\pi
f_{\mathrm{v}})^{2}\sqrt{G_{\mathrm{r}}(f_{\mathrm{v}})}\simeq 10^{-2}$
(m/s2)/Hz1/2, while having a larger noise in $f\neq f_{\mathrm{v}}$. Fig. 3
shows the simulation results, which agree with the calculation obtained using
Eq. (13). The result demonstrates the importance of broadband noise
suppression.

####  2.3.2 Common random noise for the sensor and reference

The vibration noise of the exciter $n_{x}(t)$ is commonly measured by the
sensor and reference. Ideally, such a common noise appears in the two signals
in the same way; hence, it does not affect the calibration sensitivity in Eq.
(3). However, since the sensor measures $\ddot{n}_{x}(t)$ while the reference
measures $n_{x}(t)$, the spectral leakage has different contributions for the
two signals, which results in $\hat{\ddot{n}}_{x\mathrm{,est}}\neq(2\pi
f_{\mathrm{v}})^{2}\hat{n}_{x\mathrm{,est}}$. From the similar calculations as
Section 2.3.1, the relative standard uncertainty of the sensitivity
calibration is given by

$\displaystyle\frac{u(S_{\mathrm{cal}})}{S}$ $\displaystyle=$
$\displaystyle\frac{u(|\hat{\ddot{n}}_{x\mathrm{,est}}-(2\pi
f_{\mathrm{v}})^{2}\hat{n}_{x\mathrm{,est}}|)}{(2\pi
f_{\mathrm{v}})^{2}\hat{x}_{0}}$ (18) $\displaystyle=$
$\displaystyle\frac{1}{(2\pi
f_{\mathrm{v}})^{2}\hat{x}_{0}}\sqrt{\int_{-\infty}^{\infty}\frac{|\tilde{W}(f_{\mathrm{v}}-f)|^{2}}{T^{2}}\left\\{(2\pi
f)^{2}-(2\pi f_{\mathrm{v}})^{2}\right\\}^{2}G_{x}(|f|)df}$ (19)

in the absence of the other noise sources. Here, $G_{x}(f)$ is the PSD of the
common vibration noise $n_{x}(t)$. The difference from $n_{\mathrm{s}}$ and
$n_{\mathrm{r}}$ is that the estimation errors of $\hat{V}_{\mathrm{s,est}}$
and $\hat{x}_{\mathrm{r,est}}$ are correlated. Eq. (19) indicates that the
acceleration signal amplitude estimation is affected more by the leakage from
high-frequency noise than the reference displacement signal. The same effect
can matter when using random or triangle waveform excitation to simultaneously
calibrate sensitivity at multiple frequencies.

Figure 4: (left) Amplitude spectral density of the vibration noise in the low-
frequency vibration exciter of NMIJ. The measured spectrum (green) and model
(black) used for the calculations. (right) Standard deviation of the estimated
sensitivity under the random vibration noise. The results of simulation (red
circles) and calculation using Eq. (19) (blue dashed-line).

To validate the calculation and show the order of the uncertainty
contribution, we performed a simulation similar to that shown in Section
2.3.1. The excited acceleration amplitude was fixed to $(2\pi
f_{\mathrm{v}})^{2}\hat{x}_{0}=1$ m/s2 for frequency varying from 0.1 Hz to
300 Hz. The data length was set to $T=100/f_{\mathrm{v}}$. The background
vibration spectrum in the acceleration unit, $(2\pi f)^{2}\sqrt{G_{x}(f)}$, of
the calibration system in NMIJ is shown in Fig. 4. The smoothed spectrum model
was used for the simulation. The common vibration noise $n_{x}$ was randomly
generated 20 times at each frequency, and the amplitude of the reference
displacement, $\hat{x}_{\mathrm{r,est}}$, was estimated from each
$x_{0}+n_{\mathrm{r}}$ using Eq. (2). Then, the same time series were
numerically differentiated twice to prepare the sensor signal
$S(\ddot{x}_{0}+\ddot{n}_{\mathrm{r}})$, which was used to estimate the
amplitude $\hat{V}_{\mathrm{s,est}}$. The standard deviation of the estimated
amplitude ratio $\hat{V}_{\mathrm{s,est}}/(2\pi
f_{\mathrm{v}})^{2}\hat{x}_{\mathrm{r,est}}$ relative to the true sensitivity
modulus $S$ was calculated at each frequency and compared with the theoretical
result obtained from Eq. (19). The results are shown in Fig. 4. Eq. (19)
explains the simulated standard deviation of the sensitivity modulus. The
contribution ranged from 0.1 % to 1 % around 100 Hz, which is not ignorable in
accelerometer calibration.

###  2.4 Effect of line noise and harmonics

The error components include both the line noise and the harmonics of the
input vibration having constant amplitudes and phases. The difference between
them is that the frequency of the former is fixed and independent of
$f_{\mathrm{v}}$, while the latter always appears at integer multiples of
$f_{\mathrm{v}}$. Their contributions can be treated in the same way using Eq.
(4). Here, we consider the general case of the line noise
$l(t)=\hat{l}\sin(2\pi f_{l}t+\phi_{l})$ added to a sinusoidal wave
$x_{0}(t)=\hat{x}_{0}\sin(2\pi f_{\mathrm{v}}t)$. The complex amplitude
estimated from $x(t)=x_{0}(t)+l(t)$ is calculated from Eq. (4) as

$\hat{x}_{\mathrm{est}}=i\hat{x}_{0}+i\hat{l}\left(\frac{\tilde{W}(f_{\mathrm{v}}-f_{l})}{T}e^{-i\phi_{l}}-\frac{\tilde{W}(f_{\mathrm{v}}+f_{l})}{T}e^{i\phi_{l}}\right).$
(20)

Then, the amplitude and phase estimation errors are

$\displaystyle\frac{|\hat{x}_{\mathrm{est}}|}{\hat{x}_{0}}-1$ $\displaystyle=$
$\displaystyle\frac{\hat{l}}{\hat{x}_{0}}\frac{\tilde{W}(f_{\mathrm{v}}-f_{l})-\tilde{W}(f_{\mathrm{v}}+f_{l})}{T}\cos(\phi_{l}),$
(21) $\displaystyle\arg[\hat{x}_{\mathrm{est}}]-\frac{\pi}{2}$
$\displaystyle=$
$\displaystyle\frac{\hat{l}}{\hat{x}_{0}}\frac{\tilde{W}(f_{\mathrm{v}}-f_{l})+\tilde{W}(f_{\mathrm{v}}+f_{l})}{T}\sin(\phi_{l}).$
(22)

In actual measurements, the phase $\phi_{l}$ of the line noise relative to
$x_{0}$ is random for each measurement and is uniformly distributes from 0 to
$2\pi$. Therefore, their standard uncertainties are given by

$\displaystyle\frac{u(|\hat{x}_{\mathrm{est}}|)}{\hat{x}_{0}}$
$\displaystyle=$
$\displaystyle\frac{\hat{l}}{\sqrt{2}\hat{x}_{0}}\frac{|\tilde{W}(f_{\mathrm{v}}-f_{l})-\tilde{W}(f_{\mathrm{v}}+f_{l})|}{T},$
(23) $\displaystyle u(\arg[\hat{x}_{\mathrm{est}}])$ $\displaystyle=$
$\displaystyle\frac{\hat{l}}{\sqrt{2}\hat{x}_{0}}\frac{|\tilde{W}(f_{\mathrm{v}}-f_{l})+\tilde{W}(f_{\mathrm{v}}+f_{l})|}{T}.$
(24)

The uncertainty follows a U-shaped distribution. Since $\tilde{W}(f)$ has
large value around $f\simeq 0$, the uncertainty is large when
$f_{\mathrm{v}}\simeq f_{l}$. In this case, Eqs. (23) and (24) are
approximated as

$\frac{u(|\hat{x}_{\mathrm{est}}|)}{\hat{x}_{0}}\simeq
u(\arg[\hat{x}_{\mathrm{est}}])\simeq\frac{\hat{l}}{\sqrt{2}\hat{x}_{0}}\frac{|\tilde{W}(f_{\mathrm{v}}-f_{l})|}{T}.$
(25)

Notably, the harmonics of $f_{\mathrm{v}}$ do not affect the amplitude
estimation through the spectral leakage when the rectangular window is used
because $\tilde{W}_{\mathrm{r}}(f_{\mathrm{v}}\pm Nf_{\mathrm{v}})=0$ ($N$:
integer) if the length $T$ is integer multiples of the vibration period. It is
the same for the Hanning window or some other types of windows. Note that it
does not mean that the harmonics do not become an error source in any case. If
the sensor or reference interferometer has nonlinearity, the harmonics $l_{x}$
is nonlinearly converted with the fundamental wave $x_{0}$ to the line noise
$l_{\mathrm{s}}$ or $l_{\mathrm{r}}$ at $f_{\mathrm{v}}$ and can affect the
estimated amplitude even if $T$ is properly selected. To evaluate such a
contribution, detailed information about the input harmonics $l_{x}$ and the
nonlinearity of the sensor are necessary, which is out of the scope of this
work. As far as the harmonics $l_{x}$, $l_{\mathrm{s}}$, and $l_{\mathrm{r}}$
are independent from each other, the contribution of the harmonics can be
easily eliminated by the proper choice of $T$. Therefore, we mainly discuss
the line noise in the remainder of this article.

####  2.4.1 Independent line noise for the sensor and reference

The standard uncertainties of sensitivity under the independent line noise of
the sensor and reference, $l_{\mathrm{s}}=\hat{l}_{\mathrm{s}}\sin(2\pi
f_{l,\mathrm{s}}t+\phi_{l,\mathrm{s}})$ and
$l_{\mathrm{r}}=\hat{l}_{\mathrm{r}}\sin(2\pi
f_{l,\mathrm{r}}t+\phi_{l,\mathrm{r}})$, are given by the square-root of the
square sum of their uncertainty contributions. Using Eq. (25), it is
approximated as

$\frac{u(S_{\mathrm{cal}})}{S}=\frac{1}{\sqrt{2}}\sqrt{\left(\frac{\hat{l}_{\mathrm{s}}}{(2\pi
f_{\mathrm{v}})^{2}\hat{x}_{0}}\right)^{2}\frac{|\tilde{W}(f_{\mathrm{v}}-f_{l,\mathrm{s}})|^{2}}{T^{2}}+\left(\frac{\hat{l}_{\mathrm{r}}}{\hat{x}_{0}}\right)^{2}\frac{|\tilde{W}(f_{\mathrm{v}}-f_{l,\mathrm{r}})|^{2}}{T^{2}}}.$
(26)

Here, we assumed that $\phi_{l,\mathrm{s}}$ and $\phi_{l,\mathrm{r}}$ are not
correlated to each other. Accurate expression is given by replacing
$\tilde{W}(f_{\mathrm{v}}-f_{l})$ with
$\tilde{W}(f_{\mathrm{v}}-f_{l})-\tilde{W}(f_{\mathrm{v}}+f_{l})$ for
sensitivity and
$\tilde{W}(f_{\mathrm{v}}-f_{l})-\tilde{W}(f_{\mathrm{v}}+f_{l})$ for phase
delay.

When the line noise frequency is equal to the vibration frequency,
$\tilde{W}(f_{\mathrm{v}}-f_{l})=T$. Therefore, the relative amplitude
estimation uncertainties of the sensor and reference are determined by only
the amplitude ratio:

$\frac{u(|\hat{x}_{\mathrm{r,est}}|)}{\hat{x}_{0}}=\frac{\hat{l}_{\mathrm{r}}}{\sqrt{2}\hat{x}_{0}}.$
(27)

This is the fundamental limit of uncertainty from the line noise. Although the
leakage from $f\neq f_{\mathrm{v}}$ can be reduced by proper signal
processing, the contribution at $f=f_{\mathrm{v}}$ cannot be avoided unless
the line noise amplitude is reduced.

Figure 5: Relative amplitude (left) and phase (right) estimation uncertainty
due to the line noise in the case of $\hat{l}/\hat{x}_{0}=10^{-3}$. The open
circles show the simulated result, the solid magenta and cyan lines show the
theoretical uncertainty expected from Eqs. (23) and (24), respectively, and
the black lines show the approximated expression Eq. (25).

The line noise contribution was confirmed by simulation. Fixing
$\hat{l}/\hat{x}_{0}=10^{-3}$, we estimated the amplitude and phase of
$x(t)=x_{0}(t)+l(t)$ for various $f_{\mathrm{v}}/f_{l}$ and $\phi_{l}$. The
data length was fixed to $T=100/f_{\mathrm{v}}$ at each frequency. The
standard deviations for different $\phi_{l}$ were calculated from the
simulated values and compared with Eqs. (23) and (24). Fig. 5 shows the
simulation results. Theoretical calculations explained the uncertainty well.
The approximated formula, Eq. (25), also accurately estimated the uncertainty
around the line noise frequency $f_{l}$, where the line noise contribution
becomes important. At frequencies away from $f_{l}$, Eq. (25) led to
over-/under-estimation by a few times.

####  2.4.2 Common line noise for the sensor and reference

The standard uncertainty of the sensitivity under the common line noise
$l_{x}=\hat{l}_{x}\sin(2\pi f_{l,x}t+\phi_{l,x})$ was calculated similarly as
Eq. (18). Using Eq. (25), it was approximated as

$\frac{u(S_{\mathrm{cal}})}{S}\simeq
u(\Delta\phi_{\mathrm{cal}})\simeq\frac{\hat{l}_{x}}{\sqrt{2}\hat{x}_{0}}\frac{|f_{\mathrm{v}}^{2}-f_{l,x}^{2}|}{f_{\mathrm{v}}^{2}}\frac{|\tilde{W}(f_{\mathrm{v}}-f_{l,x})|}{T}.$
(28)

Again, accurate expressions are given by replacing
$\tilde{W}(f_{\mathrm{v}}-f_{l,x})$ with
$\tilde{W}(f_{\mathrm{v}}-f_{l,x})\pm\tilde{W}(f_{\mathrm{v}}+f_{l,x})$.

Figure 6: Sensitivity (left) and phase delay (right) estimation uncertainty
due to the common line noise $l_{x}$ with
$\hat{\ddot{l}}_{x}/\hat{\ddot{x}}_{0}=10^{-3}$. The approximation formula Eq.
(28) is shown with a black line, and the accurate formulae are shown with red
and blue lines. The simulated values are plotted with open circles.

The effect of the line noise was simulated again. The common line noise
amplitude was fixed as $\hat{\ddot{l}}_{x}/\hat{\ddot{x}}_{0}=10^{-3}$, and
the data length was fixed to $T=100/f_{\mathrm{v}}$ at each frequency. The
sensor and reference signals $V_{\mathrm{s}}=S(\ddot{x}_{0}+\ddot{l}_{x})$ and
$x_{\mathrm{r}}=x_{0}+l_{x}$ were prepared for different line phase
$\phi_{l,x}$. The standard deviations of the estimated sensitivity and phase
for different $\phi_{l,x}$ were calculated. The simulation results are shown
in Fig. 6. The results agreed with the theoretical calculations using Eq.
(28). Unlike the independent line noise, the uncertainty contribution of the
common line noise is not concentrated around $f_{\mathrm{v}}\simeq f_{l,x}$,
because the peak of $\tilde{W}(f_{\mathrm{v}}-f_{l,x})$ is canceled by the
factor $f_{\mathrm{v}}^{2}-f_{l,x}^{2}$, as shown in Eq. (28), when the
rectangular window is used.

##  3 Reduction of noise contribution in calibration

The calibration uncertainty considered in this paper is explained by Eqs.
(16), (19), (26), and (28). These equations show how the spectral leakage
contributes to the calibration in the conventional acceleration calibration
with the SAM. As already mentioned, the background random/line noise of the
sensing parts (sensor and interferometer) at $f_{\mathrm{v}}$ is a fundamental
limit of signal processing. In this section, we aimed to minimize the leakage
from $f\neq f_{\mathrm{v}}$. In the following subsections, we propose three
signal processing modification methods including filtering, changing the
window function, and numerical differentiation to align the unit of the
measurand. For the line noise, the selection of the data length $T$ is also
discussed. The investigations about the three modifications are mainly focused
on the random noise, although the same methods are applicable to the line
noise. These methods are simulated for both random and line noises.

###  3.1 Filtering the signal

One of the simplest ways to reduce the leakage is filtering the signal before
amplitude estimation. The bandpass filter (BPF) centered at $f_{\mathrm{v}}$
reduces the noise PSD $G(f\neq f_{\mathrm{v}})$ in Eqs. (16) and (19) and the
line noise at $f\neq f_{\mathrm{v}}$. Depending on the background noise
spectrum, the low-pass or high-pass filters can also be used. Such filtering
is already adopted in the calibration process in NMIJ, although the shape of
the filter has been empirically determined. The filtering changes $G(f\neq
f_{\mathrm{v}})$ in Eqs. (16) and (19) to $|\tilde{F}(f)|^{2}G(f)$, where
$\tilde{F}(f)$ is the transfer function of the filter. The line noise
amplitude also decreases to $|\tilde{F}(f)|\hat{l}$ in Eqs. (26) and (28).

Figure 7: The noisy waveform with the constant sinusoidal amplitude starting
at $t=0$ (grey) and bandpass-filtered data (orange).

The disadvantage of filtering is that it takes time until the waveform becomes
stable after starting the excitation. An example is shown in Fig. 7. A sixth-
order BPF was applied to the noisy waveform to extract the frequency component
around $f_{\mathrm{v}}$ with a Q-factor ($Q$) of $\simeq 2$. Here, the
Q-factor is the ratio of the center frequency $f_{\mathrm{v}}$ over the
bandwidth of the filter. In the example, the amplitude of the filtered data
became stable after about eight vibration cycles. Consequently, the first
several vibration cycles of the data could not be used for amplitude
estimation, which increased the measurement time, especially at low
frequencies.

Figure 8: (left) Simulated settling times of bandpass filters with different
combinations of the order (2, 4, 6, and 8) and Q-factor (1, 2, 5, 10, and 20).
(right) The filter gain of the sixth-order bandpass filters used for the
calculation.

Fig. 8 shows the simulated settling time for different orders and Q-factors of
the BPF. The settling time is defined as the time length until the filtered
waveform amplitude settles within $\pm 0.1$ % of the true amplitude. The
settling time is roughly proportional to the Q-factor, and increases with the
order of the filter. Thus, a trade-off exists between the noise reduction
ability of the filter and discarded measurement time. Although it may be
possible to compensate for the amplitude change associated with settling, it
may take effort to perform a simulation according to the shape of the filter
and the input signal. In the remainder of this article, we use a sixth-order
BPF with $Q=1$ to keep the settling time within 10 vibration cycles.

###  3.2 Changing the window function

Another way to improve the vibration amplitude estimation accuracy is changing
the window function $w(t)$. The rectangular window, which is implicitly used
in the conventional SAM, is famous for large spectral leakage. The other
window functions such as a Hanning window can reduce $\tilde{W}(f)$ at $f\neq
f_{\mathrm{v}}$. Although this method is well-known for fast Fourier transform
applications, it has not been discussed for the SAM in vibration calibration.

Figure 9: The amount of the spectral leakage in vibration amplitude estimation
from $T=30/f_{\mathrm{v}}$ (upper) and $100/f_{\mathrm{v}}$ (lower) data.
Three cases using the conventional SAM (red), sixth-order bandpass filtering
(orange), and Hanning window (green), are plotted. The enlarged view around
$f/f_{\mathrm{v}}=1$ is shown in the inset.

The effect of the window function is compared with filtering in Fig. 9. The
factor $|\tilde{W}(f_{\mathrm{v}}-f)\tilde{F}(f)|/T$, which determines the
amount of the leakage in both the random noise and line noise, is plotted for
$T=30/f_{\mathrm{v}}$ and $100/f_{\mathrm{v}}$ in three cases: using the
rectangular window (conventional SAM), by applying the sixth-order BPF ($Q=1$)
with the rectangular window, and using the Hanning window without filtering.
As the figures show, the Hanning window can reduce the leakage more
effectively than the BPF around the vibration frequency. The reduction ratio
of the Hanning window in $0.7<f/f_{\mathrm{v}}<1.5$ is approximately equal to
the second-order BPF with $Q\simeq 15$ and $50$ for $T=30/f_{\mathrm{v}}$ and
$100/f_{\mathrm{v}}$, respectively. In such cases, the long measurement time
is wasted due to the settling time, as shown in Fig. 8, while changing the
window function does not require discarding the data. On the other hand, the
BPF has a better reduction ratio at frequencies away from $f_{\mathrm{v}}$,
although the details depend on the design of the filter. The window functions
also increase the leakage from a very close frequency to $f_{\mathrm{v}}$
compared to the rectangular window, as shown in the inset in Fig. 9. In other
words, the frequency resolution is compromised by changing the window
function. Therefore, when there is a large low-frequency noise or the
calibration frequency is close to a large line noise (e.g., calibration at 49
Hz when there is a large power supply line noise at 50 Hz), bandpass filtering
with the rectangular window can be a better choice. In the scope of vibration
frequency determination, the frequency resolution of signal processing is
usually not important in vibration calibration because the frequency is
accurately controlled based on reference frequency standards.

###  3.3 Numerical differentiation

For the common vibration noise $n_{x}$, the factor $((2\pi f)^{2}-(2\pi
f_{\mathrm{v}})^{2})$ is included as in Eq. (19). Similarly, the uncertainty
from the common line noise, Eq. (28), contains the factor
$|f_{\mathrm{v}}^{2}-f_{l,x}^{2}|$. These factors originate from the
difference of the measured physical quantity between the sensor and the
reference interferometer. The accelerometer measures the acceleration, while
the interferometer measures the displacement. In the conventional SAM for the
reference signal, the displacement amplitude is estimated from
$x_{\mathrm{r}}$; then, the acceleration amplitude is calculated by
multiplying $(2\pi f_{\mathrm{v}})^{2}$. Although this gives an accurate
estimation for a purely sinusoidal wave, the background noise degrades
calibration accuracy.

A natural way to avoid such an effect is converting the reference displacement
signal into acceleration by numerical differentiation. Then, the amplitude of
the reference signal is estimated in the unit of acceleration. This process
ideally eliminates the factor $((2\pi f)^{2}-(2\pi f_{\mathrm{v}})^{2})$ in
Eq. (19) and $|f_{\mathrm{v}}^{2}-f_{l,x}^{2}|$ in Eq. (28). Therefore, the
common noise sources $n_{x}$ and $l_{x}$ no longer interfere with the
calibration. A similar process is also applicable to the velocity sensor
(e.g., seismometer) by changing the number of differentiation. The important
point is to align the units of the signals.

The differentiation process of the displacement signal relatively enlarges the
high-frequency noise in the time domain. Additional processing such as low-
pass filtering or changing the window function is required to avoid increasing
noise contributions. On the other hand, the low-frequency noise, which is
often included as a drift component in the reference interferometer signal, is
reduced because the differentiation process works similarly to a second-order
high-pass filter with infinite cutoff frequency.

###  3.4 Adjusting data length

The line noise contribution can be eliminated by setting proper data length
$T$ so that both $f_{\mathrm{v}}-f_{l}$ and $f_{\mathrm{v}}+f_{l}$ are at the
zero points of $\tilde{W}(f)$. The Fourier transform of the rectangular
window, $\tilde{W}_{\mathrm{r}}(f)$, is equal to zero at $fT=N$ ($N$:
integer). Recalling that $f_{\mathrm{v}}T$ is set to an integer to eliminate
the harmonics, as mentioned in Section 2.4, $f_{l}T$ also needs to be an
integer. Consequently, if $T$ is an integer multiple of the inverse of the
greatest common divisor of $f_{\mathrm{v}}$ and $f_{l}$, both the harmonics
and line noise can be eliminated. For example, in the case of $f_{l}=50$ Hz
and $f_{\mathrm{v}}=49.2$ Hz, the greatest common divisor frequency is 0.4 Hz;
hence $T$ should be the multiple of 2.5 s. The condition is similar for the
Hanning window, although it has two fewer zero points around $f_{l}$ than the
rectangular window, as shown in Fig. 9. In usual vibration calibration,
$f_{\mathrm{v}}$ is selected from the one-third octave band specified in ISO
266:1997, and the frequencies are rounded to 0.5 Hz increment around the main
line noise frequency (50 Hz). In such a case, $T$ should be multiple of 2 s.
If there are large line noises at several frequencies, it will be good to set
$T$ so that the closest line noise to $f_{\mathrm{v}}$ is eliminated.

One of the possible drawbacks of adjusting the data length is that the
measurement time can be somewhat longer in some cases due to the limited
choice of $T$. If the elimination is impossible within a reasonable
measurement time, the line noise contribution needs to be reduced using other
methods proposed in Section 3.1 to 3.3. Note that the adjustment is not
necessary at every frequency because the line noise contribution becomes a
problem only around the line frequency. Therefore, the total measurement time
of calibration does not substantially increase over a wide frequency range.

###  3.5 Simulation of proposed methods

Figure 10: Summary of the signal processing methods for sensitivity
calibration including (0) the conventional SAM and the proposed methods; (a)
bandpass filtering, (b) multiplying window function, and (c) differentiation
of the reference displacement signal.

So far, we have proposed three signal processing methods for sensitivity
estimation, which are summarized in Fig 10. Proper selection of the data
length $T$ to remove the line noise contribution was also discussed. To
confirm their effectiveness, a simulation of amplitude estimation was
performed for the two types of random noise discussed in Section 2.3 and the
line noise mentioned in Section 2.4.

####  3.5.1 Reduction of independent random noise

Figure 11: Waveforms used in the simulation ($f_{\mathrm{v}}=1$ Hz, $(2\pi
f_{\mathrm{v}})^{2}\hat{x}_{0}=1$ m/s2) (top). Three cases using the
conventional SAM (red), sixth-order bandpass filtering (orange), and Hanning
window (green), are plotted. For the second case, the data of the first ten
vibration cycles (grey colored) is not used for amplitude estimation.
Histograms of the estimated amplitudes for the three cases (bottom).

First, the independent noise sources, $n_{\mathrm{s}}$ and $n_{\mathrm{r}}$,
were investigated. The additional signal processing methods proposed in
Section 3.1 and 3.2 were effective in this case and corresponded to (a) and
(b) in Fig. 10, respectively. Here, the amplitude estimation uncertainty of
the single signal ($V_{\mathrm{s}}$ or $x_{\mathrm{r}}$) was simulated for
simplicity. Using the results, the sensitivity calibration uncertainty can be
calculated based on Eq. (16). A simulation similar to that in Section 2.3.1
was performed; the amplitude estimation of sinusoidal waveform with $(2\pi
f_{\mathrm{v}})^{2}\hat{x}_{0}=1$ m/s2 under the background noise as the right
column case of Fig. 3. In this case, the standard deviation of the estimated
amplitude increased to 0.83 % based on the conventional SAM due to the noise
at $f\neq f_{\mathrm{v}}$, while the S/N at the vibration frequency was 0.1 %.
Here, the sixth-order Butterworth BPF with $Q=1$ and the Hanning window were
applied to the noisy waveform, and the amplitude was estimated from the
processed waveforms. The examples of the processed waveforms are shown in Fig.
11 (upper figure). The data of the first ten vibration cycles were discarded
for the filtered waveform, as explained in Section 3.1. The simulation was
repeated 300 times for the randomly generated noise, and the histogram of the
amplitude estimation error is shown in Fig. 11 (lower figure). For the
conventional SAM, BPF, and Hanning window, the simulated standard deviations
were 0.91 %, 0.13 %, and 0.15 %, and the expected deviations from Eq. (13)
were 0.83 %, 0.13 %, and 0.15 %, respectively. As expected, the BPF and
Hanning window reduce the background noise contributions at $f\neq
f_{\mathrm{v}}$. Although they do not achieve an S/N of 0.1 % at
$f_{\mathrm{v}}$ due to the leakage around $f_{\mathrm{v}}$, the difference is
sufficiently small. Consequently, the calibration uncertainty is also reduced
because Eq. (16) is the sum of the amplitude estimation uncertainties of
$V_{\mathrm{s}}$ and $x_{\mathrm{r}}$.

In any cases shown in Fig. 11, the amplitude estimation errors are distributed
around zero, which indicates that there is no bias due to the filter or
window. The correction of the filter gain is necessary, depending on the
choice of the filter. To avoid the systematic effects on sensitivity
calibration, applying the same process on both the sensor and reference
signals is recommended to cancel out the effect.

####  3.5.2 Reduction of common random noise

Figure 12: Standard deviation of the estimated sensitivity under the random
vibration noise $n_{x}$ (Fig. 4) for four cases; using the conventional SAM
(red), sixth-order bandpass filter (orange), Hanning window (green), and
numerical differentiation (blue). The simulation results are shown with the
circles, and the theoretical expectations from Eq. (19) are plotted with the
dashed lines. The theoretical expectation for the numerical differentiation is
not shown because it is zero.

Second, the common vibration noise, $n_{x}$, was investigated. The signal
processing methods proposed in Section 3.1, 3.2, and 3.3 were effective in
this case. They correspond to (a), (b), and (c) in Fig. 10, respectively. The
noise spectrum in Fig. 4 and the excitation amplitude of $(2\pi
f_{\mathrm{v}})^{2}\hat{x}_{0}=1$ m/s2 were used for the simulation. In (a)
and (b), the same BPF and Hanning window as the previous simulation (Fig. 11)
were used. The amplitudes were estimated using the processed data, and the
sensitivity was calculated from Eq. (3). In (c), the reference displacement
data $\\{x_{\mathrm{r},n}\\}$ ($n=0,1,...,N-1$) was converted to the reference
acceleration $a_{\mathrm{r},n}$:

$a_{\mathrm{r},n}=(x_{\mathrm{r},n+1}-2x_{\mathrm{r},n}+x_{\mathrm{r},n-1})f_{\mathrm{s}}^{2},$
(29)

where $f_{\mathrm{s}}$ is the sampling frequency. Then, the amplitudes were
estimated using $\\{v_{\mathrm{s},n}\\}$ and $\\{a_{\mathrm{r},n}\\}$, and the
sensitivity was calculated as
$S_{\mathrm{cal}}=|\hat{V}_{\mathrm{s,est}}/\hat{a}_{\mathrm{r,est}}|$ instead
of Eq. (3). These simulations are repeated 20 times for the randomly generated
noise at vibration frequencies
$0.1~{}\mathrm{Hz}<f_{\mathrm{v}}<300~{}\mathrm{Hz}$. The simulation results
and theoretical expectations from Eq. (19) are shown in Fig. 12. Both the BPF
and Hanning window reduced the relative standard uncertainty of sensitivity
from 0.3 % to 0.01 % and 0.001 %, respectively, around 100 Hz. The reduction
effect agreed well with the theory. Numerical differentiation is much more
effective, as shown in the figure. Therefore, it is recommended to adopt the
differentiation process in accelerometer calibration when the background
vibration noise is large.

Nevertheless, the systematic error of numerical differentiation requires
attention. The transfer function of Eq. (29) from the displacement to
acceleration is given by

$(e^{2\pi if/f_{\mathrm{s}}}-2+e^{-2\pi
if/f_{\mathrm{s}}})f_{\mathrm{s}}^{2}=-(2\pi
f)^{2}\left(\mathrm{sinc}(f/f_{\mathrm{s}})\right)^{2},$ (30)

while the second-order derivative of continuous signal is given by $-(2\pi
f)^{2}$. The factor
$\left(\mathrm{sinc}(f_{\mathrm{v}}/f_{\mathrm{s}})\right)^{2}$ needs to be
corrected when numerical differentiation is applied. Note that the factor
differs from 1 by less than 0.02 % for $f_{\mathrm{v}}<0.01f_{\mathrm{s}}$. If
the sampling frequency is sufficiently faster than the vibration frequency,
typically by 100 times, the error is ignorable.

####  3.5.3 Reduction of independent line noise

Figure 13: Uncertainty simulation for the line noise with
$\hat{l}/\hat{x}_{0}=10^{-3}$. Three cases using the conventional SAM (red
circle), sixth-order bandpass filtering (orange square), and Hanning window
(green triangle), are plotted.

The same methods applied to the random noise are also effective for the line
noise. For the independent line noise, $l_{\mathrm{s}}$ and $l_{\mathrm{r}}$,
the filter and window function are useful. They correspond to (a) and (b) in
Fig. 10, respectively. The effects of the BPF and Hanning window on the
independent line noise were simulated here. The signal $x(t)=x_{0}(t)+l(t)$
was filtered by the sixth-order BPF with $Q=1$ or windowed using the hannin
window. The amplitude estimation uncertainties were simulated using the same
process described in Section 2.4. The result is shown in Fig. 13. As Eq. (25)
indicates, the reduction effect is similar to the frequency dependence of
$|\tilde{W}(f_{\mathrm{v}}-f)\tilde{F}(f)|/T$, shown in Fig. 9. The Hanning
window is more effective than the BPF over the simulated frequency range,
where the line noise has a large uncertainty contribution.

Figure 14: Simulated amplitude estimation uncertainty at $f_{\mathrm{v}}=49.2$
Hz under the line noise at $f_{l}=50$ Hz (red) and harmonics at
$f_{l}=2f_{\mathrm{v}}$ for different data length $T$. The line amplitude is
$\hat{l}/\hat{x}_{0}=10^{-3}$ at each frequency.

To confirm that the proper choice of $T$ can eliminate the line noise
contribution, as discussed in 3.4, the dependence on the data length was also
simulated. As an example, the calibration frequency was fixed at
$f_{\mathrm{v}}=49.2$ Hz, and the data length was varied changed from 0.1 s to
6 s. The amplitude estimation uncertainty for $x(t)=x_{0}(t)+l(t)$ was
simulated using the same method described in Section 2.4.1 (using the
rectangular window withoug a filter). For $l(t)$, both the line noise at fixed
frequency $f_{l}=50$ Hz and harmonics at $f_{l}=2f_{\mathrm{v}}$ were used for
the simulation. The results are shown in Fig. 14. As discussed in Section 3.4,
both uncertainties from the line noise and harmonics reach almost zero
(limited by numerical computation error) when $T$ is multiple of 2.5 s, which
is the inverse of the greatest common divisor of $f_{\mathrm{v}}$ and $f_{l}$.

####  3.5.4 Reduction of common line noise

Figure 15: Sensitivity estimation uncertainty due to the common line noise
$l_{x}$ with $\hat{\ddot{l}}_{x}/\hat{\ddot{x}}_{0}=10^{-3}$. The simulated
uncertainties for four cases using the conventional SAM (red circle), sixth-
order bandpass filter (orange square), Hanning window (green triangle), and
numerical differentiation (blue cross), are plotted.

For the common line noise, $l_{x}$, numerical differentiation is effective.
The uncertainties with the conventional SAM, sixth-order BPF, Hanning window,
and numerical differentiation were simulated and plotted in Fig. 15.The
simulation was performed under the same conditions, as described in Section
2.4.2. The same filter, window, and differentiation process described in
Section 3.5.2 were applied. As expected, those signal process modifications
reduce the common line noise contributions. Especially, the numerical
differentiation eliminates the common line noise contribution over a wide
frequency range.

Note that the proper choice of $T$ is also effective for reducing the common
line noise. Therefore, the modification of signal processing is not necessary
for line noise reduction, if it is already eliminated by the choice of $T$.

##  4 Application to the accelerometer calibration in NMIJ

In NMIJ, the calibration system with a small excitation amplitude ($(2\pi
f_{\mathrm{v}})^{2}\hat{x}_{0}\simeq 10^{-2}$ m/s2) is under development [10].
The background noise of the calibration system becomes a significant
uncertainty source in micro vibration calibration. Here, the proposed
processing methods are applied to the actual calibration data to reduce
uncertainty. The calibration system is shown in Fig. 1. A servo accelerometer
JA-5V (Japan Aviation Electronics Industry, Ltd., $S\simeq 0.1$ V/(m/s2)
(nominal)) was used as the calibration target. Since the calibration
uncertainty discussed in this paper mainly affects repeatability, the
calibration was repeated 5 times at frequencies from 0.4 Hz to 500 Hz, and
their standard deviations were measured at each frequency. The results were
compared with the theoretical limits due to the background random noise of
$V_{\mathrm{s}}$ and $x_{\mathrm{r}}$ (as Eq. (17)) in our system. The effect
of the line noise was limited and not large compared to the random noise;
hence, we mainly discuss the random noise in this section. About the
independent background noise, $n_{\mathrm{s}}$ includes the self-noise of the
accelerometer and the noise of the signal acquisition system, and
$n_{\mathrm{r}}$ includes the self-noise of the laser interferometer, seismic
vibration noise, and signal acquisition system noise. The common background
noise $n_{x}$ is shown in Fig. 4. Relatively large low-frequency drift, which
also behaves as the common background noise, was applied above 30 Hz for
averaging the cyclic error of the interferometer.

Figure 16: Standard deviations of the calibrated sensitivities (upper) and
phase delays (lower) repeated 5 times with $(2\pi
f_{\mathrm{v}})^{2}\hat{x}_{0}\simeq 10^{-2}$ m/s2 using the calibration
system in NMIJ. The results using different signal processing methods are
plotted: the conventional SAM (red solid circle), sixth-order BPF (orange
solid square), Hanning window (green solid triangle), numerical
differentiation for $x_{\mathrm{r}}$ (blue open circle), differentiation and
BPF (brown open square), and differentiation and Hanning window (purple open
triangle). The theoretical limit due to the background noise of the system is
shown by the solid grey line.

The experimentally obtained repeatabilities with different signal processing
methods are shown in Fig. 16. The sixth-order BPF, Hanning window, numerical
differentiation (Eq. (29)), and their combinations were adopted. With the
conventional SAM without any additional processing methods, the repeatability
uncertainty was about 100 % over a broad frequency range. The calibration
results were almost meaningless in such a large uncertainty. By using the BPF,
Hanning window, or numerical differentiation, the repeatability was improved
by up to four orders of magnitude. Since numerical differentiation could
reduce uncertainty, the overall repeatability was thought to be limited by the
common vibration noise $n_{x}$ or low-frequency noise in the reference
displacement signal $n_{\mathrm{r}}$. At high frequencies above 50 Hz, the
noise reduction is insufficient only in the case of the Hanning window because
the low-frequency drift is not suppressed enough. The BPF is a better choice
in such a case. The combinations of numerical differentiation with BPF or
Hanning window were more effective than the single processing. The results
with the combined processing methods are also shown in Fig. 16, which achieved
the theoretical limit by the calibration system noise. As expected from Eq.
(14), the reduction effects are almost the same for the phase delay.

Our experiment proved that the signal processing methods proposed in this
paper are useful in actual calibration systems. The results agreed with the
theoretically expected reduction ability; therefore, the calculations shown in
Section 2 and 3 can be used for optimizing signal processing. For example, in
the case of the accelerometer calibration in NMIJ, the common vibration noise
had a dominant uncertainty contribution. The independent background noise
bottomed out near 10 Hz and increased on the low- and high-frequency side;
hence, the calibration around 10 Hz was affected by the leakage from the low-
and high-frequency ranges. Since both the common and independent noises needed
to be reduced, the combination of numerical differentiation and the BPF or
Hanning window was required. If the common vibration noise is small, numerical
differentiation is not necessary. The spectral leakage may be insignificant if
the independent background noise has a flat spectrum; then, the conventional
SAM without the BPF or Hanning window may be sufficient. Thus, the processing
methods need to be combined based on the calibration system noise.

##  5 Conclusions

Figure 17: Combination of the numerical differentiation and the filter or
window.

We proposed three signal processing methods to modify the conventional SAM
specified in ISO16063-11. The proposed methods include

  * •

filtering the signal,

  * •

changing the window function,

  * •

differentiation of the reference displacement signal.

The proper choice of the data length $T$ was also discussed to eliminate the
line noise contribution. Their effect on the reduction of calibration
uncertainty was confirmed by both the simulation and experiment. These results
showed that, in accelerometer calibration, using the combination of numerical
differentiation and the filter or window before amplitude estimation results
in a more robust calibration against the background noise compared to using
the conventional SAM.

Table 1: Summary of the effectiveness ($++$:contribution is eliminated, $+$:contribution is reduced, $-$:ineffective) and advantage/disadvantage of the proposed signal processing methods. method | effect on |  advantage |  disadvantage
---|---|---|---
$n_{x}$ | $l_{x}$ |  $n_{\mathrm{r}}$, $n_{\mathrm{s}}$ |  $l_{\mathrm{r}}$, $l_{\mathrm{s}}$
filter | $+$ | $+$ | $+$ | $+$ |  large reduction ratio at $f$ away from $f_{\mathrm{v}}$, design is flexible |  data is wasted during the settling time
window | $+$ | $+$ | $+$ | $+$ |  large reduction ratio at $f$ close to $f_{\mathrm{v}}$ |  low-frequency noise remains
differentiation | $++$ | $++$ | $-$ | $-$ |  common noise is eliminated, low-frequency reference noise is reduced |  high-frequency reference noise is amplified
adjusting $T$ | $-$ | $++$ | $-$ | $++$ |  line noise is eliminated |  not always possible within reasonable $T$

Table 1 summarizes the effectiveness, advantages, and disadvantages of the
proposed methods. The line noise and common noise contributions can be
eliminated by adjusting the data length and using numerical differentiation
for the reference displacement signal. If the independent random noise
contribution is significant, additional filtering or window function needs to
be combined with numerical differentiation. Fig. 17 shows the recommended
combinations of sensitivity calibration. The combination of numerical
differentiation and Hanning window is sufficient in standard cases. Although
filtering offers flexibility, the data of the settling time needs to be
discarded, and instability is a concern for the infinite impulse response
filter. Eqs. (16), (19), (26), and (28) are useful for the optimization based
on the system noise characterization.

This work enables accelerometer calibration with a small excitation amplitude
relative to the system background noise. Such calibration is required to
confirm the sensitivity linearity of accelerometers used for micro vibration
measurements. Additionally, for accelerometers with high sensitivity, large
vibration may not be applied for calibration because they can get saturated
with small input vibration. In the calibration at a low frequency ($<1$ Hz),
which is required for broadband seismometers, the excited acceleration
amplitude is limited due to the stroke limit of the exciter. The reduction of
the background noise is necessary in these cases, and proper signal processing
is essential to take full advantage of noise reduction.

Although we mainly discussed accelerometer calibration in this work, some of
the knowledge obtained also applies to other fields of dynamic sensor
calibration, where sinusoidal signal extraction is required. When the sensor
output signal is compared with a reference signal, their signals should be in
unit of the same physical quantity; otherwise, the common background noise
affects the estimation of the amplitude ratio. The demand for measuring small
fluctuations is increasing in not only mechanical vibration measurements but
also various fields along with the progress of industry. The proposed signal
processing methods for accelerometer calibration can contribute to improving
the reliability of those measurements.

## Acknowledgment

This work was partially based on the results obtained from a project
commissioned by the New Energy and Industrial Technology Development
Organization (NEDO), Japan.

## References

## References

  * [1] Brownjohn J M W 2007  Structural health monitoring of civil infrastructure, Philosophical Transactions A: Mathematical, Physical, and Engineering Sciences, 365, 589
  * [2] Calvi A and Roy N 2013  Spacecraft mechanical loads analysis handbook, ESA Requirements and Standards Division, Noordwijk, The Netherlands
  * [3] Deng T, Chen D, Wang J, Chen J, Sun Z and Li G 2015  Microelectromechanical systems-based electrochemical seismic sensors with insulating spacers integrated electrodes for planetary exploration, IEEE Sensors Journal, 16, 3, 650
  * [4] Isobe A, Kamada Y, Takubo C, Furubayashi Y, Oshima T, Sakuma N and Sekiguchi T 2020  Design of perforated membrane for low-noise capacitive MEMS accelerometers, IEEE Sensors Journal, 20, 1184
  * [5] Wang C et al. 2020  Micromachined accelerometers with sub-$\mu$g/$\sqrt{\rm Hz}$ noise floor: A review., Sensors , 20, 14, 4054
  * [6] Bruns T, Ripper G P and Täubner A 2014  Final report on CIPM key comparison CCAUV.V-K2, Metrologia, 51, 1A, 09002
  * [7] Qiao S et al. 2017  Final report of CCAUV.V-K3: Key comparison in the field of acceleration on the complex charge sensitivity, Metrologia, 54, 1A, 09001
  * [8] International Organization for Standardization  ISO 16063-11:1999, Methods for the calibration of vibration and shock transducers - Part 11: Primary vibration calibration by laser interferometry
  * [9] Ingerslev H, Andresen S and Winther J H 2020  Digital signal processing functions for ultra-low frequency calibrations, Acta IMEKO, 9, 5, 374
  * [10] Shimoda T, Kokuyama W and Nozato H 2020  A low-acceleration measurement using anti-vibration table with low-frequency resonance, Acta IMEKO, 9, 5, 369