On the use of Singular Spectrum Analysis
Abstract
Singular Spectrum Analysis (SSA) or Singular Value Decomposition (SVD) are often used to denoise univariate time series or to study their spectral profile. Both techniques rely on the eigendecomposition of the correlation matrix estimated after embedding the signal into its delayed coordinates. In this work we show that the eigenvectors can be used to calculate the coefficients of a set of filters which form a filter bank. The properties of these filters are derived. In particular we show that their outputs can be grouped according to their frequency response. Furthermore, the frequency at the maximum of each frequency response and the corresponding eigenvalue can provide a power spectrum estimation of the time series. Two different applications illustrate how both characteristics can be applied to analyze wideband signals in order to achieve narrowband signals or to infer their frequency occupation.
1 Introduction
Projective subspace models, applied to time series data sets, can be found in literature under various names depending on the domain of application: Singular Spectrum Analysis (SSA) (for instance in climate time series analysis) [10], [9], [24] or Singular Value Decomposition (SVD) (for instance in speech enhancement [7, 11, 15, 14] and more recently in spectrum sensing [4, 16]). The aim of SSA is to decompose a time series into a sum of a small number of interpretable components such as a slowly varying trend, oscillatory components and superimposed noise interferences. Applied to speech enhancement, SSA simply serves to eliminate noise which is usually considered additive, independent and normally distributed. In spectrum sensing applications, the eigenvalue profile is used to decide wether a communication channel is occupied or not, or simply to estimate the noise floor of the channel. Thus spectrum sensing applications are concerned mainly with the study of the spectral profile of the signals. Recently, cognitive radio applications give another perspective to the application of these singularvaluebased decomposition methods. In general, the issue is to detect the presence of a signal (user) in a given narrowband channel [27].
SSA also has been applied to process biomedical signals with different goals: eliminating highamplitude artifacts [5], suppressing noise contributions or extracting informative components [22, 8]. The elimination of highamplitude artifacts is naturally related to the components associated with the largest eigenvalues. Therefore, the number of eigenvalues is always discussed but in most of the applications eliminating the component related with the largest eigenvalue reduces significantly the artifactrelated interference without distorting the underlying signal [5, 23]. On the contrary, reducing white noise interferences usually is addressed by eliminating components related with small eigenvalues [8]. The extraction of components and their grouping according to some criterion has been addressed also by recognizing different trends on the eigenvalue spectra [8], or based on certain statistics of the extracted components [22]. Note that in most biomedical applications, as in many other fields in science as well, the power spectra show a dependence in the low frequency range which is matched by an eigenvalue distribution arranged in decreasing order of magnitude [20].
But in a more general perspective, such subspace models allow to decompose a time series into several components which reflect distinct frequency bands of the original signal. In general, SSA [10]  or SVD [7]  based subspace approaches can be summarized with the following steps :

Embedding the time series into its timedelayed coordinates resulting in a sequence of lagged multidimensional data vectors. The latter are arranged as columns of a trajectory matrix with either Toeplitz or Hankel structure.

Estimating an orthogonal basis vector matrix, using singular value decomposition or principal component analysis (PCA).

Projecting the multidimensional data vectors onto the new basis vectors.

Selecting relevant components.

Reconstructing the multidimensional embedded data using a generally lowerdimensional subspace representation.

Employing diagonal (or antidiagonal) averaging to reconstitute the Toeplitz or Hankel structure of the reconstructed trajectory matrix.

