Singular Spectrum Analysis of Tremorograms for Human Neuromotor Reaction Estimation

Abstract

Singular spectrum analysis (SSA) is a method of time series analysis and is used in various fields, including medicine. A tremorogram is a biological signal that allows evaluation of a person’s neuromotor reactions in order to infer the state of the motor parts of the central nervous system (CNS). A tremorogram has a complex structure, and its analysis requires the use of advanced methods of signal processing and intelligent analysis. The paper’s novelty lies in the application of the SSA method to extract diagnostically significant features from tremorograms with subsequent evaluation of the state of the motor parts of the CNS. The article presents the application of a method of singular spectrum decomposition, comparison of known variants of classification, and grouping of principal components for determining the components of the tremorogram corresponding to the trend, periodic components, and noise. After analyzing the results of the SSA of tremorograms, we proposed a new algorithm of grouping based on the analysis of singular values of the trajectory matrix. An example of applying the SSA method to the analysis of tremorograms is shown. Comparison of known clustering methods and the proposed algorithm showed that there is a reasonable correspondence between the proposed algorithm and the traditional methods of classification and pairing in the set of periodic components.

1. Introduction

The COVID-19 pandemic has shown that the development of remote medical care technology is important [1,2,3]. First of all, the telemetric means that provide remote monitoring of the state of human organs and systems are relevant [4]. Important parameters such as oxygen saturation, heart rate, respiration rate, and blood pressure can be transmitted at an extremely low cost. These parameters can be measured by wearable devices connected to an ordinary smartphone [5]. It is possible to implement more complex specially designed measuring devices, which allow obtaining data for monitoring of more specific characteristics.

The next level of the telemedicine development concerns the additional interaction between the patient and the physician through the available means of communication [6,7]. To optimize this, the development of specialized mathematical methods and algorithms for primary information processing are required. Primary processing can provide data analysis using artificial intelligence methods [8,9,10]. This can minimize the volume of transmitted telemetry information and facilitate the work of a physician operating remotely [11].

The control of the state of the human central nervous system is of particular importance now. The development of the pandemic has led to an increase in the number of people under stress [12,13], and as a result, the number of patients with central nervous system disorders grows. This leads to the necessity to develop software and hardware for possible detection and prediction of neurodegenerative diseases.

Recently, methods of estimation of the status of the human central nervous system based on the analysis of its motor activity have been developed. One of the manifestations of human motor activity is a tremor defined as involuntary, rhythmic, and oscillatory movements of body parts. The tremor is caused by sequential or simultaneous contractions of the agonist and antagonist muscles. There are two types of tremors: physiological and pathological [14]. The physiological tremor appears as a low-amplitude tremor with the frequency of 8–12 Hz and can be detected in any healthy person. The pathological tremor is caused by various disorders in the central or peripheral parts of the nervous system [15]. The parameters of the tremor can be used to estimate the functional state of the central nervous system (CNS) and diagnose human neurological pathologies.

There are various instrumental methods for recording and quantitative assessment of a tremor. For example, biomechanical analysis of a tremor can be based on electromyography, accelerometry, spirography, etc. [16]. However, these methods are focused on the analysis of the visible tremor, i.e., when a tremor becomes visible. If there are no external manifestations of a tremor, it is possible to stimulate its appearance. A tremor appears both in people with CNS disorders and in healthy people during formation of the isometric effort, when there is a strong muscle tension without movement. It can occur, for example, when clenching the hand into a fist or holding a heavy object [17,18].

In [19], we proposed an original method of tremor stimulation against the background of isometric effort and developed a device for measuring tremor parameters. The measurement results are represented as tremorograms—time series composed of digitized values of an involuntary force (tremor) arising when the hands are acting on the sensitive elements of the device.

After the measurement, it is necessary to analyze the obtained time series. Spectral analysis based on fast Fourier transformation can be used to solve these problems regardless of the signal type and pathology [20]. However, the Fourier transform has a number of disadvantages that narrow the area of its use in diagnostics. These disadvantages are especially manifested in the analysis of multicomponent signals containing singularities, such as jumps, changes of periods, amplitudes, and phases of harmonic components. Therefore, there is a necessity to use other methods of time series analysis, significantly different from the classical correlation-spectral analysis, which give more informative results [21].

At present, the most commonly used methods to study the principal structural components of time series are singular spectrum analysis (SSA) [22], wavelet transform [23,24], and fractal methods [25,26]. The aim of our study was to offer a method of automatic analysis, which can detect and identify the components of tremorograms, and on the basis of the data obtained, diagnose nervous system disorders.

Wavelet analysis of time series is effective for studying local features of signals and for filtering. Fractal analysis methods allow decomposition of diagnostic signals into sets of periodic signals and chaotic components. Since the purpose of examining tremorograms is to distinguish smoothly changing components, periodic components, and noise, we applied SSA in our work. The SSA method allows decomposition of stationary and nonstationary time series into principal components, and analyzing these components enables us to make conclusions about the characteristics of the original time series [27]. This method is already used when processing some biomedical signals to eliminate noise as well as to detect anomalies in quasi-periodic biosignals such as electrocardiograms, oxygen levels, blood pressure, and electroencephalograms [28,29].

The paper is organized as follows: Section 2 presents the basic method of singular spectrum decomposition and considers the known variants of classification and grouping of principal components. Section 3 is devoted to the application of the SSA method for the analysis of tremorograms. An example of singular spectrum decomposition is shown, different grouping methods are compared, and a grouping algorithm based on singular values of the trajectory matrix is proposed. It is also shown that the results of the proposed method correspond to the results of the known methods. A general algorithm for the analysis of tremorograms is formulated. Section 4 presents the interpretation of the obtained results.

2. Materials and Methods

