Estimation of the Resultant Expanded Uncertainty of the Output Quantities of the Measurement Chain Using the Discrete Wavelet Transform Algorithm

Featured Application

Estimation of expanded uncertainty in case of measurement chains containing wavelet transform algorithms without using Monte-Carlo simulations.

Abstract

This paper discusses the role of the discrete wavelet transform algorithm in processing error signals present in the input quantities of the algorithm. In considering the error model of the measurement chain, the parameters of the error signals in the input quantities of the wavelet transform algorithm are estimated. Subsequently, in accounting for the algorithm’s properties, the parameters of its output values are determined, and the resulting uncertainty values of the output quantities of the measurement chain are estimated. The interval reduction arithmetic method is employed in the calculations for estimating the expanded uncertainty. All findings were validated through measurements conducted using the implemented measurement chain.

1. Introduction

Currently, extensive applications of wavelet transform algorithms (WT) are well documented in the literature [1,2,3,4,5]. The algorithms discussed are used, for example, in medicine [6], image and sound processing [7], machine diagnostics [8], seismic vibration analysis [9], metal analysis [10], and even in the case of detecting leaks in pipelines [11]. However, there are few works that provide a universal and accessible method for measurement chain designers to estimate the uncertainty associated with the output values of these algorithms. The methods outlined in current works [12,13,14,15] are complex, requiring the measurement chain designer to possess in-depth knowledge of the algorithm being utilized.

The universal method for analyzing the metrological properties of WT algorithms was previously introduced in the [16], where the propagation of random error signals using these algorithms was discussed. Consequently, the issue of how these algorithms propagate deterministic error signals and their role in introducing self-error signals into the output values necessitates further examination. The method outlined in this paper, akin to that in the previous study, employs a matrix representation of the algorithm. Alongside the identification algorithm, which is elaborated on in [16,17,18], an analytical approach for determining the coefficients of this matrix will be presented.

The main goal of the work was to provide a quantitative description of how error signals present in the input quantities of the WT algorithm are transferred to its output and how the parameters of these signals change. This manuscript is divided into 6 sections. Section 1 contains an introduction and presents the most important assumptions of the work. Section 2 is devoted to the transmittance of the WT algorithm, explains how the quantity in question affects the error signals processed by the algorithm, and presents an analytical example of determining this quantity for selected algorithm parameters. Section 3 is devoted to the algorithm’s own errors, explaining their origins and describing the algorithm for identifying their parameters. Section 4 contains a description of the method used in this work to determine the resultant expanded uncertainty value. Section 5 presents an example of the application of the discussed analysis method and summarizes the results of a measurement experiment aimed at verifying its correctness. Section 6 contains the most important conclusions from the work.

All considerations presented in the work were verified through measurements using a previously constructed measurement chain. This article does not include a metrological analysis of the discussed measurement chain regarding the error model and metrological parameters of the input values of the WT algorithm. This aspect is addressed separately and has been detailed in another work [19]. The uncertainty budget of the signal processed by the WT algorithm has already been established and will serve as a reference for the metrological properties of this signal. Therefore, this work, in conjunction with the previous one [16,19], can offer guidance to the designers of measurement chains on how to implement the presented considerations in their own applications of measurement chains.

According to the division introduced in some works [16,19,20], the error signals can be divided according to the nature of the implementation, and we distinguish the following:

• Static signals, where subsequent values are constant within a single measurement window;
• Dynamic signals, where subsequent values change within a single measurement window, and it is possible to deterministically describe the course of these signals;
• Random signals, where subsequent values change within a single measurement window, and it is not possible to deterministically describe the course of these signals;

due to the origin of the signals, we distinguish the following as well:

• Own signals, introduced by the analyzed object, resulting from its imperfect properties;
• Propagated signals transferred from the input to the output of an object, present in the input quantities.

2. Wavelet Transform Algorithm Transmittance

As described in [16,17,18], many data processing algorithms can be represented in matrix form. In denoting the successive input quantities of the algorithm as $x\left(i\right)$ and the output quantities as $X\left(j\right)$, their relationships can be expressed in the following form:

$$\left[\begin{array}{c}X\left(0\right)\\ X\left(1\right)\\ ⋮\\ X\left(M-1\right)\end{array}\right]=\left[\begin{array}{cccc}{a}_{0,0}& a_{0,1}& \cdots & a_{0,N-1}\\ a_{1,0}& \ddots & & a_{1,N-1}\\ ⋮& & \ddots & ⋮\\ a_{M-1,0}& \cdots & \cdots & a_{M-1,N-1}\end{array}\right]\left[\begin{array}{c}x\left(0\right)\\ x\left(1\right)\\ ⋮\\ x\left(N-1\right)\end{array}\right],$$

(1)

where N is the number of input quantities, M is the number of output quantities, and the symbol $a_{i,j}$ denotes the coefficients of the transformation matrix $\mathbf{A}$ of the algorithm, obtained according to the method described in [18] or determined analytically in accordance with the assumptions of the analyzed algorithm, as in [21]. The method discussed was previously utilized in papers [16,22] that described the transfer of random errors using discrete wavelet transform algorithms.

When analyzing how the discussed algorithm transfers errors of a deterministic nature, the transmittance of this algorithm should be taken into account. This algorithm will act as a filter and thus modify the spectrum of the processed signal. This work assumed that the transmittance of the algorithm is ideal, which means that the algorithm does not introduce any deterministic, own errors, but only transfers the errors present in the signal (in case the actual transmittance of the algorithm differs from the ideal one, an additional component of the self-error signal should be considered, as proposed in work [19] in the attached example). Additionally, the rounding errors introduced by the algorithm during the multiplication and addition operations are taken into account, as presented later in this paper. In general, the analysis discussed should be carried out individually for each output quantity of the algorithm. However, as described in [16,22], the values for successive rows of the transformation matrix related to the same level of decomposition are only shifted toward each other, and therefore, their impact on the transmission of the error signal remains constant. Therefore, an analysis can be performed for each level of signal decomposition and not for each output quantity separately, as demonstrated in this paper.

2.1. Relationship between Transmittance and Matrix Form of Algorithm

Based on Equation (1), the transmittance in the $\mathbb{Z}$ domain can be determined for the selected row of the transformation matrix in the form

$$H_i\left(z\right)=a_{i,0}+a_{i,1}{z}^{-1}+\dots +a_{i,N-1}{z}^{-N+1}=\sum _{k=0}^{N-1}{a}_{i,k}{z}^{-k},$$

(2)

and after substituting $z=e^{\mathrm{j}\omega T_s}$, where $T_s=\frac{1}{f_s}$ is the sampling perioid corresponding to the sampling frequency $f_s$, the dependence describing the transmittance of the selected i-th row in the pulsation domain denoted as $G_i\left(\mathrm{j}\omega \right)=H_i\left(e^{\mathrm{j}\omega T_s}\right)$ is obtained. Based on the determined transmittance, the amplification $K_i\left(\omega \right)$ and phase shift ${\phi }_i\left(\omega \right)$ for the selected harmonic of the signal can be described in the following way:

$$\begin{array}{c}{K}_i\left(\omega \right)=\left|G_i\left(\mathrm{j}\omega \right)\right|=\sqrt{{\left(\Re \left(G_i\left(\mathrm{j}\omega \right)\right)\right)}^2+{\left(\Im \left(G_i\left(\mathrm{j}\omega \right)\right)\right)}^2},\\ {\phi }_i\left(\omega \right)=arctan\left(\frac{\Im \left(G_i\left(\mathrm{j}\omega \right)\right)}{\Re \left(G_i\left(\mathrm{j}\omega \right)\right)}\right).\end{array}$$

The transmittance of the object being discussed will impact the error signal components in a manner that alters their amplitude and phase.

2.2. Algorithm Transmittance Impact on Error Signals

Assume that the signal $x\left(t\right)$, processed by the analyzed WT algorithm, can be described in an ideal case as follows:

$$\dot{x}\left(i\right)=\sum _{j=0}^{\infty }{E}_{x,o}\left({\omega }_j\right)sin\left({\omega }_{j}i{T}_s+{\phi }_{x,o}\left({\omega }_j\right)\right),$$