Reverting the embedding to yield a univariate time series.
Thus, the SSA  or SVD  based subspace model corresponds to an orthogonal matrix whose column vectors form an eigenbasis of the multidimensional space created by the embedding. Relevant components of the signal can be obtained by projecting the data onto that eigenbasis and omitting irrelevant, for example noise related, components. Finally an improved univariate time series can be reconstructed.
The selection of relevant components is a critical issue that needs to be addressed each time a problemspecific subspace model is generated. For instance, in [1, 2], the frequency profile of the basis vectors was studied to select those vectors (and related components) which correspond to the trends in the data. In [25] it is shown that SSA can equivalently be formulated as a parallel filter bank structure where the finite impulse response (FIR) filters represent the eigenvectors of the correlation matrix. During the analysis step, the filter bank decomposes the input signal into components using the eigenvectors as causal (or anticausal) FIR filters. Afterwards, the multidimensional reconstruction, including the subsequent diagonal averaging, is correspondingly explained as an anticausal (causal) FIR filter. The latter filter coefficients simply correspond to the eigenvector entries written in reversed order. The interpretation of a signal enhancement as a linear filtering operation was discussed in [12, 11], and in [25], the frequency response of the filters is deduced in closed form. Using techniques from linear timeinvariant systems theory to compute input  output relationships [17, 18], a closedform expression for the frequency response of the filters is obtained such that the frequency profile of the output can be extracted easily. Two important issues are relevant in SSA: the value of the embedding dimension, which is an user predefined parameter, and the choice of relevant components.
In this work we show that a study of the frequency response of the equivalent FIR filters favorably complements a study of the eigenvalue spectrum and a related choice of the number of relevant eigenvalues [3]. It is also shown that the eigenvalue profile of the input signal is related to its power spectrum if the eigenvalues are ordered according to the maximal value in the absolute value of the frequency response of the corresponding filter. We illustrate the application of these ideas in two typical applications where such methods may be employed advantageously.
The first example concerns a biomedical application, more specifically an Electroencephalogram (EEG). The latter represents recordings of electrical brain activities via several sensors thus creating a multichannel biomedical signal that is often used for medical diagnosis or in brain studies. In both cases, artifacts frequently turn the visual analysis difficult and often even impossible. Therefore, preprocessing techniques like digital filtering are often applied to extract the characteristic bands of the EEG signals: delta (), theta (), alpha (), beta () and gamma (). Moreover, as nowadays the acquisition boards work with high sampling rates, the recorded input signals have a bandwidth which is larger than what is really useful in most of these studies. Hence, as is shown in this study, extracting only the problemrelated informative components of the signals represents a preprocessing technique that can be used profitably in EEG analysis.
The second application concerns a telecommunication application. Cognitive Radio (CR) is a recent technique the goal of which is to make use of at times unoccupied Radio Frequency (RF) bands through intelligent spectrum sensing [26]. Therefore, the study of spectrum analysis methods regained interest, especially studying methods which are able to cope with very low signaltonoise ratios. These algorithms need to be able to accurately detect the presence or absence of an active user. To detect the occupation of a channel, usually a decision threshold is defined based on the estimation of the noisefloor. One of the available methods to detect the presence of an active user is based on the dynamic range of the eigenvalue spectrum, i. e. the ratio between the maximal and minimal eigenvalues of the correlation matrix of the sensed signal [27].
The paper is organized as follows: in section 2 SSA and its representation as a filter bank is presented before in section 3 some illustrative applications are discussed and finally, in section 4, some conclusions are drawn.
2 Singular Spectrum Analysis
Time series analysis techniques often rely on embedding onedimensional sensor signals in the space of their timedelayed coordinates. Embedding can be regarded as a mapping that transforms a onedimensional time series into a multidimensional sequence of lagged vectors. Considering an univariate signal
with samples, its multidimensional variant is obtained by , where . These lagged vectors form the columns of the related data matrix which is called a trajectory matrix [10]. The column vectors of span an embedding space of dimension .
(1) 
Note that the trajectory matrix has identical entries along its diagonals. Such a matrix is called a Toeplitz matrix. There are other alternatives to organize a data matrix via embedding the univariate time series signal in an  dimensional space of its timedelayed coordinates. For example, an Hankel matrix is obtained if the embedding is arranged such that identical elements occur along the antidiagonals [10, 11]. Anyway, by computing twopoint time correlations between the entries of the multidimensional representation of the signal, it is possible to find an orthogonal matrix whose columns form an orthogonal basis of the dimensional space. The nonnormalized  dimensional correlation matrix is obtained as
and its eigenvalue decomposition
provides the related subspace model via the matrix of basis vectors and corresponding eigenvalues .
Furthermore, the subspace model allows to rewrite the original data matrix into a sum of rankone matrices according to
(2)  
where the denote proper weights, and each element of the th  dimensional matrix contains the dot product between the th eigenvector and the columns of the data matrix. These projections have the following properties:

They are uncorrelated, i. e.

