1 Introduction
Singular spectrum analysis (SSA) is a nonparametric and adaptive decomposition of a time series. In the SSA, the original time series is decomposed into arbitrary number of additive components and they are interpreted as the slowly varying trend, oscillatory and noise components, respectively [2,3]. In the conventional approach, the original time series is embedded into the trajectory matrix
through the lagged vectors of length
by forminglagged vectors. Then the singular value decomposition of the trajectory matrix
is carried out as(1) 
with are the rank 1 matrices. After grouping the matrices , the time series is reconstructed by projecting the matrices into Hankel form through the antidiagonal averaging. Recently, it has been shown that these SSA procedure can equivalently be interpreted as the twostep (forward and the reverse) filtering of a time series [6,8]. The eigenvector of the matrix
can be interpreted as the convolution coefficients of the noncausal FIR filter. Let us denote the discrete Fourier transform of the eigenvector
for as . In Fourier space, the first (second) step filtering can be expressed as the multiplication of () to the original sequences. Then the twostep filtering can be expressed as the multiplication of the zerophase filters [7]. The normalization and the perfectreconstruction properties can be expresses as(2) 
and
(3) 
and these are the consequences of the normalization and the completeness properties of the eigenvectors for [8].
Though the filtering interpretation of the SSA algorithm is equivalent to the conventional treatment, it has several advantages. In the conventional approach, the singular value decomposition of the trajectory matrix and the reconstruction procedure through the Hankel projection seems to be quite different procedures and this makes it difficult to extend SSA to higher dimension. With the filtering interpretation, these procedures can be seen as the forward and the reverse twostep pointsymmetric convolution with respect to the reference points. This largely simplify the algorithms and enables us to extend the SSA to the decomposition of multidimensional data with arbitrary dimension [9]. Also, we can study the properties of the filters separately and it helps us for the detailed study of the spectral meaning of the decomposition.
There have been several works in which the SSA algorithm has been applied to 2D image data [4,12] or three dimensional polygonal data [10,11]. In this paper, we focus on the SSA decomposition of the 2D image data. Under the periodic boundary condition for the image data, the matrix
becomes to be the bisymmetric type and it leads to either the symmetric or the antisymmetric eigenvectors. A single step filtering generated with these eigenvectors can be interpreted as the differential operation to the images and the symmetric (antisymmetric) eigenvalues lead to the derivative filters of the even ( odd ) order. In the image processing, various differential filters are designed to detect and enhance the edge lines in the images. In the SSA decomposition of the 2D image data, these differential filters are adaptively and optimally generated. We have examined these points in detail, focusing our attention to the characteristics of the filters in the SSA. Finally, we briefly discuss the possibility of the denoising with the SSA decomposition.
In Sec. 2, we give a brief summary of the SSA algorithm in 2D image decomposition, emphasizing the role of filters. In Sec. 3, we give an example of the SSA decomposition of the image data, and the symmetry of the eigenvector leads to the symmetry properties of the filters and these are reflected as the filter characteristics for the image data. An implication of these properties to the image denoising is briefly discussed. In Sec. 4, we summarize the results of this paper.
2 SSA Decomposition of 2D Image Data
2.1 Basic algorithm
Let us consider the 2D image data given by the matrix with and . We adopt the moving window of any shape. For simplicity we assume here the rectangular form , which is shown as the box in the following equation
(4) 
Then, by moving it from left to right and from top to bottom, for example, the trajectory matrix can be defined as
(5) 
To introduce the filter for the 2D data, we consider an arbitrary
dimensional orthonormal vector
. We arbitrarily assign the components of in the moving window,(6) 
In the SSA, we adopt the vector set which is the eigenvectors of the matrix and the filters are defined as,
(7) 
The following convolutiontype linear relation express the forward filtering
(8)  
Next, we perform the reverse filtering with the point symmetric filter
(9) 
for the reference point with the relation
(10)  
We assume the periodic boundary for the image data . In the SSA, the vectors are determined as the eigenvector of the lagcovariance matrix . It is easy to show that the completeness of the orthonormal eigenvectors leads to the perfectreproduction property of the twostep filtering irrespective of the window shape,
(11) 
The properties of the filters are clearly seen in the Fourier space. By embedding the eigenvector in the matrix
(12) 
we define the discrete Fourier transform as
(13) 
and
(14) 
Then they satisfy the normalization
(15) 
and the completeness relations (perfect reproduction in the language of digital filtering)
(16) 
The above properties Eqs. (15) and (16) are satisfied for arbitrary orthonormal vectors . The SSA algorithm employs the eigenvectors of the eigenvalue equation for
(17) 
and this criterion means that the quantities
(18) 
are maximized preserving the orthonormality of the vector
[9]. In this sense, the SSA algorithm is the principal component analysis in the space of the lagged vectors. In the next subsection, we discuss the symmetry properties of the eigenvectors
and it’s implication to the 2D image decomposition.2.2 Symmetry properties of the filters
We assign the vector components for the rectangularshaped moving window with in a natural order as in Eq. (6). We assume the periodic boundary condition for the image data. Then the lagcovariance matrix becomes to be the bisymmetric matrix, satisfying the relation
(19) 
This can be seen as follows. The moving window in Eq.(6) slides from left to right and up to bottom sweeping the whole image data . The component of the matrix is calculated by picking and summing up the matrix elements of , corresponding to the position and in the window, and this is the same for with and . The location of the components and is just point symmetric about the center of the window. This also is the same for the relation between and . Thus the relative location between and is the same as that between and . Because of this, the matrix element is equal to . This property holds for higher dimensional cases if we adopt the periodic boundary condition. For one dimensional case, a singlechannel time series for instance, the matrix is a symmetric Toeplitz matrix.
Thus the matrix commutes with the exchange matrix
(20) 
and if there is no degeneracy for the eigenvalues of , the eigenvectors are either symmetric or antisymmetric
This means that the filters generated with the eigenvectors have corresponding symmetry. In particular, for the case of square form the matrices expressing the filter
(21) 
are the centrosymmetric or skewcentrosymmetric matrices corresponding to the symmetric
or antisymmetric eigenvectors, respectively.As an example, we have carried out the SSA decomposition of the image data ’Building’ with rather small filter of square form . The nine eigenvalues of the matrix are ordered from larger to smaller ones and are shown in Fig. 1. The five eigenvectors corresponding to the larger eigenvalues are shown in Table 1. Here, the components of the eigenvectors are ordered in the filter as in Eq.(21). As seen, the eigenvectors , and are symmetric, while , are antisymmetric.
0.3327  0.3747  0.4151  0.2671  0.0899 
0.3338  0.0197  0.4336  0.4223  0.2913 
0.3326  0.4168  0.3732  0.2390  0.3057 
0.3336  0.4307  0.0198  0.1967  0.4716 
0.3348  0.0000  0.0000  0.5558  0.4272 
0.3336  0.4307  0.0198  0.1967  0.4716 
0.3326  0.4168  0.3732  0.2390  0.3057 
0.3338  0.0197  0.4336  0.4223  0.2913 
0.3327  0.3747  0.4151  0.2671  0.0899 
To see the details of the decomposition characteristics, we have analyzed the filters as the differential operation. Let us assume that the image pixels are located on the lattice points with mesh length and for both and directions, respectively. Here, the axis is taken for the horizontal direction from left to right and the axis for the vertical direction from up to bottom. The pixel value can be expressed with the function on the lattice points, and the operation of the filter can be expressed as
(22)  
By carrying out the Taylor expansion we obtain
(23)  
where the vector components are abbreviated as and the etc. The coefficients for each terms are shown in Table 2 for dominant five filters. The filters generated from the symmetric (antisymmetric) eigenvectors can be seen as the derivative filters of the even (odd) order which is simply a consequence of the centrosymmetric (skewcentrosymmetric) nature of the filters .
Filter  

