1. Introduction
The basic condition for mutual comparison of different types of accelerometers [
1,
2,
3] is that the accelerometers under test have the same frequency bandwidth and it is necessary to develop the model of special filter [
4] which is the reference for such an analysis. The frequency bandwidth of the filter must be the same as the operation frequency range of the compared accelerometers.
The reference of energies at the outputs of the accelerometer and the filter is necessary because, from the point of view of the practical implementation of measurements, only the value of energy for the operation range of the accelerometer, not for the full frequency bandwidth, is interesting [
5]. Such comparison is carried out by subtracting the energy value at the outputs of both accelerometers. The difference between these energies is defined below as the reference energy and, for the purpose of determining its maximum value, it is convenient to apply the algorithm discussed in
Section 3.
The reference energy of an accelerometer corresponds to the integral-square error [
6] of the sensor related to this type of error at the output of the special filter. Developing the procedure determines the maximum value of the integral-square error (maximum reference energy) and seems to be the justified proposal for the comparison of different accelerometers.
The basis for the implementation of the task described above is the mathematical model of the accelerometer determined through its parametric identification [
7,
8,
9,
10,
11,
12], as well as the mathematical model of the special filter. The model of the accelerometer is defined by the second-order transfer function [
1,
2,
3], whereas the model of the filter should be a high-order (a range of 10 to 15 order is sufficient). Determination of the parameters of the accelerometer model should be carried out on the basis of measurement points of both components of the frequency response (the amplitude and phase) in accordance with the procedures included in the relevant international guides [
7]. In this way (e.g., by using the weighted least squares procedure (WLSP)), the accelerometer parameters and associated uncertainties can be obtained [
3,
13,
14]. Additionally, when we want to determine the associated uncertainties, it is necessary to use the Monte Carlo procedure (MCP) [
15,
16,
17,
18], whose number of trials can be assumed in advance according to the guidelines included in the JCGM (Joint Committee for Guides in Metrology)—supplement 1 [
13]. The model of the filter is expected to meet the so-called non-distortion transformation [
19], according to which the amplitude-frequency and phase-frequency responses are flat and linear, respectively, within the frequency bandwidth of interest [
20]. However, these assumptions are ideal from a practical point of view. In practice, they are not possible to meet, and one can only strive to attain them [
20], for example, by using the above-mentioned high-order filter. The degree of achievement for the approximation therefore depends on the type of filter used and has an effect on the value of the reference energy (i.e., the higher the order of the approximation, the lower value of this energy). Hence, there is the obvious necessity to use the same model of the filter within the considered comparison of the different types of accelerometers.
The maximum value of the reference energy can be determined by using the fixed-point algorithm (FPA) [
5,
21,
22,
23] with its extension by the MCP [
13,
18]. The MCP can be implemented based on the same assumptions as in the case of the procedure intended for determining the uncertainties associated with the parameters of the accelerometer model by using the WLSP. A method based on a system of complex convolution equations could also be applied instead of the above-mentioned FPA [
4]. However, the possibilities of its application are limited due to the numerical difficulties of solving a system including more than twenty-five equations [
24]. Therefore, for the needs of this paper, the first type of approach is chosen. This choice is also supported by the fact that by employing the FPA, the maximum reference energy can be obtained quickly. The FPA also allows determination of the critical case of the input signal of the accelerometer that is limited both in magnitude and in time interval [
25], and which produces this energy. Both the number of reversals of such a signal and the times at which this reversal occurs are determined. It should be emphasized here that any other signal contained in these limits generates less energy than that produced by the critical case of the magnitude-limited signal [
4].
The FPA determines, by an application of the MCP, the numerical value of the maximum reference energy in an iterative way for the ranges of variability of the model parameter of the accelerometer by the values of associated uncertainties [
26,
27]. Taking into account the uniform distribution of parameters in the above ranges, the Monte Carlo simulation (MCS) should therefore be based on a uniform pseudo-random number generator [
28].
The assessment of the impact of the uncertainties associated with the parameters of the accelerometer on the value of the absolute error is presented in studies by Tomczyk [
26,
27]. However, this assessment does not apply to full ranges of parameter variability on the values of uncertainties associated with them. Additionally, the cited assessment concerns only the values of parameters increased or decreased by the values of the uncertainties associated with them. Therefore, it was not necessary to use the MCP in such a case.
The above-mentioned issues are developed in this paper using the MCP. The block diagrams with an exact description of both the proposed MCP as well as the algorithm for determining the maximum reference energy at the output of accelerometers are discussed here. This maximum reference energy can be the basis for mutual comparison [
4] of accelerometers with the same bandwidth but produced by different companies. This comparison relates to the maximum integral-square error, which reflects an accuracy of the accelerometer. It can therefore help one choose the proper accelerometer from many types at different prices.
The practical solutions presented here are applied to the testing of the DJB A/1800/V type of accelerometer [
29].
2. Accelerometer Modeling
The model of the accelerometer is defined by using the transfer function,
where
and
are the accelerometer voltage sensitivity and damping factor, respectively, while
where
denotes the non-damped natural frequency.
Parametric identification of the accelerometer model is directed towards determining the parameters of the transfer function given by Equation (1) together with associated uncertainties [
8,
9]. The first step in this subject is the determination of the measurement points for both frequency responses. This is achieved by applying the back-to-back method commonly used in the measurement practice. The block diagram for this type of identification method is shown in
Figure 1. The reference accelerometer is used here as the comparison signal source, while the data acquisition card together with the LabVIEW software are employed for the control of the shaker and recording of the measurement data.
The sensitivity
of the tested accelerometer is determined on the basis of the relation,
where
,
, and
are the sensitivity of the reference accelerometer, the voltage output of the system for signal conditioning at the output of the accelerometer, and the voltage output of the reference, respectively.
Based on the measurement points of the frequency responses of the accelerometer, the estimates
,
and
are associated with the parameters
,
and
which occur in the transfer function given by Equation (1), as well as the associated uncertainties
,
and
, which are determined by applying the WLSP [
3,
13]. This procedure employs the MCP which utilizes the Box–Muller and the Wichmann–Hill pseudo-random number generator with normal and uniform distributions, respectively. The minimum number
of Monte Carlo trials (MCts) should be determined based on the relation,
where
denotes the assumed coverage probability [
13].
3. Algorithm for Determining the Maximum Energy (ME)
The optimization algorithm for determining the ME is based on the integral equation,
where
while
is any signal at the input of the accelerometer,
is the Lagrange multiplier [
5],
and
are the impulse responses of the accelerometer and the filter, respectively, and
denotes the time of accelerometer testing. The impulse responses
and
are determined as the inverse Laplace transforms of the corresponding transfer functions
and
. In turn, the transfer function
is obtained by way of the parametric identification of the accelerometer, while the transfer function
is convenient to be defined by the low-pass Butterworth filter given by,
In Equation (6), L is the order of the filter, is the filter’s cut-off frequency, which is equal to the frequency bandwidth of the accelerometer, and denotes the imaginary number.
The solution to Equation (4) determines the signal
that yields the maximum value of the coefficient
for which the ratio is satisfied,
where
and
are the energies of any input and output signals
and
respectively. The error
is defined by,
where
and
are the output signals of the accelerometer and filter, respectively.
Both the signal
as well as the signal
maximizing the coefficient
are limited in time
and magnitude
. Thus, we can write,
The denominator and nominator of Equation (7) are represented by,
and
Substituting Equations (9) and (10) into Equation (13) and taking into account the relationship between the impulse responses,
we finally have
When the signal
is limited in magnitude
(i.e., it has a rectangular shape), then we can consider the condition given by,
where
denotes the signum operation.
Hence, the signal maximizing the energy
is represented by,
and
is the maximum value of
defined by Equation (7).
The procedure for determining the input signal that maximizes the energy
is based on the FPA [
5] and includes the following stages:
Calculate the function based on Equation (5);
Determine the initial input signal
where
is the initial Lagrange multiplier, assumed in advance by the formula,
Determine the maximizing signal on the basis of the iterative algorithm, given by,
where
denotes the number of assumed iterations.
It should be emphasized that for each iteration step the Lagrange multiplier must be determined based on Equation (7).
Calculate the energy by using Equation (15).
The signal
that yields the maximum value of the quotient
is the one that maximizes the energy
This ME denoted by
is obtained by substituting the signal
into Equation (15). Then, we have
The energy is the accelerometer output energy, referred to as the output energy of filter, and is valid for the accelerometer frequency bandwidth.
Figure 2 shows the block diagram of the algorithm described above applied for determining the ME at the output of accelerometers. The switch
S has the position 1 for
, and it has the position 2 when
. During the particular iterations
i = 1, 2, …,
J1 of the algorithm above, the energies
and
are compared with each other, and the value of the current ME and the corresponding input signal limited in magnitude are stored in the memory. This signal can be then processed by using specialized IT tools, e.g., neural network [
30].
The final solutions of the ME and corresponding signal are delivered for . It should be noted that this algorithm converges quickly because the maximum number J of iterations does not exceed 30.
The signal producing the ME has the property that any other signal contained in its limitations can only yield less energy than the energy value caused by the maximizing signal [
4].
Figure 3 shows examples of the maximizing signal
limited simultaneously in time and magnitude and three other signals
,
and
contained in its limitations denoted by
(magnitude) and
(time).
The limited signal
in
Figure 3 has seven time switches denoted by
,
, …,
. The number of these switches and the times at which they occur are determined by the FPA, and they constitute the shape of the signal that maximizes the energy at the output of the accelerometer.
4. Monte Carlo-Based Procedure for Determining the ME
On the basis of obtaining both estimates
,
and
of the accelerometer model and the associated uncertainties
,
and
, the values of the parameters denoted by
,
and
that produce the energy
should be determined. This can be carried out for the rectangular ranges shown in
Figure 4 (the lines slanting to the right), by applying the MCP with the number of trials determined using Equation (3), where:
and
corresponds to the estimates
,
, and
respectively, in the case of three ranges shown in
Figure 4.
The trials of the MCP executed in three ranges were carried out by employing a uniform distribution. Hence, it is recommended to use the Wichmann–Hill pseudo-random number generator [
28], according to the guidelines included in the JCGM supplement 1 [
13]. The horizontal ranges of a uniform distribution in
Figure 4 are constrained by the values of uncertainty associated with the parameter estimates of the accelerometer, while the vertical ranges of this distribution are contained in the range
where
is inversely proportional to the range of the parameter estimate variability.
Figure 5 shows the block diagram of the MCP for determining the
. This procedure includes three main sub-blocks:
Back-to-back procedure for identification and determination of the accelerometer parameters and associated uncertainties.
Selection of the special filter which is the reference for calculation of the energy at the output of accelerometer.
Execution of the FPA.
The input parameter for the MCP is the number of MCts assumed in advance according to Equation (3), while the input parameters for the FPA are: number of iterations value of time quantization step size of the calculation, values of the signal limit and initial Lagrange multiplier .
5. Results
The results of testing the accelerometer of the type DJB A/1800/V are presented below. Typical parameters for the accelerometer, taken from the corresponding datasheet [
29], are as follows:
- –
voltage sensitivity (±10%): 1.02 V/(m/s2),
- –
frequency response (accelerometer bandwidth): 0.2 Hz–1 kHz ±5%,
- –
maximum continuous g level: 4903 m/s2,
- –
resonant frequency: 4 kHz.
The B&K8305 charge mode accelerometer was applied as the reference. The parameters of this accelerometer, contained in the datasheet [
31], are:
- –
charge sensitivity (±10%): 1.08 pC/g,
- –
amplitude (±10%) and phase (±1°) range: 0.2 Hz–1 kHz.
The results of the back-to-back identification are obtained based on the measurement of the amplitude and phase responses carried out for 42 frequencies within the range 30 Hz–6 kHz. This range was selected in such a way as to exceed the value of the resonant frequency since it has a significant influence on the value of the estimates and .
The measurement points of both frequency responses, obtained for the above frequencies (the first column denoted by
F), are contained in the second and third columns (denoted by
A and
) of
Table 1.
The measurement points above are registered by the NI-6221 data acquisition card with a sampling rate equal to 100 kS/s and the computer programs developed in the LabVIEW software.
The expanded uncertainties
u(
A) and
u(
) which are contained in the fourth and fifth columns of
Table 1 were determined based on
Table 1 included in the ISO 16063-21:2003(E) standard (7). These uncertainties are valid for the coverage factor
.
The accelerometer model responses (solid lines) obtained by utilizing the WLSP are shown in
Figure 6. The covariance matrices involved in the basis kernel of the WLSP used for calculating the estimates of the accelerometer model parameters and applied for determining the uncertainties associated with these estimates are determined by using the MCP. The number
of MCts, which is equal to
was taken as the minimum acceptable value obtained based on Equation (3). The coverage probability
equal to 0.95 was assumed in advance.
The expanded uncertainties,
u(
A) and
u(
), are marked in
Figure 6 by the solid vertical lines within the measurement points of both responses. The graphs of both responses are presented by using a logarithmic scale due to the wide range of considered frequencies.
The estimates
,
and
of the parameters of the accelerometer model and the associated uncertainties
,
and
determined by using the WLSP are tabulated in
Table 2.
The model responses shown in
Figure 6 are obtained by relative approximation by using the functions
and
for the parameters tabulated in
Table 2, where
and
The estimates and associated uncertainties from
Table 2 were utilized to determine the parameters
,
and
giving the maximum value of the reference energy determined by the FPA controlled by the MCP. The input parameters for the FPA were selected as:
= 30,
= 2 ms,
= 2 μs, and
= 1, and the value of signal limit
= 1.037 V was assumed as equal to the estimate of the voltage sensitivity
of the considered accelerometer. The time
was obtained from the steady state of the impulse response
calculated based on Equation (14). This steady state was assumed to be equal to 5% of the maximum deviation (peak response) from zero of this impulse response. The tenth-order (
= 10) analog Butterworth filter was adopted and is described by Equation (6). The filter cut-off frequency
is equal to the frequency response of the considered accelerometer and is 1 kHz. The values of both the gain factor
and the signal limit
are equal to 1.03 V (i.e., the value of the estimate
of the accelerometer voltage sensitivity). The upper limit of the filter order results from the computational limitation of the MathCad program which was applied for calculation of the impulse response in the way of determining the inverse Laplace transform on the basis of the relevant transfer function
.
Figure 7 shows the impulse responses
, where the maximum deviations of these responses from the steady state value are equal to 1764. The number of MCts for the procedure controlling the FPA was
similar to the WLSP. Based on parameters and uncertainties given in
Table 2, the following values were obtained:
= 1.034 V/(m/s
2);
= 1.040 V/(m/s
2);
= 0.016 V/(m/s
2);
= 4067 Hz;
= 4079 Hz;
= 12 Hz;
= 0.3178;
= 0.3186; and
= 8 for the ranges shown in
Figure 4.
The values of the parameters = 1.038 V/(m/s2), = 4070 Hz, and = 0.3178, which give the maximum reference energy = 5.862 mV·m2, are obtained by applying the MCP for = 2 ms, = 127,548, and = 15. The obtained uncertainty of the MCP is equal to 0.008 mV·m2.
The relationship between the maximum reference energy
and the time
of the considered accelerometer is linear [
5] for the steady state of the impulse response
(i.e., for
> 2 ms). Thus, it is possible to determine the functional relationship between this energy and any time
higher than 2 ms. Hence, the value of energy
should be determined for any time higher than 2 ms, and then the relevant formula can be easily determined based on the linear regression. For example, if
= 10 ms, then we have
= 34.115 mV·m
2 and
= 0.012 mV·m
2, which are obtained for
= 93,515 and
= 9. Based on the two values of the energy
and the values of time
we have:
Figure 8 shows the signal
with two reversals producing the
at the output of the considered accelerometer obtained for
= 2 ms.
The reversals of the signal are 0.68 ms and 1.49 ms.
6. Conclusions
The paper presents a proposal for the application of the FPA controlled by the MCP to determine the maximum reference energy at the output of accelerometers. Block diagrams together with a detailed description of the proposed solutions are presented.
The reference energy of the accelerometer should be understood as its output energy decreased by the output energy of the special filter (modeled by the low-pass type) with the cut-off frequency corresponding to the upper limit of the accelerometer frequency response. The lower limit of this response is neglected in this paper due to its low value (0.2 Hz).
Additionally, the MCP was used for parametric identification of both the parameters associated with the transfer function of the accelerometer and the associated uncertainties. Then, the value of the maximum reference energy was determined along with the associated uncertainties at the output of the accelerometer of the type DJB A/1800/V, considered as an example of the accelerometer. This was performed by checking the full range of the parameter variability of the accelerometer model by value-associated uncertainties. This is possible to perform only by using the MCP based on the pseudo-random number generator with a uniform distribution. This necessity results from the fact that the energy above can be obtained for the values of parameters contained within the variability ranges and not, for example, at their ends or in their middles. Only in this way it is possible to accurately determine the maximum reference energy, which can be the basis for mutual comparison of different types of accelerometers. The presentation of the relevant procedures together with a practical example giving the possibility of such comparison is the main contribution of this paper.