A tremorogram is a time series F (f₀, f₁, f₂, …, f_N−1) of length N with the regular time step. This series consists of the components that characterize the smoothly varying component Trend(t), harmonic signals Harmonic(t) characterizing involuntary contractions of human muscles on the background of an isometric force, and noise Noise(t), which is caused by the sensor noise and interference in the primary data processing channel:

F(t) = Trend(t) + Harmonic(t) + Noise(t).

To analyze the smoothly varying and harmonic constituents of the tremor, it is necessary to extract the trend and periodic components of the series. For this purpose, we used the singular spectrum decomposition method.

2.1. Singular Spectrum Decomposition Method

A trajectory matrix X was compiled from the original time series as follows. We chose a window of length L, where 1 < L < N, and overlay the window on the time series forming K vectors X_i having the following form:

X_i = (f_i, f_i+1, f_i+2, …, f_i+L−1)^T,

where 0 ≤ i ≤ K − 1 and K = N − L+1.

As a result, we obtained the matrix:

$$\mathrm{X}=\hspace{0.17em}\left[\begin{array}{l}{\mathrm{f}}_0 {\mathrm{f}}_1 {\mathrm{f}}_2 \dots {\mathrm{f}}_{\mathrm{K}-1}\\ {\mathrm{f}}_1 {\mathrm{f}}_2 {\mathrm{f}}_3 \dots {\mathrm{f}}_{\mathrm{K}}\\ {\mathrm{f}}_2 {\mathrm{f}}_3 {\mathrm{f}}_4 \dots {\mathrm{f}}_{\mathrm{K}+1}\\ ⋮ ⋮ ⋮ \dots ⋮ \\ {\mathrm{f}}_{\mathrm{L}-1} {\mathrm{f}}_{\mathrm{L}} {\mathrm{f}}_{\mathrm{L}+1} \dots {\mathrm{f}}_{\mathrm{N}-1}\end{array}\right]$$

This is a Hankel matrix because the elements on the diagonals from bottom left to top right are equal.

Then, we performed a singular value decomposition (SVD) of the trajectory matrix:

$${\mathrm{X}=\mathrm{U}\mathsf{\Sigma }\mathrm{V}}^{\mathrm{T}}$$

where U is a unitary matrix of size L × L containing an orthonormal set of left singular vectors X as columns; Σ is a rectangular diagonal matrix of size L × K containing singular values of X, ordered from largest to smallest; and V is a unitary matrix of size K × K containing the orthonormal set of right singular vectors of X as columns.

The SVD of the trajectory matrix X can be written as

$$\mathrm{X}={\displaystyle \sum _{\mathrm{i}=0}^{\mathrm{d}-1}{\mathsf{\sigma }}_{\mathrm{i}}{\mathrm{U}}_{\mathrm{i}}{\mathrm{V}}_{\mathrm{i}}^{\mathrm{T}}}={\displaystyle \sum _{\mathrm{i}=0}^{\mathrm{d}-1}{\mathrm{X}}_{\mathrm{i}}}$$

where σ_i is the i-th singular value and d is the rank of the trajectory matrix X or the number of singular values satisfying the inequality σ_i > 0.

The columns of the trajectory matrix are a sequence of K vectors that span the space of the original time series. This space is at most L-dimensional. If the columns of the trajectory matrix are linearly dependent, the dimensionality of the space decreases, that is, the part of the matrix of singular values is zero, and the rank corresponds to the inequality d < L.

However, the appearance of linearly dependent columns in the matrix X is practically impossible because the tremorograms have a complex structure and are noisy, that is, the trajectory space has dimension L.

As a result, the original time series F is represented as a set of elementary matrices X_i, which constitute the trajectory matrix X. For the transition X_i to the time series, it is necessary to perform the Hankelization of the matrices. Hankelization is performed by diagonal averaging, which determines the values of the reconstructed time series F_i as the average values of the corresponding antidiagonals of matrices X_i.

To convert a matrix X_i of size L × K into an antidiagonal matrix ${\tilde{\mathrm{X}}}_{\mathrm{i}}$, we calculated x_m,n as the mean value of the elements of the antidiagonal, which includes x_m,n. The number of elements of the antidiagonal to be summed depends on the arrangement of x_m,n in the matrix and is taken into account in the following equation:

$${\tilde{\mathrm{x}}}_{\mathrm{m},\mathrm{n}}=\{\begin{array}{l}\frac{1}{\mathrm{s}+1}{\displaystyle \sum _{\mathrm{l}=0}^{\mathrm{s}}{\mathrm{x}}_{\mathrm{l},\mathrm{s}-\mathrm{l} 0 \le \mathrm{s} \le \mathrm{L}-1}}\\ \frac{1}{\mathrm{L}-1}{\displaystyle \sum _{\mathrm{l}=0}^{\mathrm{L}-1}{\mathrm{x}}_{\mathrm{l},\mathrm{s}-\mathrm{l} \mathrm{L} \le \mathrm{s} \le \mathrm{K}-1}}\\ \frac{1}{\mathrm{K}+\mathrm{L}-\mathrm{s}-1}{\displaystyle \sum _{\mathrm{l}=\mathrm{s}-\mathrm{K}+1}^{\mathrm{L}}{\mathrm{x}}_{\mathrm{l},\mathrm{s}-\mathrm{l} \mathrm{K} \le \mathrm{s} \le \mathrm{K}+\mathrm{L}-2}}\end{array}$$

Based on the Hankelized matrices ${\tilde{\mathrm{X}}}_{\mathrm{i}}$, we could reconstruct time series ${\tilde{\mathrm{F}}}_{\mathrm{i}}$.

2.2. Classification and Grouping

The next stage of the analysis is the classification of matrices ${\tilde{\mathrm{X}}}_{\mathrm{i}}$. It was necessary to determine how many principal components are sufficient to describe the original series, and determine the trend, periodic, and noise components.

According to the results of classification, the matrices ${\tilde{\mathrm{X}}}_{\mathrm{i}}$ could be categorized into non-intersecting sets trend, harmonic, and noise:

$$\tilde{\mathrm{X}}={\displaystyle \sum _{\mathrm{t}\in \mathrm{T}}{\tilde{\mathrm{X}}}_{\mathrm{t}}}+{\displaystyle \sum _{\mathrm{h}\in \mathrm{H}}{\tilde{\mathrm{X}}}_{\mathrm{h}}}+{\displaystyle \sum _{\mathrm{n}\in \mathrm{N}}{\tilde{\mathrm{X}}}_{\mathrm{n}}},$$

where T, H, and N are non-intersecting sets of indexes of elements of the sets trend, harmonic, and noise.

Four sets of objects were obtained according to the results of the singular spectrum decomposition. Any of them could be used to classify and group elementary matrices ${\tilde{\mathrm{X}}}_{\mathrm{i}}$ or time series ${\tilde{\mathrm{F}}}_{\mathrm{i}}$.

• The singular values of the trajectory matrix of the L-dimensional representation of the original one-dimensional series.
• A set of singular vectors of the trajectory matrix. Since their elements are ordered by the operator of the series matrix formation, they can be studied as functions of time.
• The set of principal components ${\tilde{\mathrm{X}}}_{\mathrm{i}}$ of the L-dimensional representation. They, as well as their corresponding eigenvectors, form an orthogonal system and can be considered as functions of the number i.
• Time series ${\tilde{\mathrm{F}}}_{\mathrm{i}}$ reconstructed from different sets of principal components ${\tilde{\mathrm{X}}}_{\mathrm{i}}$.

When classifying components, the components corresponding to the trend and periodic components of tremorograms are of particular interest. The classification of components can be performed visually [30], as well as with the use of methods of intelligent data analysis.

2.2.1. Visual Analysis by Singular Values

When identifying periodic components, it is necessary to take into account the fact that the singular values of the trajectory matrix corresponding to periodic constituents form pairs. The authors in [31] showed that one sinusoidal component in an ideal situation corresponds to two principal components with the same period (sine and cosine), corresponding to one singular value. In a real situation, these two principal components correspond not to one, but to two eigenvalues that are close in value. The analysis of the singular values to determine the proximity of their values can be performed visually, for example, by graph. For contrast, it is recommended to use graphs of logarithms or roots of singular values, as such graphs clearly show the segments similar to the horizontal segments. The ends of these segments correspond to a pair of principal components (${\tilde{\mathrm{X}}}_{\mathrm{i}}$, ${\tilde{\mathrm{X}}}_{\mathrm{i}+1}$) with approximately identical eigenvalues σ_i ≈ σ_i + 1. Their combination makes it possible to reveal the periodic components of the initial series.

2.2.2. Visual Analysis by Eigenvectors

Classification by eigenvectors is based on the visual analysis of two-dimensional graphs, when different pairs of eigenvectors are plotted along the x and y axes. If on the x-and y-axis we plot the values of sines with the same frequency but with different phases, we obtain an ellipse on the plane. From the orthogonality of the eigenvectors for the pair of principal components corresponding to one periodic constituent, it follows that the phase shift between such pairs is ±π/2, and the ellipse becomes a circle.

The harmonic component with an integer period is represented as a regular polygon and the number of vertices of the polygon is determined by the value of the period. When the amplitude of the harmonic component changes, the polygon transforms into a spiral. The pair of principal components (${\tilde{\mathrm{X}}}_{\mathrm{i}}$, ${\tilde{\mathrm{X}}}_{\mathrm{i}+1}$) determining the periodicity can be determined by a visual analysis: the trajectory of such a pair should be approximately represented as a spiral with a fixed center and an insignificantly varying radius.

The visual method can be used when analyzing short and simple time series based on a small number of components.

2.2.3. Classification by Principal Components and Reconstructed Time Series

Clearly, it is necessary to perform a numerical evaluation in order to provide an intelligent analysis of the tremorograms and make a decision about the grouping X_i. Various grouping methods for the singular spectral analysis have been developed [32]. When grouping by principal components or reconstructed time series, known clustering methods can be used [33].

The metrics used in determining the similarity of time series can be divided into two groups: based on a distance [34,35,36] and based on a shape [37,38]. In our research, we used the estimation of distances between objects (time series and matrices). The distance metrics p used in our study is shown in

Table 1. Metrics used to determine the distance between time series A and B.

Metric	Formula for Calculation
Euclidean distance	$\mathrm{p}(\mathrm{A}, \mathrm{B})=\sqrt{{\displaystyle \sum _{\mathrm{i}}^{\mathrm{n}}{{(\mathrm{a}}_{\mathrm{i}}-{\mathrm{b}}_{\mathrm{i}})}^2}}$
Manhattan distance	$\mathrm{p}(\mathrm{A}, \mathrm{B})={\displaystyle \sum _{\mathrm{i}}^{\mathrm{n}}{|\mathrm{a}}_{\mathrm{i}}-{\mathrm{b}}_{\mathrm{i}}|}$
Cosine measure	$\mathrm{p}(\mathrm{A}, \mathrm{B})=\arccos (\frac{\langle \mathrm{A},\mathrm{B}\rangle }{\left|\right|\mathrm{A}\left|\right| \left|\right|\mathrm{B}\left|\right|})$
Chebyshev distance	${\mathrm{p}(\mathrm{A}, \mathrm{B})=\max |\mathrm{a}}_{\mathrm{i}}-{\mathrm{b}}_{\mathrm{i}}|$

.

To group the principal components into pairs, it is necessary to construct a distance matrix ${\mathrm{D}=(\mathrm{d}}_{\mathrm{i},\mathrm{j}}{)}_{\mathrm{i},\mathrm{j}=1}^{\mathrm{d}}$. Based on the values of the distance matrix D, we could estimate the distance between pairs of matrices being compared and determine which of them can be combined into pairs to form a periodic constituent.

For comparison, we used Hankelized matrices ${\tilde{\mathrm{X}}}_{\mathrm{i}}$. To do this, we formed the matrices ${\mathrm{A}}_{\mathrm{i}}{=(\mathrm{a}}_{\mathrm{l},\mathrm{k}}{)}_{\mathrm{l},\mathrm{k}=0}^{\mathrm{L}-1,\mathrm{K}-1}$ as follows:

$${\mathrm{A}}_{\mathrm{i}}=\frac{{\tilde{\mathrm{X}}}_{\mathrm{i}}}{\left|\right|\widetilde{{\mathrm{X}}_{\mathrm{i}}}\left|\right|},$$

where ‖${\tilde{\mathrm{X}}}_{\mathrm{i}}$‖ is the Frobenius norm [39] defined by

$$\left|\right|\mathrm{X}\left|\right|=\sqrt{{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{L}}}{\displaystyle \sum _{\mathrm{j}=1}^{\mathrm{K}}{\mathrm{x}}_{\mathrm{ij}}{}^2}}.$$