1  3.000  0.000  0.000  0.999  0.000  0.999 
2  0.000  0.123  2.444  0.000  0.000  0.000 
3  0.000  2.444  0.123  0.000  0.000  0.000 
4  0.005  0.000  0.000  0.084  0.056  0.703 
5  0.004  0.000  0.000  0.687  0.432  0.076 
The dominant filter works as smoother of the pixel values. The second and the third filters and are generated from the antisymmetric eigenvectors and work as the differential filter in and directions, respectively. The filters and are generated from the symmetric eigenvectors and the dominant components are the second derivative for various directions.
The image decomposition for the ’Building’ are shown in Fig. 2. We can see that the results of the decomposition reflect the characteristics of each filter. For the decompositions (2) and (4), the horizontal lines are enhanced corresponding to the large , components as in Table 2. On the other hand, the vertical components are enhanced in the image (3) and (5). As seen, each filter has it characteristics and the SSA decompositions reflect their characteristic properties. The images correspond to smaller eigenvalues and they exhibit somewhat featureless pattern due to the large highfrequency components. These nine images are added to reproduce the original one in the SSA algorithm.
2.3 SSA decomposition of images with noise
To see how the SSA algorithm works for the decomposition of the image data with noise, we adopt an example ’Lenna’ with pixels. Each pixel takes the values . The Gaussian noise is added for every pixels. We have carried out SSA decomposition with a square filter of the size . The larger fifteen eigenvalues of the matrix are shown in Fig. 3, where the black circles (white triangles) represent the eigenvalues corresponding to the symmetric (antisymmetric) eigenvectors.
Similar to the case of the ’Building’ image, we have examined how each filter works for the image by examining the coefficients in Eq. (23) extended for the case of filter. In this case, since the filter size is large, we have replaced the mesh lengths and by and , respectively as,
(24) 
The coefficients , , etc. are shown in Table 3.
Filter  

