Evaluation of measurement uncertainty for time-dependent quantities

One of the main challenges in the analysis of dynamic measurements is the estimation of the timedependent value of the measurand and the corresponding propagation of uncertainties. In this paper we outline the design of appropriate digital compensation filters as a means of estimating the quantity of interest. For the propagation of uncertainty in the application of such digital filters we present online formulae for finite impulse response and infinite impulse response filters. We also demonstrate a recently developed efficient Monte Carlo method for uncertainty propagation in dynamic measurements which allows calculating point-wise coverage intervals in real-time.


Introduction
The "Guide to the Expression of Uncertainty in Measurement" (GUM) [1] provides a framework for harmonized evaluation and the interpretation of measurement uncertainty.The recently published Supplement 2 to the GUM extends this framework to multivariate quantities [2].However, the GUM does not address the measurement of quantities whose values depend on a continuous quantity like time, space or wavelength [3][4][5].We here mainly focus on timedependent "dynamic" measurands, although the proposed methodology is largely the same when the measurand depends on another continuous quantity such as space or wavelength.
Figure 1 Typical workflow in the analysis of a dynamic measurement where A/D denotes analogue-to-digital conversion Dynamic measurements are of growing importance in industry and metrology.Typical examples are flow measurements with a step-like change of flow, in-cylinder pressure measurements in the automotive industry, highspeed electronics or crash tests.A reliable and traceable calibration of the measurement systems employed requires a consistent evaluation of measurement uncertainty.The main challenges in the analysis of dynamic measurements are the estimation of the value of the measurand and the corresponding propagation of uncertainties.Estimation in dynamic measurements is typically an ill-posed problem and thus requires regularization to obtain a stable solution.For example, in the input estimation for the calibration of a pulse generator in [5] Tikhonov regularization in the frequency domain is advocated.However, regularization causes systematic errors which need to be accounted for in an uncertainty analysis.

Figure 2 Simulated dynamic measurement
The typical workflow in the analysis of dynamic measurements is shown in Figure 1.The continuous timedependent ܻሺ‫ݐ‬ሻǡ the value of the dynamic measurand, serves as input to the measurement device with the corresponding output ܺሺ‫ݐ‬ሻ.For instance, in an engine teststand ܻሺ‫ݐ‬ሻ could be the time-dependent torque produced by the engine whereas ܺሺ‫ݐ‬ሻ is the output of an employed torque transducer.We assume that the measurement device can be modelled as a linear time-invariant (LTI) system with the transfer function ‫ܪ‬ሺ‫ݏ‬ሻ [6].For LTI systems the relation between the system input ܻሺ‫ݐ‬ሻ and the corresponding output ܺሺ‫ݐ‬ሻ is given by a convolution [6] where ݄ሺ‫ݐ‬ሻ denotes the impulse response, i.e., the inverse Laplace transform of the transfer function.Figure 2 shows an example of a simulated dynamic measurement where the system model is an LTI system.Note that the values indicated by the measurement device, i.e., ܺሺ‫ݐ‬ሻ, are not proportional to the values of the measurand ܻሺ‫ݐ‬ሻ.This is a typical situation in dynamic measurements.Equation (1) relates the dynamic measurand ܻሺ‫ݐ‬ሻ to the dynamic observation ܺሺ‫ݐ‬ሻ.Both are continuous functions and thus the estimation of the value of the measurand is related to continuous modelling [7].Let us assume that there is a mapping ‫ܨ‬ so that the estimate ܻ ሺ‫ݐ‬ሻ of ܻሺ‫ݐ‬ሻ is calculated from the estimate ܺ ሺ‫ݐ‬ሻ by ܻ ሺ‫ݐ‬ሻ ൌ ‫ܨ‬൛ܺ ሺ‫ݐ‬ሻൟǤ Then the evaluation of uncertainties requires us to (i) assign an uncertainty to the continuous function ܺሺ‫ݐ‬ሻ and (ii) to propagate this uncertainty through the mapping ‫ܨ‬Ǥ The treatment of continuous functions is beyond the current GUM and requires the employment of a stochastic process as a model for the state of knowledge instead of a probability density function as in the GUM [7,9].An extension of the GUM methodology to continuous functions using stochastic calculus has been proposed recently, but is beyond the scope of this paper [7].Here we assume that the continuous time-dependent functions ܺሺ‫ݐ‬ሻ and ܻሺ‫ݐ‬ሻ can be represented by discrete-time sequences ൌ ሺܺ ଵ ǡ ǥ ǡ ܺ ே ሻ ் and ൌ ሺܻ ଵ ǡ ǥ ǡ ܻ ே ሻ ் with ܺ ൌ ܺሺ‫ݐ‬ ሻ and ܻ ൌ ܻሺ‫ݐ‬ ሻ for equidistant ‫ݐ‬ with a sampling interval length ܶ ௦ .This can be assumed when, for instance, the length of the sampling interval ܶ ௌ in the analogue-todigital conversion satisfies the Nyquist criterion [6].The task is then to estimate the discrete-time sequence ൌ ሺܻ ଵ ǡ ǥ ǡ ܻ ே ሻ ் , and the evaluation of uncertainties can be carried out in accordance with GUM Supplement 2 [2].Subsequently, we restrict ourselves to this case of discretetime analysis.In order to solve the ill-posed estimation problem, we consider the design of a digital filter ݃ሺ‫ݐ‬ሻ so that the filter input is the sequence of observed values , e.g., the measurement device output signal, whereas the output of the digital filter is an estimate of the (discrete-time) dynamic measurand with ઢ ൌ ሺȟ ଵ ǡ ǥ ǡ ȟ ே ሻ ் an error caused by the application of the digital filter [7].In this paper we outline the design of such digital compensation filters as a means of estimating the quantity of interest.For the propagation of uncertainty in the application of such digital filters we present sequential formulae for finite impulse response (FIR) and infinite impulse response (IIR) filters [8,9].We also demonstrate a recently developed efficient Monte Carlo method for uncertainty propagation in dynamic measurements which, in principle, allows calculating point-wise coverage intervals in real-time [10].

