A Non-Parametric Algorithm for Predicting Future Samples in Single- and Multi-Channel Time Series

1. Introduction

Difference equations have been widely used to model time series by means of exponential functions. The parameters of the exponential functions result from an optimization process that minimizes the modelling error. These techniques suffer from numerical problems that affect the accuracy of the computed parameters. The frequencies and the dumping factors of the exponential functions are the roots of polynomials that result from the optimization process. However the numerical computation of roots is rather sensitive leading to compromised accuracy. To overcome this problem, alternative approaches for the computation of the exponential parameters were proposed in literature. Among the most efficient techniques are the well known ESPRIT and HSVD algorithms that perform well when it comes to spectral estimation but continue having numerical problems to predict future samples [1,2]. In addition round-off errors add-up creating significant deviation from the correct values of the parameters which compromises the prediction based on extrapolation of the exponential functions. This paper exploits the null space of the Hankel matrix of the data allowing the prediction of future samples with better accuracy and confidence. In a more general sense the notion of null space refers to the set of eigenvectors of the data Hankel matrix which are associated with the smallest eigenvalues.

A generalized form of the new algorithm is also derived and applied to multichannel signals. Multichannel modeling has been widely used in signal processing applications. In particular, time varying multichannel signals are very common in real life applications such as video and multi source audio but also in econometrics while studying the combined evolution of various economic indices. The prediction of future vectors in a multichannel configuration has been an important problem with interesting applications. Among the most popular approaches are the multichannel autoregressive, subspace identification, deep learning and wavelet methods [8]. Here the multichannel autoregressive technique is considered and a new algorithm for the computation of future vectors is presented. It exploits the null space of the block-like Hankel matrix constructed from existing vectors of the multichannel signal. Both cases with and without cross-channel coupling are considered and different algorithms are presented.

2. The Orthogonal Approach

Let us consider the p_th order Hankel matrix H of a given set W of N samples w₁, w₂, … w_N with N>p and size px(N-p+1). The SVD of H is H=U*S*V where U is a pxp matrix of singular vectors of size p, S is a pxp diagonal matrix of singular values, arranged in ascending order and V is a px(N-p+1) matrix of singular vectors of size (N-p+1) [3]. Assuming that the size of the null space of H is q only p-q singular values in S are non-zero and q of the singular vectors of U namely u₁… u_q span the null space of H. If U_q is a pxq matrix of these vectors then the N-p+1 column vectors of H are orthogonal to U_q [4,5,6]. The (N+1)_th unknown sample of the set W that has to be estimated, namely w(N+1), along with the p-1 known samples w(N-p+2), w(N-p+3), w(N-p+4) … w(N) form a new column vector of size p, namely x₁=[w(N-p+2) w(N-p+3) w(N-p+4) … w(N) w(N+1)]^T which is part of the augmented Hankel matrix H_a1=[H x₁]. Since all the vectors of H are orthogonal to U_q the same should hold for the vectors of the augmented matrix H_a1, hence x₁ is also orthogonal to U_q. By splitting x₁ into two parts, the known part namely ${\mathit{s}}_1^{\mathit{T}}$=[w(N-p+2) w(N-p+3) … w(N)] and the unknown part y₁=w(N+1) the following equation is obtained,

$${\mathit{s}}_1^{\mathit{T}}\ast {\mathit{u}}_{1\mathit{i}}+{\mathrm{y}}_1\ast {\mathrm{u}}_{\mathrm{pi}}=0 \text{∀} \mathrm{i}\text{≤}\mathrm{q}$$

(1)

where u_1i consists of the p-1 first components of the i_th singular vector u_i and u_pi is the last component of u_i.

Equation (1) is only theoretically valid. In real life the orthogonality is only approximate and (1) should be reformulated as follows,

$${\mathit{s}}_1^{\mathit{T}}\ast {\mathit{u}}_{1\mathit{i}}+{\mathrm{y}}_1\ast {\mathrm{u}}_{\mathrm{pi}}={\mathrm{e}}_{1\mathrm{i}}\text{∀} \mathrm{i}\text{≤}\mathrm{q}$$

(1.1)

where e_1i is the approximation error. Best approximation will be achieved when the accumulated squared error is minimized,

$${\sum }_{\mathit{i}=1}^{\mathit{q}}{\mathrm{e}}_{1\mathrm{i}}^2\to \mathrm{m}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{u}\mathrm{m}$$