Their energy is equal to :
These properties can be easily verified by considering an SVD of the original trajectory matrix according to . Here, the  dimensional matrix has also orthogonal columns as it corresponds to the right eigenvector matrix of . Therefore, the above discussed projections can be expressed as , i. e. with the right eigenvector scaled by the related singular value.
The most straightforward and indeed widely used technique to select the projections is to consider binary weights, i. e. , if the th projection is selected for the reconstruction, otherwise . But in noise reduction applications, for example, the selected relevant projections, corresponding to the largest eigenvalues, are scaled by appropriate normalized weights to reduce the noise contribution. This way, the eigenspectrum of is a truncated and rescaled version of the eigenspectrum of the original data. Therefore, assigning a real value like, for instance, where denotes the variance of the th eigenvector and represents the noise variance (or energy), yields for the weighted or rescaled projections:
(3) 
Note that only eigenvectors and related eigenvalues are retained onto which the data is projected. Rescaling thus effectively subtracts out the effect of noise contributions. Thereby, the noise variance of the data can be estimated as the mean of the discarded eigenvalues [19], but more sophisticated schemes have been considered also. Though other possibilities to compute the weights can be found, all assume that the noise variance should be smaller than the selected relevant eigenvalues [11, 15]. In the case discussed here this is assured as only refers to the eigenvalues which are kept, and eigenvalues are ordered with decreasing value.
Concerning the reconstructed data matrix as well as its rankone components , they generally do not exhibit identical elements along each descending diagonal like the original trajectory matrix . In SSA or SVD applications, these distinct entries along each diagonal (or antidiagonal) are replaced by their mean to reconstitute a Toeplitz or a Hankel matrix. Finally, a reconstructed univariate time series is obtained by reverting the embedding, i.e. by building the time series, i. e. the signal, from the mean values of each diagonal (or antidiagonal) of .
After having exposed the main properties of an SSA analysis, in the following we provide a comparative exposition of a classical time domain SSA analysis and an equivalent filter bank description to highlight the complementary information one can gain from both analysis techniques.
2.1 SSA components: time domain
Several works have addressed SSA  or SVD  based subspace projection techniques using linear filtering theory [12, 6, 2, 1, 25]. Concerning signal enhancement applications, in [12] the manipulations of the input signal were described as a filter bank approach, but all operations are justified using matrix manipulations as suggested by [6]. Examples of frequency responses of the eigenfilters, i. e. the eigenvectors, are shown graphically in [12], but no analytic expressions of the filter responses were given. In [1], the frequency profile of the eigenvectors was studied by way of the corresponding periodogram. In [25], an approach based on linear invariant system theory was presented to compute transfer functions of the discrete systems related with the projection step as well as the transfer functions related with the reconstruction and diagonal averaging steps. The coefficients of both filters correspond to the components of the related eigenvectors and define a causal and anticausal filter, respectively. The transfer function of the cascaded causal and anticausal filters, the branch of the filter bank, is given by
(4)  
where the filter coefficients are the entries of the corresponding eigenvector
Fig. 1 illustrates the block diagram of the system, where the scaling coefficient modifies the contribution of each signal component to the output .
The analysis filter is described by a polynomial with negative exponents (causal), while for the synthesis filter , the polynomial has positive exponents (anticausal) and the coefficients are divided by due to the diagonal averaging. The product of the two polynomials leads to the transfer function of an equivalent system which can be written as
(5) 
The coefficients of this system of cascaded analysissynthesis filters have special properties which can be derived easily from the SSA multidimensional approach. These coefficients are the result of the sum of the diagonals of , e.g, the  th element of is the sum of the entries along the  th diagonal of divided by . Here, the main diagonal corresponds to and positive (negative) values correspond to diagonals above (below) the main diagonal. Therefore, the properties result from the following considerations:

The matrix is symmetric, i. e. after adding and dividing each diagonal by ,
the coefficients of the filter obey to the relation
(6) 
The component matrices sum up as follows: , where is the identity matrix. Therefore, the sum of the impulse responses of the filters yields the unit impulse signal:
(7)
This last property of the coefficients of the combined analysis  synthesis filter cascade implies that by applying the filters with maximal weight when decomposing any input signal into inphase components, at the output the input signal is recovered. Thus the output signal , obtained by adding up filter outputs, is equal to the input signal . The symmetry of the coefficients leads to realvalued frequency responses which can be obtained by substituting in 5, which then reads
(8) 
Note that the realvalued frequency response with normalized sampling rate, i. e. period , corresponds to a zerophase filter. This means that each extracted component is always inphase with its related input signal . Note that frequency responses of the component filters and cannot be given in closedform similar to in eqn. 8 due to lacking symmetry properties of their coefficients [18].
In [13] the correlation matrix was forced to form a symmetric Toeplitz matrix which then has symmetric or antisymmetric eigenvectors. In that case, the frequency responses of the analysis and synthesis filters can be expressed as where the represents a realvalued function related with the magnitude of the frequency responses, and represents the phase [18]. However, the resulting also follows the same characteristics as described in eqn. 6. This is because the same procedure applies to its calculation, and its final expression corresponds to eqn. 8
Naturally, eqn. 5 can be similarly expressed in the time domain, leading to the implementation equation
(9) 
which corresponds to the time domain convolution of the noncausal filtering operation as each sample depends on future samples as it is defined in the second term of the eqn. 9.
2.2 SSA components: frequency domain
The power spectrum of a times series is defined as the Fourier Transform of its autocorrelation function
The autocorrelation function is an even function, therefore the power spectrum is a realvalued function with positive values. In an SSA approach, the entries of the matrix represent estimates of the nonnormalized autocorrelation function, computed for while using the considered segment of data encompassing samples. The diagonal of the matrix has the values , where represents the main diagonal and are upper diagonals starting from the main diagonal.
However notice that the in the different entries of the diagonal are numerically different, as the timedelayed correlations, estimated with subsegments of the data segment, are slightly different. This issue can be overcome if the correlation matrix is filled with the entries of the correlation function , instead of computing it with the embedding data matrix .
Note that estimating the power spectrum with samples of the autocorrelation function using a Discrete Fourier Transform (DFT), the resolution in frequency is the sampling rate divided by .
Furthermore, in linear invariant systems, the power spectrum of the output is related to the power spectrum of the input by the transfer function according to
2.2.1 Power Spectrum
Denoting by and the extreme values of the power spectrum of the input signal , it is possible to show that the outputs of the analysis filters have energies . Computing the energy at the output of the filter yields
Then as the filter has unit energy, it follows from Parseval’s theorem that
Hence, the upper bound is . Similarly it is possible to show that energy is also lower bounded by . Therefore, it can be concluded that the eigenvalues of the matrix must be also bounded by as those values represent the energy of the extracted components. Remember that by eqn. 5 the eigenvectors of the correlation matrix correspond to the filters. Then one can allocate the related eigenvalues to the frequency of the corresponding maximum of the frequency responses of the filters. In this way, the distribution of the eigenvalues matches the power spectrum of the input signal.
The synthesis filters show the same frequency profiles of the magnitude of the frequency response as the corresponding analysis filters. The frequency range (passband) of both, the analysis and the synthesis filter, is identical, but the total energy of the equivalent cascaded filter is different from the total energy of the individual filters. However, by computing the Fourier Transform of the equality eqn. 7, it follows that . Therefore, the power spectrum of is similar to the input if all components are summed up.
Illustrative Example
: A segment of the discrete signal
where represents Gaussian random noise (zero mean and unit variance), was analyzed with an SSA. Ideally, the power spectrum of the signal would consist of two pairs of delta Dirac functions at the frequencies of the sinusoids. However, due to the use of a finite number of samples, the analysis shows the socalled spectral leakage. Furthermore, the number of samples determines the frequency resolution of the estimated spectra. The numerical simulations were conducted for SSA with and and are compared with the estimation of the power spectrum using the DFT of the autocorrelation function with autocorrelation values. Fig. 2 illustrates the different results. On top, the frequency response of the filters are shown, and the maximal value of each filter is marked. On the bottom, the three plots correspond to the scree plot of the eigenvalues, and to the plots of the power spectrum estimates. Note that the scree plot shows four significant eigenvalues, whatever is the value of . But also notice that for a scree plot, the values on the abscissa have no meaning. The power spectrum estimates correspond to the DFT of the autocorrelation function. Similarly, a power spectrum estimate is obtained by plotting the eigenvalues in accord with the frequency of the maxima of the corresponding filters (eigenvectors). It has to be noticed that the frequency localization provided by the maximal values of the filters automatically identifies the frequencies of the related sinusoids. And note that, as the maxima of the filters are not equidistant, the localization can be more precise in cases when the frequency is not a multiple of the spectral resolution. The latter happens for the normalized frequency equal to as shown in the fig. 2.
3 Numerical Simulations
The goal of the two following examples is