Compensation Filters
In this section we apply concepts of digital signal processing (DSP) to correct the output signal of the measurement device for the dynamic effects caused by the measurement device employed, see Figure 2. We assume that the device can be modelled as a linear time-invariant system [6].The goal is to design a digital compensation filter that approximates the inverse of this LTI system up to a certain frequency and attenuates measurement noise beyond that frequency [11].Therefore, we consider a cascade of an inverse filter ݃ ௩ ሺ‫ݐ‬ሻ and a low-pass filter ݃ ௪ ሺ‫ݐ‬ሻ with frequency response ‫ܩ‬ ௩ ሺ݁ ఠ் ೞ ሻ and ‫ܩ‬ ௪ ሺ݁ ఠ் ೞ ሻ, respectively [7].The compensation filter is then the cascade of the inverse filter and the low-pass filter, and its frequency response is given by with ߱ ൌ ʹߨ݂ and satisfies the following relation where ‫ܪ‬ሺ݆߱ሻ denotes the frequency response of the measurement device LTI system model and ߱ ଵ തതതത ൏ ߱ ଶ തതതത are chosen frequencies [9].The inverse filter ݃ ௩ ሺ‫ݐ‬ሻ is designed from knowledge about the measurement device.This may be available in terms of a parametric model, a sequence of frequency response values or a sequence of discrete-time impulse response values [6,11].Let us consider, for instance, that a set of frequency response values ൌ ൫‫ܪ‬ሺ݆߱ ଵ ሻǡ ǥ ǡ ‫ܪ‬ሺ݆߱ ெ ሻ൯ ் is available in terms of real and imaginary parts ൌ ሺܴ݁ሾ ் ሿǡ ‫݉ܫ‬ሾ ሿሻ ் and their associated uncertainties.This is the case when, for instance, a sensor is calibrated using sinusoidal excitations over a certain range of frequencies.Consider furthermore, that a time delay of ݊ ௗ samples is introduced for in order to improve the estimation quality [11][12][13].Hence, the time-shifted response of the measurement device with frequency response values ‫ܪ‬ ෩ ሺ߱ሻ ൌ ‫ܪ‬ሺ݆߱ሻ݁ ఠ ் ೞ is considered.A finite-impulse response (FIR) digital filter with coefficients and a z-domain transfer function can be designed as an approximation to ෩ ି ൌ ቀ‫ܪ‬ ෩ ିଵ ሺ߱ ଵ ሻǡ ǥ ǡ ‫ܪ‬ ෩ ିଵ ሺ߱ ெ ሻቁ ் , the reciprocal of the shifted frequency response values, using linear least-squares [9].Therefore, the coefficients (5) of the inverse filter are calculated as the solution of , and design matrix ൌ ሺܴ݁ሾ ் ሿǡ ‫݉ܫ‬ሾ ் ሿሻ ் , where the matrix is given by ǡ‫ז‬ ൌ ሺͳǡ ݁ ିఠ ೖ ் ೞ ǡ ǥ ǡ ݁ ିே ್ ఠ ೖ ் ೞ ሻǤ The weighting matrix may be chosen to reflect incomplete knowledge of the frequency response Ǥ The estimate of the inverse filter coefficients is thus given by provided that the applied weighting is chosen as the inverse of the uncertainty associated with the real and imaginary parts of the shifted reciprocal frequency response, i.e., ൌ ෩ ష .Note that for this approach a parametric identification of the measurement device itself is not necessary since the least-squares fit of the inverse filter is applied directly to the reciprocal of the calibrated frequency response.The low-pass filter ݃ ௪ ሺ‫ݐ‬ሻ is designed so that highfrequency noise is attenuated while low-frequency components of the estimated measurand are kept unchanged.A very small cut-off frequency reduces measurement noise strongly, but may result in a large estimation error ઢ, whereas a large cut-off results in a small error, but also in a larger noise amplification [7].The choice of the low-pass filter cut-off frequency is thus always a trade-off between noise attenuation and a small estimation error ઢ.Note that this choice depends on the spectral content of the measurand and is thus carried out for each measurement individually.