(2)

which leads to the solution of the following minimization problem,

$${\sum }_{\mathit{i}=1}^{\mathit{q}}{\left({\mathit{s}}_1^{\mathit{T}}\ast {\mathit{u}}_{1\mathit{i}}+{\mathrm{y}}_1\ast {\mathrm{u}}_{\mathrm{pi}}\right)}^2\to \mathrm{m}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{u}\mathrm{m}$$

(3)

An equivalent minimization problem results after elaborating the above expression (3),

$${\sum }_{\mathit{i}=1}^{\mathit{q}}\left({\mathit{s}}_1^{\mathit{T}}\ast {\mathit{u}}_{1\mathit{i}}\ast {\mathit{u}}_{1\mathit{i}}^{\mathit{T}}\ast {\mathit{s}}_1+2\ast {\mathrm{y}}_1\ast {\mathrm{u}}_{\mathrm{pi}}\ast {\mathit{u}}_{1\mathit{i}}^{\mathit{T}}\ast {\mathit{s}}_1+{\mathrm{y}}_1^2\ast {\mathrm{u}}_{\mathrm{pi}}^2\right)\to \mathrm{m}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{u}\mathrm{m}$$

(4)

By setting to zero the first derivative of (4) with respect to y₁ the following solution is obtained,

$${\mathrm{y}}_1=-{\sum }_{\mathit{i}=1}^{\mathit{q}}{\mathrm{u}}_{\mathrm{pi}}\ast {\mathit{u}}_{1\mathit{i}}^{\mathit{T}}\ast {\mathit{s}}_1/{\sum }_{\mathit{i}=1}^{\mathit{q}}{\mathrm{u}}_{\mathrm{pi}}^2$$

(5)

3. The Case of Two Unknown Samples

To estimate w(N+1) and w(N+2) simultaneously a similar approach is used. Besides the vector x₁ a second column vector x₂ of size p, consisting of p-2 known samples w(N-p+3), w(N-p+4),.. w(N) and two unknown samples, namely w(N+1) and w(N+2) is formed. The augmented Hankel matrix H_a1 becomes H_a2=[H x₁ x₂] and both vectors x₁ and x₂ have to be orthogonal to U_q. Consequently next to equation (1) an additional equation to guarantee the orthogonality of x₂ is needed. Vector x₂ is split into two parts, the known part s₂ and the unknown components y₁ and y₂ with ${\mathit{s}}_2^{\mathit{T}}$=[w(N-p+3) w(N-p+4) … w(N)] and [y₁ y₂]=[w(N+1) w(N+2)]. Due to the orthogonality property the following must hold,

$${\mathit{s}}_2^{\mathit{T}}\ast {\mathit{u}}_{2\mathit{i}}+{\mathrm{y}}_1\ast {\mathrm{u}}_{\mathrm{p}-1,\mathrm{i}}+{\mathrm{y}}_2\ast {\mathrm{u}}_{\mathrm{pi}}=0 \text{∀} \mathrm{i}\text{≤}\mathrm{q}$$

(6)

where u_2i consists of the p-2 first components of the i_th singular vector u_i and u_p-1,i and u_pi are the last two components of u_i.

As it was mentioned before, equation (6) is only theoretically valid. In real life the orthogonality is only approximate and (6) should be reformulated as follows,

$${\mathit{s}}_2^{\mathit{T}}\ast {\mathit{u}}_{2\mathit{i}}+{\mathrm{y}}_1\ast {\mathrm{u}}_{\mathrm{p}-1,\mathrm{i}}+{\mathrm{y}}_2\ast {\mathrm{u}}_{\mathrm{pi}}={\mathrm{e}}_{2\mathrm{i}}\text{∀} \mathrm{i}\text{≤}\mathrm{q}$$

(7)

where e_2i is the approximation error. For best estimates of both w(N+1) and w(N+2) the following combined optimization problem must be solved,

$${\sum }_{\mathit{i}=1}^{\mathit{q}}{\mathrm{e}}_{1\mathrm{i}}^2+{\sum }_{\mathit{i}=1}^{\mathit{q}}{\mathrm{e}}_{2\mathrm{i}}^2\to \mathrm{m}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{u}\mathrm{m}$$

(8)

By using (1.1) and (7) the above minimization problem can be expressed as,

