Algorithm

Singular Spectrum Decomposition (SSD) algorithm is highly based on the Singular Spectrum Analysis (SSA) algorithm, which has been shown suitable for short and noisy time series analysis. SSA is a principal component analysis-based nonparametric spectral estimation method, which uses no prior knowledge about the biology and physics of the system. The main advantage of SSA is that its components are not necessarily harmonic functions and can capture highly non-harmonic oscillatory shapes. One main disadvantage is that both the window length, or embedding dimension, and the principal components to choose for reconstructing a specific component series have to be selected manually. Those two processes have been automated in the SSD algorithm.

First step of SSD algorithm is finding the trajectory matrix X. With respect to standard way of building it, a new approach is proposed.
Given a time series x(n) of length N and an embedding dimension M, an ( M x N ) matrix X is generated such that its i-th row (where i = 1,…,M) is obtained as:
$x_{i} = (x(i),\:.\:.\:.\:,x(N),\:x(1),\:.\:.\:.\:,\:x(i-1))$ Hence: $X= [ x_{{1}}^{T}, x_{{2}}^{T},\:.\:.\:.\:, x_{{M}}^{T}]^{T}$

For example, given the time series $x(n)= \left \{ 1,2,3,4,5 \right \}$ , and an embedding dimension M=3, then the matrix X is:

$\ X = \left| \begin{array}{ccccc} \ 1 & 2 & 3&4&5 \\ 2&3&4&5&1 \\ 3&4&5&1&2 \end{array} \right|$

In the SSA algorithm matrix X would correspond to the 3 x 3 left block of the above matrix:

$\ X = \left| \begin{array}{ccc} \ 1 & 2 & 3 \\ 2&3&4 \\ 3&4&5 \end{array} \right|$

Thanks to this new definition of the trajectory matrix in the SSD method, the oscillatory content of X is enhanced and useful properties for the decrease of energy of the residual are provided.

Next, the diagonal averaging procedure needs to be adapted to this new definition of trajectory matrix. In order to carry out the average along the i-th cross-diagonal of X, the wrapped part of the right hand block must be correctly appended to the top right of the left hand block. With this approach every cross-diagonal has length M. An example is shown below:

$\ X = \left| \begin{array}{ccccc} \ &&1&& \\ &1&2&& \\ 1 & 2 & 3&4&5 \\ 2&3&4&5& \\ 3&4&5&& \end{array} \right|$

However, care has to be taken in case of a sizable trend in the series. To avoid that situation at first iteration, a test for sizable trend is performed. In case the sizable trend is detected, then a standard trajectory matrix is generated to guarantee reliable estimation of the trend. The X found by that method is used in all the following iterations.

Another vital step is choosing the embedding dimension. The following criterion has been adopted for the automated choice of M at iteration j:

the Power Spectral Density (PSD) of the residual time series at iteration j, $v_{j}(n)=x(n) - \sum^{j-1}_{k=1} g_{k}(n)$ , is computed, where $g_k(n)$ is the k-th component series estimated by SSD at iteration k. The frequency in its PSD associated with the dominant peak, $f_{{max}}$ is then estimated.
at the first iteration, if the normalized frequency $f_{{max}}/F_{{s}}$ (With $F_{{s}}$ being the sampling frequency), is less than a given threshold (default equal to $10^{{-3}}$ in the implementation), a sizable trend is assumed to characterize the residual. M is then set to $\frac{N}{3}$ in case of a sizable trend present in the time series.
Otherwise and for iterations j > 1, the window length is defined as $M=1.2*F_{s}}/{f_{max}}$ . The factor 1.2 allows the window length to cover a time span 20% larger than the average period of the wanted component series, to improve its identification by SSA.

The adaptive choice for M guarantees SSD to provide adaptive spectral filters associated with the dominant oscillations in the residual series, exploiting the property of SSA of concentrating the transfer functions related to the dominant eigenfunction within regions where sharp peaks in the power spectrum occur, facilitating the detection of such oscillations.

Further step is reconstruction of the j-th component series, $g_j(n)$ , which is carried out as follows: at the first iteration, if a sizable trend has been detected, only the first left and right eigenvectors are used to obtain $g_1(n)$ , such that $\b X_{{1}}$ $= \sigma_{{1}}u_{1}{v_{1}}^{T}$ and $g_1(n)$ is obtained from diagonal averaging of $\b X_{{1}}$ . Indeed, as a sizable trend is the component with the highest energy in the analyzed signal, is should be mainly reflected in the first eigentriple provided by SSA.

