Fanny Godet
email: fanny.godet@math.univ-nantes.fr

Linear Prediction of Long-Memory Processes: Asymptotic Results on Mean-squared Errors

Keywords: long-memory, linear model, autoregressive process, forecast error ARMA

We present two approaches for linear prediction of long-memory time series. The first approach consists in truncating the Wiener-Kolmogorov predictor by restricting the observations to the last k terms, which are the only available values in practice. We derive the asymptotic behaviour of the mean-squared error as k tends to +∞. By contrast, the second approach is non-parametric. An AR(k) model is fitted to the long-memory time series and we study the error that arises in this misspecified model.

In Section 3 we discuss the asymptotic properties of the forecast error if we fit a misspecified AR(k) model to a long-memory time series. This approach has been proposed by [Ray, 1993] for fractional noise series F(d). His simulations show that high-order AR-models predict fractional integrated noise very well.

Finally in Section 4 we compare the two previous methods for h-step prediction. We give some asymptotic properties of the mean-squared error of the linear least-squares predictor as h tends to +∞ in the particular case of long-memory processes. Then we study our k-th order predictors order as k tends to +∞.

Model

Let (X n ) n∈Z be a discrete-time (weakly) stationary process in L 2 with mean 0 and σ its autocovariance function. We assume that the process (X n ) n∈Z is a long-memory process i.e.:

∞ k=-∞ |σ(k)| = ∞.
The process (X n ) n∈Z admits an infinite moving average representation as follows:

X n = ∞ j=0 b j ε n-j (1)
where (ε n ) n∈Z is a white-noise series consisting of uncorrelated random variables, each with mean zero and variance σ 2 ε and (b j ) j∈N are square-summable. We shall further assume that (X n ) n∈Z admits an infinite autoregressive representation:

ε n = ∞ j=0 a j X n-j , (2)
where the (a j ) j∈N are absolutely summable. We assume also that (a j ) j∈N and (b j ) j∈N , occurring respectively in (2) and (1), satisfy the following conditions for all δ > 0:

|a j | ≤ C 1 j -d-1+δ
(3)

