{'fraction_non_alphanumeric': 0.05141688577271399, 'fraction_numerical': 0.014052001168565587, 'mean_word_length': 4.26711801815664, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stationarity of Time-Series on Graph via Bivariate Translation Invariance

Stationarity is a cornerstone in classical signal processing (CSP) for modeling and characterizing various stochastic signals for the ensuing analysis. However, in many complex real world scenarios, where the stochastic process lies over an irregular graph structure, CSP discards the underlying structure in analyzing such structured data. Then it is essential to establish a new framework to analyze the high-dimensional graph structured stochastic signals by taking the underlying structure into account. To this end, looking through the lens of operator theory, we first propose a new bivariate isometric joint translation operator (JTO) consistent with the structural characteristic of translation operators in other signal domains. Moreover, we characterize time-vertex filtering based on the proposed JTO. Thereupon, we put forth a new definition of joint wide-sense stationary (JWSS) signals in time-vertex domain using the proposed isometric JTO together with its spectral characterization. Then a new joint power spectral density (JPSD) estimator, called generalized Welch method (GWM), is presented. Simulation results are provided to show the efficacy of this JPSD estimator. Furthermore, to show the usefulness of JWSS modeling, we focus on the classification of time-series on graph. To that end, by modeling the brain Electroencephalography (EEG) signals as JWSS processes, we use JPSD as the feature for the Emotion and Alzheimer's disease (AD) recognition. Experimental results demonstrate that JPSD yields superior Emotion and AD recognition accuracy in comparison with the classical power spectral density (PSD) and graph PSD (GPSD) as the feature set for both applications. Eventually, we provide some concluding remarks.


Introduction
Beyond doubt, we are in the era of big data in which massive amount of information is generated at a fast pace and this poses new challenges for the data science. Often, the big structured data lies over an irregular structure, but the classical signal processing (CSP) disregards the underlying topological structure. Graphs, as the powerful mathematical tool, have been widely used in biology, physics, and computer science to model the underlying complex structure of non-Euclidean data through pairwise relations between objects. Connecting the concepts from algebraic graph theory to the CSP gave birth to the field of graph signal processing (GSP) [2][3][4][5][6][7] as a theoretical discipline for analyzing the structured data. In the recent years, GSP emerged via numerous theoretical research works for tackling challenging problems in modern signal processing and data science [8][9][10][11][12][13][14][15]. In particular, [15] proposed a new framework for analyzing time-varying graph signals through a meaningful representations of timeseries on graph. Stationarity and its important subclass wide-sense stationarity play an essential role in statistical signal processing and time-series analysis. By the classical definition, a signal is temporal wide-sense stationary (TWSS) if its mean and autocorrelation functions are translation invariant. Likewise, this concept is of paramount importance concerning time-series on irregular graph structures.
Stationarity on graph is first defined by Girault et al. [16,17] via isometric graph translation operator. Afterwards, Perraudin and Vandergheynst [18] defined stationarity on graph via a novel localization operator. Moreover, Marques et al. [19] defined weak stationarity of random graph signals using the adjacency matrix (or graph Laplacian matrix) as the graph shift operator. These approaches lead to almost the same definition of stationarity on graph in spite of their different initial ideas [20]. Then, Segarra et al. [21] defined stationarity in time-vertex domain based on the weighted adjacency matrix of joint/product graph. By generalizing a classical "filtering interpretation" of stationarity from Euclidean space, Loukas and Perraudin [22] defined joint wide-sense stationary (JWSS) signals via time-vertex filtering. Noticeably, they have shown that stationarity in the time-vertex domain is more general than defining the stationarity on joint/product graph. Using this general notion of joint wide-sense stationarity, Isufi et al. [23] extended the classical vector autoregressive and vector autoregressive moving average recursions for modeling and forecasting time-series on graph.
In this paper, a novel approach, beyond [21][22][23], is proposed for characterizing the stationarity of time-series on graph via bivariate translation invariance in the time-vertex domain. First, we design a bivariate isometric joint translation operator (JTO) in the time-vertex domain. Then we put forth a new definition of wide-sense stationarity of time-series on graph based on the proposed JTO followed by its spectral characterization. Next, the joint power spectral density (JPSD) estimation of JWSS processes is proposed using a generalized Welch method (GWM) followed by some simulations to show its effectiveness. Finally, we demonstrate the applicability of proposed framework for the classification of time-series on irregular graphs by providing experimental results concerning Emotion and Alzheimer's disease (AD) recognition from Electroencephalography (EEG) data.
Notations. Matrices and vectors are denoted by uppercase and lowercase boldface letters, A and a, respectively. The n-th element of a vector is indexed by a[n], and the entry in n-th row and m-th column of a matrix is denoted by A[n, m]. R = [R (m,n) ] is a block matrix where R (m,n) is its submatrix in the m-th row and n-th column partition. Other notations are as follows: A T , A, and A * = (A) T stand for the transpose, conjugate, and complex conjugate transpose of the matrix A, respectively. Moreover, vec(A) stands for the column vector by stacking all the columns of A sequentially, Diag(a) represents a diagonal matrix by placing the elements of vector a on the main diagonal, and Diag(A) is equivalent to Diag(vec(A)), while diag(A) represents the column vector containing the diagonal elements of matrix A. We use a p and |a| as the p-norm of a and absolute value of a, respectively. On the other hand, |A| is a matrix with (n, m)-th element equal to |A[n, m]|. row k (A) stands for the k-th row of matrix A. Also, I, 0, and 1 denote the identity matrix, matrix/column vector of all zeros, and matrix/column vector of all ones (their dimensions may be indicated by their subscript for some emphasized cases), respectively. Symbols ⊗ and represent the Kronecker product and Hadamard (element-wise) product, respectively. Besides, i = √ −1 and a, b represents the set of integers between a and b inclusive. The Kronecker sum, denoted by ⊕, is defined as: On the other hand, the direct sum of a set of matrices is the set of N × 1 complex (real) vectors. R + (Z + ) accounts for the set of nonnegative real (integer) numbers. Note that E[·] accounts for the statistical expectation. Moreover, C, D, G, and J represent the continuous-time, discrete-time, graph, and joint time-vertex domains, respectively. Note 1. Let L be an N × N Hermitian matrix with eigenvalue decomposition L = ΦΛΦ * where Φ is a unitary eigenbasis matrix ΦΦ * = Φ * Φ = I and Λ = Diag ([λ 0 , . . . , λ N −1 ]) is the diagonal eigenvalue matrix. Then, the univariate matrix function is defined as h(L) [24]. Note that h can also be a multivariate matrix function.

Background
Vertex Harmonic Analysis. Let G := (V, E, W ) denote a fixed graph with finite vertex set V with the cardinality |V| = N , E = {(n 1 , n 2 )| n 1 , n 2 ∈ V, n 1 ∼ n 2 } ⊆ V × V is the edge set and W : V × V → R + is a weight function. This function yields the weighted adjacency matrix as W G = [w n 1 ,n 2 ] ∈ R N ×N + where w n 1 ,n 2 represents the strength of the connection between nodes n 1 and n 2 . The neighborhood of vertex v n is the set of vertices connected to v n by an edge. Throughout this paper, we assume that the graph is finite, weighted, connected, and undirected. A graph signal, represented in a compact form by the vector f ∈ C N , is defined by the function f G : V → C where f [n] is the function value at the vertex n. Then the graph Laplacian matrix is defined as L G := Diag (W G 1) − W G . Since L G is a symmetric positive semi-definite matrix, the eigenvalue decomposition (EVD) of L G can be written as where Φ G = [ϕ G,0 , . . . , ϕ G,N −1 ] ∈ C N ×N is a unitary matrix consisting of a complete set of orthonormal eigenvectors B G := {ϕ G, , ∀ ∈ 0, N − 1 }, such that ϕ G, is equal to the -th graph Fourier mode [2], and contains the eigenvalues of L G where without loss of generality, we assume: 0 = λ G,0 < λ G,1 ≤ · · · ≤ λ G,N −1 = λ G,max . Shuman et al. [2] specified that λ for ∈ 0, N − 1 carries the frequency notion in graph setting. The reduced graph frequencies (one may consider as the angular frequencies on graph) is defined as [17] W G := ω G, := π λ G, /ρ G , where ρ G is an upper bound of λ G, defined as where d i = j∈N i w i,j . On the other hand, the graph Fourier transform (GFT) and its inverse can be expressed as [2] GFT: where F G and F −1 G account for the GFT and inverse GFT (IGFT) operator, respectively.
Discrete−Time Harmonic Analysis. Let x(n), n ∈ 1, M be a discrete-time signal of finite length M . The discrete Fourier transform (DFT) operator F D and its inverse F −1 D can be represented in a matrix form as [25] DFT: where x := [x(1), x(2), . . . , x(M )] T is the discrete-time signal in vector form. Here, Φ * D and Φ D account for the DFT and inverse DFT (IDFT) matrices, respectively, such that and the set of angular frequencies is defined as By a classic interpretation, discrete-time domain can be modeled as a M -Cycle graph D with all edge weights equal to unity. Moreover, the symmetric time Laplacian matrix L D , corresponding to the graph D, is defined as the second-order derivative in discrete-time domain up to a negative sign and, as a circulant matrix, it can be diagonalized as [15] where the corresponding eigenvalue matrix is such that λ D,k = 2(1 − cos(ω D,k )) for all k ∈ 0, M − 1 . is defined as X := Φ * G XΦ D where Φ * G and Φ D are the GFT and IDFT matrices, respectively [15,26]. The two-dimensional spectrum of time-vertex signal X, namely joint spectrum, can be described as containing frequency pairs (ω G, , ω D,k ) where ω G, and ω D,k are the -th and k-th angular frequencies on graphs G and D, respectively. In a compact form, JFT and its inverse can be rewritten as [26] JFT: where x = vec(X) is called the "time-vertex signal" and is a unitary matrix. This consists of a complete set of orthonormal eigenvectors such that ϕ J, is equivalent to the -th joint Fourier mode. One can easily see that ϕ J,0 = 1 N M / √ N M . Likewise the graph setting and discrete-time domain, joint time-vertex domain can also be modeled by a joint graph (or product graph) which is equal to the Cartesian product of undirected graph G and M -Cycle graph D (see [26,Figure 2]). Moreover, the joint Laplacian matrix L J , corresponding to the joint graph, can be described as [15] where the corresponding eigenvalue matrix is contains the eigenvalues of L J corresponding to the sum of all eigenvalue pairs of L G and L D (cf. (2), (9)). One can write the -th eigenvalue of L J as follows where

Joint Time-Vertex Translation Operator
First, we briefly discuss the translation operator in continuous-time, discrete-time, and graph domains followed by designing bivariate isometric translation operator in time-vertex domain.
Remark 1. Let x(t) be a continuous-time signal and x(ξ) = (F C x)(t) be its Fourier transform where F C is the Fourier transform operator. The continuous translation operator is defined as (T τ C x)(t) = x(t − τ ) with τ as the translation value. This operator can be formulated in the abstract form as where P τ C = exp (−2πiτ ξ) is the linear phase shift operator. Proof: See Appendix A. where is the angular frequency matrix in discrete-time domain.

Definition 1. (Graph Shift
Operator [17]) The isometric graph shift operator on graph G is defined as The isometric graph shift operator can be written as where P G = exp (−iM G ) such that is a diagonal matrix containing the reduced graph frequencies, and Φ * G is the GFT matrix.
Definition 2. (Joint Translation Operator ) Let T D and T G be the translation operators in discretetime and graph domains, respectively. We define T where ϑ, υ ∈ Z + account for the translation value in graph and discrete-time domains, respectively.

Remark 3.
To represent the bivariate joint translation operation in matrix-vector multiplication form, one can write where x (υ,ϑ) = vec X (υ,ϑ) , x = vec(X), and thanks to the nice property of Kronecker product, the matrix representation of bivariate isometric JTO can be obtained as It should be underlined that the proposed bivariate JTO is genuinely defined in the time-vertex domain as a two-dimensional operation on the matrix X. However, for notational simplicity, in the ensuing presentation, we may use the compact form of this operation given by (24) instead.
Definition 2) is a unitary operator, and hence isometric.
Proof: It is sufficient to prove it for the unit joint time-vertex translation simply denoted by T J . Then we have where the third equality holds since T G and T D are unitary matrices. Similarly, it can be shown that can be written in matrix form as follows where the bivariate phase shift operator is and M which consists of all the combinations of frequencies in discrete-time and graph domains as Proof: From Definition 2, one can write  (17)) where the phase shift matrix is defined as Using (19) and (22), one can easily write the elements of this matrix as (29).
Proposition 2. Some properties associated with the isometric JTO given by (30) are as follows: (i) It is a linear convolutive operator since P (ii) The power spectrum of time-varying signal X is invariant under the operator T where X[n, m] denotes the JFT coefficient of time-varying graph signal X corresponding to the frequency pair (ω G,n−1 , ω D,m−1 ) for all n ∈ 1, N , m ∈ 1, M .

