1. Introduction
Online monitoring of the technical condition of technical facilities in operation in maritime transport is associated with digital processing of complex signals. The frequency spectrum of the analyzed signals, as a rule, contains the harmonic components of high frequencies, which are estimated on the Fourier basis. Evaluations are performed using hardware, including harmonic analyzers of various designs. Signals of complex configuration include voltages and currents of systems and elements of power complexes containing nonlinear semiconductor converters, sources of noise and vibrations during the operation of the main and auxiliary engines of a ship’s power systems, nonstationary vibrations of the ship’s hull structures in waves, shallow waterways, etc. They also include resonance vibrations that can occur on ships at certain speed modes of powerful energy converters. Signals of a complex configuration, along with the deterministic component, may contain interference, individual pulses with a weak energy spectrum. A wide variety of monitored signals creates conditions for improving the quality of diagnostics of the ship’s operating parameters, using digital technologies under the conditions of simultaneous use of machines and mechanisms which must be supplemented with tools for rapid signal analysis, since they may contain components containing information with signs of violations of the normal operating mode. These features include increased fuel consumption, changes in the temperature regimes of ship systems and mechanisms, the occurrence of vibration modes due to misalignment of the shafts as well as insufficient rigidity of the mountings of the units, and many more. In the modes of monitoring the technical condition, it is required to make decisions promptly, with retrieving the most complete information about the state of the object and assessing the type and location as well as the cause of a defect or malfunction in an automatic mode [
1,
2,
3].
The task, here, is to create such algorithms and systems of technical diagnostics, which should, by individual signs and their combination, determine the likelihood of a malfunction in advance [
1]. This will allow taking preventive measures to ensure troublefree operation of a ship’s technical systems [
1,
2,
3].
Maritime diesel installations (MDIs) are complex complexes consisting of many subsystems and units. During operation, complex chemical–physical processes and energy transformations take place in them. On the other hand, marine diesel installations during operation can be exposed to vibrations, shock loads, moisture, salt, changes in ambient temperature in a wide range, etc. Such influences accelerate the process of degradation of parts and assemblies of the MDIs and reduce their reliability. Maintaining the reliability of MDI operation during autonomous navigation of ships not only ensures the safety of the ship’s crew and cargo but also increases the economy and efficiency of the operation of ships. These goals can be achieved by monitoring the performance of installations and early detection of symptoms of a violation of their normal operation. This is the task of technical diagnostics of the MDI during operation. In general, technical diagnostics can include one or more of the following tasks:
 
Determination of the technical condition (performance) of the object;
 
Searching for a defect that has arisen;
 