to illustrate the application of SSA to gather relevant information from signals by

extracting narrowband components from a wideband biomedical (EEG) signal or

detecting the presence of a signal in communication (LTE) channels


to show that the frequency responses of the filters can be used to complement the eigenvalue spectrum in the conjugated Fourier domain.
3.1 Extraction of NarrowBand Signals
Typically in biomedical signal analysis, one acquires broadband signals but is interested only in some narrowband signal components. As an example, figure 3 shows the power spectrum of one segment of a single EEG channel recording, with a duration of and a sampling rate of . The power spectrum was estimated, employing Welch’s method, in subsegments of samples with of overlap. Figure 3 clearly exhibits, outside of the band of interest, i. e. , interferences representing higher harmonics of the power line. The same EEG segment was also analyzed with an SSA methodology, using an embedding dimension . Next, the eigenvalues are plotted semilogarithmically (see Fig. 3) versus frequency in the range . The eigenvalues are ordered in accord with the corresponding maximum of frequency response of the FIR filter. In this way, a spectral profile of the eigenvalues is generated which resembles the power spectrum of the input signal. Note that the signal is decomposed into components in order to assure a frequency resolution close to .
In practice, one never is interested in all components resulting from the decomposition, rather only the sum of few components within certain frequency ranges are of interest mostly. Fig. 4 shows an example of signal components representing characteristics EEG bands. Here the top row signal represents the original signal, and the bottom trace shows signal components related to nonneuronal interferences like electrooculograms, head/body movements and power line harmonics. The corresponding power spectra of the extracted signals are shown in figure 5. The signal components, resulting in traces , are obtained by grouping together synthesized components, i. e. components that have been obtained after the last stage of the SSA processing chain as shown in fig. 1, in accord with the maxima of the frequency response of the corresponding FIR filter .
Therefore the frequencies of the maxima of the filters are used to group the components into four clusters: three related with different characteristic EEG bands and the fourth is considered noise related. Then, the components forming one cluster are added up to result in a narrow band signal containing its frequency contents in one of the characteristic EEG bands. More specifically, referring to Fig. 4, six SSA components were combined to result in traces 2 and 3, each, while trace 4 encompasses SSA components and , finally, trace 5 represents a combination of the remaining SSA components. The related power spectra are show in Fig. 5.
3.2 Detection of radio frequency signals
The next example is related with a universal mobile telecommunication system (UMTS) signal, more specifically a signal from the Physical Uplink Control Channel (PUCCH) of the Long Term Evolution (LTE) uplink [21]. The data was acquired with an USRP B200 Software Defined Radio board ^{1}^{1}1https://www.ettus.com/product/details/UB200KIT, at a sampling rate of and with a central frequency of . Figure 6 illustrates the spectrogram of the acquired signal from DC to . It can be verified that the signals are present in short segments of time (). In this example, the LTE signal comprises three resource blocks (RB) of the LTE, each with bandwidth of which are assigned to PUCCH channels. The PUCCH is only occupied when a User Equipment (UE) needs to send control information.
For each segment of samples of the SSA model with an embedding dimension , an eigenvalue ratio is estimated for the subsegments of the signal shown in the spectrogram. The figure 7 shows the eigenvalue ratio, in semilogarithm scale, for all subsegments . It can be concluded that it is possible to define a threshold to decide if a signal is present or not, in each of the subsegments. However, note that with this detection it is not possible to know which of the bands is occupied. However if the eigenvalues are plotted according to the frequency of the maxima of the frequency response of the corresponding filters the frequency localization is possible. Figure 8 shows that in the second segment () the low power signal occupies the highest frequency band (). The figure also shows the eigenvalue spectrum when there is no signal (). The next two plots concern the occupancy of the lowpass band and the bandpass, respectively.
4 Conclusions
The interpretation of singular spectrum analysis (SSA) as a bank of filters can be useful to attain a more clearcut insight into the outcomes of the method. By the frequency responses of the filter bank, corresponding to the eigenvectors of the correlation matrix, the frequency content of the different component signals can be easily attained. SSA filters are dataadaptive, and the relevance of one component to the energy of the input signal is deduced from the corresponding eigenvalue. In this study, the relation between the eigenvalues and the power spectrum was addressed. It was shown that it is advantageous to reorder the eigenvalues in accord with the occurrence of the maxima of the frequency responses of the related filters rather than simply arranging them by decreasing value as is common practice. In summary, exploiting the information from both analysis strategies, i. e. the eigenvalue spectrum and the related frequency response of the filter, provides additional insight compared to a canonical analysis.
Two typical applications illustrated the advantage of this analysis technique. Analyzing a typical biomedical signal, more specifically an EEG, often amounts to extracting narrowband signal components from a broadband acquired signal. It is known, and corroborated in fig. 3 that the power spectra of most biomedical signals show a spectral distribution [20], which naturally results in a standard ordering of the eigenvalues by decreasing value. But the peaks, belonging to the higher order harmonics of the line noise spectral contribution, can only be seen if the eigenvalues are properly matched to the maxima of the frequency responses of the filters. Thus, to assure that the extracted components fall in a chosen narrowband range, namely the classical EEG frequency bands, the frequency profiles of the related filters have to be studied to identify those signal components which need to be combined to result in the desired narrowband signals. A second prototypical application is channel identification in spectrum sensing. The goal here was to identify at times unoccupied transmission channels which then can be used for additional signal transmission. These simple applications illustrate the usefulness of a combined analysis of the eigenspectrum of the correlation matrix in the time domain as well as the related frequency response. It was shown that reordering the eigenvalues in accord with the maxima of the frequency responses of the filters allows to identify the most relevant signal components for the problem at hand.
Acknowledgements
This work was supported by FCT: grant to A. R. Teixeira (Ref:SFRH/BPD/101112/2014) and the project EXCL/EEITEL/0067/2012 (CREaTION).
References
 Alexandrov et al. [2008] T. Alexandrov, N. Golyandina, and A. Spirov. Singular spectrum analysis of gene expression profiles of early drosophila embryo: exponentialindistance patterns. Res. Let. Signal Proc., 2008:1–5, 2008.
 Alexandrov [2009] Theodore Alexandrov. A method of trend extraction using singular spectrum analysis. REVSTAT, 7(1):1–22, 2009.
 Alharbi and Hassani [2016] N. Alharbi and H. Hassani. A new approach for selecting the number of the eigenvalues in singular spectrum analysis. Journal of the Franklin Institute, 353(1):1 – 16, 2016. ISSN 00160032. doi: 10.1016/j.jfranklin.2015.10.015.
 Axell et al. [2012] E. Axell, G. Leus, E. G. Larsson, and H. V. Poor. Spectrum sensing for cognitive radio: Stateoftheart and recent advances. IEEE Signal Processing Magazine, 29(3):101–116, May 2012. ISSN 10535888. doi: 10.1109/MSP.2012.2183771.
 De Sanctis et al. [2011] Silvia De Sanctis, Wilhelm M. Malloni, Werner Kremer, Ana M. Tomé, Elmar W. Lang, KlausPeter Neidig, and Hans Robert Kalbitzer. Singular spectrum analysis for an automated solvent artifact removal and baseline correction of 1 D NMR spectra. Journal of Magnetic Resonance, 210(2):177–183, 2011.
 Dologlou and Carayannis [1991] I. Dologlou and G. Carayannis. Physical interpretation of signal reconstruction from reduced rank matrices. IEEE Transactions on Signal Processing, 39(7):1681–1682, 1991.
 Ephraim and Trees [1995] Yariv Ephraim and Harry L. Van Trees. A signal subspace approach for speech enhancement. IEEE Transactions on Acoustic, Speech and Signal Processing,, 3(4):251–266, 1995.
 Ghaderi et al. [2011] F. Ghaderi, H. R. Mohseni, and S. Sanei. Localizing heart sounds in respiratory signals using singular spectrum analysis. IEEE Transactions on Biomedical Engineering, 58(12):3360–3367, Dec 2011. ISSN 00189294. doi: 10.1109/TBME.2011.2162728.
 Ghil et al. [2002] M. Ghil, M.R. Allen, M. D. Dettinger, K. Ide, and et al. Advanced spectral methods for climatic time series. Reviews of Geophysics, 40(1):3.1–3.41, 2002.
 Golyandina et al. [2001] Nina Golyandina, Vladimir Nekrutkin, and Anatoly Zhigljavsky. Analysis of Time Series Structure: SSA and Related Techniques. Chapman & HALL/CRC, 2001.
 Hansen and Jensen [2007] Per Christian Hansen and S. Holdt Jensen. Subspacebased noise reduction for speech signals via diagonal and triangular matrix decompositions: Survey and analysis. Eurasip Journal on Advances in Signal Processing, Vol 2007(1):092953, 2007. ISSN 16876180. doi: 10.1155/2007/92953.
 Hansen and Jensen [1998] Per Christian Hansen and Soren Holdt Jensen. Fir filter representations of reducedrank noise reduction. IEEE Transactions on Signal Processing, 46(6):1737–1741, 1998.
 Harris and Yuan [2010] T.J. Harris and Hui Yuan. Filtering and frequency interpretations of singular spectrum analysis. Physica D: Nonlinear Phenomena, 239(20 – 22):1958 – 1967, 2010. ISSN 01672789. doi: 10.1016/j.physd.2010.07.005.
 Hassanpour et al. [2012] Hamid Hassanpour, Amin Zehtabian, and S.J. Sadati. Time domain signal enhancement based on an optimized singular vector denoising algorithm. Digital Signal Processing, 22(5):786 – 794, 2012. ISSN 10512004. doi: 10.1016/j.dsp.2012.03.009.
 Hermus et al. [2006] Kris Hermus, Patrick Wambacq, and Hugo Van hamme. A review of signal subspace speech enhancement and its application to noise robust speech recognition. Eurasip Journal on Advances in Signal Processing, 2007(1), 2006. ISSN 16876180. doi: 10.1155/2007/45821.
 Huang et al. [2014] Lei Huang, H.C. So, and Cheng Qian. Volumebased method for spectrum sensing. Digital Signal Processing, 28:48 – 56, 2014. ISSN 10512004. doi: http://dx.doi.org/10.1016/j.dsp.2014.02.003.
 Jackson [1991] Leland B. Jackson. Signals, Systems and Transforms. AddisonWesley, 1991.
 Jackson [1996] Leland B. Jackson. Digital Filters and Signal Processing. Kluwer Academic Publishers, 1996.
 Liavas and Regalia [2001] Athanasios P. Liavas and Philip A. Regalia. On the behavior of information theoretic criteria for model order selection. IEEE Transactions on Signal Processing, 49(8):1689–1695, 2001.
 M. et al. [2012] Li M., Zhao W, and Chen B. Heavytailed prediction error: A difficulty in predicting biomedical signals of 1/f noise type. Computational and Mathematical Methods in Medicine, page 5, 2012. doi: 10.1155/2012/291510.
 Perez [2015] Andre Perez. LTE and LTE Advanced: 4G Network Radio Interface. WileyISTE, 2015.
 Sanei et al. [2011] Saeid Sanei, Mansoureh Ghodsi, and Hossein Hassani. An adaptive singular spectrum analysis approach to murmur detection from heart sounds. Medical Engineering and Physics, 33(3):362 – 367, 2011. ISSN 13504533. doi: http://dx.doi.org/10.1016/j.medengphy.2010.11.004.
 Santos et al. [2016] I.M. Santos, A.R. Teixeira, A.M. Tomé, A.T. Pereira, P. Rodrigues, P. Vagos, J. Costa, M.L. Carrito, B. Oliveira, N.A. DeFilippis, and C.F. Silva. ERP correlates of error processing during performance on the Halstead Category Test. International Journal of Psychophysiology, 106:97 – 105, 2016. doi: /10.1016/j.ijpsycho.2016.06.010.
 Teixeira et al. [2006] A. R. Teixeira, A. M. Tomé, E.W. Lang, P. Gruber, and A. Martins da Silva. Automatic removal of highamplitude artifacts from singlechannnel electroencephalograms. Computer Methods and Programs in Biomedicine, 83(2):125–138, 2006.
 Tomé et al. [2010] A. M. Tomé, A. R. Teixeira, N. Figueiredo, I. M. Santos, P. Georgieva, and E. W. Lang. SSA of biomedical signals: A linear invariant systems approach. Statistics and Its Interface, 3(3):345–355, 2010.
 Z. Idrees et al. [2015] Z. Z. Idrees, F. A. Bhatti, and A. Rashdi. Spectrum sensing using lowcomplexity principal components for cognitive radios. EURASIP Journal on Wireless Communications and Networking, 2015(1):1–8, 2015. ISSN 16871499. doi: 10.1186/s1363801504124.
 Zeng and Liang [2009] Y. Zeng and Y. C. Liang. Eigenvaluebased spectrum sensing algorithms for cognitive radio. IEEE Transactions on Communications, 57(6):1784–1793, June 2009. ISSN 00906778. doi: 10.1109/TCOMM.2009.06.070402.