where $T_s=\frac{1}{f_s}$ is the sampling period corresponding to sampling frequency $f_s$, $E_{x,o}\left(\omega \right)$ is the amplitude, and ${\phi }_{x,o}\left(\omega \right)$ is the phase shift of the signal harmonic of pulsation $\omega$. In the real case, the quantity $x\left(i\right)$ can be described as follows:

$$\tilde{x}\left(i\right)=\dot{x}\left(i\right)+e_{x,r}\left(i\right)+\sum _{j=0}^{\infty }{E}_{x,e}\left({\omega }_j\right)sin\left(i{T}_s{\omega }_j+{\phi }_{x,e}\left({\omega }_j\right)\right),$$

(3)

where $e_{x,r}\left(i\right)$ is random error (non-deterministic), $E_{x,e}\left(\omega \right)$ is the amplitude, and ${\phi }_{x,e}\left(\omega \right)$ is the phase shift of the deterministic error.

According to Equation (3), it is possible to define static error signal $e_{x,s}\left(i\right)$ (for which subsequent values do not change for a single WT algorithm run) as follows:

$$e_{x,s}\left(i\right)=E_{x,e}\left(0\right)sin\left({\phi }_{x,e}\left(0\right)\right),$$

where in case of ${\phi }_{x,e}\left(0\right)=\pi \mathrm{rad}$, the equation in question can be written as $e_{x,s}\left(i\right)=E_{x,e}\left(0\right)$. In case of harmonics with non-zero pulsation, the dynamic error signal $e_{x,d}\left(i\right)$ can be defined as follows:

$$e_{x,d}\left(i\right)=\sum _{j=1}^{\infty }{E}_{x,e}\left({\omega }_j\right)sin\left(i{T}_s{\omega }_j+{\phi }_{x,e}\left({\omega }_j\right)\right).$$

The random error signal $e_{x,r}\left(i\right)$ cannot be described in a deterministic form, but statistical parameters such as variance, expected value, and the shape of the probability density function can describe its properties.

According to [5,23,24], a single row of the transmission matrix, as indicated in Equation (1), and the related transmittance described in Equation (2), can be analyzed in a manner analogous to a finite impulse response (FIR) filter [25]. For each output quantity, it is possible to analyze how the algorithm introduces error signals present at its input to this quantity. It should also be noted that in the case of the discussed family of algorithms, the transfer function associated with the selected quantity is linear and time invariant. Therefore, the deterministic error component of the input quantity, as described in Equation (3), is transferred to the i-th output of the WT algorithm according to the following relationship:

$$\begin{array}{c}{e}_{X,s,i}\left(j\right)=K_i\left(0\right)E_{x,e}\left(0\right)sin\left({\phi }_{x,e}\left(0\right)+{\phi }_i\left(0\right)\right),\\ e_{x,d,i}\left(j\right)=\sum _{k=1}^{\infty }{K}_i\left({\omega }_k\right)E_{x,e}\left({\omega }_k\right)sin\left(j{T}_s{\omega }_k+{\phi }_{x,e}\left({\omega }_k\right)+{\phi }_i\left({\omega }_k\right)\right).\end{array}$$

The discussed transmittance also affects the variance of the analyzed error signals. In the case of both deterministic and random signals, the following relationship can be written [23,26]:

$${\sigma }_{X,i}^2\left(\omega \right)=K_i^2\left(\omega \right){\sigma }_x^2\left(\omega \right)={\sigma }_x^2\left(\omega \right){\left|H_i\left(e^{\mathrm{j}\omega T_s}\right)\right|}^2,$$

(4)

where for analyzed error signal ${\sigma }_x^2$ is the variance in the WT algorithm input, and ${\sigma }_{X,i}^2$ is the variance in the i-th WT algorithm output. Given the above equation, the single harmonic variance of the dynamic error signal $e_{X,d,i}\left(j\right)$ in the algorithm output can be determined according to the relationship [24]

$${\sigma }_{X,d,i}^2\left(\omega \right)=\frac{1}{2}{E}_{x,e}^2\left(\omega \right){\left|H_i\left(e^{\mathrm{j}\omega T_s}\right)\right|}^2.$$

(5)

In the case of static error signals, the presented equation simplifies (for $\omega =0 \mathrm{rad/s}$, it becomes $e^{\mathrm{j}\omega T_s}=1$), so it can be written as follows:

$${\sigma }_{X,s,i}^2={\left|H_i\left(1\right)\right|}^2{\sigma }_{x,s}^2={\sigma }_{x,s}^2{\left(\sum _{j=0}^{N-1}{a}_{i,j}\right)}^2.$$

In the case of random error signals, it can be observed that the algorithm in question processes N subsequent realizations of these signals. If these realizations are uncorrelated with each other, or the autocorrelation mentioned is small, and the assumption of the same power spectral density of these signals in the frequency range $\widehat{f}\in [0; \frac{1}{2}{f}_s]$ holds, the variance of the random error signal at the algorithm’s output can be described by the following equation:

$${\sigma }_{X,r,i}^2=a_{i,0}^2{\sigma }_{x,r}^2+a_{i,1}^2{\sigma }_{x,r}^2+\dots +a_{i,N-1}^2{\sigma }_{x,r}^2={\sigma }_{x,r}^2\sum _{j=0}^{N-1}{a}_{i,j}^2.$$

(6)

If the specified conditions are not satisfied for the analyzed signal, its variance should be calculated using Equation (4). In this scenario, it is also feasible to calculate the average variance of the analyzed random error signal:

$${\sigma }_{X,r,i}^2=\frac{1}{\pi }{\int }_0^{\pi }{\sigma }_{x,r}^2\left(\frac{{\omega }_n}{T_s}\right){\left|H_i\left(e^{\mathrm{j}{\omega }_n}\right)\right|}^2\mathrm{d}{\omega }_n,$$

(7)

where ${\omega }_n=\omega T_s$ is the normalized pulsation [24]. Note that Equations (6) and (7) are equivalent when the $e_{x,r}\left(i\right)$ signal has a constant power spectral density (in case where ${\sigma }_{x,r}^2\left(\omega \right)=\mathrm{const}$), which occurs, i.e., in case of white noise or in case of a quantization error signal [27,28,29,30]. In the case of non-constant power spectral density, only Equations (4) and (7) can be used to calculate the analyzed signal variance correctly.

2.3. Identification of Algorithm Transmittance

The determination of the algorithm’s transmittance concerning the subsequent output quantities can be carried out using Equation (2) or by understanding the properties of the wavelet utilized by the analyzed algorithm. This process can also be executed for an existing implementation of the algorithm (e.g., implemented in MATLAB [31], GNU Octave [32], or the PyWavelets package [33]), employing a suitable identification algorithm [18]. Since the mentioned identification algorithm has been utilized in previous studies [16,22], and its application is straightforward, the analytical approach will be elaborated below. The illustration pertains to the “Daubechies” [21] wavelet family, assuming that the WT algorithm will process $N=8$ input quantities to produce $M=8$ output quantities, utilizing the “db2” wavelet and undergoing $K=2$ iterations of the signal decomposition process.

In the case under analysis, the vector of input quantities, previously mentioned in Equation (1), can be represented as [5]

$${\mathbf{X}}^T=\left[\begin{array}{cccccccc}{S}_{2,0}& S_{2,1}& T_{2,0}& T_{2,1}& T_{1,0}& T_{1,1}& T_{1,2}& T_{1,3}\end{array}\right],$$

(8)

where the symbol $S_{m,n}$ denotes the approximations, and the symbol $T_{m,n}$ denotes the details of the signal for the scale number m and the time shift number n, while the vector of output quantities can be described in the form

$${\mathbf{x}}^T=\left[\begin{array}{cccccccc}{S}_{0,0}& S_{0,1}& S_{0,2}& S_{0,3}& S_{0,4}& S_{0,5}& S_{0,6}& S_{0,7}\end{array}\right].$$

(9)

The signal decomposition process, enabling the calculation of the values of the quantities referenced in Equation (8) from the values of the quantities specified in Equation (9), can be recursively defined using the equations outlined in [5]:

$$\begin{array}{c}{S}_{m+1,n}=\frac{1}{\sqrt{2}}\sum _{k=0}^{N_k-1}{c}_k{S}_{m,2n+k},\end{array}$$