Proof:
We have which clearly reduces to (33) since P ϑ G = exp(iϑM G ) and P υ D = exp(iυM D ).
(iii) It can be expressed as where for all j ∈ 0, N M − 1 denotes the 'j + 1'-th entry in the main diagonal of P (29)). Note that γ is a bivariate function of ω G, and ω D,k . Table 1 summarizes the characteristics of isometric translation/transition operators in different signal domains including continuous-time, discrete-time, graph, and time-vertex domains where they share similar structural characteristics.

Joint Filtering via Joint Translation Operator
In discrete-time domain, a linear shift invariant (LSI) filter is equivalent to the circular convolution operator [25]. Then the filtering operation in this domain can be represented in a compact form as y = H D x where x and y are the input and output signal vectors, respectively, and the filter matrix can be expressed as follows where h D,0 , . . . , h D,L 1 −1 are the L 1 filter coefficients. In the spectral domain, the dual of H D can be expressed as Filtering in the graph setting is defined as x where x and y are the input and output of the filter, respectively, and the graph filter H G ∈ C N ×N in the form of matrix function is defined as [2] where h(Λ G ) is a diagonal matrix and h : [0, λ G,max ] → C (cf. Note 1). Sandryhaila and Moura [5,27] considered weighted adjacency matrix W G as the graph shift operator and showed that any LSI graph filter can be written as a polynomial of W G [5, Theorem 1]. In the same spirit, the graph filter can be written as a polynomial expansion of GTO as follows where h G,q ∈ C is the q-th tap of the filter. In the graph spectral domain, the dual of H G can be written as Filtering in the time-vertex domain, or joint filtering, is defined as y = H J x = h(L G , L D )x where x and y are the input and output of the filter, respectively, and the joint filter is expressed as [15,26] .
be the matrix representation of isometric JTO. A joint filter H J is linear translation invariant (i.e., T Suppose h(x), with x := (x 1 , x 2 ), is the bivariate polynomial of degrees L 1 − 1 and L 2 − 1 such that h(γ J,j ) = O J [j + 1, j + 1] (cf. (36)). In other words, where h J,(p,q) for all p, q are the polynomial coefficients and γ (p,q) J,j is given by (36). Therefore, which clearly reduces to (43) by (35). Next, we prove the sufficiency of the theorem. Since (43) holds true, we have since O J and P J are diagonal, implying that the joint filter H J is linear translation invariant.

Stationarity in Joint Time-Vertex Domain
Let x be a TWSS process in discrete-time domain and R x = E[xx * ] be its autocorrelation matrix. It is known that isometric translation operator in discrete-time domain T D and autocorrelation matrix R x are simultaneously diagonalizable by the DFT matrix.
Via different approaches for defining wide-sense stationarity on graph, it is shown [16] or assumed [18,19] that the GFT matrix diagonalizes the autocorrelation matrix as R and R G,(m,n) is a diagonal matrix.
The bivariate JTO (Definition 2) enables us to generalize wide-sense stationarity to the time-vertex domain as follows.

Definition 6. (JWSS via Bivariate translation Invariance)
A zero-mean stochastic process x = vec(X) defined on graph G, where X ∈ R N ×M , is called JWSS under T is its submatrix in m-th row and n-th column partition, and m, n ∈ 1, M . Then R x can be reexpressed as the following block matrix (2,2) . . . Ξ (2,M ) . . . . . . . . . . . . where and implying that Ξ (m,n) depends on the discrete time difference m − n. Clearly, R x is a block circulant matrix and hence x is a MTWSS process (cf. (ii) in Remark 4). Furthermore, by (59), the submatrices of R x , namely Ξ (m,n) for all n, m are simultaneously diagonalizable with the graph Laplacian. Therefore, it is MVWSS (cf. (ii) in Definition 5). On the other hand, for the proof of sufficiency, assuming that the process x is simultaneously MTWSS and MVWSS, by the reverse implications above, it can be easily shown that x is JWSS. . Then, x can be written as the output of a joint time-vertex filter H J (cf. (49)) in response to a white noise z with zero-mean and autocorrelation matrix R z = I N M . One can deduce that S x = | H J | 2 (cf. (56)). For the simplicity of notation, let θ j := S x [j + 1, j + 1] be the j-th JPSD component where j ∈ 0, N M − 1 . Then by (50), we can write where h : (ω G, , ω D,k ) → C is a bivariate function and γ (p,q) J,j = exp(−i(qω G, + pω D,k )) (cf. (36), (16)). It implies that θ j for all j is a bivariate function of ω G, and ω D,k . On the other hand, if one chooses the isometric GTO -using joint Laplacian matrix, as the translation (or translation) operator on joint graph (cf. (15), (16)) -then the j-th component of JPSD obtained via this approach is where g : λ J,j → C. Clearly, by defining the stationarity on joint graph using isometric GTO, the JPSD is a univariate function of eigenvalues of joint Laplacian matrix, i.e., λ J,j (cf. (15)). From (61), (62), and Theorem 3 one can infer that the JWSS process, in the second-order sense, is more general than wide-sense stationarity on joint graph where the latter is actually a special case of former. The reader is referred to [20] for a further detailed discussion about the difference between the general notion of stationarity in time-vertex domain and stationarity on joint graph.

Joint Power Spectral Density Estimation
Analogous to the stochastic processes in Euclidean space, a reliable JPSD estimator is important for analyzing time-series on graph. In this section, given a data X := {x q , q ∈ 1, Q } where Q is the number of realizations, we present the estimation of JPSD vector denoted by of a JWSS process x = vec(X). The generalized Welch method (GWM) is defined as the average of windowed periodograms (cf. (66) below). So, we begin with the definition of windowing in joint time-vertex domain.
Definition 7. (Joint Windowing) Let x = vec (X) be the given time-varying graph signal. Let A D = Diag(a D ) and A G = Diag(a G ) be the windowing matrices corresponding to the windows a D and a G in discrete-time and graph domains, respectively. The time-vertex windowing is defined as In the vector form, is the joint window matrix and x w = vec(X w ).
By this definition, we have is the dual joint windowing matrix in spectral domain. Following this definition, the generalized Welch JPSD estimator • θ x,GWM is defined as: It is worth noting that by assuming A J = I N M in (66), GWM reduces to the generalized Bartlett method (GBM) -also called sample estimator [20], denoted by • θ x,GBM in Section 7, as an unbiased estimator. The next theorem provides the bias and variance of generalized Welch JPSD estimator.
Remark 6. For a given single realization of a JWSS process denoted by x = vec(X) ∈ C N M , we exploit a bank of joint windows for JPSD estimation. In discrete-time domain, the time-series data of length M is split up into overlapping segments of length L where is the number of windows, · is the floor function, and ∆τ is the length of overlap. By this, we have a set of discrete-time windows A D := {A D,k 1 : k 1 ∈ 1, K 1 }. Moreover, following the same concept of local windowing [19], we obtain a set of graph windows A G := {A G,k 2 : k 2 ∈ 1, K 2 }. Then we come up with a bank of joint windows as follows Then, we calculate where ∈ 1, N M .

Simulation and Experimental Results
This section presents simulation and experimental results to demonstrate the effectiveness of the proposed JPSD estimator  (50)). Likewise [19], we also use Erdős-Rényi graph [31] and Watts-Strogatz small-world graph [32] with N ∈ {100, 200} nodes in our simulations. In discrete-time window, we use the Hamming window with 50% overlapping (cf. (67), (68)). Then, with the obtained bank of joint windows stated in Remark 6, we calculate • θ x,GWM (cf. (69)). Then, for the estimated JPSD • θ x (via GBM or GWM), we compute the normalized mean-squared error (NMSE), bias, and standard deviation (Std) as follows: where E[·] is the average over all the realizations. Figure 1 depicts the NMSE performance of JPSD estimator • θ x,GBM of Q ∈ 1, 10 independent realizations of JWSS process. The Erdős-Rényi and Watts-Strogatz graphs of N ∈ {100, 200} nodes are used, and meanwhile the true JPSD is obtained from a JFIR filter of degrees L 1 = 7 and L 2 = 4. For all the considered cases shown in this figure, NMSE performance of GBM is better for larger Q. However, this estimator suffers from high NMSE when only a small number of realizations of the process are available. Next, we will focus on the case where there exists only a single realization of the process. Figure 2 presents the NMSE performance of • θ x,GWM versus the degree pairs of JFIR filters (L 1 , L 2 ). Here, there is only a single realization of JWSS process available. The Erdős-Rényi and Watts-Strogatz graphs with N ∈ {100, 200} vertices are used. Moreover, the simulation is performed for window number pairs (K 1 , K 2 ) ∈ {(7, 5), (7, 10)}. It can be seen that, for all scenarios, GWM performs significantly better than GBM for Q = 1 as shown in Figure 1. The trend of NMSE linearly and mildly increases with L 1 and L 2 for all the considered cases indicating its efficacy for both graphs. Figure 3 illustrates the performance of GWM w.r.t. the length of discrete-time window L and number of windows in graph setting K 2 (cf. (67), (69)). In this simulation, there is only a single realization of JWSS process and the Erdős-Rényi graph with N = 200 nodes is used. Overall, from Figure 3a, 3b and 3c, one can observe that the larger the L, the better the NMSE performance along with lower bias and higher Std, whereas the variation of NMSE versus K 2 is mild indicating its low sensitivity to K 2 . Furthermore, it can be observed from Figure 3(d), that it is also computationally efficient for smaller K 2 and larger L.
Experimental Results (ER). To show the practicality of JWSS modeling, [20] focused on the recovery of JWSS processes while [23] concentrated on forecasting of time-series on graphs. In our experiment, we present the usefulness of JWSS modeling for the classification of "time-series on graph" via JPSD features, using real EEG data, which can potentially be a new application for graph signal processing. To this end, we model EEG signals as JWSS processes for the purpose of: (1) Emotion recognition and (2) Alzheimer's disease (AD) recognition from EEG data. To the best of our knowledge, Dong [14] is the first who addressed the classification of "graph signals", using features derived from graph spectral information, where the principal component analysis (PCA) and support vector machine (SVM) are used for feature reduction and classification, respectively. Likewise, it is sufficient to exploit the PCA and SVM as the fundamental machine learning tools (rather than advanced tools) in our experiment as in [14].
In our experiments, we divide EEG signals into frames of length 150 ms as the realizations of JWSS processes. Then, estimated JPSD via (66) with A J = I N M is used for the purpose of classification on EEG time-series on graph. It should be noted that our experiments using machine learning approach performed under the same experimental conditions. The common procedure between the two experiments are as follows. We first normalize the EEG dataset D = {x k = vec(X k ) : k ∈ 1, N s } to the range [0, 1]. After feature extraction, the feature matrix is normalized by exploiting the z-score technique which makes the values of each feature have zero-mean (µ z = 0) and unit-variance (σ 2 z = 1); the standardized feature matrix is denoted by Z = [z 1 , . . . , z N M ] where z = (z − µ z 1)/σ z , and z is the -th column of the original feature matrix Z. Note that µ z and σ z are the mean and variance of z , respectively. ER 1 : Emotion Recognition From EEG Data. Recent attempts for emotion recognition using Electroencephalography (EEG) signals have demonstrated its effectiveness in human-machine interactions [33]. In this section, we apply the concept of joint wide-sense stationarity for the challenging task of emotion recognition from brain EEG signals. The SEED-IV [34] is a publicly available EEG dataset obtained from 15 subjects each participating in 3 sessions, each session including 24 trials. In each trial, every participant watched one out of 72 movie clips while his/her EEG signals are collected via 62-channel ESI NeuroScan System. The corresponding EEG channels are illustrated in Figure 4a  which is borrowed from [34]. The samples are categorized into four emotions as fear, happy, sad, and neutral. Our experiments are based on all the N s = 1080 available samples in this dataset. As an initial preprocessing step, we filter out the frequencies less than 2 Hz (also called slow-wave EEG activity or Delta wave [35]) and high frequency noise using a bandpass Butterworth filter (2-150 Hz). Moreover, a notch filter is applied to get rid of the 60 Hz artifact which is a fundamental filtering step in EEG data analysis. Some studies have shown that the asymmetry in neuronal activities between the left and right hemispheres is useful for emotion recognition [36][37][38]. Zhong et al. [33] exploited this differential asymmetry information to initialize the adjacency matrix for developing the graph convolution network for emotion recognition. It is shown experimentally that the following set of channel pairs, denoted by E glb , balances the wiring cost and global efficiency, having particular importance in Emotion recognition [33,39]: where γ = 5.1 is a scaling parameter, κ = 2 if (i, j) ∈ E glb and κ = 1 if (i, j) ∈ E \ E glb , and is the Manhattan distance between two EEG channels i and j with coordinate vectors v i and v j , respectively. Note that the values for γ and κ are chosen empirically. In this modeling, we set κ = 1 for the local inter-channel relations, however, for the global connections we employ κ = 2 due to abovementioned differential asymmetry information between right and left brain hemispheres in emotion recognition. Figure 5 shows the considered Emotion recognition system from EEG signals using JPSD features. Table 2 presents the correct classification rate (CCR) of emotion recognition obtained from three modelings of EEG signals namely TWSS, VWSS, and JWSS using the corresponding PSD estimates in each case as the feature sets. The methods have been tested on the SEED-IV dataset using SVM with radial basis kernel and stratified 10-fold cross validation. Features obtained through TWSS, VWSS, and JWSS modeling reach up to 55.4%, 35.4%, and 59.5% recognition accuracy (on average over all emotions), all using cross-validated PCA. Moreover, the obtained average of recognition accuracy over 30 Monte Carlo runs for the JWSS modeling is 57.8 ± 0.8 which is superior to the accuracy obtained   from TWSS and VWSS modeling as 53.2 ± 1.0 and 32.4 ± 0.9, respectively 5 . Clearly, the feature set obtained from JWSS modeling (cf. Definition 6 and Theorem 3) performs better than VWSS and TWSS modeling in terms of recognition accuracy for emotion recognition from EEG signals.
ER 2 : Alzheimer s Disease Recognition From EEG Data. Alzheimer is a chronic neurodegenaritive disorder which appears in the form of progressive dementia with no specific cure. Due to its major public health challenge, analysis of EEG signals associated with AD is particularly important. Here, we aim at exploiting the PSD and JPSD as the features for AD recognition from EEG signals. For that, we use the AD dataset [40] containing EEG signals taken from 12 AD patients (5 men and 7 women) and 11 healthy control subjects with average age 72.8 ± 8.0 recruited from the Alzheimer's Patients' Relatives Association of Valladolid. Note that there are N s = 663 samples of EEG signals collected via the international 10-20 system (cf. Figure 4b). The N = 16 channels are as follows: F3, F4, F7, F8, FP1, FP2, T3, T4, T5, T6, C3, C4, P3, P4, O1, and O2. To make this experiment more challenging and to show the effectiveness of JPSD in noisy data, we add additive white Gaussian noise to the raw data and do the classification afterwards. We define the weighted adjacency matrix based upon the locations of EEG channels via using (70) with γ = 6 and κ = 1. The methods have been tested on the AD dataset by SVM classifier with radial basis kernel and stratified 10-fold cross validation. The system for AD recognition is the same as emotion recognition