Figure 3
Amplitude of simulated measurement device frequency response (dotted), FIR compensation filter (solid) and cascade of system model and compensation filter (dashed) Figure 3 shows the frequency response amplitude of the simulated system employed for Figure 2, together with the frequency response amplitude of an FIR-type compensation filter.The low-pass filter has been designed as an FIR filter using the window technique with a Kaiser window.The cut-off frequency of the low-pass filter is at ͳʹǤͷ kHz.Note that the compensation filter is of FIR-type provided that the inverse filter and the low-pass filter are both of FIR-type.As a consequence, the compensation filter has a linear phase.For instance, the time-delay of the compensation filter in Figure 3 is ݊ ௗ ଶ ൌ ͳͳ samples for ܶ ௦ ൌ ͲǤͲͲʹ ms, where ‫ܮ‬ denotes the order of the lowpass filter and ݊ ௗ the time delay introduced for the design of the inverse filter [9].When the estimation error ઢ, due to an imperfect inverse filter and the chosen cut-off of the low-pass filter, is significant compared to the other sources of uncertainty, it has to be taken into account in the uncertainty budget.To this end, an upper bound ȁȟ ȁ ߛ may be derived from prior knowledge about the measurand.According to GUM Supplement 2 [2], the uncertainty associated with the estimation error is then given as the variance of the corresponding rectangular probability distribution, i.e., ‫ݑ‬ሺȟ ሻ ൌ ߛȀξ͵ for all time instants ‫ݐ‬ .When, for instance, prior knowledge about the amplitude of the frequency spectrum of the measurand is available in terms of a frequency function ȁܻሺ݆߱ሻȁ ൏ ‫ܤ‬ሺ߱ሻ, then [7] with ‫ܨ‬ ௌ ൌ ͳȀܶ ௌ the sampling frequency.
Uncertainties associated with the filter coefficients are then typically evaluated using the Monte Carlo method of Supplement 2 to the GUM [2].