Appendix
SSA MATLAB implementation and two of the examples of the report.
a=[2 4] f=[0.1 0.4] x=zeros(1,1024) n=0:length(x)1; for k=1:2 x=x+a(k)*sin(2*pi*f(k)*n); end x=x+randn(1,length(x)); M=9 %Correlation matrix symmetric model=model=SSAglobal(x,M,’toeplitz’,’eigenvectors’) %Correlation matrix nonsymmetric model1=SSAglobal(x,M,’embedding’,’eigenvectors’) % Grouping components according to model.peaks % EEG sampled 1000Hz,embedding dimension 201 %model= SSAglobal(s,200,’toeplitz’,’eigenvectors’,1000) F=[2 8;8 15;15 32; 32 100] L=size(F,1) % Components group=zeros(L+1,9000) contador=zeros(1,L+1); for c=1:L a= (model.peaks>= F(c,1) & model.peaks <F(c,2))’; switch sum(a) case 0 % no peak group(c,:)=zeros(1,M) case 1 group(c,:)=model.comp(a,:); otherwise group(c,:)=sum(model.comp(a,:)); end contador(c)=sum(a); end %the a=(model.peaks< F(1,1)  model.peaks>F(L,2))’; contador(L+1)=sum(a) group(L+1,:)=sum(model.comp(a,:)) %%% %%%%%% function [model]=SSAglobal(x,M,TYPE,ORDERING,Fs) % model=SSAglobal(x,M,TYPE,ORDERING,Fs) % model=SSAglobal(x,M) % Default: TYPE=’toeplitz’,ORDERING=’eigenvalues’,Fs=1 (normalized % frequencyradians) %model=SSAglobal(x,M,’toeplitz’,’eigenvectors’) %OUTPUT: % model.eigval : eigenvalues (ordering by decreasing order or according to % the peaks of TF) % model.all: filter bank (analysis and synthesis) % model.TF: frequency responses of the analysissynthesis % model.peaks: frequency of the peaks (if Fs=1, 0<values<0.5) % model.f: f=(0:N1)*(Fs/N); % model.comp : components outputs of all filters of model.all % AMT October 2015 N=4000; %N if nargin <2  nargin>5 error (’invalid number of arguments’); else switch nargin case 2 TYPE=’toeplitz’; ORDERING=’eigenvalues’; Fs=1; case 3 ORDERING=’eigenvalues’; Fs=1; case 4 Fs=1; end end if strcmp(TYPE,’toeplitz’) % %compute autocorrelation function (from delay M up to M) rx=xcorr(double(x),M); R=toeplitz(rx(M+1:end)); elseif strcmp(TYPE,’embedding’) fim=0; for k=1:M+1 X(k,:)=x(M+1k+1:endfim); fim=fim+1; end R=X*X’; else error(’valid options: embedding or toeplitz’); end % order of eigenvalues by decreasing [U,S,dummy]=svd(R); % % % % filter bank for k=1:M+1 u=U(:,k); %normalization t(k,:)=conv(u,u(end:1:1))./(M+1); end % % %%%ORDERING BY eigenvalues (default) % %filter bank and eigenvalues by dcreasing order model.all=t; model.eigval=diag(S); model.U=U; % f=(0:N1)*(Fs/N); Tf=abs(fft(t,N,2)); % %Only half... 0 to Fs/2 Tf(:,N/2+2:end)=[]; model.TF=Tf; model.peaks=[]; model.f=f(1:size(Tf,2)); if strcmp(ORDERING,’eigenvectors’) [dummy,inx]=max(Tf,[],2); [dummy,inx]=sort(f(inx)); model.eigval=model.eigval(inx); model.U=U(:,inx); model.all=model.all(inx,:); model.TF=Tf(inx,:); model.peaks=dummy; end %%The components for k=1:M+1 y=conv(model.all(k,:),x); model.comp(k,:)=y(M+1:endM); end