The distance between the two components X_i and X_j, can be assessed through the distance between the matrices A_i and A_j. The elements of the distance matrix D can be defined, for example, as ${\mathrm{d}}_{\mathrm{i},\mathrm{j}}=\left|\right|{\mathrm{A}}_{\mathrm{i}}-{\mathrm{A}}_{\mathrm{j}}\left|\right|$.

Another method of similarity analysis of Hankelized matrices was proposed by the authors in [34]. The method is based on the use of w-correlation matrix. This matrix consists of weighted correlations between the reconstructed components of the time series ${\tilde{\mathrm{F}}}_{\mathrm{i}}$. The weights w reflect the number of occurrences of the components of the time series in its trajectory matrix. Well-separated components have a small correlation, while poorly-separated components have a large correlation. Therefore, considering the correlation matrix, we could find pairs of correlated reconstructed series and use this information for grouping.

For two reconstructed time series ${\tilde{\mathrm{F}}}_{\mathrm{i}}$ and ${\tilde{\mathrm{F}}}_{\mathrm{j}}$ of length N and window length L, we could determine the weighted inner product:

$${\left(\widetilde{{\mathrm{F}}_{\mathrm{i}}}, \widetilde{{\mathrm{F}}_{\mathrm{j}}}\right)}_{\mathrm{w}}={\displaystyle \sum _{\mathrm{k}=0}^{\mathrm{N}-1}{\mathrm{w}}_{\mathrm{k}}{\tilde{\mathrm{f}}}_{\mathrm{i},\mathrm{k}}{\tilde{\mathrm{f}}}_{\mathrm{j},\mathrm{k}}},$$

where ${\tilde{\mathrm{f}}}_{\mathrm{i},\mathrm{k}}$ and ${\tilde{\mathrm{f}}}_{\mathrm{j},\mathrm{k}}$ are k-th value of ${\tilde{\mathrm{F}}}_{\mathrm{i}}$ and ${\tilde{\mathrm{F}}}_{\mathrm{j}}$, respectively, and w_k is determined as follows:

$${\mathrm{w}}_{\mathrm{k}}=\{\begin{array}{l}\mathrm{k}+1 0 \le \mathrm{k} \le \mathrm{L}-1\\ \mathrm{L} \mathrm{L} \le \mathrm{k} \le \mathrm{K}-1\\ \mathrm{N}-\mathrm{k} \mathrm{K} \le \mathrm{k} \le \mathrm{N}-1\end{array}.$$

If (${\tilde{\mathrm{F}}}_{\mathrm{i}}$, ${\tilde{\mathrm{F}}}_{\mathrm{j}}$)_w = 0, then ${\tilde{\mathrm{F}}}_{\mathrm{i}}$ and ${\tilde{\mathrm{F}}}_{\mathrm{j}}$ are w-orthogonal, and the components of the time series are separable.

Clearly, the full w-orthogonality does not appear in the analysis of real signals, so we could use the weighted correlation matrix Wcorr to analyze the matrices corresponding to real signals. This matrix of dimension d × d determines the deviation of components F_i and F_j from w-orthogonality. The elements of Wcorr are specified as follows:

$${\mathrm{W}}_{\mathrm{i},\mathrm{j}}=\frac{{\left(\widetilde{{\mathrm{F}}_{\mathrm{i}}}, \widetilde{{\mathrm{F}}_{\mathrm{j}}}\right)}_{\mathrm{w}}}{\left|\right|\widetilde{{\mathrm{F}}_{\mathrm{i}}}{\left|\right|}_{\mathrm{w}} \left|\right|\widetilde{{\mathrm{F}}_{\mathrm{j}}}|{|}_{\mathrm{w}}},$$

where $\left|\right|\widetilde{{\mathrm{F}}_{\mathrm{k}}}{\left|\right|}_{\mathrm{w}}=\sqrt{{\left(\widetilde{{\mathrm{F}}_{\mathrm{k}}}, \widetilde{{\mathrm{F}}_{\mathrm{k}}}\right)}_{\mathrm{w}}}$ for k = i, j.

The interpretation of W_i,j is as follows: if ${\tilde{\mathrm{F}}}_{\mathrm{i}}$ and ${\tilde{\mathrm{F}}}_{\mathrm{j}}$ are close to each other (not identical), then (${\tilde{\mathrm{F}}}_{\mathrm{i}}$, ${\tilde{\mathrm{F}}}_{\mathrm{j}}$)_w → ‖${\tilde{\mathrm{F}}}_{\mathrm{i}}$‖_w ‖${\tilde{\mathrm{F}}}_{\mathrm{j}}$‖_w and therefore, w_i,j → 1. If ${\tilde{\mathrm{F}}}_{\mathrm{i}}$ and ${\tilde{\mathrm{F}}}_{\mathrm{j}}$ are w-orthogonal, then w_i,j = 0. Depending on the data set, some threshold d has to be defined, and if w_i,j ≥ d, the components ${\tilde{\mathrm{F}}}_{\mathrm{i}}$ and ${\tilde{\mathrm{F}}}_{\mathrm{j}}$ can be grouped.

3. Result

In this section, we apply SSA to the analysis of tremorograms, review the grouping algorithm we developed, and compare the results of our algorithm with the results of grouping by known methods.

3.1. Data Acquisition

To obtain data, we applied registration of an involuntary effort on the background of a voluntary isometric effort [19]. The person under the test sat in front of the monitor at the table, and on the table, a device with two strain gauges for the left and right arms’ effort measurement was also placed. With the fingers of the arms straightened forward, the person acted on the corresponding support elements of the meters and observed the value of the generated effort by the arms on displacement of the marks on the monitor. The marks moved along the vertical axis of the screen in proportion to the applied effort. The task of the test person was to hold a certain effort for a given time, combining the marks for the right and left arms at the same level.

The effort of each arm was recorded during the continuous test duration of 10 s. The effort was converted into a voltage, which was transmitted to a computer through an analog-to-digital converter. The voluntary and involuntary force components were stored for subsequent processing. The involuntary force was also used to control the displacement of the marks on the monitor screen in real time. The quantization frequency of the analog signal was 320 Hz. Such measurement settings allowed us to obtain a sufficient number of measurements for statistically reliable estimation of the time series parameters.

Healthy people aged 25–70 years as well as patients with CNS disorders took part in the study. All participants held different forces from 0.5 kg to 5 kg for several minutes with control of efforts via visual feedback.

Figure 1 shows 2 s fragments of tremorograms of two persons under the test. The graph in Figure 1a was obtained when studying the person without detected pathologies, while the graph in Figure 1b is the result of measuring tremor in the patient with a detected CNS pathology. The graphs show the complex structure of isometric force oscillations, which can be decomposed into components correlated with different levels of motor system control.

3.2. Singular Spectrum Decomposition of the Original Series

The singular spectrum decomposition given in Section 2.1 and various variants of principal component classification were implemented using the Python language and relevant libraries: numpy for matrix processing, matplotlib for visualization of arrays used, csv for reading and writing data in the standard csv format.

For processing, we used an original time series consisting of 640 elements corresponding to 2 s recordings. It was shown in [22] that when converting a series into a matrix, the applied window should not be too small, since a more detailed decomposition is obtained when the window is longer. It was shown in [40] that the window should not exceed half of the time series under study. In our research, we used the window length L = 310.

Figure 2 shows the obtained trajectory matrices X1–X4. For the convenience of visual analysis, the matrices are represented by the color equivalents of the corresponding element values.

The ranks d of the trajectory matrices are equal to the length of the window L:

Rank X1 = Rank X2 = Rank X3 = Rank X4 = 310.

Next, consider the process of analyzing the matrix X1. The results of similar operations with matrices X2–X4 are shown in Appendix A.

Figure 3 shows the first 14 matrices X1_i, obtained as a result of decomposition of the trajectory matrix X1. For the convenience of visual analysis, the matrices are represented by color equivalents of the corresponding element values.

From visual inspection of the matrices X1_i, it is obvious that the elementary matrices were not antidiagonal. However, the nature of the components could be determined by their appearance. For example, the values of X1₀, X1₁, and X1₂ changed slowly throughout the matrix, indicating that these components can be related to slowly changing constituents of the time series. The patterns of matrices X1₃–X1₆ showed periodic changes. X1₇–X1₁₀ also characterized periodic components, but their frequencies were different from the frequencies of the components corresponding to X1₃–X1₆. The remaining components refer to noise.

Let us estimate the relative contribution (RC) of the decomposition matrices X1_i to the trajectory matrix X1 using the values of the singular values σ_i as follows:

$${\mathrm{RC}}_{\mathrm{i}}=\frac{{\mathsf{\sigma }}_{\mathrm{i}}^2}{{\displaystyle \sum _{\mathrm{k}=0}^{\mathrm{d}-1}{\mathsf{\sigma }}_{\mathrm{k}}^2}}.$$

The relative contribution for the first 20 values of σ is shown in

Table 2. The relative contribution of the first 20 decomposition matrices X1_i.

i	RC	i	RC	i	RC	i	RC
0	0.5028	5	0.0064	10	0.0047	15	0.0015
1	0.3859	6	0.0063	11	0.0045	16	0.0013
2	0.0211	7	0.0052	12	0.0038	17	0.0010
3	0.0168	8	0.0051	13	0.0030	18	0.0010
4	0.0164	9	0.0048	14	0.0016	19	0.0008

.

Analysis of the data in the Table 2 showed that elementary matrices X1₀ and X1₁ contributed 50.28% and 38.59% to matrix X1, respectively. Together, the first 20 elementary matrices contributed 99.39%:

$${\displaystyle \sum _{\mathrm{i}=0}^{19}{\mathrm{RC}}_{\mathrm{i}}}=0.9939.$$

For the following analysis, we used only the first 20 components of the series, considering the contribution of the others insignificant.

3.3. Classification of Principal Components

Elementary matrices with close singular values, i.e., when σ_i ≈ σ_i + 1, contributed equally to the matrix X1 and could be grouped together when reconstructing of the time series. For example, Table 2 shows that X1₃ and X1₄, as well as X1₅ and X1₆, should be grouped together. The pairs X1₇ and X1₈ as well as X1₉ and X1₁₀ should be investigated further, since they can refer to both periodic and noise. Similar data were obtained when analyzing the graph of the singular values shown in Figure 4.

Figure 4 shows that the singular values corresponding to the principal components $\tilde{\mathrm{X}}{1}_3$ and $\tilde{\mathrm{X}}{1}_4$, $\tilde{\mathrm{X}}{1}_5$ and $\tilde{\mathrm{X}}{1}_6$, $\tilde{\mathrm{X}}{1}_7$ and $\tilde{\mathrm{X}}{1}_8$, and $\tilde{\mathrm{X}}{1}_9$ and $\tilde{\mathrm{X}}{1}_{10}$ formed segments parallel to the x-axis. These pairs of singular values probably determine harmonic components of the series.

Let us focus on the problem of separating the signal components from the noise components, e.g., the pair $\tilde{\mathrm{X}}{1}_7$ and $\tilde{\mathrm{X}}{1}_8$, and $\tilde{\mathrm{X}}{1}_9$ and $\tilde{\mathrm{X}}{1}_{10}$. Firstly, the irregular behavior of singular vectors can indicate that they belong to the set generated by a noise component (these cases should not be confused with mixing of components due to the absence of strong separability of series). Secondly, a slow decrease of the singular values starting from a certain number almost without jumps also points to this. Finally, the large set of singular values generating correlated reconstructed components most likely refers to noise.

Starting from the value σ₉, the graph decreased smoothly, which indicates the beginning of the appearance of the noise components of the series, and therefore, the components from σ₉ should be referred to noise.

It is known from the SSA theory that singular values belonging to slowly varying components have the highest weights, i.e., they are in the beginning of the list of singular values (since the singular decomposition sorts by decreasing their weights). Harmonics correspond to pairs of adjacent vectors that are on the same “step”, i.e., have similar eigenvectors. It is also noted that the distance between pairs of eigenvalues of interdependent harmonics is many times greater than the distance in the pairs themselves between the eigenvalues. Noise is defined as a slow decrease of singular values.

Such behavior of singular values is typical for tremorograms of healthy people. SVD decomposition of tremorograms of patients with CNS pathologies shows a different pattern. Single singular values are observed between explicit pairs corresponding to harmonic components. In this case, we can attribute it to slowly changing components. There may appear grouped pairs in the first singular values. Such cases are shown in Figure A2b and Figure A3b in Appendix A.

Using these observations, we could perform classification into trend, harmonic, and noise components according to the following Algorithm 1.

Algorithm 1 Principal component classification algorithm

Input: Set of singular values σi Output: Trend, harmonic, and noise sets STEP 1: Calculate the set of logarithms of the singular values σi.

STEP 2: Determine the threshold value dd for the distance between logarithms of neighboring singular values considered as pairs, and the scaling factor c for determining the distances between pairs.

STEP 3: Analyze sequentially the pairs of neighboring values from the set of logarithms of singular values. If the distance between values in the pair is greater than dd, then the first point is added to the trend set. Else, we assume that the points can form a pair; to identify the pair, go to Step 4. STEP 4: Calculate the distance between the second point of a pair and the one next to it. If the distance is greater than dd*c, then the pair is found. Return to step 3. If the distance is less than dd*c, we assume that the components corresponding to the analyzed singular values and all subsequent ones are noisy.

When dd = 0.05 and c = 2 the proposed algorithm classified the singular values into three groups:

• Trend: σ₀, σ₁, σ₂;
• Harmonic: σ₃, σ₄; σ₅, σ₆; σ₇, σ₈;
• Noise: σ₉, and so on.

Another confirmation of the validity of the grouping are the two-dimensional scattergrams shown in Figure 5. A scattergram is a graphical representation of selected pairs of eigenvectors U_i and U_i+1 obtained as a result of SVD. When the eigenvectors are orthogonal, i.e., correspond to the same harmonic, the scattergram is a spiral (possibly fictitious).

Clearly, the presented scattergrams of the principal components $\tilde{\mathrm{X}}{1}_3-\tilde{\mathrm{X}}{1}_8$ are the most informative. The subsequent components probably belong to the noise components of the series and were excluded from the analysis. The first three components $\tilde{\mathrm{X}}{1}_0-\tilde{\mathrm{X}}{1}_2$ were also excluded from the analysis, as they belong to the smoothly changing components.

Thus, for the analysis of time series for its periodic components will be used three pairs of components: $\tilde{\mathrm{X}}{1}_3-\tilde{\mathrm{X}}{1}_8$.

A similar result was obtained by visual analysis of the Hankelized matrices (principal components) ${\tilde{\mathrm{X}}}_{\mathrm{i}}$. According to the appearance of the color representation of the principal components shown in Figure 6, the frequency components corresponded to matrices $\tilde{\mathrm{X}}{1}_3-\tilde{\mathrm{X}}{1}_8$. They had a periodic antidiagonal pattern. The pairing was performed on the basis of frequency.

Generally, there are no recommendations on the number of components determining the reconstructed series. In specific tasks, this number is chosen from compromise considerations, based on the desire to filter out the series from the noise components as best as possible and the need to preserve the informative components.

It is important to note that elementary matrices represent an optimal (perhaps not unique) separation of components in the trajectory space: by definition, the rows and columns of one elementary matrix are orthogonal to the rows and columns of other elementary matrices. However, this separation is not always interpretable. There are restrictions on the types of components of the time series, which are exactly separated according to the chosen formal representation.

Figure 7 shows the results of estimating the similarity of the principal components based on the distances between the time series calculated using different metrics.

The correlation between time series is clearly visible on the matrix formed on the basis of the cosine distance. The following pairs $\{\tilde{\mathrm{X}}{1}_3,\tilde{\mathrm{X}}{1}_4\}$, $\{\tilde{\mathrm{X}}{1}_5,\tilde{\mathrm{X}}{1}_6\}$, $\{\tilde{\mathrm{X}}{1}_7,\tilde{\mathrm{X}}{1}_8\}$, and $\{\tilde{\mathrm{X}}{1}_9,\tilde{\mathrm{X}}{1}_{10}\}$ were distinguished.

Figure 8 shows the distances between the principal component matrices calculated from the normalized matrices as well as from the w-correlation matrix. The analysis of the matrices in Figure 8 showed a similar grouping result. We also observed correlations between components of the time series, especially in the range for $\tilde{\mathrm{X}}{1}_9-\tilde{\mathrm{X}}{1}_{10}$. These components were classified as noise by the proposed algorithm, so it is not surprising that there were some correlations between them, but they were smaller than the correlations between pairs of periodic components.

From the results above, we could see that there is a reasonable correspondence between the proposed algorithm and the traditional methods of classification and pairing in the set of periodic components, at least for the first most significant principal components. This fact is illustrated in

Table 3. Comparison of the results of component grouping by the considered methods.

Method	Component Grouping Result
Visual analysis of singular values (Figure 4)	Trend: {σ0, σ1, σ2} Harmonic: {σ3, σ4; σ5, σ6; σ7, σ8; σ9, σ10} Noise: {σ11–σ19}
Visual analysis of scattergrams of pairs of eigenvectors (Figure 5)	Harmonic: {σ3, σ4; σ5, σ6; σ7, σ8}
Classification based on reconstructed time series using cosine measure (Figure 7)	Trend: {σ0, σ1, σ2} Harmonic: {σ3, σ4; σ5, σ6; σ7, σ8; σ9, σ10} Noise: {σ11–σ19}
Classification based on the w-correlation matrices of the principal components (Figure 8a)	Trend: {σ0, σ1, σ2} Harmonic: {σ3, σ4; σ5, σ6; σ7, σ8; σ9, σ10} Noise: {σ11–σ19}
Classification based on the correlation of the principal component matrices using the Frobenius norm (Figure 8b)	Trend: {σ0, σ1, σ2} Harmonic: {σ3, σ4; σ5, σ6; σ7, σ8; σ9, σ10} Noise: {σ11–σ19}
Proposed algorithm	Trend: {σ0, σ1, σ2} Harmonic: {σ3, σ4; σ5, σ6; σ7, σ8; σ9, σ10} Noise: {σ11–σ19}

, which shows the results of the time series analysis obtained by analyzing a tremorogram of the left arm of a patient without detectable pathology.

When performing a grouping of components based on distances between time series using metrics, we used the cosine measure, because it is difficult to conclude confidently using other variants (Euclidean distance, Manhattan distance, and Chebyshev distance).

Table 4. Comparison of the results of component grouping by the considered methods for different tests.

Method	Hankelized Principal Components Matrices X1	Hankelized Principal Components Matrices X2 (Figure A1)	Hankelized Principal Components Matrices X3 (Figure A2)	Hankelized Principal Components Matrices X4 (Figure A3)
Visual analysis of singular values	{σ0, σ1, σ2} {σ3, σ4; σ5, σ6; σ7, σ8; σ9, σ10} {σ11–σ19}	{σ0, σ1, σ2, σ3} {σ4, σ5} {σ6–σ19}	{σ6, σ7} {σ0, σ1; σ2, σ3; σ4, σ5; σ8, σ9; σ10, σ11} {σ12–σ19}	{σ4} {σ0, σ1; σ2, σ3; σ5, σ6} {σ7–σ19}
Classification based on the w-correlation matrices of the principal components	{σ0, σ1, σ2} {σ3, σ4; σ5, σ6; σ7, σ8; σ9, σ10} {σ11–σ19}	{σ0, σ1, σ2, σ3, σ6, σ7} {σ4, σ5; σ8, σ9} {σ10–σ19}	{σ6, σ7} {σ0, σ1; σ2, σ3; σ4, σ5; σ8, σ9} {σ10–σ19}	{σ4} {σ0, σ1; σ2, σ3; σ5, σ6} {σ7–σ19}
Classification based on the correlation of the principal component matrices using the Frobenius norm	{σ0, σ1, σ2} {σ3, σ4; σ5, σ6; σ7, σ8; σ9, σ10} {σ11–σ19}	{σ0, σ1, σ2, σ3, σ6, σ7} {σ4, σ5; σ8, σ9} {σ10–σ19}	{σ6, σ7} {σ0, σ1; σ2, σ3; σ4, σ5; σ8, σ9} {σ10–σ19}	{σ4} {σ0, σ1; σ2, σ3; σ5, σ6} {σ7–σ19}
Proposed algorithm	{σ0, σ1, σ2} {σ3, σ4; σ5, σ6; σ7, σ8; σ9, σ10} {σ11–σ19}	{σ0, σ1, σ2, σ3} {σ4, σ5} {σ6–σ19}	{σ6, σ7} {σ0, σ1; σ2, σ3; σ4, σ5; σ8, σ9; σ10, σ11} {σ12–σ19}	{σ4} {σ0, σ1; σ2, σ3; σ5, σ6} {σ7–σ19}

shows the results of comparison of the proposed algorithm and the most representative methods: visual analysis of singular values [31]; classification using the w-correlation principal component matrices [32]; and classification based on correlation by principal component matrices using the Frobenius norm [34]. In the analysis of matrix X1, we used Figure 4 and Figure 8a,b, and in the analysis of matrices X2, X3, and X4 we analyzed the Figure A1, Figure A2 and Figure A3 of Appendix A, respectively.

3.4. Generalized Algorithm for the Analysis of Tremorograms

Based on the results of the studies, we propose a general algorithm for analyzing tremorograms, and its block diagram is shown in Figure 9.

The input data for the algorithm are tremorograms. They are a time series of amplitudes of involuntary human efforts formed against the background of isometric efforts under the control of biofeedback.

At the first stage, a trajectory matrix is formed. We know [39] that the SSA method is stable in relation to changes of the window length, and the decomposition effect appears not so much in a numerical as in a qualitative sense. To obtain the most accurate decomposition, we recommend choosing a window length equal to half of the time series length.

At the second stage, the trajectory matrix is decomposed by singular values and a vector of singular values is formed. The singular values are ordered from the largest to the smallest.

At the third stage, the singular values are analyzed using the proposed algorithm. As a result, the values related to trend, harmonic components, and noise are determined.

As a result of the analysis, we obtained three sets trend, harmonic and noise, which contain indexes of singular values. Hankelized matrices with such indices correspond to the principal components that determine the trend and periodic components of the time series.

4. Discussion

The classification of singular spectrum decomposition components based on singular values of the trajectory matrix gave results comparable to those obtained by other methods. It was determined that the classification by singular numbers is easy to implement, the analysis algorithm has low complexity, and it gives the correct result for the interpretation. An important advantage is that the algorithm does not require prior information about the number of pairs of singular numbers belonging to the harmonic components.

To test the proposed method, we used 60 tremorograms obtained from healthy people and patients with CNS pathologies using the device [19], as well as sets of specially prepared synthetic data. The testing was performed for two purposes. The first purpose was checking the quality of selection of components in the initial time series (robustness of the proposed algorithm). The second purpose was to assess the quality of medical diagnosis (application of SSA method in combination with the proposed algorithm). In both cases, synthetic data were additionally used for more complete testing. Such data are often applied when testing software products and it gives good results when test suites are adequately formed [41].

To test the robustness of the proposed algorithm, values corresponding to decreasing singular numbers were used as input data. Both sets with paired values and sets consisting only of single points were generated. Testing on such sets showed the performance of the proposed algorithm.

When assessing the quality of diagnosis, the test data were generated in the form of initial time series as a combination of smoothly changing components corresponding to isometric effort, frequency components with certain amplitude–frequency characteristics of functional or pathological tremor according to [16,17], and noise components. In this case, the possibility of detecting the components indicative of CNS pathology was evaluated and the diagnostic value of the proposed analytical method was verified.

Analysis of real tremorograms obtained from healthy people and patients with CNS pathologies revealed some peculiarities. The tremorograms of healthy persons were characterized by the following distribution of singular values, arranged in descending order: several numbers, corresponding to slowly changing components; some number of pairs, corresponding to harmonic signals; and singular values of noise components. The appearance of unpaired singular values between pairs corresponding to harmonic components, as well as the localization of harmonic components in the initial positions of the set of singular values arranged in descending order, indicated a CNS pathology in the test person.

Such a result is medically justified. In healthy people, involuntary oscillations formed against the background of isometric effort have a small amplitude and do not interfere with movement performance [42]. Functional tremor occurring in healthy persons makes a contribution to the tremorogram no greater than that of the smoothly varying force. At the same time, the contribution of periodic components decreases with increasing frequency. In patients with CNS pathologies, the contribution of frequency components is significant; they begin to manifest themselves more actively than the smoothly varying effort. At the same time, the singular values of trend components and pairs corresponding to periodic components are intermixed. Thus, based on the results of the classification of singular values, we can determine the pathological states of the human CNS.

As the analysis of publications related to the application of the SSA method has shown, SSA-based methods of time series analysis in most cases involve visual analysis at one or more stages [30,32,33,34,35]. A variant of automation was shown in [34], where agglomerative hierarchical clustering methods were used to perform operations on matrices. In our algorithm, the analysis was performed on vectors, as it is simpler and requires less computational resources.

In some cases, automatic identification of pairs of singular values gives worse results than visual identification performed by a specialist. Moreover, the proposed algorithm, as well as other known methods of automatic identification, is not fully automated, because they require the definition of thresholds. A number of problems arise in selecting the threshold. First of all, we must be sure that there exist such optimal thresholds that provide component recovery accuracy no worse than visual analysis. Secondly, we need to develop an approach for automatic threshold determination that allows us to clearly classify pairs of singular values belonging to the same periodic component.

The coefficients, which are constants for this algorithm, were selected by hand, based on the results of manual testing. These coefficients will not change significantly during the study of new results, and their possible insignificant correction will only improve the classification accuracy. We plan to solve these problems in future studies. Additionally, in the future, we plan to form the rules of differential diagnosis, when not only the presence of a functional disorder of the central nervous system is detected, but also the type of pathology is determined.

5. Conclusions

The paper showed that the method of singular spectral analysis allows us to qualitatively determine the composition of the tremorogram for the analysis of its components.

The practical result is demonstration of an opportunity for the use of singular spectral analysis for the investigation of various elements of tremorograms and detection of a number of physiologically significant phenomena by tremorography. Their further study is important in the analysis of human neuromotor reactions.

At the current stage, the prepared general algorithm of the tremorogram analysis can be used together with some additional method of visual identification, for example, with w-correlation matrix. On the one hand, this allows additional control of the result of the algorithm to have maximal confidence in the correctness of the decision about grouping. On the other hand, we obtain sets of singular values and the results of their analysis. These data can be further used as marked sets, for example, in clustering singular values with neural networks.