Propagation of Uncertainties
The propagation of uncertainties for the application of a digital filter can be carried out using the GUM law of propagation of uncertainty [7][8][9] or the Monte Carlo method of GUM Supplement 2 [2,10].
The measurand is the sequence of discrete-time values ܻሺ‫ݐ‬ ሻ.The dimension of is often very high.Typical dimensions are in the order of ܰ ൌ ʹͲͲͲ.The values to be estimated evolve over time and thus methods are preferred which can operate in real-time, i.e., during the measurement.This allows, for instance, an evaluation of uncertainties in monitoring tasks or when the dimension of is too high to be stored in computer memory and postprocessed.Let us consider the measurement model with input quantities being the error ઢ, the filter coefficients ൌ ൫ܾ ǡ ǥ ǡ ܾ ே ್ ǡ ܽ ଵ ǡ ǥ ǡ ܽ ே ೌ ൯ ் and the sequence of observed values ൌ ሺܺ ଵ ǡ ǥ ǡ ܺ ே ሻ ் .Let us assume that knowledge about ȟ is modelled by the rectangular distribution ȟ ̱ܷሺെߛǡ ߛሻ for a known ߛ independent from ‫ݐ‬ .The knowledge about the filter coefficients is assumed to be modelled by the probability distribution with density function ‫‬ .The observational model is assumed to be with ߳ measurement noise modelled by a stationary zeromean auto-regressive moving average stochastic process with known covariance matrix ࣕ .Hence, according to [2] the state of knowledge about is modelled by a multivariate normal distribution ܰሺǡ ࣕ ሻ.Knowledge about the input quantities is assumed to be obtained independently.Hence, the joint distribution of the input quantities ǡ ઢ and to the measurement model ( 8) is the product of the individual distributions.Note that equation ( 8) and the methodologies outlined in this section are not limited to compensation filters, but are applicable for any kind of digital filtering.

Figure 4
Uncertainties associated with an estimate of the measurand shown in Figure 2, obtained utilizing the compensation filter with frequency response shown in Figure 3 Figure 4 shows the uncertainty associated with the estimation result for the simulated example shown in Figure 2 using the compensation filter with the frequency response amplitude shown in Figure 3.The error ȟ can be assumed to be negligible, and the noise process ߳ is modelled as zero-mean Gaussian noise with known standard deviation ߪ ൌ ͲǤͲͲͳ.The relative uncertainty associated with the inverse filter coefficients is assumed to be about ʹΨ whereas the uncertainty associated with the low-pass filter coefficients is assumed to be zero.Note that the uncertainty shown in Figure 4 shows a significant timedependence.This is typical of dynamic measurements.

Formula for FIR filters
Let denote the covariance matrix of the probability distribution which models knowledge about the compensation filter coefficients and let denote the local (auto-)covariance matrix of the discretetime observations at the time instants ‫ݐ‬ ିே ್ ǡ ǥ ǡ ‫ݐ‬ with If the digital compensation filter is of the FIR-type, i.e., ܽ ൌ Ͳ for all ݉, the uncertainty associated with the estimate ‫ݕ‬ ො of ܻሺ‫ݐ‬ ሻ is given by [7,9] where ‫ݎܶ‬ denotes the trace of a matrix.The trace term accounts for the non-linearity of the measurement model and is typically small compared to the other components.

Formula for IIR filters
For the evaluation of uncertainty for the application of an IIR filter, equation ( 8) is transformed into the corresponding state-space equation system [8] EPJ Web of Conferences where ǡ and are calculated from the IIR filter coefficients and ࢠሾ݊ሿ denotes the vector of system state variables at a time instant ‫ݐ‬ .
The uncertainty associated with the estimate ‫ݕ‬ ො obtained by the application of an IIR-type filter for the case of uncorrelated measurement noise, i.e., ܸ ఢ diagonal, is given by [7,8] with Note that the formulae ( 10) and ( 12) for the evaluation of uncertainty allow a real-time calculation of uncertainties.
That is, the uncertainty can, in principle, be calculated during the measurement.