1  11.000  0.000  0.000  6.079  0.021  6.101 
2  0.000  11.337  1.672  0.000  0.000  0.000 
3  0.000  1.669  11.133  0.000  0.000  0.000 
4  0.067  0.000  0.000  5.092  2.654  0.555 
5  0.005  0.000  0.000  1.245  11.140  0.667 
For the ’Lenna’ image with noise, the filter characteristics are quite similar to the case of the ’Building’ photo image. The first filter works as the smoother for the original image. The smoothing operation works stronger than the case of ’Building’, since the filter size is larger in this case. The second and the third filters work as the differentiation to the and directions, respectively. For the fourth and the fifth filters and , the expansion coefficients and are large as seen in Table 3.
For the filters corresponding to the smaller eigenvalues, the derivative components are fragmented and the filtered images gradually become to be structureless including the the higher frequency components.
In Fig. 4, we have shown the ’Lenna’ image with and without noise together with the dominant five SSA decomposition images with noise. Corresponding to the filter characteristics as seen in the Table 3, we can see that the edge lines in various directions are enhanced and are extracted from the original image. For the images (2) and (3), the vertical and the horizontal lines are enhanced. The filters and have large and components, and hence the finer details of the edges in the vertical or the diagonal lines are enhanced for the images (4) and (5).
In order to see the relation between the SSA decomposition and the noise components, we show all the 121 eigenvalues for the ’Lenna’ image with and without the Gaussian noise in Fig.5. The added noise enhances the smaller eigenvalues coming from the highfrequency components.
We have also calculated the RMS distance between the original noiseless image and the partial sum of the SSA decomposition [RodriguezAragon],
(25) 
The results are shown in Fig. 6. Obviously, for the noiseless case, the distance converges to zero as , which merely means that the sum of the SSA decomposition reproduce the original image. For the case with noise term, the RMS distance decreases as , while it increases as . The RMS distance takes minimum value at . We could discard the components in order to obtain better images. The image corresponding to the partial sum is also shown in Fig. 4. The RMS distance between noiseless and noisy images is
which is just the standard deviation of the added noise. The RMS distances for the partial sum of the SSA decomposition
with and without noise are and , respectively.So far, much efforts have been devoted to the image restoration. Various linear or nonlinear filtering operations have been examined to improve the noisy images [1,5]. The performance of the filtering operation depends on the properties of the image itself or the nature of the noise. In the SSA, the noisy images are exactly decomposed into arbitrary number of components and each image decomposition has its own characteristics. Thus, by combining the conventional nonlinear denoising algorithms and SSA decomposition, we expect to develop better denoising algorithms. Since the main purpose of this paper is to show the characteristic feature of the SSA decomposition of the 2D image data, we do not pursue these points further.
3 Summary and Conclusions
The filtering interpretation of the SSA algorithm enabled us to study the characteristics of the image data decomposition in detail. The periodic boundary condition and the rectangular moving window leads to the bisymmetric lagcovariance matrix and its eigenvectors are either symmetric or antisymmetric. For the symmetric (antisymmetric) eigenvectors, the corresponding filter are centrosymmetric (skewcentrosymmetric) matrices and they are interpreted as the derivative filter of even (odd) order. The SSA algorithm can be understood as the optimal and adaptive generation of the filters
and the twostep pointsymmetric operation of these filters to the reference points of the original image. The twostep filtering ensures the reality of the filters in the Fourier space and it introduces no distortion or shift to the data. Also the completeness of the eigenvectors ensures the perfect reproductivity of the SSA decomposition. Each image components in the SSA decomposition exhibit characteristic features, which comes from the properties of the adaptive filters generated from the eigenvectors of the lagcovariance matrix. We briefly discussed the relation of the SSA decomposition and the noise reduction of the images. Much efforts have been devoted so far for the image denoising. We expect that the filtering approach for SSA algorithm enables us to develop an novel denoising algorithm by combining the SSA decomposition with conventional nonlinear filtering algorithms. The study in this direction is now under way.
Acknowledgment
This work is supported by GrantinAid for Scientific Research from the Japan Society for the Promotion of Science; Grant Numbers (24650150) and (26330237).
References