Feature Extraction
Training Samples Figure 5: Emotion and AD recognition system from EEG signals using the JPSD features system which is shown in Figure 5. Figure A.6 depicts the CCR of AD recognition for a range of signal to noise ratio (SNR) associated with features obtained from TWSS, VWSS, and JWSS modeling. Here, the CCR is obtained by averaging recognition accuracy over 50 Monte Carlo simulation runs and the corresponding Std in each case is shown by the bounded vertical bar. Clearly, JWSS modeling has superior performance to the TWSS and VWSS for all the tested SNRs and demonstrates higher robustness of JPSD features against noise. It is worth noting that impact of temporal information for AD recognition is much more significant than the spatial information. Also note that as the SNR increases, the Std of recognition accuracy also decreases.

Additional Declarations
Availability of Data. The data that we used in our experiments are as follows: • SEED IV [34]: The dataset concerning emotion recognition from EEG signals is available at: http://bcmi.sjtu.edu.cn/~seed/seed-iv.html.
• AD dataset [40]: The dataset concerning AD recognition from EEG signals is available at: https://osf.io/jbysn/. Note 2. We need to specify that the datasets concerning with neurological disorders are in general not publicly available, however, this framework can potentially be used in the case of availability of challenging datasets concerning with patients in different stages of such diseases.