$${\sum }_{\mathit{i}=1}^{\mathit{q}}{(\mathit{s}}_1^{\mathit{T}}\ast {\mathit{u}}_{1\mathit{i}}+{\mathrm{y}}_1\ast {\mathrm{u}}_{\mathrm{pi}}{)}^2+{\sum }_{\mathit{i}=1}^{\mathit{q}}{\left({\mathit{s}}_2^{\mathit{T}}\ast {\mathit{u}}_{2\mathit{i}}+{\mathrm{y}}_1\ast {\mathrm{u}}_{\mathrm{p}-1,\mathrm{i}}+{\mathrm{y}}_2\ast {\mathrm{u}}_{\mathrm{pi}}\right)}^2\to \mathrm{m}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{u}\mathrm{m}$$

(9)

To solve the above the first derivatives of expression (9) with respect to y₁ and y₂ are set to zero. This leads to the following two equations,

$${\mathrm{y}}_1\ast {\sum }_{\mathit{i}=1}^{\mathit{q}}\left({\mathrm{u}}_{\mathrm{pi}}^2+{\mathrm{u}}_{\mathrm{p}-1,\mathrm{i}}^2\right)+{\mathrm{y}}_2\ast {\sum }_{\mathit{i}=1}^{\mathit{q}}({\mathrm{u}}_{\mathrm{p}-1,\mathrm{i}}\ast {\mathrm{u}}_{\mathrm{pi}})=-{\sum }_{\mathit{i}=1}^{\mathit{q}}({\mathrm{u}}_{\mathrm{pi}}\ast {\mathit{s}}_1^{\mathit{T}}\ast {\mathit{u}}_{1\mathit{i}}+{\mathrm{u}}_{\mathrm{p}-1,\mathrm{i}}\ast {\mathit{s}}_2^{\mathit{T}}\ast {\mathit{u}}_{2i})$$

(10.1)

$${\mathrm{y}}_1\ast {\sum }_{\mathit{i}=1}^{\mathit{q}}{\mathrm{u}}_{\mathrm{p}-1,\mathrm{i}}\ast {\mathrm{u}}_{\mathrm{pi}}+{\mathrm{y}}_2\ast {\sum }_{\mathit{i}=1}^{\mathit{q}}{\mathrm{u}}_{\mathrm{pi}}^2=-{\sum }_{\mathit{i}=1}^{\mathit{q}}{\mathrm{u}}_{\mathrm{pi}}\ast {\mathit{s}}_2^{\mathit{T}}\ast {\mathit{u}}_{2\mathit{i}}$$

(10.2)

The above linear system of equations provides best estimates for y ₁ and y₂.

4. The General Case of Many Unknown Samples

When three future samples need to be estimated, a third term in the optimization problem (8) is added [7,9]. By following the same type of notation equation (8) becomes,

$${\sum }_{\mathit{i}=1}^{\mathit{q}}{\mathrm{e}}_{1\mathrm{i}}^2+{\sum }_{\mathit{i}=1}^{\mathit{q}}{\mathrm{e}}_{2\mathrm{i}}^2+{\sum }_{\mathit{i}=1}^{\mathit{q}}{\mathrm{e}}_{3\mathrm{i}}^2\to \mathrm{m}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{u}\mathrm{m}$$

(11)

where e_3i is the error due to the third unknown sample w(N+3)=y₃

By following the same notation as in the case of two unknown samples the error e_3i can be expressed as follows,

$${\mathit{s}}_3^{\mathit{T}}\ast {\mathit{u}}_{3\mathit{i}}+{\mathrm{y}}_1\ast {\mathrm{u}}_{\mathrm{p}-2,\mathrm{i} }+{\mathrm{y}}_2\ast {\mathrm{u}}_{\mathrm{p}-1,\mathrm{i}}+{\mathrm{y}}_3\ast {\mathrm{u}}_{\mathrm{pi}}={\mathrm{e}}_{3\mathrm{i}}\text{∀} \mathrm{i}\text{≤}\mathrm{q}$$

(12)

where ${\mathit{s}}_3^{\mathit{T}}$=[w(N-p+4) w(N-p+4) … w(N)], u_3i consists of the p-3 first components of the i_th singular vector u_i and u_p-2,i , u_p-1,i and u_pi are the last three components of u_i.