The GUM-S2 Monte Carlo method
The measurement model ( 8) is non-linear in its input quantities which holds in particular for the IIR filter case.Hence, a linearized propagation of uncertainty as in equation ( 12) may be unreliable.Moreover, the uncertainty formulae ( 10) and ( 12) do not allow a direct calculation of (point-wise) coverage intervals as the probability distribution associated with the estimate of the measurand is unknown [2].To this end, the Monte Carlo method for the propagation of probability density functions as described in GUM Supplement 2 can be employed [10].
Realizations of the input quantities are thus drawn from their associated probability distributions and propagated through the measurement model ( 8) to obtain realizations of the value of the measurand.In more detail, the following steps are carried out for ݇ ൌ ͳǡ ǥ ǡ ‫,ܭ‬ where ‫ܭ‬ denotes the number of Monte Carlo trials.
1. Draw realizations ሺሻ ǡ ሺሻ ǡ ઢ ሺሻ of the input quantities from their associated probability distributions.2. Calculate the corresponding realization ሺሻ of the estimate of the measurand from model (8).
The set ൛ ሺሻ ൟ then represents samples from the probability distribution with density ‫‬ , which models knowledge about the value of the measurand.The estimate of the measurand is then calculated as the mean over all ‫ܭ‬ realizations and its associated uncertainty as the covariance In addition, point-wise or simultaneous coverage intervals ‫ܫ‬ ൌ ቂ‫ݕ‬ ሺ௪ሻ ǡ ‫ݕ‬ ሺሻ ቃ to achieve some prescribed coverage probability ܲ (e.g.95%) ሺೢሻ can be calculated as described in GUM-S2 [2].
As mentioned above, the dimension of and is often very high -in the order of several thousand samples.This limits the straightforward application of the above Monte Carlo procedure as it requires the storage of the complete set ൛ ሺሻ ൟ of realizations, which has dimensions ‫ܭ‬ ൈ ܰ.
The value of ‫,ܭ‬ the number of Monte Carlo trials, depends on the sought accuracy of the Monte Carlo result.That means, since the Monte Carlo method is random in nature, the results differ when the whole procedure is repeated.
The number of digits that remain constant from one repetition to another depends on the number of Monte Carlo trials in each repetition.Hence, a large ‫ܭ‬ is required in order to achieve a reliable estimate and uncertainty.However, a large value of ‫ܭ‬ is not always achievable using the above straightforward Monte Carlo implementation.For instance, when ܰ, the length of , is about ʹͲͲͲ, the number of Monte Carlo trials ‫ܭ‬ on a PC with a 32-bit operating system is limited to about ʹ ൈ ͳͲ ସ .Even on a workstation computer with 24 cores and 96 GB of RAM running a 64-bit Linux operating system with MATLAB, the number of trials for ܰ ൌ ͵ͲͲͲ is limited to about ͳǤʹ ൈ ͳͲ .Hence, there is a need for more efficient implementations.
When only the mean (14) and the covariance (15) are sought, multivariate update formulae can be employed.Let ‫ܭ‬ denote an initial number of Monte Carlo trials, and let ෝ and ܷ ෝ బ denote the corresponding mean and covariance, respectively.After an additional number of ‫ܭ‬ ଵ Monte Carlo trials, the mean and covariance are then calculated as [10] ෝ respectively.This can then be repeated for an additional ‫ܭ‬ ଶ ǡ ǥ ‫ܭ‬ ெ number of Monte Carlo trials.The numbers ‫ܭ‬ can be kept small, and hence, the calculation of mean and covariance can be carried out for arbitrarily large ‫ܭ‬ ൌ σ‫ܭ‬ on standard PCs.
For the calculation of point-wise coverage intervals we propose a recently developed sequential Monte Carlo method [10].The main idea is to utilize the sequential character of the measurement model (8).That is, the estimate of the value of the measurand at one time instant is calculated from the previous estimates and previous observations.Thus, estimation evolves sequentially over time.In the same way, knowledge about the (point-wise) knowledge of the measurand evolves sequentially.For the sequential Monte Carlo method the following steps are carried out.to be stored which are required to calculate the estimate at one time instant from the previous ones.This reduces the requirement for computer memory significantly and allows highly accurate Monte Carlo results to be achieved even on standard PCs [10].Moreover, the sequential implementation allows an on-line evaluation of mean, variance and point-wise coverage intervals.However, the drawback of this approach is that covariances are only available in a certain time window.