[1]
Buades, A., Coll, B. and Morel, JM. (2008). Nonlocal image and movie denoising.
International Journal of Computer Vision
,762:123–139, and references therein.  [2]
 [3] Elsner, J. B. and Tsonis, A. A. (1996). Singular Spectrum Analysis. A new tool in time series analysis. Plenum Press.
 [4]
 [5] Golyandina, N., Nekrutkin, V. and Zhigljavsky, A. (2001). Analysis of time series structure: SSA and related techniques. Chapman and Hall/CRC.
 [6]
 [7] Golyandina, N. and Usevich K. (2010). 2Dextension of singular spectrum analysis: algorithm and elements of theory. Matrix methods: Theory, algorithms and applications (Eds. Olshevsky, V and Tyrtyshnikov, E.) ; World Scientific Publishing, 2010, 449–473.
 [8]
 [9] Gupta, K and Gupta, S. K. (2013). Image Denoising Techniques A Review Paper. International Journal of Innovative Technology and Exploring Engineering,24:6–9., and references therein.
 [10]
 [11] Harris, T. J. and Yuan, H. (2010). Filtering and frequency interpretations of singular spectrum analysis. Physica D, 239: 1958–1967, and references therein.
 [12]
 [13] Kormylo, J. and Jain, V. (1974). Twopass recursive digital filter with zero phase shift, IEEE Trans. Acoustics, Speech and Signal Process., 22: 384–387.
 [14]
 [15] Kume, K. (2012). Interpretation of singular spectrum analysis as complete eigenfilter decomposition. Advances in Adaptive Data Analysis, 44: 1250023 1–18.
 [16]
 [17] Kume, K. and NoseTogawa, N. (2014). Multidimensional extension of singular spectrum analysis based on filtering interpretation. Advances in Adaptive Data Analysis, 61: 1450005 1–17.
 [18]
 [19] Murotani, K. and Sugihara, K. (2005). New spectral decomposition method for threedimensional shape models and its application, Journal of computing and information science in engineering, 5: 277–282.
 [20]
 [21] Murotani, K. and Yagawa, G.(2008). Noise filtering of images using generalized singular spectrum analysis, The 16th International Conference in Central Europe on Computer Graphics, Visualization and Computer Vision 2008. University of West Bohemia, Plzen, Czech Republic, 47–54.
 [22]
 [23] RodriguezAragon, Licesio J. and Zhigljavsky, A. (2010). Singular spectrum analysis for image processing. Statistics and Its Interface, 3: 419–426.
Comments
There are no comments yet.