By setting the first derivatives of expression (11) with respect to y₁, y₂ and y₃ to zero, a linear system of three equations is formed and its solution provides the estimates of the three unknown samples y₁, y₂ and y₃ .

By replacing (1), (7) and (12) into expression (11) and by elaborating on the algebraic manipulations, the following set of equations is obtained,

$${\mathrm{y}}_1\ast {\sum }_{\mathit{i}=1}^{\mathit{q}}\left({\mathrm{u}}_{\mathrm{pi}}^2+{\mathrm{u}}_{\mathrm{p}-\mathrm{i},1}^2+{\mathrm{u}}_{\mathrm{p}-2,\mathrm{i}}^2\right)+{\mathrm{y}}_2\ast {\sum }_{\mathit{i}=1}^{\mathit{q}}({\mathrm{u}}_{\mathrm{pi}}\ast {\mathrm{u}}_{\mathrm{p}-1,\mathrm{i}}+{\mathrm{u}}_{\mathrm{p}-1,\mathrm{i}}\ast {\mathrm{u}}_{\mathrm{p}-2,\mathrm{i}})+{\mathrm{y}}_3\ast {\sum }_{\mathit{i}=1}^{\mathit{q}}{\mathrm{u}}_{\mathrm{pi}}\ast {\mathrm{u}}_{\mathrm{p}-2,\mathrm{i}}=-{\sum }_{\mathit{i}=1}^{\mathit{q}}({\mathrm{u}}_{\mathrm{pi}}\ast {\mathit{s}}_1^{\mathit{T}}\ast {\mathit{u}}_{1\mathit{i}}+{\mathrm{u}}_{\mathrm{p}-1,\mathrm{i}}\ast {\mathit{s}}_2^{\mathit{T}}\ast {\mathit{u}}_{2\mathit{i}}+{\mathrm{u}}_{\mathrm{p}-2,\mathrm{i}}\ast {\mathit{s}}_3^{\mathit{T}}\ast {\mathit{u}}_{3\mathit{i}})$$

(13.1)

$${\mathrm{y}}_1\ast {\sum }_{\mathit{i}=1}^{\mathit{q}}\left({\mathrm{u}}_{\mathrm{pi}}\ast {\mathrm{u}}_{\mathrm{p}-1,\mathrm{i}}+{\mathrm{u}}_{\mathrm{p}-1,\mathrm{i}}\ast {\mathrm{u}}_{\mathrm{p}-2,\mathrm{i}}\right)+{\mathrm{y}}_2\ast {\sum }_{\mathit{i}=1}^{\mathit{q}}\left({\mathrm{u}}_{\mathrm{pi}}^2+{\mathrm{u}}_{\mathrm{p}-1,\mathrm{i}}^2\right)+{\mathrm{y}}_3\ast {\sum }_{\mathit{i}=1}^{\mathit{q}}{\mathrm{u}}_{\mathrm{pi}}\ast {\mathrm{u}}_{\mathrm{p}-1,\mathrm{i}}=-{\sum }_{\mathit{i}=1}^{\mathit{q}}{(\mathrm{u}}_{\mathrm{pi}}\ast {\mathit{s}}_2^{\mathit{T}}\ast {\mathit{u}}_{2\mathit{i}}+{\mathrm{u}}_{\mathrm{p}-1\mathrm{i}}\ast {\mathit{s}}_3^{\mathit{T}}\ast {\mathit{u}}_{3\mathit{i}})$$

(13.2)

$${\mathrm{y}}_1\ast {\sum }_{\mathit{i}=1}^{\mathit{q}}{\mathrm{u}}_{\mathrm{pi}}\ast {\mathrm{u}}_{\mathrm{p}-2,\mathrm{i}}+{\mathrm{y}}_2\ast {\sum }_{\mathit{i}=1}^{\mathit{q}}{\mathrm{u}}_{\mathrm{pi}}\ast {\mathrm{u}}_{\mathrm{p}-1,\mathrm{i}}+{\mathrm{y}}_3\ast {\sum }_{\mathit{i}=1}^{\mathit{q}}{\mathrm{u}}_{\mathrm{pi}}^2=-{\sum }_{\mathit{i}=1}^{\mathit{q}}{\mathrm{u}}_{\mathrm{pi}}\ast {\mathit{s}}_3^{\mathit{T}}\ast {\mathit{u}}_{3\mathit{i}}$$