|b j | ≤ C 2 j d-1+δ . ( 4
)
where C 1 and C 2 are constants and d is a parameter verifying d ∈]0, 1/2[. For example, a FARIMA process (X n ) n∈Z is the stationary solution to the difference equations:

φ(B)(1 -B) d X n = θ(B)ε n
where (ε n ) n∈Z is a white noise series, B is the backward shift operator and φ and θ are polynomials with no zeroes on the unit disk. Its coefficients verify

|a j | ∼ +∞ C 1 j -d-1 |b j | ∼ +∞ C 2 j d-1
and thus (3) and (4) hold. When φ = θ = 1, the process (X n ) n∈Z is called fractionally integrated noise and denoted F(d). More generally, (a j ) j∈N and (b j ) j∈N verify conditions (3) and (4) if:

|a j | ∼ +∞ L(j)j -d-1 |b j | ∼ +∞ L ′ (j)j d-1
where L and L ′ are slowly varying functions. A positive function L is a slowly varying function in the sense of [Zygmund, 1968] if, for any δ > 0, x → x -δ L(x) is decreasing and x → x δ L(x) is increasing.

The condition (4) implies that the autocovariance function σ of the process (X n ) n∈Z verifies:

∀δ > 0, ∃C 3 ∈ R, |σ(j)| ≤ C 3 j 2d-1+δ . (5)
Notice that it suffices to prove (5) for δ near 0 in order to verify (5) for δ > 0 arbitrarily chosen. So we prove (5) for δ < 1-2d 2 :

σ(k) = +∞ j=0 b j b j+k |σ(k)| ≤ +∞ j=1 |b j b j+k | + |b 0 b k | ≤ C 2 2 +∞ j=1 j d-1+δ (k + j) d-1+δ + |b 0 b k | ≤ C 2 2 +∞ 0 j d-1+δ (k + j) d-1+δ dj + |b 0 b k | ≤ C 2 2 k 2d-1+2δ +∞ 0 j d-1+δ (1 + j) d-1+δ dj + C 2 k d-1+δ ≤ C 3 k 2d-1+2δ
More accurately, [Inoue, 1997] has proved than if:

b j ∼ L (j) j d-1 then σ(j) ∼ j 2d-1 [L (j)] 2 β(1 -2d, d)
where L is a slowly varying function and β is the beta function. The converse is not true, we must have more assumptions about the series (b j ) j∈N in order to get an asymptotic equivalent for (σ(j)) j∈N (see [Inoue, 2000]).

2 Wiener-Kolmogorov Next Step Prediction Theory

Wiener-Kolmogorov Predictor

The aim of this part is to compute the best linear one-step predictor (with minimum mean-square distance from the true random variable) knowing all the past {X k+1-j , j 1}. Our predictor is therefore an infinite linear combination of the infinite past:

X k (1) = ∞ j=0 λ(j)X k-j
where (λ(j)) j∈N are chosen to ensure that the mean squared prediction error:

E X k (1) -X k+1 2
is as small as possible. Following [Whittle, 1963], and in view of the moving average representation of (X n ) n∈Z , we may rewrite our predictor X k (1) as:

X k (1) = ∞ j=0 φ(j)ε k-j .
where (φ(j)) j∈N depends only on (λ(j)) j∈N and (ǫ n ) n∈Z and (a j ) j∈N are defined in (2). From the infinite moving average representation of (X n ) n∈Z given below in (1), we can rewrite the meansquared prediction error as:

E X k (1) -X k+1 2 = E     ∞ j=0 φ(j)ε k-j - ∞ j=0 b(j)ε k+1-j   2   = E     ε k+1 - ∞ j=0 (φ(j) -b(j + 1)) ε k-j   2   =   1 + ∞ j=0 b j+1 -φ(j) 2   σ 2 ε
since the random variables (ε n ) n∈Z are uncorrelated with variance σ 2 ε . The smallest mean-squared prediction error is obtained when setting φ(j) = b j+1 for j ≥ 0.

The smallest prediction error of (X n ) n∈Z is σ 2 ε within the class of linear predictors. Furthermore, if

A(z) = +∞ j=0 a j z j ,
denotes the characteristic polynomial of the (a(j)) j∈Z and

B(z) = +∞ j=0 b j z j ,
that of the (a(j)) j∈Z , then in view of the identity, A(z) = B(z) -1 , |z| ≤ 1, we may write:

X k (1) = - ∞ j=1 a j X k+1-j .
(6)

Mean Squared Prediction Error when the Predictor is Truncated

In practice, we only know a finite subset of the past, the one which we have observed. So the predictor should only depend on the observations. Assume that we only know the set {X 1 , . . . , X k } and that we replace the unknown values by 0, then we have the following new predictor:

X ′ k (1) = - k j=1 a j X k+1-j . (7)
It is equivalent to say that we have truncated the infinite series (6) to k terms. The following proposition provides us the asymptotic properties of the mean squared prediction error as a function of k.

Proposition 2.2.1. Let (X n ) n∈Z be a linear stationary process defined by (1), ( 2) and verifying conditions (3) and (4). We can approximate the mean-squared prediction error of X ′ k (1) by:

∀δ > 0, E X k+1 -X ′ k (1) 2 = σ 2 ε + O(k -1+δ
). Furthermore, this rate of convergence O(k -1 ) is optimal since for fractionally integrated noise, we have the following asymptotic equivalent:

E X k+1 -X ′ k (1) 2 = σ 2 ε + Ck -1 + o k -1
. Note that the prediction error is the sum of σ 2 ε , the error of Wiener-Kolmogorov model and the error due to the truncation to k terms which is bounded by O(k -1+δ ) for all δ > 0.

Proof.

X k+1 -X ′ k (1) = X k+1 -X k (1) + X k (1) -X ′ k (1) = X k+1 - +∞ j=0 b j+1 ε k-j - +∞ j=k+1 a j X k+1-j = ε k+1 - +∞ j=k+1 a j X k+1-j . (8)
The two parts of the sum (8) are orthogonal for the inner product associated with the mean square norm. Consequently:

E X k+1 -X ′ k (1) 2 = σ 2 ε + ∞ j=k+1 ∞ l=k+1 a j a l σ(l -j).
For the second term of the sum we have:

+∞ j=k+1 +∞ l=k+1 a j a l σ(l -j) = 2 +∞ j=k+1 a j +∞ l=j+1 a l σ(l -j) + +∞ j=k+1 a 2 j σ(0) ≤ 2 +∞ j=k+1 |a j | |a j+1 | |σ(1)| + +∞ j=k+1 a 2 j σ(0) +2 +∞ j=k+1 |a j | +∞ l=j+2 |a l ||σ(l -j)|
from the triangle inequality, it follows that:

+∞ j=k+1 +∞ l=k+1 a j a l σ(l -j) ≤ C 2 1 C 3   2 +∞ j=k+1 j -d-1+δ (j + 1) -d-1+δ + +∞ j=k+1 j -d-1+δ 2   (9) + 2C 2 1 C 3 +∞ j=k+1 j -d-1+δ +∞ l=j+2 l -d-1+δ |l -j| 2d-1+δ (10)
for all δ > 0 from inequalities (3) and ( 5). Assume now that δ < 1/2 -d. For the terms (9), since j → j -d-1+δ (j + 1) -d-1+δ is a positive and decreasing function on R + , we have the following approximations:

2C 2 1 C 3 +∞ j=k+1 j -d-1+δ (j + 1) -d-1+δ ∼ 2C 2 1 C 3 +∞ k j -d-1+δ (j + 1) -d-1+δ dj ∼ 2C 2 1 C 3 1 + 2d -2δ k -2d-1+2δ
Since the function j → j -d-1+δ 2 is also positive and decreasing, we can establish in a similar way that:

C 2 1 C 3 +∞ j=k+1 j -d-1+δ 2 ∼ C 2 1 C 3 +∞ k j -d-1+δ 2 dj ∼ C 2 1 C 3 1 + 2d -2δ k -2d-1+2δ .
For the infinite double series (10), we will similarly compare the series with an integral. In the next Lemma, we establish the necessary result for this comparison:

Lemma 2.2.1. Let g the function (l, j) → j -d-1+δ l -d-1+δ |l -j| 2d-1+δ
. Let m and n be two positive integers. We assume that δ < 1 -2d and m ≥ δ-d-1 δ+2d-1 for all δ ∈ 0, δ-d-1 δ+2d-1 . We will call

A n,m the square [n, n + 1] × [m, m + 1]. If n ≥ m + 1 then An,m g(l, j) dj dl ≥ g(n + 1, m).
Assume now that δ < 1 -2d without loss of generality. Thanks to the previous Lemma and the asymptotic equivalents of (9), there exists

K ∈ N such that if k > K: +∞ j=k+1 +∞ l=k+1 a j a l σ(l -j) ≤ C +∞ k+1 j -d-1+δ +∞ j l -d-1+δ (l -j) 2d-1+δ dl dj + O k -2d-1+2δ
By using the substitution jl ′ = l in the integral over l we obtain:

+∞ j=k+1 +∞ l=k+1 a j a l σ(l -j) ≤ C ′ +∞ k+1 j -2+3δ +∞ 1 l -d-1+δ (l -1) 2d-1+δ dldj + O k -2d-1 . Since if δ < (1 -d)/2 +∞ 1 l -d-1+δ (l -1) 2d-1+δ dl < +∞, it follows: +∞ j=k+1 +∞ l=k+1 a j a l σ(l -j) ≤ O k -1+3δ + O k -2d-1 ≤ O k -1+3δ . (11) If δ > 0, δ < 1 -2d and δ < (1 -d)/2, we have: +∞ j=k+1 +∞ l=k+1 a j a l σ(l -j) = O k -1+3δ .
Notice that if the equality is true under the assumptions δ > 0, δ < 1 -2d and δ < (1 -d)/2, it is also true for any δ > 0. Therefore we have proven the first part of the theorem.

We prove now that there exists long-memory processes whose prediction error attains the rate of convergence k -1 . Assume now that (X n ) n∈Z is fractionally integrated noise F(d), which is the stationary solution of the difference equation:

X n = (1 -B) -d ε n (12)
with B the usual backward shift operator, (ε n ) n∈Z is a white-noise series and d ∈ ]0, 1/2[ (see for example [START_REF] Brockwell | Time Series : Theory and Methods[END_REF]). We can compute the coefficients and obtain that:

∀j > 0, a j = Γ(j -d) Γ(j + 1)Γ(-d) and ∀j ≥ 0, σ(j) = (-1) j Γ(1 -2d) Γ(j -d + 1)Γ(1 -j -d)
σ 2 ε then we have: ∀j > 0, a j < 0 and ∀j ≥ 0, σ(j) > 0 and

a j ∼ j -d-1 Γ(-d) and σ(j) ∼ j 2d-1 Γ(1 -2d) Γ(d)Γ(1 -d) when j → ∞.
In this particular case, we can estimate the prediction error more precisely:

+∞ k+1 +∞ k+1 a j a l σ(l -j) = +∞ k+1 |a j | +∞ j+1 |a l ||σ(l -j)| + +∞ k+1 a 2 j σ(0) ∼ Γ(1 -2d) Γ(-d) 2 Γ(d)Γ(1 -d) +∞ k+1 j -2 +∞ 1/j+1 l -d-1 (l -1) 2d-1 dldj + O k -2d-1 +∞ k+1 +∞ k+1 a j a l σ(l -j) ∼ Γ(1 -2d)Γ(2d) Γ(-d) 2 Γ(d)Γ(1 + d) k -1 (13)
The asymptotic bound O(k -1 ) is therefore as small as possible.

In the specific case of fractionally integrated noise, we may write the prediction error as:

E X k+1 -X ′ k (1) 2 = σ 2 ε + C(d)k -1 + o k -1
and we can express C(d) as a function of d:

C(d) = Γ(1 -2d)Γ(2d) Γ(-d) 2 Γ(d)Γ(1 + d) . ( 14
)
It is easy to prove that C(d) → +∞ as d → 1/2 and we may write the following asymptotic equivalent as d → 1/2:

C(d) ∼ 1 (1 -2d)Γ(-1/2) 2 Γ(1/2)Γ(3/2) . ( 15
)
As d → 0, C(d) → 0 and we have the following equivalent as d → 0:

C(d) ∼ d 2 . Figure 2.1: Behaviour of constant C(d), d ∈ [0, 1/2[, defined in (14) 0.0 0.1 0.2 0.3 0.4 0.0 0.1 0.2 0.3 0.4 d C(d)
As the figure 2.1 suggests and the asymptotic equivalent given in (15) proves, the mean-squared error tends to +∞ as d → 1/2. By contrast, the constant C(d) takes small values for d in a large interval of [0, 1/2[. Although the rate of convergence has a constant order k -1 , the forecast error is bigger when d → 1/2. This result is not surprising since the correlation between the random variable, which we want to predict, and the random variables, which we take equal to 0, increases when d → 1/2.

Truncating to k terms the series which defines the Wiener-Kolmogorov predictor amounts to using an AR(k) model for predicting. Therefore in the following section we look for the AR(k) which minimizes the forecast error.

The Autoregressive Models Fitting Approach

In this section we develop a generalisation of the "autoregressive model fitting" approach developed by [Ray, 1993] in the case of fractionally integrated noise F(d) (defined in ( 12)). We study the asymptotic properties of the forecast mean-squared error when we fit a misspecified AR(k) model to the long-memory time series (X n ) n∈Z .

Rationale

Let Φ a k th degree polynomial defined by:

Φ(z) = 1 -a 1,k z -. . . -a k,k z k .
We assume that Φ has no zeroes on the unit disk. We define the process (η n ) n∈Z by:

∀n ∈ Z, η n = Φ(B)X n
where B is the backward shift operator. Note that (η n ) n∈Z is not a white noise series because (X n ) n∈Z is a long-memory process and hence does not belong to the class of autoregressive processes. Since Φ has no root on the unit disk, (X n ) n∈Z admits a moving-average representation as the fitted AR(k) model in terms of (η n ) n∈Z :

X n = ∞ j=0 c(j)η n-j .
If (X n ) n∈Z was an AR(k) associated with the polynomial Φ, the best next step linear predictor would be:

X n (1) = ∞ j=1 c(i)η t+1-i = a 1,k X n + . . . + a k,k X n+1-k si n k.
Here (X n ) n∈Z is a long-memory process which verifies the assumptions of Section 1. Our goal is to derive a closed formula for the polynomial Φ which minimizes the forecast error and to estimate this error.

Mean-Squared Error

There exists two approaches in order to define the coefficients of the k th degree polynomial Φ: the spectral approach and the time approach.

In the time approach, we choose to define the predictor as the projection mapping on to the closed span of the subset {X n , . . . , X n+1-k } of the Hilbert space L 2 (Ω, F, P) with inner product < X, Y >= E(XY ). Consequently the coefficients of Φ verify the equations, which are called the k th order Yule-Walker equations:

∀j ∈ 1, k , k i=1 a i,k σ(i -j) = σ(j) (16)
The mean-squared prediction error is:

E X n (1) -X n+1 2 = c(0) 2 E(η 2 n+1 ) = E(η 2 n+1 ).
We may write the moving average representation of (η n ) n∈N in terms of (ε n ) n∈N :

η n = ∞ j=0 min(j,p) k=0 Φ k b(j -k)ε n-j = ∞ j=0 t(j)ε n-j with ∀j ∈ N, t(j) = min(j,p) k=0 Φ k b(j -k).
Finally we obtain:

E X n (1) -X n+1 2 = ∞ j=0 t(j) 2 σ 2 ε .
In the spectral approach, minimizing the prediction error is equivalent to minimizing a contrast between two spectral densities:

π -π f (λ) g(λ, Φ) dλ
where f is the spectral density of X n and g(., Φ) is the spectral density of the AR(p) process defined by the polynomial Φ (see for example [Yajima, 1993]),so:

π -π f (λ) g(λ, Φ) dλ = π -π ∞ j=0 b(j)e -ijλ 2 Φ(e -iλ ) 2 dλ = π -π | ∞ j=0 t(j)e -ijλ | 2 dλ = 2π ∞ j=0 t(j) 2 .
In both approaches we need to minimize ∞ j=0 t(j).

Rate of Convergence of the Error by AR(k) Model Fitting

In the next theorem we derive an asymptotic expression for the prediction error by fitting autoregressive models to the series:

Theorem 3.3.1. Assume that (X n ) n∈Z is a long-memory process which verifies the assumptions of Section 1. If

0 < d < 1 2 : E X k (1) -X k+1 2 -σ 2 ε = O(k -1 )
Proof. Since fitting an AR(k) model minimizes the forecast error using k observations, the error by using truncation is bigger. Since the truncation method involves an error bounded by O k -1 , we obtain:

E X k (1) -X k+1 2 -σ 2 ε = O(k -1
). Consequently we only need to prove that this rate of convergence is attained . This is the case for the fractionally integrated processes defined in (12). We want the error made when fitting an AR(k) model in terms of the Wiener-Kolmogorov truncation error. Note first that the variance of the white noise series is equal to:

σ 2 ε = π -π f (λ) +∞ j=0 a j e ijλ 2 dλ.
Therefore in the case of a fractionally integrated process F(d) we need only show that:

π -π f (λ) +∞ j=0 a j e ijλ 2 dλ - σ 2 ε 2π π -π f (λ) g(λ, Φ k ) dλ ∼ C(k -1 ). π -π f (λ) +∞ j=0 a j e ijλ 2 dλ - σ 2 ε 2π π -π f (λ) g(λ, Φ k ) dλ = π -π f (λ)   +∞ j=0 a j e ijλ 2 - k j=0 a j,k e ijλ 2   dλ = +∞ j=0 +∞ l=0 (a j a l -a j,k a l,k ) σ(j -l) we set a j,k = 0 if j > k. +∞ j=0 +∞ l=0 (a j a l -a j,k a l,k ) σ(j -l) (17) = +∞ j=0 +∞ l=0 (a j a l -a j,k a l )σ(j -l) + +∞ j=0 +∞ l=0 (a j,k a l -a j,k a l,k )σ(j -l) = +∞ j=0 (a j -a j,k ) +∞ l=0 a l σ(l -j) + k j=0 a j,k +∞ l=0 (a l -a l,k )σ(j -l) (18)
We first study the first term of the sum (18). For any j > 0 , we have +∞ l=0 a l σ(l -j) = 0:

ε n = ∞ j=0 a l X n-l X n-j ε n = ∞ l=0 a l X n-l X n-j E (X n-j ε n ) = ∞ l=0 a l σ(l -j) E ∞ l=0 b l ε n-j-l ε n = ∞ l=0 a l σ(l -j)
and we conclude that +∞ l=0 a l σ(l -j) = 0 because (ε n ) n∈Z is an uncorrelated white noise. We can thus rewrite the first term of (18) like:

+∞ j=0 (a j -a j,k ) +∞ l=0 a l σ(l -j) = (a 0 -a 0,k ) +∞ l=0 a l σ(l) = 0
since a 0 = a 0,k = 1 according to definition. Next we study the second term of the sum (18):

k j=0 a j,k +∞ l=0 (a l -a l,k )σ(j -l).
And we obtain that:

k j=0 a j,k +∞ l=0 (a l -a l,k )σ(j -l) = k j=1 (a j,k -a j ) k l=1 (a l -a l,k )σ(j -l) + k j=1 (a j,k -a j ) +∞ l=k+1 a l σ(j -l) (19) + k j=0 a j k l=1 (a l -a l,k )σ(j -l) (20) + k j=0 a j +∞ l=k+1 a l σ(j -l)
Similarly we rewrite the term (19) using the Yule-Walker equations:

k j=1 (a j,k -a j ) +∞ l=k+1 a l σ(j -l) = - k j=1 (a j,k -a j ) k l=0 a l σ(j -l)
We then remark that this is equal to (20). Hence it follows that:

k j=0 a j,k +∞ l=0 (a l -a l,k )σ(j -l) = k j=1 (a j,k -a j ) k l=1 (a l -a l,k )σ(j -l) +2 k j=1 (a j,k -a j ) +∞ l=k+1 a l σ(j -l) + k j=0 a j +∞ l=k+1 a l σ(j -l) (21)
On a similar way we can rewrite the third term of the sum (21) using Fubini Theorem:

k j=0 a j +∞ l=k+1 a l σ(j -l) = - +∞ j=k+1 +∞ l=k+1 a j a l σ(j -l).
This third term is therefore equal to the forecast error in the method of prediction by truncation.

In order to compare the prediction error by truncating the Wiener-Kolmogorov predictor and by fitting an autoregressive model to a fractionally integrated process F(d), we need the sign of all the components of the sum (21). For a fractionally integrated noise, we know the explicit formula for a j and σ(j):

∀j > 0, a j = Γ(j -d) Γ(j + 1)Γ(-d) < 0 and ∀j ≥ 0, σ(j) = (-1) j Γ(1 -2d) Γ(j -d + 1)Γ(1 -j -d) σ 2 ε > 0.
In order to get the sign of a j,k -a j we use the explicit formula given in [START_REF] Brockwell | Simple consistent estimation of the coefficients of a linear filter[END_REF] and we easily obtain that a j,k -a j is negative for all j ∈ 1, k .

a j -a j,k = Γ(j -d) Γ(j + 1)Γ(-d) - Γ(k + 1)Γ(j -d)Γ(k -d -j + 1) Γ(k -j + 1)Γ(j + 1)Γ(-d)Γ(k -d + 1) = -a j -1 + Γ(k + 1)Γ(k -d -j + 1) Γ(k -j + 1)Γ(k -d + 1) = -a j k...(k -j + 1) (k -d)...(k -d -j + 1) -1 > 0 since ∀j ∈ N * a j < 0.
To give an asymptotic equivalent for the prediction error, we use the sum given in ( 21). We have the sign of the three terms: the first is negative, the second is positive and the last is negative. Moreover the third is equal to the forecast error by truncation and we have proved that this asymptotic equivalent has order O(k -1 ). The prediction error by fitting an autoregressive model converges faster to 0 than the error by truncation only if the second term is equivalent to Ck -1 , with C constant. Consequently, we search for a bound for a j -a j,k given the explicit formula for these coefficients (see for example [START_REF] Brockwell | Simple consistent estimation of the coefficients of a linear filter[END_REF]):

a j -a j,k = Γ(j -d) Γ(j + 1)Γ(-d) - Γ(k + 1)Γ(j -d)Γ(k -d -j + 1) Γ(k -j + 1)Γ(j + 1)Γ(-d)Γ(k -d + 1) = -a j -1 + Γ(k + 1)Γ(k -d -j + 1) Γ(k -j + 1)Γ(k -d + 1) = -a j k...(k -j + 1) (k -d)...(k -d -j + 1) -1 = -a j j-1 m=0 1 -l k 1 -l+d k -1 = -a j j-1 m=0 1 + d k 1 -d+l k -1 .
Then we use the following inequality:

∀x ∈ R, 1 + x ≤ exp(x)
which gives us:

a j -a j,k ≤ -a j exp j-1 m=0 d k 1 -d+l k -1 ≤ -a j exp d j-1 m=0 1 k -d -l -1 ≤ -a j exp d j-1 m=0 1 k -d -l
According to the previous inequality, we have:

k j=1 (a j -a j,k ) +∞ l=k+1 -a l σ(j -l) = k-1 j=1 (a j -a j,k ) +∞ l=k+1 -a l σ(j -l) +(a k -a k,k ) +∞ l=k+1 -a l σ(k -l) ≤ k-1 j=1 -a j exp d j-1 m=0 1 k -d -m +∞ l=k+1 -a l σ(j -l) +(-a k ) exp d k-1 m=0 1 k -d -m +∞ l=k+1 -a l σ(k -l) ≤ k-1 j=1 -a j exp d j 0 1 k -d -m dm +∞ l=k+1 -a l σ(j -l) +(-a k )k 3 2 d +∞ l=k+1 -a l σ(k -l)
As the function x → 1 k-d-x is increasing, we use the Integral Test Theorem. The inequality on the second term follows from:

k-1 m=0 1 k -d -m ∼ ln(k) ≤ 3 2 ln(k)
for k large enough. Therefore there exists K such that for all k ≥ K:

k j=1 (a j -a j,k ) +∞ l=k+1 -a l σ(j -l) ≤ k-1 j=1 -a j exp d ln k -d k -d -j +∞ l=k+1 -a l σ(j -l) +(-a k )k 3 2 d +∞ l=k+1 -a l σ(0) ≤ C(k -d) d k-1 j=1 j -d-1 (k -d -j) -d +∞ l=k+1 l -d-1 (l -j) 2d-1 +Ck -d-1 k 3 2 d k -d ≤ C (k -d) 2 1 1/(k-d) j -d-1 (1 -j) -d +∞ 1 l -d-1 (l -1) 2d-1 dldj +Ck -1 2 d-1 ≤ C ′ (k -d) -2+d + Ck -1 2 d-1
and so the positive term has a smaller asymptotic order than the forecast error made by truncating. Therefore we have proved that in the particular case of F(d) processes, the two prediction errors are equivalent to Ck -1 with C constant.

The two approaches to next-step prediction, by truncation to k terms or by fitting an autoregressive model AR(k) have consequently a prediction error with the same rate of convergence k -1 . So it is interesting to study how the second approach improves the prediction. The following quotient:

r(k) := k j=1 (a j,k -a j ) k l=1 (a l -a l,k )σ(j -l) + 2 k j=1 (a j,k -a j ) +∞ l=k+1 a l σ(j -l) k j=0 a j +∞ l=k+1 a l σ(j -l) (22)
is the ratio of the difference between the two prediction errors and the prediction error by truncating in the particular case of a fractionally integrated noise F(d). The figure 3.1 shows that the prediction by truncation incurs a larger performance loss when d → 1/2. The improvement reaches 50 per cent when d > 0.3 and k > 20.

After obtaining asymptotic equivalent for next step predictor, we will generalize the two methods of h-step prediction and aim to obtain their asymptotic behaviour as k → +∞ but also as h → +∞.

The h-Step Predictors

Since we assume that the process (X n ) n∈Z has an autoregressive representation (2) and moving average representation (1), the linear least-squares predictor, X k+h , of X k+h based on the infinite past (X j , j ≤ k) is given by:

X k (h) = - +∞ j=1 a j X k (h -j) = +∞ j=h b j ε k+h-j
(see for example Theorem 5.5.1 of [START_REF] Brockwell | Time Series : Theory and Methods[END_REF]). The corresponding mean squared error of prediction is:

E X k (h) -X k+h 2 = σ 2 ε h-1 j=0 b 2 j .
As the prediction step h tends to infinity, The mean-squared prediction error converges to σ 2 ε +∞ j=0 , which is the the variance of the process (X n ) n∈Z . But if the mean-squared prediction error is equal to σ(0), we have no more interest in the prediction method since its error is equal to the error of predicting the future by 0. Remark that the mean-squared error increases more slowly to σ(0) in the long-memory case than in the short-memory case since the sequence b j decays more slowly to 0. More precisely in the case of a long-memory process, if we assume that:

b j ∼ +∞ j d-1 L(j)
where L is a slowly varying function, we can express the asymptotic behaviour of the prediction error. As j → L 2 (j) is also a slowly varying function according to the definition of [Zygmund, 1968], b 2 j = j 2d-2 L 2 (j) is ultimately decreasing. The rest of the series and the integral are then equivalent and we may write:

σ(0) -E X k (h) -X k+h 2 = σ 2 ε +∞ j=h b 2 j ∼ +∞ j=h j 2d-2 L 2 (j) ∼ +∞ h j 2d-2 L 2 (j)dj
According to Proposition 1.5.10 of [START_REF] Bingham | Regular variation[END_REF]:

σ(0) -E X k (h) -X k+h 2 ∼ +∞ h j 2d-2 L 2 (j)dj ∼ h→+∞ 1 1 -2d h 2d-1 L 2 (h) (23)
In the case of a long-memory process with parameter d which verifies b j ∼ j d-1 L(j), the convergence of the mean-squared error to σ(0) is slow as h tends to infinity. On the contrary, for a moving average process of order q, the sequence σ(0

) -E X k (h) -X k+h 2
is constant and equal to 0 as soon as h > q. More generally, we can study the case of an ARMA process, which canonical representation is given by: Φ(X t ) = Θ(ε t )

where Φ and Θ are two coprime polynomials with coefficients of degree 0 are equal to 1 and ε t is a white noise. Φ has no root in the unit disk |z| ≤ 1 and Θ has no root in the open disk |z| < 1. b j is bounded by:

|b j | ≤ Cj m-1 ρ -j
where ρ is the smallest absolute value of the roots of Φ and m the multiplicity of the corresponding root (see for example [START_REF] Brockwell | Time Series : Theory and Methods[END_REF] p92). Thus the mean-squared prediction error is bounded by:

σ(0) -E X k (h) -X k+h 2 = σ 2 ε +∞ j=h b 2 j ≤ σ 2 ε C 2 +∞ j=h j 2m-2 ρ -2j ≤ σ 2 ε C 2 +∞ j=h j 2m-2 exp (-2j log(ρ)) ≤ σ 2 ε C 2 +∞ h j 2m-2 exp (-2j log(ρ)) dj
By using the substitution t = 2 log(ρ)j ,

σ(0) -E X k (h) -X k+h 2 ≤ σ 2 ε C 2 (2 log(ρ)) 1-2m +∞ 2 log(ρ)h t 2m-2 exp (t) dt ≤ σ 2 ε C 2 (2 log(ρ)) 1-2m Γ(2m -1, 2 log(ρ)h)
where Γ(., .) is the incomplete Gamma function defined in equation 6.5.3 of [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF].

We know an equivalent of this function:

Γ(2m -1, 2 log(ρ)h) ∼ h→+∞ (2 log(ρ)h) 2m-2 exp (2 log(ρ)h)
We conclude that the rate of convergence is exponential. The mean-squared prediction error goes faster to σ(0) when the predicting process is ARMA than when the process is a long-memory process.

The h-step prediction is then more interesting for the long-memory process than for the shortmemory process, having observed the infinite past. We consider the truncating effect next.

Truncated Wiener-Kolmogorov predictor

In practice, we only observe a finite number of samples. We assume now that we only know k observations (X 1 , . . . , X k ). We then define the h-step truncated Wiener-Kolmogorov of order k as:

X ′ k (h) = - h-1 j=1 a j X ′ k (h -j) - k l=1 a h-1+j X k+1-j (24)
We now describe the asymptotic behaviour of the mean-squared error of the predictor (24). First we write the difference between the predicting random variable and its predictor:

X ′ k (h) -X k+h = - h-1 j=1 a j X ′ k (h -j) - k l=1 a h-1+j X k+1-j -ε k+h + +∞ j=1 a j X k+h-j = -ε k+h + h-1 j=1 a j X k+h-j -X ′ k (h -j) + k j=1 a h-1+j (X k+1-j -X k+1-j ) + +∞ j=k+1 a h-1+j X k+1-j = -ε k+h + h-1 j=1 a j X k+h-j -X ′ k (h -j) + +∞ j=k+1 a h-1+j X k+1-j
We will use the process of induction on h to show that

X ′ k (h) -X k+h = - h-1 l=0   j 1 +j 2 +...+j h =l (-1) card({j,j =0}) a j 1 a j 2 . . . a j h   ε k+h-l + +∞ j=k+1   i 1 +i 2 +...+i h =h-1 (-1) card({i l ,i l =0,l>1}) a j+i 1 a i 2 . . . a i h   X k+1-j .
For h = 2, we have for example

X ′ k (2) -X k+2 = -(a 0 ε k+2 -a 1 ε k+1 ) + +∞ j=k+1
(-a 1 a j + a j+1 )X k+1-j .

Let A(z) and B(z) denote A(z) = 1 + +∞ j=1 a j z j and B(z) = 1 + +∞ j=1 b j z j . Since we have A(z) = B(z) -1 , we obtain the following conditions on the coefficients:

b 1 = -a 1 b 2 = -a 2 + a 2 1 b 3 = -a 3 + 2a 1 a 2 -a 3 1 . . .
So we obtain:

X ′ k (h) -X k+h = - h-1 l=0 b l ε k+h-l + +∞ j=k+1 h-1 m=0 a j+m b h-1-m X k+1-j . (25)
Since the process (ε n ) n∈Z is uncorrelated and then the two terms of the sum ( 25) are orthogonal, we can rewrite the mean-squared error:

E X ′ k (h) -X k+h 2 = h-1 l=0 b 2 l σ 2 ε (26) +E   +∞ j=k+1 h-1 m=0 a j+h-1-m b m X k+1-j   2 . ( 27
)
The first part of the error ( 26) is due to the prediction method and the second (27) due to the truncating of the predictor. We now approximate the error term (27) by using (3) and (4). We obtain the following upper bound:

∀δ > 0, h-1 m=0 a j+h-1-m b m ≤ h-1 m=1 |a j+h-1-m b m | + |b 0 a j+h-1 | ≤ C 1 C 2 h 0 (j + h -1 -l) -d-1+δ l d-1+δ dl + C 1 (j + h) -d-1 ≤ C 1 C 2 h -1+2δ 1 0 j h + 1 -l -d-1+δ l d-1+δ dl + C 1 (j + h) -d-1 ≤ C 1 C 2 h -1+2δ j -d-1+δ 1 0 1 h + 1 -l j -d-1+δ l d-1+δ dl + C 1 (j + h) -d-1 ≤ C 1 C 2 h d+2δ j -d-1+δ 1 0 l d-1+δ dl + C 1 (j + h) -d-1 h m=0 a j+h-m b m ≤ C 1 C 2 h d+2δ d j -d-1+δ (28)
This bound is in fact an asymptotic equivalent for the fractionally integrated process F(d) because, in that case, the sequences a j and b j have a constant signs. Using Proposition 2.2.1 for the one-step prediction and we have:

Proposition 4.1.1. Let (X n ) n∈Z be a linear stationary process defined by (1), (2) and possessing the features (3) and (4). We can approximate the mean-squared prediction error of X ′ k (1) by:

∀δ > 0, E X ′ k (h) -X k+h 2 = h-1 l=0 b 2 l σ 2 ε + O h 2d+δ k -1+δ . ( 29
)
Having k observations, we search for the step h for which the variance of the predictor has for upper bound σ(0). Then the prediction error have for asymptotic bound O h 2d k -1 . We want to choose h to have the prediction error negligible with respect to the information given by the linear least-squares predictor given the infinite past (see ( 23)) and we obtain:

h 2d k -1 = o(h 2d-1 )
and then h = o(k). With the truncated Wiener-Kolmogorov predictor, it is interesting to compute the h-step predictor if we have k observations h = o(k).

The k-th Order Linear Least-Squares Predictor

For next step predictor, when we fitted an autoregressive process, we search the linear least-squares predictor knowing the finite past (X 1 , . . . , X k ) and the predictor is then the projection of the random variable onto the past. Let X k (h) denote the projection of X k+h onto the span of (X 1 , . . . , X k ). X k (h) verifies the recurrence relationship

X k (h) = - k j=1 a j,k X k (h -j)
where X k (h -j) is the direct linear least-squares predictor of X k+h-j based on the finite past (X 1 , . . . , X k ). By induction, we obtain the predictor as a function of (X 1 , . . . , X k ): For next step prediction by fitting an autoregressive process, the best linear least-squares predictor knowing the finite past is a projection of the random variable X k+1 onto the past.

X k (h) = - k j=1 c j,k X k+1-j .
Since X k (h) is the projection of X k+h onto (X 1 , . . . , X k ) in L 2 , the vector (c j,k ) 1≤j≤k minimizes the mean-squared error:

E X k (h) -X k+h 2 = π -π f (λ) exp(iλ(h -1)) + k j=1 c j,k exp(-iλj) 2 dλ
The vector (c j,k ) 1≤j≤k is a solution of the equation:

∇ c E X k (h) -X k+h 2 = 0

Figure 3 . 1 :

 31 Figure 3.1: Ratio r(k), d ∈]0, 1/2[ defined in (22)

Figure 4

 4 Figure 4.1: Mean-squared error of X k (h) (MMSE), X ′ k (h) (TPMSE) and X k (h) (LLSPE) for d = 0.4 and k = 80

where ∇ c is the gradient. The vector (c j,k ) 1≤j≤k is then equal to:

The corresponding mean squared error of prediction is given by:

The matrix Σ -1 k is symmetric positive definite and the prediction error of this method is always lower than σ(0).

As X k (h) is the projection of X k+h onto (X 1 , . . . , X k ), the mean-squared prediction error is also lower than the prediction error of the truncated Wiener-Kolmogorov predictor (see figure 4.1). The mean-squared error of prediction due to the projection onto the span of (X 1 , . . . , X k ) tends at least as fast to zero as the mean-squared due to truncation of the least-squares predictor. For one-step predictor,we have shown that the two methods can have the same rate of convergence.


Bruno Torr Ésani

SPREADING FUNCTION REPRESENTATION OF OPERATORS AND GABOR MULTIPLIER APPROXIMATION

Modification of signals in the time-frequency domain are used in many applications. However, the modification is often restricted to be purely multiplicative. In this paper, it is shown that, in the continuous case, a quite general class of operators can be represented by a twisted convolution in the short-time Fourier transform domain. The discrete case of Gabor transforms turns out to be more intricate. A similar representation will however be derived by means of a special form for the operator's spreading function (twisted spline type function). The connection between STFT-and Gabor-multipliers, their spreading function and the twisted convolution representation will be investigated. A precise characterization of the best approximation and its existence is given for both cases. Finally, the concept of Gabor multipliers is generalized to better approximate "overspread" operators.

Introduction

Time-frequency expansions have been shown to provide extremely efficient representations for functions and distributions, and found many relevant applications, notably in signal processing. Efficient representation of linear operators in the time-frequency domain still remains a difficult challenge. Many theoretical results have been obtained, but these have only had little practical impact so far. We here address this problem from a slightly different point of view, emphasizing the geometrical aspects of short time Fourier and Gabor transforms.

It is well known that twisted convolution plays a central role as soon as functions and operators are studied directly on the phase space. We first show that the spreading function representation of Hilbert-Schmidt operators is equivalent to a representation by left twisted convolution in the STFT domain. Unfortunately, such a remark does not hold true any more when the STFT is discretized, i.e. replaced with a Gabor transform.

Gabor multipliers are often proposed as an efficient alternative for time-frequency operator representation. It is well known that for "nice operators", i.e. operators that do not involve time-frequency shifts of large magnitude, Gabor multipliers can provide accurate approximations. We show that for more general operators suitable generalizations of Gabor multipliers may still provide good approximation. We mainly propose two possible extensions. In the first one, the analyzed operator is approximated by a linear combination of adapted Gabor multipliers, and when the latter are given on a suitably chosen lattice in the time-frequency plane, the question of optimal approximation (in Hilbert-Schmidt sense) may be well formulated.

In the second extension, an adaptation of the synthesis window to the analyzed operator is proposed, while the time-frequency transfer function is fixed. In both situations, the final result reveals interesting connections to twisted convolution.

The Weyl-Heisenberg group, Time-frequency representations

Classical time-frequency analysis is based on translation and modulation operators, which generate the Weyl-Heisenberg group. Even though the group theoretical language may be completely avoided in most applications and theoretical developments, it is nevertheless interesting to trace back some tools which play a central role in Gabor theory (including covariance issues, and twisted convolution) to their algebraic origins.

2.1.

Preliminaries: the Weyl-Heisenberg group and its unitary representations. We start by briefly reviewing the main results concerning the Weyl-Heisenberg group and its representation theory that will be of interest here. We refer the interested reader to [START_REF] Führ | Abstract harmonic analysis of continuous wavelet transforms[END_REF][START_REF] Schempp | Harmonic analysis on the Heisenberg nilpotent Lie group[END_REF] for a more thorough analysis.

Classical time-frequency analysis rests on translation operators T b and modulation operators M ν , defined by

T b f (t) = f (t -b) , M ν f (t) = e 2iπνt f (t) , f ∈ L 2 (R) ,
and the corresponding time-frequency shifts π(b, ν) = M ν T b , which satisfy the canonical commutation relations

M ν T b = e 2iπνb T b M ν .
These operators generate a three-parameter Lie group, called the (reduced) Weyl-Heisenberg group 1 H.

(1)

H = {(b, ν, ϕ) ∈ R × R × [0, 1]} , with group multiplication (2) (b, ν, ϕ)(b ′ , ν ′ , ϕ ′ ) = (b + b ′ , ν + ν ′ , ϕ + ϕ ′ -ν ′ b) .
The zero element is clearly e = (0, 0, 0), and the inverse operation is (b, ν, ϕ) -1 = (-b, -ν, -ϕ -bν). The group H is unimodular, i.e. its two Haar measures coincide, and read

(3) dµ(b, ν, ϕ) = db dν dϕ .

The left regular representation of H on L 2 (H) takes the following form:

(4) L(b ′ , ν ′ , ϕ ′ )F (b, ν, ϕ) = F (b ′ , ν ′ , ϕ ′ ) -1 (b, ν, ϕ) = F (b-b ′ , ν-ν ′ , ϕ-ϕ ′ +b ′ (ν-ν ′ )) .
The unitary dual H of H (i.e. the set of unitary equivalence classes of irreducible unitary representations of H) has the following form

H ≃ R 2 ∪ Z * , with Z * = Z\{0}. The R 2 part consists of the unitary characters χ α,β , α, β ∈ R χ α,β (b, ν, ϕ) = e 2iπαb e 2iπβν , while the Z * part consists in a series of representations π λ on H λ ≃ L 2 (R d ): (5) π λ (b, ν, ϕ)f (t) = e 2iπλ[ϕ+νt] f (t -b) .
In other terms, π λ (b, ν, ϕ) = e 2iπλϕ M λν T b ,

1 Notice that there exist other versions of the Weyl-Heisenberg group, including a more symmetric one generated by time-frequency shifts M ν/2 T b M ν/2 , and another one often used in signal processing, generated by time-frequency shifts T b Mν .

The representations π λ , which form the so-called discrete series, may be shown to be square-integrable, i.e. such that there exist g ∈ H λ such that

0 < H | π λ (h)g, g | 2 dµ(h) < ∞ .
It follows (see [START_REF] Grossmann | Transforms associated to square integrable group representations I: General results[END_REF][START_REF] Grossmann | Transforms associated to square integrable group representations II: Examples[END_REF] and references therein) that the integral transform V g , defined by

(6) f ∈ H λ -→ V g f = 1 g f, π λ (•)g ∈ L 2 (H) maps isometrically H λ to L 2 (H).
2.2. More details on H and V g . The dual H of the Weyl-Heisenberg group H and the integral transform V g turn out to be closely connected, as we outline next. Classical L 2 Fourier series theory yields the following Lemma 1. For F ∈ L 2 (H), the mapping F → {F ℓ , ℓ ∈ Z} defined by

(7) F ℓ (b, ν, ϕ) = e 2iπℓϕ F ℓ (b, ν, 0) = e 2iπℓϕ 1 0
F (b, ν, α)e -2iπℓα dα establishes a bijective isometry

L 2 (H) ≃ ℓ∈Z + E ℓ where (8) E ℓ = F ∈ L 2 (H), F (b, ν, ϕ) = e 2iπℓϕ F ℓ (b, ν, 0) ≃ L 2 (R 2 ) .
Hence, we have for all

F ∈ L 2 (H) (9) F (b, ν, ϕ) = ∞ ℓ=-∞ F ℓ (b, ν, ϕ) = ∞ ℓ=-∞ F ℓ (b, ν, 0)e 2iπℓϕ and G |F (b, ν, ϕ)| 2 dµ L (b, ν, ϕ) = ∞ ℓ=-∞ F ℓ 2 L 2 (R 2 ) .
In addition, we immediately get Lemma 2. Let π λ be an irreducible unitary representation in the discrete series of H, g ∈ H and let V g be defined as above in [START_REF] Gröchenig | Foundations of Time-Frequency Analysis[END_REF]. Then V g x ∈ E -λ .

From now on, we shall limit ourselves to the case λ = 1, i.e. V g x ∈ E -1 for all x ∈ L 2 (R). We denote by π = π 1 the corresponding irreducible unitary representation of H on H -1 , and identify the latter with L 2 (R).

. Let F, G ∈ E -1 . Then F * G(b, ν, ϕ) = F (b ′ , ν ′ , ϕ ′ )G(b -b ′ , ν -ν ′ , ϕ -ϕ ′ + b ′ (ν -ν ′ ))db ′ dν ′ dϕ ′ = e -2iπ[ϕ ′ +b ′ (ν-ν ′ )] R 2 F (b ′ , ν ′ , 0)G(b -b ′ , ν -ν ′ , 0)e -2iπb ′ (ν-ν ′ ) db ′ dν ′ = e -2iπϕ ′ F * G(b, ν, 0) .
Therefore, F * G ∈ E -1 . In addition, identifying E -1 with L2 (R 2 ) yields the following twisted convolution product on L 2 (R 2 ) (see [START_REF] Folland | Harmonic Analysis in Phase Space[END_REF] for more details)

(10) (F ♮G)(b, ν) = R 2 F (b ′ , ν ′ )G(b -b ′ , ν -ν ′ )e -2iπb ′ (ν-ν ′ ) db ′ dν ′ .
The Plancherel transform may be seen as a generalization of usual Fourier transform to more general groups. In the case of the Heisenberg group, the Plancherel transform associates with any function F ∈ L 2 (H) the operator valued function on H defined by [START_REF] Schempp | Harmonic analysis on the Heisenberg nilpotent Lie group[END_REF] F :

σ ∈ H -→ F (σ) = H f (h)σ(h) dµ(h) ,
the equality being understood in the weak operator sense. In the particular case σ = π, identifying again E -1 with L 2 (R 2 ) one may associate with any F ∈ L 2 (R 2 ) the operator valued function F (π) defined as

(12) F (π)f = ∞ -∞ ∞ -∞ F (b, ν)M ν T b f dbdν , f ∈ L 2 (R)
2.4. Connection to short time Fourier transform. Given the representation of the Weyl-Heisenberg group π = π 1 (see Eq. ( 5)), and g ∈ L 2 (R), g = 1, define the short time Fourier transform (STFT, also called continuous Gabor transform, CGT)

V g x ∈ L 2 (R 2 ) of x ∈ L 2 (R) by (13) V g x(b, ν) = x, π(b, ν, 0)g = x, M ν T b g ,
where by an abuse of notation we have identified

E -1 with L 2 (R 2 ). 2 Then it is a well known fact that the image of L 2 (R 2 ) by V g is a proper subspace K g of L 2 (R 2 )
, consisting of solutions of the kernel equation

K g = F ∈ L 2 (R 2 ), F (b, ν) = F (b ′ , ν ′ ) π(b ′ , ν ′ )ψ, π(b, ν)ψ db ′ dν ′ , ∀b, ν ∈ R .
The following result is proved by direct calculations. We denote by H the set of Hilbert-Schmidt operators on L 2 (R).

Proposition 1. The short time Fourier transform intertwines twisted convolution and action of Hilbert-Schmidt operators on L 2 (R). More precisely, let F ∈ L 2 (R 2 ).

(1) There exists a linear operator

H F = F (π 1 ) ∈ H such that (14) F ♮V g x = V g H F x . Therefore, F ♮V g x ∈ K g . (2)
The linear operator H F ∈ H also satisfies

V g x♮F = V H * F g x . Remark 1.
The consequence of this result is the fact that the image of the STFT of a signal x ∈ L 2 (R d ) by a twisted convolution with some function of two variables is itself the STFT of some signal H F x: K g is invariant under left twisted convolution. Similarly, right twisted convolution of V g x with some function yields another STFT of x with respect to a different window: right twisted convolution maps the reproducing kernel space K g to the reproducing kernel space K H * F g .

The Gabor transform (see for example [START_REF] Gröchenig | Foundations of Time-Frequency Analysis[END_REF] for a review) is defined as a sampled version of the STFT. Given two constants a 0 , ν 0 ∈ R + , the corresponding Gabor transform associates with any f ∈ L 2 (R) the sequence of Gabor coefficients (15)

V g f (mb 0 , nν 0 ) = f, M nν0 T mb0 g = f, g mn ,
the functions g mn = M nν0 T mb0 g being the Gabor atoms associated with g and the lattice constants b 0 , ν 0 . Under suitable assumptions, the Gabor transform is left invertible, and there exists h ∈ L 2 (R) such that any f may be expanded as

(16) f = m,n V g f (mb 0 , nν 0 )h mn .
Proposition 1 does not have any simple analogue if STFT is replaced with Gabor transform. However, connections between Gabor representations of Hilbert-Schmidt operators and twisted convolutions will appear below.

Spreading function and twisted convolution operator representation

3.1. The continuous case. Let H ∈ H , the class of Hilbert-Schmidt operator on L 2 (R). Then there exists a function η = η H ∈ L 2 (R 2 ), called the spreading function (see for example [START_REF] Gröchenig | Foundations of Time-Frequency Analysis[END_REF] and references therein), such that 3

(17) H = ∞ -∞ ∞ -∞ η(b, ν)π(b, ν) dbdν .
Such a representation turns out to be closely related to the twisted convolution, defined in (10) above. Our first remark is the fact that the spreading function representation of the operator H ∈ H actually takes the form of a twisted convolution in the (continuous) time-frequency domain. η be its spreading function. Indeed, the following result [START_REF] Dörfler | On the time-frequency representation of operators and generalized Gabor multiplier approximations[END_REF] is a direct consequence of Proposition 1, that follows from the comparison of the spreading function representation with the Plancherel transform [START_REF] Strohmer | Pseudodifferential operators and Banach algebras in mobile communications[END_REF]. Theorem 1. Let H ∈ H , with spreading function η H . Let g, h ∈ L 2 (R) be such that g, h = 1. Then H may be realized as a left twisted convolution in the timefrequency domain: for all f ∈ L 2 (R),

(18) Hf = ∞ -∞ ∞ -∞ (η H ♮V g f ) (b, ν)M ν T b h dbdν .
An immediate consequence is the fact, that the range of V g is invariant under left twisted convolution.

Remark 2. Again, let us notice that no simple analogue of this expression is available when the STFT is replaced with the Gabor transform. Replacing the integrals with discrete sums does not yield simple operators. 3 Here, this representation must be interpreted in the weak sense. More general, the spreading representation can be given an interpretation for η ∈ S ′ , i.e. , if the spreading function is a tempered distribution. For more regular spreading functions, e.g., η ∈ L 1 , the integral is absolutely convergent in operator norm.

3.2.

The discrete Heisenberg group. Although the focus of the present article is not on discrete groups, let us shortly mention that the above discussion applies almost directly to the case of the discrete Weyl-Heisenberg group

H N = Z N × Z N × Z N , with group law (m, n, ϕ)(m ′ , n ′ , ϕ ′ ) = (m + m ′ , n + n ′ , ϕ + ϕ ′ -n ′ b) ,
all operations being understood modulo N .

The unitary dual H N has a structure similar to that of H:

H N = Z 2 N ∪ Z * N and the Z * N component of the dual consists in representations π λ on H λ = C N , given by π λ (m, n, ϕ) = e 2iπλϕ/N M n T m ,
with the following definition for periodic discrete translation and modulation operators

T m x(k) = x((k -m)[mod N ]) , M n x(k) = e 2iπkn/N x(k) .
Given g ∈ C N , the corresponding short time Fourier transform maps any

x ∈ C N to V g x ∈ C 2N , defined by V g x(m, n) = x, π λ (m, n, 0)g .
Focusing on the case λ = 1, and following the lines developed in the continuous case, one is naturally led to the following discrete version of the twisted convolution: for

F, G ∈ C N , ( 19
) F ♮G(m, n) = N -1 m ′ ,n ′ =0 F (m ′ , n ′ )G(m -m ′ , n -n ′ )e -2iπm(n-n ′ )/N .
The spreading function representation of operators takes here a particularly simple form: the family of operators

{M n T m , m, n ∈ Z N } is an orthonormal basis of Hom(C N ). Every H ∈ Hom(C N ) is characterized by a spreading function η H ∈ C N 2 , such that H = N -1 m,n=0 η H (m, n)M n T m .
The consequence of this result is the same as in the continuous theory: let H ∈ Hom(C N ), and let η H be its spreading function. Let ψ ∈ C N be a unit norm vector. Then H may be realized as a twisted convolution in the time-frequency domain: for all x ∈ C N ,

Hx = m,n (η H ♮V ψ x) (m, n)M n T n ψ .
Actually, a similar theory may be developed for Weyl-Heisenberg groups associated with any locally compact abelian group. Examples. As a first example, a simple transposition has been realized by a twisted convolution with a "dirac pulse" in the TF-domain. See Figure 1. The reconstruction is performed with the dual of the original Gabor frame 4 and yields a perfect result. This twisted convolution corresponds, of course, to a pure timeor frequency shift. However, adding several dirac pulses, i.e. setting the spreading function to 1 in several places, will yield an amalgam of the signal components present in the original signal.

Figure 2 shows a nicely concentrated spreading function and the images of a random signal under the corresponding operator. The second pair of images shows the spreading function of a convolution operator and again the result of this operator being applied to a random signal. Both operators have a strong smoothing character, however, the convolution performs infinite time-shifts. 4 If synthesis is desired, we need to assume that the Gabor system is a frame, meaning that there exist constants 0 < A, B < ∞, such that (20)

A f 2 ≤ m n | f, gmn | 2 ≤ B f 2 .

STFT and Gabor multipliers

Twisted convolution operator representations such as the ones described above are not suitable for practical implementations, mainly because of the difficulty of sampling such representations, as stressed in Remark 2. In addition, even if such twisted convolutions could be treated practically, they would hardly lead to efficient algorithms, as a twisted convolution involves quite a large number of terms.

STFT and Gabor multipliers have been proposed as a valuable alternative for the approximation of operators in the time-frequency domain. In the present section we analyse the relationship between such multipliers and the spreading function representations, referring to [START_REF] Feichtinger | A first survey of Gabor multipliers[END_REF] for a survey. 4.1. Time-frequency multipliers. Let g, h ∈ L 2 (R 2 ) be such that g, h = 1, and m ∈ L ∞ (R 2 ), and define the STFT multiplier M m;g,h by

(21) M m;g,h f = ∞ -∞ ∞ -∞ m(b, ν)V g f (b, ν) π(b, ν)h dbdν, f ∈ L 2 (R d ) .
This clearly defines a bounded operator on L 2 (R d ). We shall call m the timefrequency transfer function, in analogy with signal processing notation, where timeinvariant linear filters (convolution operator) are defined by multiplication in the Fourier domain with a bounded function called transfer function (which is the Fourier transform of the filter's impulse response). We refer to [START_REF] Hlawatsch | Linear time-frequency filters[END_REF] and references therein for a presentation in a signal processing context. In the mathematical literature, m is also called the symbol of the multiplier, but we prefer to reserve this terminology for the Kohn-Nirenberg symbol or the Weyl symbol of operators (see [START_REF] Folland | Harmonic Analysis in Phase Space[END_REF][START_REF] Gröchenig | Foundations of Time-Frequency Analysis[END_REF] for example). m is sometimes also called a mask.

STFT multipliers may also be extended to various functional settings. For example, suitable assumptions on the transfer function m ensure that the corresponding multiplier is compact. We refer to [START_REF] Feichtinger | A first survey of Gabor multipliers[END_REF] for a thorough review of these results.

In a similar way, given lattice constants b 0 , ν 0 ∈ R + , set π mn = π(mb 0 , nν 0 ) = M nν0 T mb0 . Then, for m ∈ ℓ ∞ (Z 2 ), the corresponding Gabor multiplier is defined as

(22) M G m;g,h f = ∞ -∞ ∞ -∞ m(m, n)V g f (mb 0 , nν 0 ) π mn h .
In both cases, the spreading function may be computed explicitely, and yields the following result [START_REF] Dörfler | On the time-frequency representation of operators and generalized Gabor multiplier approximations[END_REF].

Lemma 3.

(1) The spreading function of the STFT multiplier M m;g,h is given by

(23) η Mm (b, ν) = M (b, ν)V g h(b, ν) ,
where M is the symplectic Fourier transform of the transfer function m

M (t, ξ) = ∞ -∞ ∞ -∞
m(b, ν)e 2iπ(νt-ξb) dbdν .

(2) The spreading function of the Gabor multiplier M G m;g,h is given by d) is the symplectic Fourier transform of the transfer function m

(24) η G Mm (b, ν) = M (d) (b, ν)V g h(b, ν) , where the (ν 0 -1 , b 0 -1 )-periodic function M (
M (d) (t, ξ) = ∞ m=-∞ ∞ n=-∞ m(m, n)e 2iπ(nν0t-mb0ξ) .
Therefore, we can see that the spreading function of a STFT or a Gabor multiplier has a very specific form. In the STFT case, the following result is easily proved. In a nutshell, when the spreading function of H decays faster than the ambiguity function, H admits a STFT multiplier representation. However, when the decay is not fast, the transfer function m may lack smoothness. A simple way of controlling such a property is to assume compact support for the spreading function [START_REF] Hlawatsch | Linear time-frequency filters[END_REF][START_REF] Kozek | Matched Weyl-Heisenberg Expansions of Nonstationary Environments[END_REF].

In the case of Gabor multiplier approximations, the periodicity of the function M (d) imposes strong constraints, which we analyze in more detail below. 4.2. Further Examples. We now give some examples for illustration of the relation between Gabor Multipliers as well as more general operators and their spreading functions. Figure 3 shows the STFT of a synthezised signal with two components, the spreading function of a Gabor multiplier, which removes one of the components, and the STFT of the resulting signal. Figure 4 shows the spreading function of an operator that again filters out one of the components, but simultaneously moves it in the TF-plane, as well as, again, the STFT of the resulting signal. Note that, while still being well-concentrated, this operator's spreading function is no longer centered about the origin. Finally, Figure 5 shows the spreading function of an operator that filters the same component as before, but "duplicates" it: once, it is shifted in time and once in frequency. Quite obviously, the spreading function has two peaks in this case. The two latter operators will not be well-approximated by a simple Gabor multiplier and require a more general treatment, as suggested in Section 5.1.

4.3.

Approximation by Gabor multipliers. It follows from (24) that not any η ∈ L 2 (R 2 ) can be the spreading function of a STFT or Gabor multiplier. In particular, for Gabor multipliers, the function η/V g h has to be periodic, which is generally not true. In such cases, the optimal (in Hilbert-Schmidt sense) approximation may be seeked, provided that the family of rank one operators

P mb0,nν0 : f ∈ L 2 (R) → P mb0,nν0 f = f, g mb0,nν0 h mb0,nν0 , m, n ∈ Z
is a Riesz sequence. This question has been addressed by several authors (see [START_REF] Feichtinger | A first survey of Gabor multipliers[END_REF] and references therein). The answer turns out to be positive [START_REF] Dörfler | On the time-frequency representation of operators and generalized Gabor multiplier approximations[END_REF] if and only if there exist real constants 0 < A ≤ B < ∞ such that for all (t, ξ)

(25) 0 < A ≤ ∞ k,ℓ=-∞ V g h t + k ν 0 , ξ + ℓ b 0 2 ≤ B < ∞ .
The latter condition is obtained from the more classical one found in [START_REF] Feichtinger | A first survey of Gabor multipliers[END_REF] by applying the Poisson summation formula. Now, if the above condition is fulfilled, the best Gabor multiplier approximation (in Hilbert-Schmidt sense) of H ∈ H with spreading function η H is defined by the time-frequency transfer function m whose discrete symplectic Fourier transform reads

(26) M (b, ν) = ∞ k,ℓ=-∞ V g h (b + k/ν 0 , ν + ℓ/b 0 ) η H (b + k/ν 0 , ν + ℓ/b 0 ) ∞ k,ℓ=-∞ |V g h (b + k/ν 0 , ν + ℓ/b 0 )| 2
Therefore, assuming that the windows g, h and the lattice constants b 0 , ν 0 have been chosen so that condition (25) holds, the above expression may be used to compute the best Gabor multiplier approximation for any linear operator H ∈ H . For numerical realisation of approximation by Gabor multipliers, see, for example, [START_REF] Feichtinger | Approximation of matrices by Gabor multipliers[END_REF].

Generalized Gabor Multipliers

Gabor multipliers have the advantage of having a diagonal matrix in a suitable Gabor representation. Unfortunately, approximating a Hilbert-Schmidt operator by a Gabor multiplier may turn out to result in poor approximations when the spreading function of the operator is not "concentrated" enough, i.e. when the operators are not underspread [START_REF] Hlawatsch | Linear time-frequency filters[END_REF][START_REF] Kozek | Matched Weyl-Heisenberg Expansions of Nonstationary Environments[END_REF].

Two strategies may be developped in such a situation. A first strategy amounts to "split" the operator into a finite sum of operators, whose spreading function is much more concentrated in the time-frequency domain, and may thus be conveniently approximated by a Gabor multiplier. When the construction preserves the action of the Weyl-Heisenberg group, the corresponding optimal time-frequency transfer functions may be given a closed form. This leads to multiple Gabor multipliers, whose contruction is briefly outlined below.

The second possible strategy rests on a spline type approximation of the spreading function, expressed as a sum of building blocks with smaller support. Again, when this is done in a H-covariant way, and if the building blocks are spreading functions of Gabor multipliers, the resulting approximation turns out to take an A simple instance of this construction is described below, and we refer again to [START_REF] Dörfler | On the time-frequency representation of operators and generalized Gabor multiplier approximations[END_REF] for a thorough analysis of the general case. 5.1. Multiple Gabor multipliers. The strategy to restrict the modification of the Gabor transform to be multiplicative (i.e. the operator is represented by a diagonal matrix in the Gabor system) may be generalized to allow for banded matrices in this modification. This is the equivalent to allowing generalized Gabor multipliers, in the following sense.

Definition 1 (Multiple Gabor Multiplier). Let a lattice Λ = b 0 Z × ν 0 Z, g, h k ∈ L 2 (R 2
) and bounded sequences m k = (m k (λ)) λ∈Λ , with k ∈ Z, be given. Then the multiple Gabor multiplier (MGM) G g,h k ,Λ,m k is defined by

(27) G g,h k ,Λ,m k (f ) = k λ∈Λ m k (λ) f, π(λ)g π(λ)h k . Remark 3. For h k = π(µ k )h with µ k ∈ Λ, |µ k | ≤ K,
this extends Gabor multipliers by replacing their diagonal matrix by a banded one. This is a generalization of ample relevance in applications, for example as prominent ones as mobile communication, see [START_REF] Strohmer | Pseudodifferential operators and Banach algebras in mobile communications[END_REF] for an overview and justification of banded matrices in mobile The operator described in the last example of Section 4.2 would actually be a (simple) example of a Multiple Gabor Multiplier.

Given such a scheme, the same questions as before may be asked. In particular, for a given Hilbert-Schmidt operator H ∈ H , is it possible to find time-frequency transfer functions such that the corresponding multiple Gabor multiplier minimizes the quadratic approximation error. To our knowledge, this question hasn't received an answer so far.

However, such a very general formulation turns out to assume a much simpler form when the synthesis windows h k are constructed through the action of the Weyl-Heisenberg group on the time-frequency plane. Indeed, it may be shown [START_REF] Dörfler | On the time-frequency representation of operators and generalized Gabor multiplier approximations[END_REF] that when the synthesis windows are taken from a suitable lattice in the time-frequency plane, i.e.

h k (t) = π(b k , ν k )h(t) = e 2iπν k t h(t -b k ) ,
the existence of optimal transfer functions m k is actually equivalent to the invertibility of some (right) discrete twisted convolution operator. When this twisted convolution operator is invertible, then a closed form may be obtained for the timefrequency transfer functions m k .

5.2. Twisted spline type spreading functions. Let us finally discuss a different form of operators, which allow for a discrete twisted convolution representation in the Gabor transform domain. Let φ be a square-integrable function, and assume that the spreading function η of H ∈ H may be written as a twisted spline type (TST) function Proof: Simply compute

(28) η(b, ν) = k,ℓ α kℓ φ(b -kb 0 , ν -ℓν 0 )e -2iπ(ν-ℓν0)kb0 = k,ℓ α kℓ φ kℓ (bν) , i.e.
H = k,ℓ α kℓ φ(b ′ , ν ′ )e -2iπν ′ kb0 M ν ′ +ℓν0 T b ′ +kb0 db ′ dν ′ = k,ℓ α kℓ φ(b ′ , ν ′ )e -2iπν ′ kb0 M ℓν0 M ν ′ T kb0 T b ′ db ′ dν ′ = k,ℓ α kℓ π(kb 0 , ℓν 0 ) φ(b ′ , ν ′ )M ν ′ T b ′ db ′ dν ′ ,
which proves the result. ♠ Now, suppose that H φ admits a Gabor multiplier representation, with timefrequency transfer function m; then 5 for f = m,n V g f (mb 0 , nν 0 )h mn , we have

H φ f = m,n m(m, n)V g f (mb 0 , nν 0 )h mn ,
for some h ∈ L 2 (R d ), with h mn = π mn h, and some bounded sequence of weights m(m, n). Plugging this expression into the above lemma, we obtain

Hf = k,ℓ α kℓ π kℓ m,n m(m, n)V g f (mb 0 , nν 0 )h mn = m,n m(m, n)V g f (mb 0 , nν 0 ) k,ℓ α kℓ π(kb 0 , ℓν 0 )π mn h .
This yields the following result 5 Here, once more, we implicitly assume, that the system of time-frequency translated functions gmn constitutes a frame and thus allows for reconstruction by means of a dual frame generated by a dual window h, see [START_REF] Gröchenig | Foundations of Time-Frequency Analysis[END_REF] for a review of Gabor frames.

Theorem 2. Assume that the spreading function η of H ∈ H is of the form (28), with φ the spreading function of a Gabor multiplier M m;g,h . Then H may be realized as follows (30) Hf = p,q

V g f (p, q)h pq , f ∈ L 2 (R d ) , where V g f (p, q) = (α♮W ) pq (31)

W mn = m(m, n)V g f (mb 0 , nν 0 ) . ( 32
)
In addition, if the right twisted convolution operator α → α♮G with the sequence G defined by

G mn = ∞ -∞ ∞ -∞
φ(b, ν)φ(b -mb 0 , ν -nν 0 )e -2iπmb0(ν-nν0) dbdν is invertible, then the coefficients α in (28) may be recovered from the scalar products η, φ kℓ . Remark 4. Notice that the expression provided in the theorem above is not a Gabor multiplier. It indeed expresses Hf as a sum of elementary atoms h pq , but the latter are not time-frequency shifts of a unique window function.

From a practical point of view, Gabor coefficients of Hf may be obtained by the following two-steps procedure: given f ∈ L 2 (R d ), and a (discrete) Gabor transform V g f (1) Weight coefficients V g f (mb 0 , nν 0 ) using the time-frequency transfer function m of the Gabor multiplier M m;g,h . (2) Evaluate the twisted convolution of the so-obtained weighted coefficients with the coefficients α of the TST expansion of the spreading function η of H. This results in a fairly simple algorithm. Corresponding numerical results will be reported elsewhere.

Let us point out that this scheme may be modified in several respects. In particular, it may be shown that using a different sampling lattice for the TST functions, one ends up with approximations involving true Gabor multipliers. These results are described in [START_REF] Dörfler | On the time-frequency representation of operators and generalized Gabor multiplier approximations[END_REF].

Conclusions and perspectives

In this paper, we have mainly focused on the approximation of 1D continuous time Hilbert Schmidt operators by Gabor multipliers and generalizations. Generalizing Gabor multipliers may be performed in various ways. We have described mainly two approaches, in which the synthesis windows or the time-frequency transfer function are varied. More general results are described in [START_REF] Dörfler | On the time-frequency representation of operators and generalized Gabor multiplier approximations[END_REF].

Let us also point out that this approach generalizes mutatis mutandis to more complex situations, for example Gabor transforms in higher dimensions, or Gabor transforms on non-rectangular (regular) lattices. It may also be formulated in a completely discrete setting, which makes it suitable for numerical applications.

This was actually one of our main motivations for developing such methods in a purely group theoretical setting. The latter opens interesting perspectives for generalizations, including wavelet multiplier approximations, or even more general situations.

2. 3 .

 3 Twisted convolution and Plancherel transform. The left convolution on the Weyl-Heisenberg group induces a twisted convolution product on E -1 , as follows

Figure 1 .

 1 Figure 1. Transposition performed by twisted convolution in the TF-domain

Figure 2 .

 2 Figure 2. Spreading functions and resulting operators applied to random signal

Lemma 4 .

 4 Let H ∈ H , with spreading function η. Assume there exists a non-negative function µ ∈ L 2 (R 2 ) such that |η(b, ν)| ≤ µ(b, ν) |V g h(b, ν)| ,then H admits a STFT multiplier representation.

Figure 3 .

 3 Figure 3. Spreading function of a Gabor multiplier

Figure 4 .

 4 Figure 4. Spreading function of a modifying Gabor multiplier

Figure 5 .

 5 Figure 5. Spreading function of a modifying Gabor multiplier

  as a sum of time-frequency translates of φ, using the H-covariant time-frequency translations. Then the following result expresses that covariance at the operator level Lemma 5. Let H be a Hilbert-Schmidt operator associated with a TST spreading function. Let H φ ∈ H denote the operator with spreading function φ. Then

	(29)	H =
		k,ℓ

α kℓ π(kb 0 , ℓν 0 )H φ .

Notice that this transform is only modulation covariant, i.e.V ψ Mν 0 x(b, ν) = V ψ x(b, ν + ν 0 ), but not translation-covariant, i.e. V ψ T b 0 x(b, ν) = V ψ (b + b 0 , ν).The other non-symmetric version of the Weyl-Heisenberg group yields a waveform transform that is shift-covariant, but not modulation covariant.

Acknowledgements

We wish to thank P. Balazs, H.G. Feichtinger, A.Klotz, Y.B. Kouagou and D. Marelli for extremely stimulating discussions. This work received support from the European Commission funded Research Training Network HASSIP (HPRN-CT-2002-00285).