Conclusion
We have presented a generalized framework for modeling the stochastic time-varying graph signal as a JWSS process via the proposed bivariate isometric JTO. It was shown that the resulting notion of wide-sense stationarity via bivariate JTO in the time-vertex domain is more general than defining stationarity via isometric GTO on joint graph. This approach can be applied to the JPSD estimation of time-series on graph, from which one can use the resulting JPSD as the features for machine learning based applications. Then we presented the generalized Welch estimator for JPSD estimation, supported by some simulation results. Eventually, we tested the applicability of JPSD features for the classification of time-series on graph including Emotion and AD recognition from EEG signals by modeling them as JWSS processes. The experimental results demonstrated that the proposed JWSS modeling with JPSD features outperforms the classical TWSS and VWSS modelings with PSD and GPSD as features, respectively. In spite of a simple graph learning method used for the classification of EEG data, more effective methods will be our future study for upgrading emotion recognition in the GSP framework. Moreover, as another future study, it is interesting to exploit the JPSD features to classify more challenging datasets with multiple classes associated with different stages of neurological disorders like AD, Parkinson, and Epilepsy. Finally, the JPSD features can potentially be used for other real world applications associated with the classification of time-series on graph, where the JWSS is a suitable modeling.
Appendix A. Proof of Remark 1 Let D C := ∂ t denote the differentiation w.r.t. t and x(ξ) be the Fourier transform of x(t). Using the Taylor series expansion From the theory of Fourier transform, it is easy to verify the following property as F C D C x(t) = i2πξ x(ξ) [41]. In a compact notation one can write is an explicit multiplier. Then one can obtain