(13.3)

The above system of equations (13) maybe expressed as,

$$\mathit{y}\ast \mathit{A}=\mathit{b}$$

(14)

where y=[y₁ y₂ y₃], A is a symmetric matrix and b is a vector.

In a more general notation, the elements of the symmetric matrix A for all j ≤ k and k≤n can be expressed as,

$${\mathrm{a}}_{\mathrm{jk}}={\sum }_{\mathit{i}=1}^{\mathit{q}}{\sum }_{\mathit{m}=1}^{\mathit{n}-\mathit{k}+1}{\mathrm{u}}_{\mathrm{p}+1-\mathrm{m},\mathrm{i}}\ast {\mathrm{u}}_{\mathrm{p}+\mathrm{j}-\mathrm{k}+1-\mathrm{m},\mathrm{i}}$$

(15)

where n stands for the number of samples to be estimated, namely n=3 in this case.

Similarly, the elements of b can be expressed as,

$${\mathit{b}}_{\mathit{j}}=-{\sum }_{\mathit{i}=1}^{\mathit{q}}{\sum }_{\mathit{k}=\mathit{j}}^{\mathit{n}}{\mathrm{u}}_{\mathrm{p}+\mathrm{j}-\mathrm{k},\mathrm{i}}\ast {\mathit{s}}_{\mathit{k}}^{\mathit{T}}\ast {\mathit{u}}_{\mathit{k}\mathit{i}}$$

(16)

With A and b known, the estimation of y is straightforward based on equation (14),

$$\mathit{y}=\mathit{b}\ast {\mathit{A}}^{-1}$$

(17)

where y is the vector of n unknown future samples of the set W.

5. The Multichannel Solution Without Cross Channel Coupling

The above technique can be generalized and applied to the multichannel configuration. In particular for the case without cross channel coupling each vector is a linear combination of previous vectors according to the following equation,

$${\mathit{x}}_{\mathit{n}}=a_1\ast {\mathit{x}}_{\mathit{n}-1}+a_2\ast {\mathit{x}}_{\mathit{n}-2}+\dots a_{p-1}\ast {\mathit{x}}_{\mathit{n}-\mathit{p}+1}+{\mathit{e}}_{\mathit{n}}$$

(18)

where the coefficients $a_i, i=1\dots p-1$ are scalars and ${\mathit{x}}_{\mathit{n}} , n=1,\dots N$ are the N vectors of size m of the multichannel signal [8].

If matrix H_b in (19) stands for the block Hankel matrix of the above multichannel signal, the error ${\mathit{e}}_n$ for all $n=p\dots N,$ will be zero when the vector ${\mathit{a}}^{\mathit{T}}=\left[ a_{p-1}\dots a_1-1\right]$ is part of the null space of matrix H_b

$$\mathbf{H}\mathbf{b}=\left(\begin{array}{ccc}{\mathit{x}}_1^T& \dots & {\mathit{x}}_{N-p+1}^T\\ ⋮& \ddots & ⋮\\ {\mathit{x}}_p^T& \dots & {\mathit{x}}_N^T\end{array} \right)$$

(19)

Assuming that the augmented block Hankel matrix H_ba=[H_b h₁] shares the same null space as H_b, where h₁ is a pxm matrix with its last row the unknown vector ${\mathit{x}}_{N+1}$ as shown below,

$${\mathbf{h}1=\left({\mathit{x}}_{\mathit{N}-\mathit{p}+2}\dots {\mathit{x}}_{\mathit{N}}{\mathit{x}}_{\mathit{N}+1}\right)}^T$$

(20)

then all column vectors of matrix h₁ should also be orthogonal to the null space of matrix H_b

The SVD of H_b is H_b=U_b*S_b*V_b where U_b is a pxp matrix of singular vectors of size p, S_b is a pxp diagonal matrix of singular values, arranged in ascending order and V_b is a px[(N-p+1)*m] matrix of singular vectors of size (N-p+1)*m. Assuming that the size of the null space of H_b is q only p-q singular values in S_b are non-zero and q of the singular vectors of U_b namely u_b1… u_bq span the null space of H_b. If U_bq is a pxq matrix of these vectors then the (N-p+1)*m column vectors of H_b are orthogonal to U_bq. The m column vectors of h₁ which should also be orthogonal to U_bq, consist of (p-1) known values and one unknown, their last element. They can be treated independently with null space the vectors in U_bq. Therefore to determine each of the m elements of vector ${\mathit{x}}_{N+1}$ , equation (5) can be used separately with each one of the m column vectors of matrix ${\mathit{h}}_1$, after replacing the corresponding null space vectors u_b1… u_bq.