Forecasting changes in the technical state of the object.
The theory of technical diagnostics systems in general, and for MDIs in particular, has arisen and developed since the 1970s in the works of worldrenowned researchers. However, insufficient attention is paid to the development of technical means for diagnosing MDIs, which can assist a ship’s crew in assessing the technical condition of the power plant during autonomous navigation of ships [
4,
5,
6].
Increasing requirements for the accuracy of diagnosing objects when processing signals with local impulses and other features leads to the necessity of improving the theory, methods, and techniques for assessing the technical state of a ship’s systems. Traditional methods based on the use of orthogonal series become ineffective for a profound and comprehensive analysis of signals with complex configurations, including cases when they are expanding into an orthogonal series using Fourier approximation. In practical calculations, an infinite number of terms in the Fourier series is not used, and the use of a limited number of harmonics leads to large errors [
7]. In such cases, wavelet transform technologies can be considered very promising solutions [
6]. A large package of wavelets, presented in modern computing systems, as well as a set of functions designed to process signals coming from information sensors, makes it possible to use them effectively to solve practical problems of monitoring and diagnosing the technical state of objects in steadystate and transient modes.
Wavelet technologies are developed on basic functions that preserve localization properties in the frequency and time domains [
4]. They provide the implementation of highprecision algorithms for the approximation of onedimensional and multidimensional characteristics.
2. Methods and Materials
Highprecision approximation procedures in quasistatic modes during the working cycle are required to be used when measuring the pressure, temperature, weight, and composition of exhaust gases in a diesel cylinder. Highprecision technologies are required when studying the appearance of local zones while nitrogen oxides are being generated at various speed modes, researching the correlation properties of fuel combustion products by indirect methods, for producing data arrays obtained using indicator diagrams, and creating “reference” characteristics of technological processes. In this case, information can be stored in the form of arrays of wavelet transform coefficients, and cluster models can constitute a tool for invariants. Parameters beyond the boundaries set by the cluster can serve as a sign of deviation of the diagnosable parameter from the standard one.
Studies have shown that the frequency properties of indicator devices can significantly affect the frequency characteristics of the indicator channel, which leads to distortion of signals and increasing the error in assessing the indicators of the working process [
8].
The digitalization of processes in the Fourier basis has shown that with Fourier series, it is difficult to detect local signals in the form of short pulses of a small amplitude, which may indicate the development of deviations of controlled signals from the standard ones against the background of measurements [
9,
10].
Wavelet analysis is a tool for studying signals with local features. Wavelets and wavelet transforms allow for obtaining spectrograms of signals with local features. It is easy to determine their location in the amplitudetime domain, which is practically impossible to perform in the Fourier basis [
5]. Wavelet transforms are based on modern digital technologies, the implementation of which requires the use of tools of modern computing systems containing a set of functions with established syntactic rules.
The MATLAB Wavelet Toolbox contains a list of basic wavelets, such as Haar, Morlet, Meyer, Shannon, Daubechies, frequency Bspline, etc. The Haar and Daubechies wavelets are most widely used in practice, having a lot in common in information coding [
2]. In this regard, the Haar wavelet in the packet is denoted by the symbol db1, and the Daubechies wavelets can be selected from the db2, db3, …, db10 family, represented by nine wavelets. Modification of Daubechies wavelets allows for an increase in the symmetry properties of filters while maintaining the simplicity of working with them [
11].
Wavelet transforms are designed to process large amounts of information. At the same time, their technologies and algorithms are suitable for working with simple signals that allow one to observe, step by step, the progress of transformations carried out by the algorithm. Let us consider the following example of this mode.
Suppose that the measurements of the hourly fuel consumption during the progress of a vessel with a displacement of 2000 tons during 8 running hours were (kg/h) as follows:
By delta coding, we will obtain an approximation of the fuel consumption vector y using the “db1” wavelet. For this purpose, according to the Haar transformation, we split the elements of the vector
y into pairs and determine the halfsum and halfdifference of values in each of them. We define the halfsums and halfdifferences using the matrix R of the following form:
Vector f = R × y’ = (153.0 1.0 156.5 0.5 156.5 0.5 157.00 −1.0)
^{T} contains mean values of groups (halfsum) and halfdifference. Next, from the mean values, we form the row vector Fcp = (f (1) f (1) f (3) f (3) f (5) f (5) f (7) f (7)) and use it to obtain deviations for the group elements from mean values:
The vectors Fcp and r can also be obtained using the Haar wavelet transform for onestep signal decomposition using the “db1” wavelet, for which we use the following operators contained in the Wavelet Toolbox:
[cA1,cD1] = dwt(y,‘db1’).
As a result, we will have
A1 = upcoef(‘a’,cA1,‘db1’,1,ls); D1 = upcoef(‘d’,cD1,‘db1’,1,ls).
The upcoef operator is used to approximate A1 and refine D1 of the first wavelet transform level by the obtained vectors cA1 and cD1. The identity of the estimates is determined by the fact that A1 = Fcp and D1 = r. According to the syntax of the operator, the following ones are sequentially given in parentheses: the symbols ‘a’ and ‘d’, specifying which kind of decomposition should be performed; vector cA1 is introduced for approximation, and vector cD1 for refining; Daubechies wavelet with the corresponding number ‘db1’ is selected; the number of the operation level is entered—for example, 1; the number of elements of the vector of the processed signal ls is indicated (for the vector y, the value ls = 8).
Now, we will turn our attention to the operator [cA1, cD1] = dwt (y, ‘db1’). As applied to the example under consideration, the operator performs a linear transformation of each group of two elements of the vector yc using the rotation matrix M consisting of the following elements:
As a result, the elements of vectors cA1 and cD1 are determined. In particular, we have, element by element, for four groups:
The Haar transform is normalized using the matrix M, since the determinant of the matrix is det (M) = 1.0000. Therefore, when a pair of numbers rotates, the area of the normalized object does not change, i.e., information is not lost, which makes it possible to process the vector y by the formed pairs of numbers independently of each other. In the MATLAB environment, such a transformation is carried out using a block matrix consisting of blocks M in the main diagonal. Such a matrix can be obtained if the matrix R is multiplied, element by element, by sqrt (2):
Note that det(R) = 0.0625 with the determinant of every block 8 × det(R) = 0.5, while
The rows of the matrix M have orthogonality properties.
As a result, inv(W) = W^{T}, and the inverse of the matrix can be obtained using the transpose operation.
Thus, vectors A1 and D1, obtained by the Haar transform, are two filters that separate the signal into lowfrequency and highfrequency components (
Figure 1). Lowfrequency resolution characterizes signal power and highfrequency resolution characterizes its frequency associated with short pulses (bursts).
A person at the modern level of development of technology and its automation, due to the limited psychophysiological capabilities of a human individual in difficult conditions, cannot always find a managerial action adequate for a problematic or other operational situation. This is especially pronounced on modern ships with a high level of automation, when the navigator must almost instantly make decisions aimed at meeting the many requirements put forward by the operational situation (safety of navigation, energy savings, environmental safety, etc.). Another example is the impossibility of a representative of the organization for the classification and survey of ships for a relatively short survey period to completely check the technical condition of a ship and its elements, or at least all elements of the power plant (EPP) of the ship (in reality, the method of selective control, as a result of which the technical condition of a significant part of the EPP facilities turns out to be untested) (
Figure 1). It is here that digital technologies based on knowledge in various fields of technology and the experience of specialists—in other words, artificial intelligence—should replace the intelligence of a human operator, which cannot always cover and evaluate a multilevel problem in its entire multidimensional difficulties.
Now, let us restore y. We perform signal regeneration using the wavelet transform inversion:
3. Mathematical and Software Implementation
With multilevel processing of signals of a complex configuration based on the ‘db1’ wavelet, it is possible to guarantee an increase in the conversion quality due to a larger volume of digitalization of the process [
8]. For definiteness, we perform the decomposition of y at detail level three (
Figure 2). Then,
[C,L] = wavedec(y,3,‘db1’)
C = [440.5275 −2.8284 −3.5000 −0.5000 −1.4142 −0.7071 −0.7071 1.4142],
L = (1 1 2 4 8).
Using the obtained vectors C and L, using the appcoef function, we estimate the coefficient cA3 of the third level of approximation and detail. Observing the syntax, we obtain:
and for refining at levels 3, 2, and 1, we will have:
cD3 = detcoef(C,L,3) = −2.8284,
cD2 = detcoef(C,L,2) = [ −3.5000 −0.5000], cD1 = detcoef(C,L,1) = [−1.4142 −0.7071 −0.7071 1.4142].
Daubechies wavelets of the ‘db2’, …, ‘db10’ families in the Wavelet Toolbox are formed according to a different principle of forming groups of number sequences, which significantly improves the signal conversion process. Instead of groups with two points, for example, four points are taken with the shift of the elements of the subsequent group at each step by two elements to the left. For the vector y, such groups are as follows:
Now, when realizing the lowpass and highpass filters, each group of four numbers will be replaced by two numbers. Let the values of the elements of the vector y in the group be equal to x, v, z, and t. Then, we write the first filter as follows:
The second filter d is represented by a linear model
Such a permutation of elements provides the orthogonality property of the block Daubechies transformation matrix, which has the form:
To calculate the numerical values of the elements of the DB matrix, four equations are required. The first equation can be determined by the requirement of orthogonality
It must be true for all linear pairs. The other three equations are determined by the requirements of the Daubechies transform.
The sum of the squares of the coefficients must be equal to one
since for a single determinant, DB will also be equal to one.
Transform
for any equal k—for example, k = 1, it must take zero values.
For linearly increasing values of m, for example, m = (1 2 3 4),
Transform
must also take zero values.
The solution of the system of nonlinear Equations (2)–(5) is possible only by numerical methods. For this purpose, we will use the fsolve function in the Optimization Toolbox of the MATLAB package. Here is a fragment of the solution, with the account of the function syntax.
For this purpose, we introduce the following notation:
Let us compose the file function sah568hh.m with Equations (2)–(5):
function F = sah568hh(x);
F = [x(1)^2 + x(2)^2 + x(3)^2 + x(4)^2 − 1;x(1)*x(3) + x(2)*x(4);−x(1) + x(2) − x(3) + x(4);−4*x(1) + 3*x(2) − 2*x(3) + x(4)].
Note that the ultimate goal of solving F is to turn the rows of the matrix equation back to “zero”, which creates certain convenience while working with the function. The function is called at each iteration from the main file sah568h.m, a fragment of which contains the initial approximation vector x0, options for displaying the solution on the display, and the fsolve function. Upon completion of the search for a solution, the vector x and the criterion for evaluating fval are displayed (
Table 1):
% sah568h.m
x0 = [0.2 0.2 0.2 −0.2]’;
options = optimset(‘Display’,‘iter’);
[x,fval] = fsolve(‘sah568hh’,x0,options)
Iteration process
Equation is solved.
Substituting the numerical values of the elements of the vector x into Equation (1), we obtain the determinant of the Daubechies matrix equal to det (DB) = 1.0000, as well as the condition DB^{−1} = DB^{T}, which indicates the accuracy of the solutions performed.
To produce Daubechies wavelets, a greedy algorithm is used, which makes it possible to establish a scaling (refinable) function elementbyelement if known coefficients are available and to select local features of the approximated signals. It is advisable to work with wavelet transforms with a high dimension of the analyzed signals according to the following Algorithm 1:
Algorithm 1 Approximation of the Input Signal 
1. Signal loading. 
2. Spline—an approximation of the input signal. 
3. Assignment of the wavelet decomposition level. 
4. Wavelet selection. 
5. Decomposition of the signal with the determination of the approximating and refining components. 
6. Signal recovery by wavelet inversiontransform. 
7. Estimation of the firstlevel approximation error. 
8. Performing multilevel wavelet decomposition with the selection of approximation and refining coefficients. 
9. Reconstruction of the thirdlevel approximation. 
10. Reconstruction of details of 1, 2, and 3 levels. 
11. Display of the produced components on the display. 
12. Analysis of the results of wavelet transforms and ensuring the noise immunity of signals. 
According to the algorithm, we will perform a wavelet analysis of the characteristics of the pressure change in the cylinder of a marine diesel engine at the nominal mode, obtained experimentally by taking the indicator diagram. To indicate a diesel engine, an indicator was used, which made it possible to record a diagram expanded along the angle of rotation of the crankshaft, which was used for wavelet analysis.
The analysis and synthesis of the diagram were performed using the sah568f.m file in MATLAB codes, where the points of the algorithm being implemented are also indicated by numbers. The file is compiled according to the given Algorithm 2. The data of the indicator diagram after appropriate processing is represented by a matrix of Z dimensions (73 × 2). The increase in the accuracy of the experimental data input is provided by the use of splines with a discreteness step of 0.5 deg. As a result, the pressure vector P of dimension (1 × 721). In our experiment, we will use the popular pressure sensor (OWEN PD100DI sensors models 111, 171, 181) [
12,
13,
14,
15]. The characteristics of the pressure sensor used to obtain experimental data, such as resolution and sampling rate, are presented in the program listing and
Figure 3 and
Figure 4.
Algorithm 2 Wavelet—Process Diagnostics 
% sah568f.m 
% Wavelet—process diagnostics. 
%======================================================= 
% The first step of the algorithm 
load Z 
% The second step of the algorithm. 
fi1 = −160:0.5:200; 