(10)

$$\begin{array}{c}{T}_{m+1,n}=\frac{1}{\sqrt{2}}\sum _{k=0}^{N_k-1}{b}_k{S}_{m,2n+k}.\end{array}$$

(11)

where the symbol $c_k$ represents consecutive non-zero scaling factors, which will be elaborated on later in this paper, while the symbol $b_k$ signifies the coefficients computed in accordance with the following equation:

$$b_k={\left(-1\right)}^k{c}_{N_k-k-1},$$

(12)

where $N_k$ is the number of non-zero scaling factors resulting from the properties of the wavelet used.

According to Equations (10) and (11), to determine the values of the quantities outlined in Equation (8), one must have knowledge of the number and values of consecutive scaling factors. In considering the assumptions regarding the analyzed wavelet “db2”, this wavelet possesses $N_k=4$ non-zero scaling factors: $c_0$, $c_1$, $c_2$, and $c_3$ [5,21]. Furthermore, it is characterized by the following assumptions [5,21]:

$$\begin{array}{c}\sum _{k=0}^{N_k-1}{c}_k=2,\end{array}$$

(13)

$$\begin{array}{c}\sum _{k=0}^{N_k-1}{\left(-1\right)}^k{c}_k{k}^m=0,\end{array}$$

(14)

$$\begin{array}{c}\sum _{k=0}^{N_k-1}{c}_k{c}_{k+2k^{\prime }}=\left\{\begin{array}{cc}\hfill 2\hfill & \mathrm{where} k^{\prime }=0\hfill \\ \hfill 0\hfill & \mathrm{in} \mathrm{other} \mathrm{cases}\hfill \end{array}\right.,\end{array}$$

(15)

where in Equation (14), $m\in [0; \frac{N_k}{2}-1]$, and in Equation (15), $k^{\prime }\in \mathbb{N}$. Based on the above assumptions, the system of equations can be written as follows:

$$\left\{\begin{array}{cc}\hfill c_0+c_1+c_2+c_3=2\hfill & \mathrm{according} \mathrm{to}\left(13\right)\hfill \\ \hfill c_0-c_1+c_2-c_3=0\hfill & \mathrm{according} \mathrm{to}\left(14\right) \mathrm{for} m=0\hfill \\ \hfill -1c_1+2c_2-3c_3=0\hfill & \mathrm{according} \mathrm{to}\left(14\right) \mathrm{for} m=1\hfill \\ \hfill c_0^2+c_1^2+c_2^2+c_3^2=2\hfill & \mathrm{according} \mathrm{to}\left(15\right)\hfill \end{array}\right.,$$

and solving this system allows the values of the analyzed scaling factors to be determined:

$$c_0=\frac{1+\sqrt{3}}{4}, c_1=\frac{3+\sqrt{3}}{4}, c_2=\frac{3-\sqrt{3}}{4}, c_3=\frac{1-\sqrt{3}}{4}.$$

(16)

Based on Equations (10) and (11), relationships describing the output values of the analyzed algorithm, indicated in Equation (8), are obtained:

$$\begin{array}{c}\begin{array}{cc}\hfill S_{2,0}=& \frac{1}{\sqrt{2}}\left(c_0{S}_{1,0}+c_1{S}_{1,1}+c_2{S}_{1,2}+c_3{S}_{1,3}\right)=\hfill \\ & \frac{1}{\sqrt{2}}\frac{c_0}{\sqrt{2}}\left(c_0{S}_{0,0}+c_1{S}_{0,1}+c_2{S}_{0,2}+c_3{S}_{0,3}\right)+\hfill \\ & \frac{1}{\sqrt{2}}\frac{c_1}{\sqrt{2}}\left(c_0{S}_{0,2}+c_1{S}_{0,3}+c_2{S}_{0,4}+c_3{S}_{0,5}\right)+\hfill \\ & \frac{1}{\sqrt{2}}\frac{c_2}{\sqrt{2}}\left(c_0{S}_{0,4}+c_1{S}_{0,5}+c_2{S}_{0,6}+c_3{S}_{0,7}\right)+\hfill \\ & \frac{1}{\sqrt{2}}\frac{c_3}{\sqrt{2}}\left(c_0{S}_{0,6}+c_1{S}_{0,7}+c_2{S}_{0,0}+c_3{S}_{0,1}\right)\hfill \end{array},\end{array}$$

(17)

$$\begin{array}{c}\begin{array}{cc}\hfill S_{2,1}=& \frac{1}{\sqrt{2}}\left(c_0{S}_{1,2}+c_1{S}_{1,3}+c_2{S}_{1,4}+c_3{S}_{1,5}\right)=\hfill \\ & \frac{1}{\sqrt{2}}\frac{c_0}{\sqrt{2}}\left(c_0{S}_{0,4}+c_1{S}_{0,5}+c_2{S}_{0,6}+c_3{S}_{0,7}\right)+\hfill \\ & \frac{1}{\sqrt{2}}\frac{c_1}{\sqrt{2}}\left(c_0{S}_{0,6}+c_1{S}_{0,7}+c_2{S}_{0,0}+c_3{S}_{0,1}\right)+\hfill \\ & \frac{1}{\sqrt{2}}\frac{c_2}{\sqrt{2}}\left(c_0{S}_{0,0}+c_1{S}_{0,1}+c_2{S}_{0,2}+c_3{S}_{0,3}\right)+\hfill \\ & \frac{1}{\sqrt{2}}\frac{c_3}{\sqrt{2}}\left(c_0{S}_{0,2}+c_1{S}_{0,3}+c_2{S}_{0,4}+c_3{S}_{0,5}\right)\hfill \end{array},\end{array}$$

(18)

$$\begin{array}{c}\begin{array}{cc}\hfill T_{2,0}=& \frac{1}{\sqrt{2}}\left(c_3{S}_{1,0}-c_2{S}_{1,1}+c_1{S}_{1,2}-c_0{S}_{1,3}\right)=\hfill \\ & \frac{1}{\sqrt{2}}\frac{c_3}{\sqrt{2}}\left(c_0{S}_{0,0}+c_1{S}_{0,1}+c_2{S}_{0,2}+c_3{S}_{0,3}\right)-\hfill \\ & \frac{1}{\sqrt{2}}\frac{c_2}{\sqrt{2}}\left(c_0{S}_{0,2}+c_1{S}_{0,3}+c_2{S}_{0,4}+c_3{S}_{0,5}\right)+\hfill \\ & \frac{1}{\sqrt{2}}\frac{c_1}{\sqrt{2}}\left(c_0{S}_{0,4}+c_1{S}_{0,5}+c_2{S}_{0,6}+c_3{S}_{0,7}\right)-\hfill \\ & \frac{1}{\sqrt{2}}\frac{c_0}{\sqrt{2}}\left(c_0{S}_{0,6}+c_1{S}_{0,7}+c_2{S}_{0,0}+c_3{S}_{0,1}\right)\hfill \end{array},\end{array}$$

(19)

$$\begin{array}{c}\begin{array}{cc}\hfill T_{2,1}=& \frac{1}{\sqrt{2}}\left(c_3{S}_{1,2}-c_2{S}_{1,3}+c_1{S}_{1,4}-c_0{S}_{1,5}\right)=\hfill \\ & \frac{1}{\sqrt{2}}\frac{c_3}{\sqrt{2}}\left(c_0{S}_{0,4}+c_1{S}_{0,5}+c_2{S}_{0,6}+c_3{S}_{0,7}\right)-\hfill \\ & \frac{1}{\sqrt{2}}\frac{c_2}{\sqrt{2}}\left(c_0{S}_{0,6}+c_1{S}_{0,7}+c_2{S}_{0,0}+c_3{S}_{0,1}\right)+\hfill \\ & \frac{1}{\sqrt{2}}\frac{c_1}{\sqrt{2}}\left(c_0{S}_{0,0}+c_1{S}_{0,1}+c_2{S}_{0,2}+c_3{S}_{0,3}\right)-\hfill \\ & \frac{1}{\sqrt{2}}\frac{c_0}{\sqrt{2}}\left(c_0{S}_{0,2}+c_1{S}_{0,3}+c_2{S}_{0,4}+c_3{S}_{0,5}\right)\hfill \end{array},\end{array}$$

(20)

$$\begin{array}{c}{T}_{1,0}=\frac{1}{\sqrt{2}}\left(c_3{S}_{0,0}-c_2{S}_{0,1}+c_1{S}_{0,2}-c_0{S}_{0,3}\right),\end{array}$$

(21)

$$\begin{array}{c}{T}_{1,1}=\frac{1}{\sqrt{2}}\left(c_3{S}_{0,2}-c_2{S}_{0,3}+c_1{S}_{0,4}-c_0{S}_{0,5}\right),\end{array}$$

(22)

$$\begin{array}{c}{T}_{1,2}=\frac{1}{\sqrt{2}}\left(c_3{S}_{0,4}-c_2{S}_{0,5}+c_1{S}_{0,6}-c_0{S}_{0,7}\right),\end{array}$$

(23)

$$\begin{array}{c}{T}_{1,3}=\frac{1}{\sqrt{2}}\left(c_3{S}_{0,6}-c_2{S}_{0,7}+c_1{S}_{0,0}-c_0{S}_{0,1}\right).\end{array}$$

(24)

With Equation (1) and taking into account the order of the elements of the vector of the output quantities consistent with that assumed in Equation (8) and substituting on the basis of Equation (9) $S_{0,i}=x\left(i\right)=x_i$, we have

$$\begin{array}{c}{S}_{2,0}=a_{0,0}{x}_0+a_{0,1}{x}_1+a_{0,2}{x}_2+a_{0,3}{x}_3+a_{0,4}{x}_4+a_{0,5}{x}_5+a_{0,6}{x}_6+a_{0,7}{x}_7,\end{array}$$

(25)

$$\begin{array}{c}{S}_{2,1}=a_{1,0}{x}_0+a_{1,1}{x}_1+a_{1,2}{x}_2+a_{1,3}{x}_3+a_{1,4}{x}_4+a_{1,5}{x}_5+a_{1,6}{x}_6+a_{1,7}{x}_7,\end{array}$$

(26)

$$\begin{array}{c}{T}_{2,0}=a_{2,0}{x}_0+a_{2,1}{x}_1+a_{2,2}{x}_2+a_{2,3}{x}_3+a_{2,4}{x}_4+a_{2,5}{x}_5+a_{2,6}{x}_6+a_{2,7}{x}_7,\end{array}$$

(27)

$$\begin{array}{c}{T}_{2,1}=a_{3,0}{x}_0+a_{3,1}{x}_1+a_{3,2}{x}_2+a_{3,3}{x}_3+a_{3,4}{x}_4+a_{3,5}{x}_5+a_{3,6}{x}_6+a_{3,7}{x}_7,\end{array}$$

(28)

$$\begin{array}{c}{T}_{1,0}=a_{4,0}{x}_0+a_{4,1}{x}_1+a_{4,2}{x}_2+a_{4,3}{x}_3+a_{4,4}{x}_4+a_{4,5}{x}_5+a_{4,6}{x}_6+a_{4,7}{x}_7,\end{array}$$

(29)

$$\begin{array}{c}{T}_{1,1}=a_{5,0}{x}_0+a_{5,1}{x}_1+a_{5,2}{x}_2+a_{5,3}{x}_3+a_{5,4}{x}_4+a_{5,5}{x}_5+a_{5,6}{x}_6+a_{5,7}{x}_7,\end{array}$$

(30)

$$\begin{array}{c}{T}_{1,2}=a_{6,0}{x}_0+a_{6,1}{x}_1+a_{6,2}{x}_2+a_{6,3}{x}_3+a_{6,4}{x}_4+a_{6,5}{x}_5+a_{6,6}{x}_6+a_{6,7}{x}_7,\end{array}$$

(31)

$$\begin{array}{c}{T}_{1,3}=a_{7,0}{x}_0+a_{7,1}{x}_1+a_{7,2}{x}_2+a_{7,3}{x}_3+a_{7,4}{x}_4+a_{7,5}{x}_5+a_{7,6}{x}_6+a_{7,7}{x}_7,\end{array}$$

(32)

Therefore, based on the relationships indicated so far, the transformation matrix $\mathbf{A}$ described in Equation (1), which is appropriate for the analyzed WT algorithm, takes the form

$$\mathbf{A}=\left[\begin{array}{cccccccc}\frac{5-\sqrt{3}}{16}& \frac{5+\sqrt{3}}{16}& \frac{3+3\sqrt{3}}{16}& \frac{5+3\sqrt{3}}{16}& \frac{3+\sqrt{3}}{16}& \frac{3-\sqrt{3}}{16}& \frac{5-3\sqrt{3}}{16}& \frac{3-3\sqrt{3}}{16}\\ \frac{3+\sqrt{3}}{16}& \frac{3-\sqrt{3}}{16}& \frac{5-3\sqrt{3}}{16}& \frac{3-3\sqrt{3}}{16}& \frac{5-\sqrt{3}}{16}& \frac{5+\sqrt{3}}{16}& \frac{3+3\sqrt{3}}{16}& \frac{5+3\sqrt{3}}{16}\\ -\frac{1+\sqrt{3}}{16}& \frac{1-\sqrt{3}}{16}& \frac{3-3\sqrt{3}}{16}& -\frac{1+\sqrt{3}}{16}& -\frac{3-5\sqrt{3}}{16}& \frac{3+5\sqrt{3}}{16}& \frac{1-\sqrt{3}}{16}& -\frac{3+3\sqrt{3}}{16}\\ -\frac{3-5\sqrt{3}}{16}& \frac{3+5\sqrt{3}}{16}& \frac{1-\sqrt{3}}{16}& -\frac{3+3\sqrt{3}}{16}& -\frac{1+\sqrt{3}}{16}& \frac{1-\sqrt{3}}{16}& \frac{3-3\sqrt{3}}{16}& -\frac{1+\sqrt{3}}{16}\\ \frac{1-\sqrt{3}}{4\sqrt{2}}& -\frac{3-\sqrt{3}}{4\sqrt{2}}& \frac{3+\sqrt{3}}{4\sqrt{2}}& -\frac{1+\sqrt{3}}{4\sqrt{2}}& 0& 0& 0& 0\\ 0& 0& \frac{1-\sqrt{3}}{4\sqrt{2}}& -\frac{3-\sqrt{3}}{4\sqrt{2}}& \frac{3+\sqrt{3}}{4\sqrt{2}}& -\frac{1+\sqrt{3}}{4\sqrt{2}}& 0& 0\\ 0& 0& 0& 0& \frac{1-\sqrt{3}}{4\sqrt{2}}& -\frac{3-\sqrt{3}}{4\sqrt{2}}& \frac{3+\sqrt{3}}{4\sqrt{2}}& -\frac{1+\sqrt{3}}{4\sqrt{2}}\\ \frac{3+\sqrt{3}}{4\sqrt{2}}& -\frac{1+\sqrt{3}}{4\sqrt{2}}& 0& 0& 0& 0& \frac{1-\sqrt{3}}{4\sqrt{2}}& -\frac{3-\sqrt{3}}{4\sqrt{2}}\end{array}\right].$$

(33)

Knowing the values of the subsequent transformation coefficients $a_{i,j}$ allows us, in accordance with Equation (2), to determine the transfer functions related to the subsequent output quantities of the analyzed algorithm. A certain regularity can be noticed here—the transmittances related to the output quantities with the same scale number will be identical, regardless of the time shift number (this applies separately to details and signal approximation). Hence,

$$\begin{array}{c}{H}_{S_{m,n}}\left(z\right)=H_{S_{m,n^{\prime }}}\left(z\right),\end{array}$$

(34)

$$\begin{array}{c}{H}_{T_{m,n}}\left(z\right)=H_{T_{m,n^{\prime }}}\left(z\right),\end{array}$$

(35)

for $n^{\prime }\in \mathbb{N}$. This property significantly simplifies the analysis of the metrological characteristics of WT algorithms, especially considering that the quantity of output values generated by these algorithms typically exceeds the number of scales analyzed, stemming from the iterations in the input signal decomposition process [3,16].

It is essential to highlight that the analysis presented can also be conducted for other wavelet families, varying numbers of decomposition process iterations, and different quantities of algorithm input. In instances involving wavelet families with explicitly defined mother wavelet equations (instead of utilizing appropriate relations as shown in the example), it becomes feasible to determine the algorithm’s transfer function through suitable transformations [26]. Nonetheless, the connection between the algorithm’s transfer function and the transformation matrix coefficient values is consistently described by Equation (2).

Another scenario that warrants discussion is when adjustments can be made to the transformation matrix coefficients. This situation may involve modifications such as the incorporation of an additional window function $w\left(n\right)$, where the algorithm’s output values are determined based on the following relationship:

$$X\left(i\right)=a_{i,0}w\left(0\right)x\left(0\right)+a_{i,1}w\left(1\right)x\left(1\right)+\cdots +a_{i,N-1}w\left(N-1\right)x\left(N-1\right).$$

(36)

In the discussed case, the coefficients of the transformation matrix of the algorithm are modified according to the equation

$$a_{i,j}^{\prime }=w\left(j\right)a_{i,j},$$

(37)

where $a_{i,j}^{\prime }$ is the new value of the transformation coefficient $a_{i,j}$. In this case, a different transfer function is associated with each output quantity, and Equations (34) and (35) are no longer valid.

3. Own Errors of the Wavelet Transform Algorithm

In analyzing the values of the subsequent coefficients of the transformation matrix given in Equation (33), it can be noticed that most of the coefficients are irrational. Based on the analysis of Equation (1), it can be noticed that in order to determine a single value for the implementation of the selected output quantity of this algorithm, it is necessary to perform N multiplications and N additions. These calculations are most often performed by a microprocessor and are therefore associated with errors related to the limited precision of real numbers [34]. It can therefore be noticed that, in addition to processing error signals present in the algorithm’s input quantities, this algorithm will introduce its own errors resulting from the presented phenomena into the output quantities.

In order to determine the parameters of the own error signal $e_{X,z,i}\left(j\right)$ for the subsequent output values of the algorithm, it is proposed to perform an appropriate experiment using the Monte-Carlo method [35]. During the experiment, random values of input quantities from the range of values that the algorithm processes in real conditions will be provided as input for the actual implementation of the algorithm. Simultaneously, with the result ${\tilde{X}}_i\left(j\right)$ for the real algorithm, the result ${\dot{X}}_i\left(j\right)$ should be determined in the case of the ideal algorithm, implemented according to Equation (1). In assuming that the input quantities $x\left(i\right)$ of the algorithm are not subject to any error signal (i.e., $\tilde{x}\left(i\right)=\dot{x}\left(i\right)$), the error signal $e_{X,z,i}\left(j\right)$ can be defined as follows:

$$e_{X,z,i}\left(j\right)={\tilde{X}}_i\left(j\right)-{\dot{X}}_i\left(j\right).$$

(38)

Based on the values of the subsequent realizations of the error signal $e_{X,z,i}\left(j\right)$, its variances, expected value, and associated expanded uncertainty for a given confidence level can be determined. Unfortunately, it is actually impossible to carry out the presented experiment. The process of determining the value of the realization of the quantity ${\dot{X}}_i\left(j\right)$ in the case of an ideal algorithm must also be carried out using a microprocessor, and this, therefore, also introduces errors related to the phenomena discussed earlier in this quantity. However, it is proposed to estimate the parameters of the error signal $e_{X,z,i}\left(j\right)$ by carrying out the process of determining the value of ${\dot{X}}_i\left(j\right)$ using numbers with much greater precision than in the case of the value ${\tilde{X}}_i\left(j\right)$. Since modern microcontrollers usually use real numbers with word lengths of 16 or 32 bits [36,37], the values of ${\dot{X}}_i\left(j\right)$ may be determined using a word length of 128 bits [34,38].

The discussed experiment was performed for the algorithm described in the previous section, for which its values for the transformation matrix coefficients are presented in Equation (33). In each experiment, 100,000 random values of the $x\left(i\right)$ signal realization from the selected range were fed as input to the algorithm. Then, based on Equation (38), subsequent realizations of the algorithm’s own error signal were determined. The obtained results allowed for the estimation of the variance, expanded uncertainty, and expansion coefficient for the analyzed cases. The experiment was conducted for an implementation of the algorithm using 16- and 32-bit floating-point numbers. The algorithm was implemented in C using the GNU GCC [38] compiler to generate the machine code. This compiler is also utilized for the “ARM” and “AVR” platforms [37,38,39]. The experimental results are summarized in

Table 1. Summary of the simulation-obtained values of the rounding error signal variance of subsequent output quantities of the discrete wavelet transform algorithm for the wavelet “db2” with two iterations of the decomposition process, for numbers with a length of 16 bits, depending on the range of possible values of the implementation of the input quantities.

Quantity	WT Algorithm Input Values Range
	[−1; 1]	[−2; 2]	[−3; 3]	[0; 2]	[0; 4]	[3; 9]
$S_{2,0}$	1.03 × 10−7	4.12 × 10−7	9.74 × 10−7	7.08 × 10−7	2.83 × 10−6	2.16 × 10−5
$S_{2,1}$	1.11 × 10−7	4.43 × 10−7	1.01 × 10−6	9.55 × 10−7	3.82 × 10−6	2.75 × 10−5
$T_{2,0}$	1.31 × 10−7	5.24 × 10−7	1.21 × 10−6	2.87 × 10−7	1.15 × 10−6	7.06 × 10−6
$T_{2,1}$	1.07 × 10−7	4.29 × 10−7	9.73 × 10−7	3.36 × 10−7	1.34 × 10−6	9.30 × 10−6
$T_{1,0}$	7.09 × 10−8	2.85 × 10−7	6.99 × 10−7	1.65 × 10−7	3.56 × 10−7	4.29 × 10−6
$T_{1,1}$	5.83 × 10−8	2.33 × 10−7	5.87 × 10−7	1.49 × 10−7	5.95 × 10−7	4.05 × 10−6
$T_{1,2}$	5.84 × 10−8	2.34 × 10−7	5.85 × 10−7	1.49 × 10−7	5.97 × 10−7	4.05 × 10−6
$T_{1,3}$	5.53 × 10−8	2.21 × 10−7	5.41 × 10−7	1.78 × 10−7	7.08 × 10−7	5.26 × 10−6

and

Table 2. Summary of the simulation-obtained values of the rounding error signal variance of subsequent output quantities of the discrete wavelet transform algorithm for the wavelet “db2” with two iterations of the decomposition process, for numbers with a length of 32 bits, depending on the range of possible values of the implementation of the input quantities.

Quantity	WT Algorithm Input Values Range
	[−1; 1]	[−2; 2]	[−3; 3]	[0; 2]	[0; 4]	[3; 9]
$S_{2,0}$	1.40 × 10−15	5.57 × 10−15	1.30 × 10−14	9.29 × 10−15	3.71 × 10−14	2.78 × 10−13
$S_{2,1}$	1.40 × 10−15	5.61 × 10−15	1.29 × 10−14	1.25 × 10−14	5.00 × 10−14	3.54 × 10−13
$T_{2,0}$	1.71 × 10−15	6.83 × 10−15	1.58 × 10−14	3.33 × 10−15	1.33 × 10−14	7.78 × 10−14
$T_{2,1}$	1.39 × 10−15	5.56 × 10−15	1.29 × 10−14	4.27 × 10−15	1.71 × 10−14	1.16 × 10−13
$T_{1,0}$	8.28 × 10−16	3.54 × 10−15	8.38 × 10−15	1.70 × 10−15	6.84 × 10−15	4.25 × 10−14
$T_{1,1}$	6.82 × 10−16	2.72 × 10−15	6.67 × 10−15	1.46 × 10−15	5.84 × 10−15	4.02 × 10−14
$T_{1,2}$	6.82 × 10−16	2.72 × 10−15	6.66 × 10−15	1.47 × 10−15	5.86 × 10−15	4.02 × 10−14
$T_{1,3}$	6.68 × 10−16	2.67 × 10−15	6.48 × 10−15	2.02 × 10−15	8.10 × 10−15	6.15 × 10−14

. Based on the results obtained, the expansion coefficient for the distributions of the analyzed signals was estimated to be $c_z=2.15$ on average.

In analyzing the results presented in Table 1 and Table 2, several of the most important relationships can be noticed. The value of the variance of the self-error signal has the following characteristics:

• It depends on the number of arithmetic operations performed on non-zero coefficients of the transformation matrix and increases with the number of these operations;
• It depends on the range of possible implementation values of the algorithm’s input quantities and increases as this range is extended;
• It depends on the length of the word used by the algorithm and decreases as this length increases.

Therefore, from the perspective of the measuring chain design, this value will rely on the number of the input quantities of the WT algorithm, the iterations of the signal decomposition process, and the order and type of the mother wavelet used. It is also noticeable that for the same stage of signal decomposition, the error signal parameters are nearly identical for each output quantity. Hence, for each algorithm implementation, the parameters of the self-error signals should be determined for the subsequent decomposition stages. The experiment conditions should closely resemble the actual operational conditions of the analyzed algorithm.

4. Method for Determining the Resultant Expanded Uncertainty

In many cases, merely knowing the value of the standard uncertainty parameter or the error signal variance may not suffice [40]. The most universal method used to determine the resulting expanded uncertainty value in a general situation is the Monte-Carlo method. However, this approach necessitates numerous iterations to obtain the final value of the analyzed error signal and may not be suitable for the real-time evaluation of the metrological properties of the measurement chain due to the time required to repeat the experiment.

An alternative to the Monte-Carlo method could be an analytical method, such as the propagation of distribution functions method described in sources like [41] or the extended rule of the combination of uncertainties method outlined in [42]. However, these methods are intricate, and their application to the error model proposed in this study might be less effective compared to the reduction interval arithmetic method, as discussed in works like [43,44,45]. Another alternative could be a method based on fuzzy logic, as detailed in [46].

Due to the advantages of the interval arithmetic reduction method, such as its simplicity of application, its ability to determine new values of the resulting expanded uncertainty when the error model parameters change without the need for a Monte-Carlo experiment, and its low computational complexity, the reduction interval arithmetic method was employed to ascertain the resulting value of the expanded uncertainty in this study.

According to the method of reduction interval arithmetic, the expanded uncertainty in case of an error signal $e_{\Sigma }\left(t\right)=e_0+\cdots +e_{N-1}$ can be determined according to the following relation [43,44]:

$$U_{\Sigma }=\sqrt{{\left[\begin{array}{c}{U}_0\\ U_1\\ ⋮\\ U_{N-1}\end{array}\right]}^T\left[\begin{array}{cccc}1& h_{0,1}& \cdots & h_{0,N-1}\\ h_{1,0}& 1& & h_{1,N-1}\\ ⋮& & \ddots & ⋮\\ h_{N-1,0}& \cdots & \cdots & 1\end{array}\right]\left[\begin{array}{c}{U}_0\\ U_1\\ ⋮\\ U_{N-1}\end{array}\right]},$$

(39)

in which the successive values of the expanded uncertainties $U_i$ for the analyzed component $e_i\left(t\right)$ of the error signal $e_{\Sigma }\left(t\right)$ are determined according to the relation [40]

$$U_i=c_i{\sigma }_i,$$

where $c_i$ is the coverage factor for the i-th $e_{\Sigma }\left(t\right)$ signal component. The successive coherence coefficients $h_{i,j}$ are determined according to the following equation [43]:

$$h_{i,j}=h_{j,i}=s_{i,j}\left(\frac{U_i^2+U_j^2}{{\sum }_{k=0}^{N-1}{U}_k^2}\right),$$

(40)

where symbol $s_{i,j}$ denotes the shape factor, which is determined for a pair of signals $e_a\left(t\right)$ and $e_b\left(t\right)$ has a specific distribution shape of the following form [43]:

$$s_{a,b}=s_{b,a}=\frac{U_{a,b}^2-U_a^2-U_b^2}{2U_a{U}_b}=\frac{U_{a,b}^2}{2U^2}-1,$$

assuming that signals $e_a\left(t\right)$ and $e_b\left(t\right)$ are uncorrelated and have the same value of expanded uncertainty $U_a=U_b$ for the same confidence level $1-\alpha$ [40].

The coherence coefficients $h_{i,j}$ determined consider the relationships between the shapes of the provided error distributions and their correlations, incorporating corrections arising from the central limit theorem [40]. The shape factor coefficients $s_{i,j}$ need to be computed only once, which can be achieved through a simulation experiment (Monte-Carlo method) or analytically. Various works such as [43,44] describe the process of calculating these coefficients in different manners.

Table 3. Summary of shape factor values for pairs of signals with typical probability density functions for the confidence level 95%, where the following symbols denote a distribution: $\left(n\right)$ normal, $\left(u\right)$ uniform, $\left(t\right)$ triangular, $\left(d\right)$ u-shape (sine function distribution), $\left(z\right)$ WT self-error signal distribution.

${\mathit{s}}_{\mathit{a},\mathit{b}}$	n	u	t	d	z
$\mathit{n}$	$0.0000$	$0.1561$	$0.0250$	$0.2988$	$-0.0091$
$\mathit{u}$	$0.1561$	$0.3356$	$0.1773$	$0.5337$	$0.0662$
$\mathit{t}$	$0.0250$	$0.1773$	$0.0419$	$0.3504$	$-0.0104$
$\mathit{d}$	$0.2988$	$0.5337$	$0.3504$	$0.7136$	$0.1971$
$\mathit{z}$	$-0.0091$	$0.0662$	$-0.0104$	$0.1971$	$0.0273$

displays values for typical scenarios at a 95% confidence level for $1-\alpha$.

5. Application of the Proposed Analysis Method

The effectiveness of the proposed analysis method was verified using the measurement chain detailed in a previous study [19]. This measurement chain processes a time-varying voltage signal $s\left(t\right)$, with realization values falling within the range $\widehat{s}\left(t\right)\in [0; 1]\mathrm{V}$ on discrete representations of $x\left(i\right)$ of this magnitude, where $x\left(i\right)=s\left(i{T}_s\right)$, and $f_s=\frac{1}{T_s}=48 \mathrm{kHZ}$ denotes the sampling period. In a single measurement series, with $N=8$, consecutive samples of quantities $x\left(i\right)$ are inputted into the WT algorithm described in the previous section, leading to the determination of a vector $M=8$ of output quantities $X\left(j\right)$. The algorithm’s implementation utilizes floating-point numbers with a word length of 32 bits.

The metrological properties of the section of the measurement chain responsible for converting the quantity $s\left(t\right)$ into $x\left(i\right)$ were previously discussed in [19]. The uncertainty budget related to $x\left(i\right)$ encompasses the random error signal $e_{x,r}\left(i\right)$ and the dynamic error signal $e_{x,d}\left(i\right)$. For the random error signal $e_{x,r}\left(i\right)$, the variance is ${\sigma }_{x,r}^2=0.13 \mu \mathrm{V}$, and the realization values of this signal are distributed close to normal. Hence, for a confidence level of $1-\alpha =95\%$, the coverage factor $c_{x,r}=c_n=1.96$, and the expanded uncertainty is $U_{x,r}=0.71 \mathrm{mV}$.

Regarding the dynamic error signal $e_{x,d}\left(i\right)$, its characteristics are influenced by the spectrum of the processed signal $s\left(t\right)$. If this signal is a sinusoidally varying signal with pulsation ${\omega }_{s,o}$ and amplitude $E_{s,o}$, the following conditions apply [19]:

$${\sigma }_{x,d,sin}^2=\frac{{\left(E_{s,o}cos\left({\tilde{\phi }}_y\left({\omega }_{s,o}\right)\right)-E_{s,o}\right)}^2+{\left(E_{s,o}sin\left({\tilde{\phi }}_y\left({\omega }_{s,o}\right)\right)\right)}^2}{2},$$

however, if the signal is a triangular signal with the given parameters, then [19]:

$${\sigma }_{x,d,tri}^2=\sum _{i=1}^{\infty }{\sigma }_{x,d,tri,i}^2,$$

$${\sigma }_{x,d,tri,i}^2=\frac{{\left(E_{s,o,i}cos\left({\tilde{\phi }}_y\left({\omega }_{s,o,i}\right)\right)-E_{s,o,i}\right)}^2+{\left(E_{s,o,i}sin\left({\tilde{\phi }}_y\left({\omega }_{s,o,i}\right)\right)\right)}^2}{2},$$

(41)

$$E_{s,o,i}=\frac{\pi }{8}{\left(2i-1\right)}^{-2}{E}_{s,o},$$

where in Equation (41), only those harmonics of the signal $s\left(t\right)$ for which $k{f}_{s,o}\le \frac{1}{2}{f}_s$ holds are taken into account, where $k=(2i-1)$. The introduced phase shift ${\tilde{\phi }}_y$, which is the origin of the dynamic error signal, is estimated according to the relationship [19]

$${\tilde{\phi }}_y\left(\omega \right)\approx -6.26\times {10}^{-13}{\omega }^2-5.73\times {10}^{-7}\omega .$$

As the ambient conditions did not change during the experiments, it is assumed that the static error signal $e_{x,s}\left(i\right)$ does not occur.

In considering the aforementioned relationships, in the scenario of a sinusoidally varying signal $s\left(t\right)$, the resulting error signal $e_{x,\Sigma ,sin}\left(i\right)$ will comprise the signal $e_{x,r}\left(i\right)$ and the signal $e_{x,d,sin}\left(i\right)$ consisting of a single harmonic with pulsation ${\omega }_{s,o}$. In the instance of a triangular signal, the error signal $e_{x,\Sigma ,tri}\left(i\right)$ will also encompass a random error component $e_{x,r}\left(i\right)$, and depending on the signal’s pulsation $s\left(t\right)$, it will include a specific number of harmonics of the random error signal $e_{x,d,tri}\left(i\right)$, where the i-th harmonic is characterized by a pulsation of ${\omega }_{x,e,i}=(2i-1){\omega }_{s,o}$.

These error signals discussed will be propagated to the algorithm output as per Equation (1). The parameters at the algorithm output can be determined for subsequent output quantities following Equation (6) for the random error signals and Equation (5) for the successive harmonics of the resultant dynamic error signal. Moreover, in the scenarios discussed, the WT algorithm will introduce its own error signals $e_{X,z,i}\left(j\right)$ to the output values, linked to the previously mentioned roundings.

The vector of expanded uncertainties related to the resultant error signal for subsequent output quantities of the algorithm, as required in Equation (39), can thus be delineated as

$${\mathbf{U}}_{\ast ,sin}=\left[\begin{array}{ccc}{U}_{\ast ,z}& U_{\ast ,r}& U_{\ast ,d,sin}\end{array}\right],$$

in the case of a sinusoidal signal, and in the case of a triangular signal, in the following form:

$${\mathbf{U}}_{\ast ,tri}=\left[\begin{array}{cccccc}{U}_{\ast ,z}& U_{\ast ,r}& U_{\ast ,d,tri,1}& U_{\ast ,d,tri,2}& \dots & U_{\ast ,d,tri,N}\end{array}\right].$$

The symbol “∗” represents the output quantity number, assigned following the symbols introduced earlier in Equation (8). Based on the previously addressed characteristics of the WT algorithm, for random error signals, the uncertainty linked to the propagation of these signals to the WT algorithm’s output can be calculated using the following relationship:

$$U_{\ast ,r}=c_n{\sigma }_{x,r}\sqrt{\sum _{j=0}^{N-1}{a}_{\ast ,j}^2},$$

resulting from Equation (6), while in the case of the subsequent harmonics of the resultant dynamic error signal,

$$U_{\ast ,d}\left(\omega \right)=c_d{\sigma }_{x,d}\left(\omega \right)\left|H_{\ast }\left(e^{\mathrm{j}\omega T_s}\right)\right|=c_d\frac{1}{\sqrt{2}}{E}_{x,e}\left(\omega \right)\left|H_{\ast }\left(e^{\mathrm{j}\omega T_s}\right)\right|,$$

which results from Equation (5). The resulting expanded uncertainty $U_{\ast ,\Sigma }$ related to the subsequent output quantities can be determined according to Equation (39), with the values of the coherence coefficients estimated based on relationship (40). For the confidence level $1-\alpha =95\%$ $c_n=1.96$ and $c_d=1.41$ [40], $c_z=2.15$.

To validate the indicated relationships, a Monte-Carlo measurement experiment was conducted, gathering 30,000 values of the $X_{\ast }\left(j\right)$ signal realization each time. Throughout the experiment, the signal source $s\left(t\right)$ originated from the RIGOL DG1011 arbitrary waveform generator [47]. The initial phase of this signal was randomized within the interval $[-\pi ; \pi ]$, which was determined using the generator’s synchronizing output. The signal’s frequency $s\left(t\right)$ ranged from ${\widehat{f}}_{s,o}\in [1; 20]\mathrm{kHz}$ for a monoharmonic signal and from ${\widehat{f}}_{s,o}\in [1; 5]\mathrm{kHz}$ for a polyharmonic signal. The signal parameters remained constant at $D_{s,o}=0.5\mathrm{V}$ and $E_{s,o}=0.475\mathrm{V}$.

Based on the collected values of the $X_{\ast }\left(j\right)$ quantity realization, we have the following the equation:

$$e_{\ast ,\Sigma }\left(j\right)={\tilde{X}}_{\ast }\left(j\right)-{\dot{X}}_{\ast }\left(j\right),$$

The error signal values $e_{\ast ,\Sigma }\left(j\right)$ were determined, and their variances were calculated along with the expanded uncertainty. The measured expanded uncertainty value $U_m$ was compared with the value $U_c$ determined using Equation (39), and the relative error in estimating this value was computed. The outcomes for specific values of signal $s\left(t\right)$ pulsation are outlined in

Table 4. Summary of the values of the relative error of estimating the expanded uncertainty values obtained by means of a measurement experiment (the case of a monoharmonic signal).

${\mathit{f}}_{\mathit{s},\mathit{o}}$, Hz	Relative Error Value $\mathit{\delta }$ for Expanded Uncertainty Estimation, %
	$S_{2,0}$	$S_{2,1}$	$T_{2,0}$	$T_{2,1}$	$T_{1,0}$	$T_{1,1}$	$T_{1,2}$	$T_{1,3}$
1000	$-3.52$	$-4.27$	$+9.29$	$+7.09$	$+30.49$	$+24.32$	$+24.32$	$+9.49$
2000	$+4.66$	$+0.49$	$+6.57$	$+4.77$	$+15.34$	$+27.58$	$+27.58$	$+5.04$
3000	$+10.24$	$+5.10$	$+7.81$	$+6.53$	$+11.67$	$+21.96$	$+21.96$	$+6.19$
4000	$-0.38$	$+9.21$	$+10.38$	$+9.09$	$+13.80$	$+16.49$	$+16.49$	$+10.12$
5000	$+3.99$	$+2.69$	$+3.51$	$+2.45$	$+5.71$	$+5.62$	$+5.62$	$+4.02$
6000	$+6.21$	$+5.83$	$+5.94$	$+5.84$	$+8.78$	$+9.08$	$+9.08$	$+7.98$
7000	$+9.29$	$+10.38$	$+8.84$	$+8.15$	$+9.86$	$+10.86$	$+10.86$	$+9.87$
8000	$+0.83$	$+1.51$	$+0.65$	$+1.82$	$+1.77$	$+1.63$	$+1.63$	$-0.17$
9000	$+2.99$	$+2.06$	$+1.74$	$+4.68$	$+1.82$	$+1.30$	$+1.30$	$+1.67$
10,000	$+6.33$	$+6.23$	$+5.04$	$+5.55$	$+5.46$	$+4.36$	$+4.36$	$+5.05$
11,000	$-1.64$	$-0.06$	$+0.49$	$+0.75$	$-0.41$	$-0.61$	$-0.61$	$+0.05$
12,000	$-11.33$	$-11.98$	$-1.06$	$-1.05$	$-1.68$	$-1.68$	$-1.68$	$-1.72$
13,000	$+5.20$	$+3.38$	$+4.78$	$+5.54$	$+5.23$	$+4.46$	$+4.46$	$+5.46$
14,000	$+2.83$	$+2.93$	$+2.50$	$+2.20$	$+3.40$	$+2.39$	$+2.39$	$+3.60$
15,000	$+2.50$	$+5.49$	$+2.42$	$+2.89$	$+3.04$	$+2.10$	$+2.17$	$+2.53$
16,000	$-0.34$	$-1.63$	$-2.09$	$-2.94$	$-2.44$	$-2.76$	$-2.76$	$-2.34$
17,000	$+3.20$	$+5.51$	$+4.40$	$+0.97$	$+3.18$	$+3.05$	$+3.05$	$+3.91$
18,000	$+6.19$	$+7.02$	$+1.95$	$+2.66$	$+4.57$	$+4.56$	$+4.56$	$+4.50$
19,000	$-4.08$	$-3.37$	$-13.04$	$-3.45$	$-3.70$	$-4.00$	$-4.00$	$-3.70$
20,000	$+5.87$	$+2.91$	$-15.64$	$+6.81$	$+8.00$	$+7.68$	$+7.68$	$+7.34$
Mean 1	$4.58$	$4.60$	$5.41$	$4.26$	$7.02$	$7.82$	$7.83$	$4.74$

¹ Mean of absolute $\delta$ values for selected output quantity.

and

Table 5. Summary of the values of the relative error of estimating the expanded uncertainty values obtained by means of a measurement experiment (the case of a polyharmonic signal).

${\mathit{f}}_{\mathit{s},\mathit{o}}$, Hz	Relative Error Value $\mathit{\delta }$ for Expanded Uncertainty Estimation, %
	$S_{2,0}$	$S_{2,1}$	$T_{2,0}$	$T_{2,1}$	$T_{1,0}$	$T_{1,1}$	$T_{1,2}$	$T_{1,3}$
1000	$+20.52$	$+18.96$	$-4.26$	$+4.16$	$+21.57$	$+10.47$	$+10.47$	$+9.63$
2000	$+8.99$	$+5.90$	$-5.19$	$-0.01$	$-22.21$	$-22.25$	$-22.25$	$-7.15$
3000	$-12.69$	$+16.92$	$+8.45$	$+9.16$	$-15.20$	$-13.56$	$-13.56$	$+7.15$
4000	$-19.09$	$-2.30$	$+9.05$	$+9.93$	$-8.68$	$-12.27$	$-12.27$	$+15.61$
5000	$-15.90$	$-0.52$	$+9.88$	$+0.00$	$-12.82$	$-20.08$	$-20.08$	$+22.54$
Mean 1	$15.44$	$8.92$	$7.37$	$4.65$	$16.10$	$15.73$	$15.73$	$12.42$

¹ Mean of absolute $\delta$ values for selected output quantity.

.

For example, in the case of a monoharmonic signal with a frequency $f_{s,o}=\frac{{\omega }_{s,o}}{2\pi }=5 \mathrm{kHz}$ for the output quantity $T_{2,1}$, the following occurs:

$$\begin{gathered} U_{T_{2,1},z}=c_z{\sigma }_{T_{2,1}}=2.15·3.73\times {10}^{-8}=8.02\times {10}^{-5} \mathrm{mV}, \\ {U}_{T_{2,1},r}=c_n{\sigma }_{x,r}\sqrt{\sum _{j=0}^{N-1}{a}_{T_{2,1},j}^2}=1.96·0.36\times {10}^{-4}·1.0=0.70 \mathrm{mV}, \\ {U}_{T_{2,1},d,sin}=c_d{\sigma }_{x,d,sin}\left({\omega }_{s,o}\right)\left|H_{T_{2,1}}\left(e^{\mathrm{j}{\omega }_{s,o}{T}_s}\right)\right|=1.41·6.31\times {10}^{-4}·1.59=14.09 \mathrm{mV}, \\ {\mathbf{U}}_{T_{2,1},sin}=\left[\begin{array}{ccc}{U}_{T_{2,1},z}& U_{T_{2,1},r}& U_{T_{2,1},d,sin}\end{array}\right]=\left[\begin{array}{ccc}8.02\times {10}^{-8}& 0.70\times {10}^{-3}& 14.09\times {10}^{-3}\end{array}\right]\mathrm{V}, \\ {U}_{T_{2,1},\Sigma }=\sqrt{{\left[\begin{array}{c}8.02\times {10}^{-8}\\ 0.70\times {10}^{-3}\\ 14.09\times {10}^{-3}\end{array}\right]}^T\left[\begin{array}{ccc}1.000& 0.000& 0.197\\ 0.000& 1.000& 0.534\\ 0.197& 0.534& 1.000\end{array}\right]\left[\begin{array}{c}8.02\times {10}^{-8}\\ 0.70\times {10}^{-3}\\ 14.09\times {10}^{-3}\end{array}\right]}=14.48 \mathrm{mV}, \end{gathered}$$

and for the indicated case, 14.13 mV was measured. Expanded uncertainty values can be determined similarly in other cases, but due to the large number of components of the indicated equations, further examples were not published in this paper.

6. Conclusions

The expanded uncertainty values determined for a monoharmonic signal for subsequent output values of the measurement chain consistently aligned with those obtained experimentally. The typical deviation between the actual and estimated expanded uncertainty fell within the range of ±5%. However, in the case of a polyharmonic signal, the resultant expanded uncertainty values were determined with a larger margin of error. This discrepancy is primarily due to the imprecise determination of the error model parameter of the analyzed measurement chain. These parameters were established in a previous study [19], incorporating several simplifications to ensure that the presented analysis example maintained a minimal level of complexity.

Upon analyzing the instances of a polyharmonic signal where there is a notable disparity between the estimated and measured expanded uncertainty values, a specific pattern emerges. Each output quantity of the algorithm corresponds to a particular transmittance, which effectively dampens or enhances error signals with a specific spectrum. Consequently, if there is an inaccurate estimation of the parameter of a particular error signal, especially if this signal emerges as the predominant signal in the algorithm’s output, the estimation will be erroneous.

This observation highlights a crucial insight—even with highly accurate estimations of the resultant error signal parameters of the input quantities of the WT algorithm, there remains a possibility that the parameters of the error signal at the algorithm’s output are inaccurately estimated. The dominant error signal (with accurately determined parameters) could be suppressed, while the less significant signal (with imprecisely determined parameters) could undergo significant amplification.

It can be seen that the transfer function of the WT algorithm can be determined in two ways. The first method, described in [16,17,18], is the most accessible to the designer of the measurement track and does not require any knowledge about the algorithm used; however, its use requires a ready implementation of this algorithm. The second method, presented in this paper, requires knowledge of the assumptions of the wavelet family used and their transformation to the form described in Equation (2), by following the steps below:

• Indicate the form of the vector of the output quantities based on the number of input quantities, the type of wavelet, and the number of iterations of the signal decomposition process (as shown in Equations (8) and (9) in the case of the analyzed example).
• Determine the values of the scaling factors based on the assumptions of the selected family and the order of the selected mother wavelet (as shown in Equations (13)–(16) in the case of the analyzed example).
• Determine the equations describing the output quantities indicated in step 1, based on Equations (10) and (11) and substitute the values determined in step 2 into them (as shown in Equations (17)–(32) in the case of the analyzed example).

This procedure is much more complex due to the necessary calculations, but it does not require the implementation of the algorithm used.

It should be noted that the use of the interval arithmetic reduction method to estimate the resulting expanded uncertainty value comes down to the following:

• Determining the value of the vector of partial uncertainties;
• Determining the values of the coherence coefficients;
• The application of Equation (39) for the obtained data.

This implies that in the event of a change in the error model parameter (e.g., alterations in the spectrum of the processed signal, the emergence of an additional error signal, a modification in the parameters of existing error signals), there is no need to conduct Monte-Carlo simulations or any other intricate procedures to ascertain the current resulting value of the expanded uncertainty. The sole task that requires additional effort is determining the values of shape coefficients, which occurs once during the analysis preparation stage. However, this feature indicates that the obtained results may be inaccurate if the analyzed error signals exhibit a realization distribution shape different from the one assumed in the calculations.

In conclusion, the proposed analysis method and the previously suggested error model can be deemed suitable. It is essential to highlight that the accuracy of the proposed analysis method hinges on the precision of determining the parameters of the proposed error model. The authors of this article anticipate that the proposed method will be beneficial for the designers of measurement chains utilizing WT algorithms or other types of linear measurement data processing algorithms.