The elements of future vectors ${\mathit{x}}_{N+k}$, k=2,…n can be determined by applying the equations (15), (16) and (17) and by using the corresponding null space vectors u_b1… u_bq.

More specifically, by introducing the matrices ${{\mathit{h}}_k=\left({\mathit{x}}_{\mathit{N}-\mathit{p}+\mathit{k}+1}\dots {\mathit{x}}_{\mathit{N}}{\mathit{x}}_{\mathit{N}+\mathit{k}}\right)}^{T}k=2,\dots n$ and assuming that they share the same null space as H_b, each one of the m column vectors of these matrices can be processed separately using equations (15), (16) and (17) to provide n future estimates for each individual channel.

6. The Multichannel Solution with Cross Channel Coupling

When cross channel coupling is assumed, each vector is a linear combination of previous vectors according to the following equation,

$${\mathit{x}}_{\mathit{n}}={\mathit{A}}_1\ast {\mathit{x}}_{\mathit{n}-1}+{\mathit{A}}_2\ast {\mathit{x}}_{\mathit{n}-2}+\dots {\mathit{A}}_{p-1}\ast {\mathit{x}}_{\mathit{n}-\mathit{p}+1}+{\mathit{e}}_{\mathit{n}}$$

(21)

where the coefficients ${\mathit{A}}_i, i=1\dots p-1$ are mxm matrices, ${\mathit{x}}_{\mathit{n}} , n=1,\dots N$ are the N vectors of size m and ${\mathit{e}}_n$ is the modelling error. By expanding equation (21) the k_th element of vector ${\mathit{x}}_n$ , denoted as $x_{nk}$, is equal to

$$x_{nk}={\mathit{a}}_{1k}*{\mathit{x}}_{n-1}+{\mathit{a}}_{2k}*{\mathit{x}}_{n-2}+\dots {\mathit{a}}_{\left(p-1\right)k}*{\mathit{x}}_{n-p+1}+e_{nk}$$

(22)

where ${\mathit{a}}_{1k} , {\mathit{a}}_{2k} , \dots {\mathit{a}}_{\left(p-1\right)k}$are the k_th rows of the matrices, ${\mathit{A}}_1 , {\mathit{A}}_2 , \dots {\mathit{A}}_{p-1}$ respectively and $e_{nk}$is the k_th element of the modelling error.

If matrix H_bk in (23) stands for the block-like Hankel matrix which is associated to the k_th element of the multichannel signal, the error $e_{nk}$ for all $n=p\dots N,$ will be zero when the vector ${\mathit{a}}^{\mathit{T}}=\left[ {\mathit{a}}_{(p-1)k}\dots {\mathit{a}}_{1k}-1\right]$ is part of the null space of matrix H_bk,

$${\mathit{H}}_{bk}=\left(\begin{array}{ccc}{\mathit{x}}_1& \begin{array}{c} \\ {\mathit{x}}_2\\ \end{array} \dots & {\mathit{x}}_{N-p+1}\\ \begin{array}{c}\begin{array}{c}⋮\\ {\mathit{x}}_{p-1}\end{array}\\ \end{array}\begin{array}{c} \\ \end{array}& \begin{array}{c}\begin{array}{cc}⋮& \end{array}\\ \begin{array}{cc}{\mathit{x}}_p& \dots \end{array}\end{array}& \begin{array}{c}⋮\\ \begin{array}{ccc} & {\mathit{x}}_{N-1}& \end{array}\end{array}\\ x_{pk}& \begin{array}{cc}{x}_{\left(p+1\right)k}& \end{array}\dots & x_{Nk}\end{array}\right)$$

(23)

Assuming that the augmented block-like Hankel matrix H_bka=[H_bk h_k] shares the same null space as H_bk, where h_k is a (p-1)*m+1 vector with its last element the k_th unknown sample $x_{\left(N+1\right)k}$ of the vector ${\mathit{x}}_{N+1}$, as shown below,