fi = −160:0.5:200; 
% Generation of the vector of Pdimension (1 × 721) 

% Overlay of the interference f on P signal as two spikes. 
a = 0.1; b = 0.3; 

P = P + f; 
v = size(P); 
s = P(1:v(2)); 
ls = length(s); 
% The third and fourth steps in the algorithm 
% Level 1, Daubechies wavelet ‘db1’. Transform algorithm Sm1 = A1 + D1 
[cA1,cD1] = dwt(s,‘db1’); 
% The fifth step in the algorithm 
A1 = upcoef(‘a’,cA1,‘db1’,1,l_s); 
D1 = upcoef(‘d’,cD1,‘db1’,1,l_s); 
% The sixth step in the algorithm 
Sm1 = idwt(cA1,cD1,‘db1’,l_s); 
% The seventh step in the algorithm 
err = max(abs(sSm1)); 
% The eighth step in the algorithm 
[C,L] = wavedec(s,3,‘db1’); 
cA3 = appcoef(C,L,‘db1’,3); 
cD3 = detcoef(C,L,3); 
cD2 = detcoef(C,L,2); 
cD1 = detcoef(C,L,1); 
[cD1,cD2,cD3] = detcoef(C,L,(1,2,3)); 
% The ninth step in the algorithm 
A3 = wrcoef(‘a’,C,L,‘db1’,3); 
D1 = wrcoef(‘d’,C,L,‘db1’,1); 
D2 = wrcoef(‘d’,C,L,‘db1’,2); 
D3 = wrcoef(‘d’,C,L,‘db1’,3); 
Sm3 = waverec(C,L,‘db1’); 
err = max(abs(sSm3)) 
% Error: 
% err = 4.9738 × 10^{−14}. 
% The tenth and eleventh steps in the algorithm 
subplot(2,2,1); plot(fi,A3);grid; title(‘Approximation level 3’); 
subplot(2,2,2); plot(fi,D1);grid; title(‘Refining level 1’); 
subplot(2,2,3); plot(fi,D2);grid; title(‘‘Refining level 2’); 
subplot(2,2,4); plot(fi,D3);grid; title(‘Refining level 3’); 
% The twelfth step in the algorithm 
[thr,sorh,keepapp] = ddencmp(‘den’,‘wv’,s); 
Sm = wdencmp(‘gbl’,C,L,‘db1’,3,thr,sorh,keepapp); 
subplot(2,1,1); plot(fi,P),grid; title(‘Original’) 
subplot(2,1,2); plot(fi,Sm),grid; title(Waveletmodel’) 
With a threestep approximation, as per the scheme
the maximum error was estimated by the formula
The graphs of the coefficients for the threelevel decomposition are shown in
Figure 3.
The inverse wavelet transform is implemented using the signal reconstruction functions provided in the file.
Figure 4 shows the original and recovered signal characteristics. It can be seen that two “spikes” in the original characteristic are also displayed in the model generated by the inverse wavelet transform.