Otherwise and for iterations j > 1, a component series $g_j(n)$ has to be retrieved such as to describe a well defined time scale. In this sense, its frequency content is to be concentrated in the frequency band [ $f_{{max}}$ $- \delta f$ , $f_{{max}}$ $+ \delta f$ ] with $\delta f$ representing half the width of dominant peak in the PSD of the residual time series $v_{{j}}(n)$ . Hence, a subset $I_{j}$ is created from all eigentriples whose left eigenvectors show a dominant frequency in their spectra in the range [ $f_{{max}}$ $- \delta f$ , $f_{{max}}$ $+ \delta f$ ] and the one eigentriple which contributes the most to the energy of the dominant peak of the selected mode. The corresponding component series is then reconstructed by diagonal averaging of the matrix $X_{Ij}$ along the cross-diagonals.

Next the width of the dominant peak $\delta f$ has to be estimated. A spectral model defined by a superposition of Gaussian functions is used for describing PSD profile, The model is defined as a sum of three Gaussian functions, each function representing a spectral peak:

$\gamma(f,\theta)=\sum_{i=1}^{3} A_ie^{-\frac{(f-\mu_i)^2}{2\sigma_i^2}}$

where $A_i$ is the magnitude of the $i^{th}$ Gaussian $\sigma_i$ its width $\mu_i$ its location and $\theta=$ $\left [ A\sigma \right ]^T$ is the parameter vector, with A=[ $A_1$ $,A_2$ $,A_3$ ] and $\sigma=$ $[\sigma_1,\sigma_2,\sigma_3]$ . Then the optimization procedure fits the model to the entire PSD using the following initial parameter values (indexed as $(\cdot)^{(0)}$ ):

$A _{1}^{(0)}=\frac{1}{2}PSD(f_{max})\:\:\:\:\:\:\:\:\:\:\:\:\sigma{_{1}}^{(0)}= f:PSD(f)=\frac{2}{3}PSD(f_{max})$ $A _{2}^{(0)}=\frac{1}{2}PSD(f_{2})\:\:\:\:\:\:\:\:\:\:\:\:\sigma{_{2}}^{(0)}= f:PSD(f)=\frac{2}{3}PSD(f_{2})$ $A _{3}^{(0)}=\frac{1}{4}PSD(f_{3})\:\:\:\:\:\:\:\:\:\:\:\:\sigma{_{3}}^{(0)}= 4|f_{max}-f_{2})|$

Next, the optimal values are determined using nonlinear minimization algorithm, such as the Levenberg-Marquardt algorithm. Given the estimated value for $\sigma_1$ , a value for $\delta f$ is then obtained as $\delta f=2.5 \sigma_1$ ,

The choice for $X_{Ij}$ avoids the need of identifying the leading eigenvalues with respect to the noise floor. Hence, components describing different scales are automatically discarded as iteration j and left to be identified at future iterations. The underlying hypothesis is that at each iteration SSD picks up the oscillation describing most of the energy in the residual series, under the assumption that SSA tries to find spectral bands with a high percentage of explained energy.

Moreover, a second run of the algorithm is carried out on the j-th component series just retrieved $g_j^{[1]}(n)$ (for j = 1, the second run is performed in case of no sizable trend detected). This allows to achieve a better estimate of $g_j^{[1]}(n)$ by polishing up the result of the first run.

Finally, a scaling factor a is applied to $g_j(n)$ such as to adjust its variance to the residual $v_j(n)$ , s.t.:

$\hat{a}=\min_{a}||v_j(n)-ag_j(n)||^2_{2}\:,\:giving\: \hat{a}=\frac{g^Tv}{g^Tg}\: and \:\overset{\sim}{g}_j}=\hat{a} g_j(n)$

The process is stopped when the stopping criterion is met. The normalized mean square error (NMSE) between the residual and the original signal is computed and the decomposition process is stopped when the NMSE is less than a given threshold (default th = 0.1%).

$NMSE^{(j)}=\frac{\sum_{i=1}^{N}(v^{(j+1)}(i))^2}{\sum_{i=1}^{N}(x(i))^2} \:where \:\:\:v^{(j+1)}=v^{(j)}-\overset{\sim}{g}^{(j)}(n)$

Then, the final decomposition gives:

$x(n)=\sum_{k=1}^{m}\overset{\sim}g_k(n)+v_{m+1}(n)$

where m is the number of component series

In order to download the article visit this link of the Maastricht University Library:
SSD taverne version