$${\mathit{h}}_k={\left({\mathit{x}}_{\mathit{N}-\mathit{p}+2}^{\mathit{T}}\dots {\mathit{x}}_{\mathit{N}}^{\mathit{T}}{x}_{\left(N+1\right)k}\right)}^T$$

(24)

then vector h_k should also be orthogonal to the null space of matrix H_bk.

The SVD of H_bk is H_bk=U_bk*S_bk*V_bk where U_bk is a [(p-1)*m+1]x[(p-1)*m+1] matrix of singular vectors of size (p-1)*m+1, S_bk is a [(p-1)*m+1]x[(p-1)*m+1] diagonal matrix of singular values, arranged in ascending order and V_bk is a [(p-1)*m+1]x(N-p+1) matrix of singular vectors of size (N-p+1). Assuming that the size of the null space of H_bk is q only [(p-1)*m+1]-q singular values in S_bk are non-zero and q of the singular vectors of U_bk namely u_bk1… u_bkq span the null space of H_bk. If U_bkq is a [(p-1)*m+1]xq matrix of these vectors then the (N-p+1) column vectors of H_bk along with vector h_k are orthogonal to U_bkq. The vector h_k consist of (p-1)*m known values and one unknown, its last element. The estimation of the unknown element $x_{\left(N+1\right)k}$ can be done independently with corresponding null space the vectors in U_bkq. Hence for the unknown element of vector ${\mathit{h}}_k$ equation (5) can be used with corresponding null space the q vectors u_bk1… u_bkq.

In order to estimate all the elements of the unknown vector ${\mathit{x}}_{N+1}$ the same procedure is repeated for all m possible values of k. For each value of k the SVD of the block-like Hankel matrix in (23) has to be computed and the new null space vectors have to be used with equation (5).

7. Experimentation – Testing

The proposed algorithm was tested against the well known method ESPRIT that models the signal by means of exponential sinusoids. A large set of signals, four thousands in total, with variable complexity in terms of the rank of their covariance matrix was considered. The comparison was based on the accuracy of predicting future samples as well as on the overall consistency of the two algorithms. The term consistency refers to the homogeneity of the prediction errors. Each signal consists of 2000 to 2500 samples and it was generated from white noise that was filtered using random low pass FIR filters with hundreds of taps. Four sets of future samples were selected for the comparison. The latter is based on the RMSE among the predicted values and the actual values of the signals. The four sets refer to predicting the first, tenth, twentieth and thirtieth future samples using both the proposed and the ESPRIT algorithms and comparing the RMSEs of the two methods. To evaluate the algorithms with respect to their consistency the variance of the errors is used. In order to eliminate variations due to the differences of the signal amplitudes, all errors were normalized with respect to the actual value of the predicted sample.

For the ESPRIT method the prediction for each one of the 4000 signals is based on twenty different models with orders equally spaced between 140 and 900. By averaging these twenty different results the overall prediction is computed.

For the new method the size p of the Hankel matrix H is set to 600. Alternative choices for p were also tested showing no significant change in the results as long as p is not too small compared to the length N of the signals. Regarding the null space not only one single size for each signal is used. Instead, an optimum set of null spaces for each signal is considered and the overall prediction is computed by averaging individual results. There are two constraints while choosing the optimum set of null spaces, a) the sizes of the null spaces of the optimum set are consecutive numbers and b) the variance of the results of the set is minimal. The second constraint is imposed to guarantee maximum consistency of the chosen set.

Table 1 shows the obtained results for the four different setups in predicting future samples. In all cases the new method performs better than ESPRIT in terms of both RMSE and variance.

Figure 1, Figure 2, Figure 3 and Figure 4 show the shape of representative signals among the 4000 used. The rank of the covariance matrix of all 4000 signals varies considerably, ranging from low to high rank. This way the results are not biased with respect to the degree of predictability of the signals.

8. Conclusion

A new method to predict future samples is presented and it is compared against the ESPRIT method. It was proven to be more accurate and more consistent than ESPRIT, based on results from four thousand time series of variable spectral complexity. The new method exploits the properties of the eigenvectors of the Hankel matrix of the time series that are associated with the low eigenvalues. A generalization of the algorithm to multichannel signals was also presented. Both cases, with and without cross-channel coupling, were considered and different algorithms were proposed. The algorithms are fast and stable and can predict any set of future samples/vectors. Further research will focus on both the optimal selection of the null space of the time series Hankel matrix as well as the choice of the model size.