Example
To illustrate the proposed framework for the analysis of dynamic measurements we applied input estimation and evaluation of uncertainties for the analysis of an acceleration measurement with impulse-like input [15].The employed sensor was of type 2270 from manufacturer Endevco®.Its dynamic characterization has been carried out in terms of sinusoidal calibration experiments at selected frequencies in the range from 500 Hz to 20 kHz, see Figure 5.In order to carry out input estimation in terms of deconvolution, an FIR-type filter of order 12 with additional time delay of 6 samples was determined via ordinary least-squares (7), see Figure 6.Note that the deconvolution filter was determined without knowledge of a parametric characterization of the transducer [11].
The uncertainties associated with the FIR filter coefficients given by their covariance matrix are shown in Figure 7.The estimate of the sensor input was then determined by cascaded application of the low-pass filter and the FIR deconvolution filter.Uncertainties were evaluated using equation (10), see Figure 9.In addition, 95% point-wise coverage intervals were calculated using the sequential Monte Carlo approach, see Figure 10.

EPJ Web of Conferences
As can be seen from the Figures 9 and 10, the deconvolution results in a good estimate of the actual input acceleration.

Summary
The analysis of dynamic measurements and in particular the corresponding evaluation of uncertainties is a topic of growing importance in metrology.However, the current guidelines, GUM and its supplements, do not address such measurement tasks.Strictly speaking, an extension of the GUM to include the treatment of continuous functions by using stochastic calculus would be required.Here we assumed that the continuous functions can be represented by discrete-time sequences.This allows an evaluation of uncertainties in accordance with GUM Supplement 2. We outlined several recently developed approaches to deal with the corresponding discrete-time dynamic measurements when the measurement device can be modelled by a linear time-invariant system.We addressed the design of digital filters as a means to estimate the value of the dynamic measurand.For the evaluation of measurement uncertainty for measurement models that involve digital filtering, we presented several approaches which are in accordance with GUM.The proposed procedures for uncertainty evaluation can, in principle, be applied in real-time during the measurement.

1 .
Draw ‫ܭ‬ realizations of the filter coefficients ሺሻ , the observations ሺሻ at time instant ‫ݐ‬ ଵ and the error ઢ ଵ ሺሻ at time instant ‫ݐ‬ ଵ .2. Calculate the ‫ܭ‬ realizations ෝ ଵ ሺሻ of the estimate at time instant ‫ݐ‬ ଵ using equation (8). 3. Draw ‫ܭ‬ realizations of the observations ሺሻ and the error ઢ ଶ ሺሻ at time instant ‫ݐ‬ ଶ .4. Calculate the ‫ܭ‬ realizations ෝ ଶ ሺሻ of the estimate at time instant ‫ݐ‬ ଶ using equation (8).Steps 3 and 4 are then repeated for the remaining time instants.Note that realizations of the filter coefficients are drawn only once at the initial time instant.The calculation of the ‫ܭ‬ realizations of the estimate in steps 2 and 4 using model (8) can be implemented very efficiently using matrix-vector operations.An implementation in MATLAB can be found in [14].At each time instant ‫ݐ‬ the estimate and its associated uncertainty are calculated as ‫ݕ‬ ො ൌ ݉݁ܽ݊ቄ‫ݕ‬ ሺሻ ቅ and ‫ݑ‬ ଶ ሺ‫ݕ‬ ො ሻ ൌ ‫ݕ‪ቄ‬ݎܽݒ‬ ሺሻ ቅǡ respectively.In addition, point-wise coverage intervals ‫ܫ‬ can be calculated from the set of realizations ቄ‫ݕ‬ ሺሻ ቅ.The calculation of these local characteristics (mean, variance and point-wise coverage intervals) does not require the storage of the whole ‫ܭ‬ ൈ ܰ dimensional set of multivariate estimation results ൛ ሺሻ ൟ.Instead, only those values have

Figure 5
Figure 5 Calibrated amplitude and phase response of considered accelerometer Endevco®, type 2270

Figure 6
Figure 6 Reciprocal of calibrated real and imaginary part of sensor's frequency response with associated uncertainties and corresponding FIR filter fit

Figure 7
Figure 7 Uncertainties associated with the FIR deconvolution filter coefficients

Figure 8
Figure 8 Sensor output signal for measured input acceleration

Figure 9
Figure 9 Input acceleration, deconvolution result and associated standard uncertainties

Figure 10
Figure 10 Point-wise 95% coverage intervals of estimated input acceleration and measured input acceleration