Atte Aalto
email: atte.ej.aalto@gmail.com

CONVERGENCE OF DISCRETE-TIME KALMAN FILTER ESTIMATE TO CONTINUOUS-TIME ESTIMATE FOR SYSTEMS WITH UNBOUNDED OBSERVATION*

Keywords: Kalman filter, Infinite-dimensional systems, Boundary control systems, temporal discretization, sampled data 2010 AMS subject classification: 93E11, 47D06, 93C05, 60G15

Convergence of discrete-time Kalman filter estimate to continuous-time estimate for systems with unbounded observation

Introduction

The minimum variance state estimate for linear systems with Gaussian noise processes is given by the continuous-time Kalman filter. However, for obvious reasons, in a practical implementation the continuous-time system is often first discretized, and then the discrete-time Kalman filter is used on the discretized system. The objective of this article is to expand the recent results presented in [START_REF] Aalto | Convergence of discrete time Kalman filter estimate to continuous time estimate[END_REF] by the author on the convergence of the state estimate given by the discrete-time Kalman filter on the sampled system to the continuous-time estimate. There convergence results were shown for finite-dimensional systems and infinite-dimensional systems with bounded observation operators. The expansion in this paper covers systems with unbounded observation operators and systems whose dynamics are governed by an analytic semigroup. In particular, we shall show convergence rate estimates for the variance of the discrepancy between the discrete-and continuous-time estimates.

We study systems whose dynamics are given by (1.1)      dz(t) = Az(t) dt + Bdu(t), t ∈ R + , dy(t) = Cz(t) dt + dw(t), z(0) = x where A : X → X, B : U → X, and C : X → Y. The Hilbert spaces X, U = R q , and Y = R r are called the state space, the input space, and the output space, respectively. The mapping A is the generator of a C 0 -semigroup e At on X with domain D(A), B : R q → X is the control operator, and C : X → R r is called the observation operator. The dynamics equations (1.1) are given in the form of stochastic differential equations, see [START_REF]Stochastic differential equations: An introduction with applications[END_REF] by Øksendal for background. The input and output noise processes u and w are assumed to be q-and r-dimensional Brownian motions with incremental covariance matrices Q ≥ 0 and R > 0, respectively. Without loss of generality, we assume that there is no deterministic input, as it can always be removed by the usual techniques. The initial state x ∈ X is assumed to be a Gaussian random variable with mean m and covariance P 0 , denoted x ∼ N (m, P 0 ), and u, w, and x are assumed to be mutually independent.

The purpose of this paper is to study the discrepancy of the discrete-and continuous-time state estimates, defined by and ẑ(T ) := E z(T ) y(s), s ≤ T , respectively, and in particular, find convergence rate estimates for the variance E ||ẑ T,n -ẑ(T )|| 2 X as n → ∞ when the observation operator C is not bounded, which typically occurs when we get a pointwise or a boundary measurement from the computational domain of a system whose dynamics are governed through a partial differential equation. However, we do assume that C ∈ L(D(A), Y).

The state estimates ẑT,n and ẑ(T ) are obtained by the Kalman(-Bucy) filter -provided that the continuous-time Kalman filter equations are solvable. The Kalman filter was originally presented by Kalman in [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF] for discrete-time systems and by Kalman and Bucy in [START_REF] Kalman | New results in linear filtering and prediction theory[END_REF] for continuous-time systems. The infinite-dimensional generalization has been treated for example by Falb in [START_REF] Falb | Infinite-dimensional filtering: the Kalman-Bucy filter in Hilbert space[END_REF], by Bensoussan in [START_REF] Bensoussan | Filtrage Optimal des Systèmes Linéaires[END_REF], by Curtain and Pritchard in [START_REF] Curtain | Infinite Dimensional Linear Systems Theory[END_REF], and by Horowitz in [START_REF] Horowitz | Optimal filtering for gyroscopic noise[END_REF]. Of course the infinite-dimensional setting gives rise to many technical issues, such as unbounded control and observation operators and the solvability of the corresponding Riccati equations. These problems are tackled for example by Da Prato and Ichikawa in [START_REF] Da Prato | Riccati equations with unbounded coefficients[END_REF] and by Flandoli in [START_REF] Flandoli | Direct solution of a Riccati equation arising in a stochastic control problem with control and observation on the boundary[END_REF].

In the results of this paper we assume that the temporal discretization can be done perfectly, so that the only error source is the sampling of the continuous-time output signal. We refer to the review article [START_REF] Goodwin | Sampling and sampled-data models[END_REF] by Goodwin et al. for a discussion on the sampling of continuous-time systems and in particular [START_REF] Salgado | Connection between continuous and discrete Riccati equations with applications to Kalman filtering[END_REF] by Salgado et al. for a study on the sampled data Riccati equations and the Kalman filter. In practice, approximative numerical schemes are used for solving both the state estimate ẑT,n and the corresponding error covariance. For a discussion on this topic, see [START_REF] Axelsson | Discrete-time solutions to the continuous-time differential Lyapunov equation with applications to Kalman filtering[END_REF] by Axelsson and Gustafsson and [10] by Frogerais et al. treating nonlinear systems.

In Section 2, we shall introduce the ingredients for the proofs of our results. The main idea is to apply the discrete-time Kalman filter on a dense, numerable subset of the interval [0, T ]. This way we obtain a martingale that starts from the discrete-time estimate ẑT,n and converges almost surely to the continuous-time estimate ẑ(T ). We shall then find bounds for the increments of this martingale. These bounds are obtained by studying the regularity of the semigroup e At and in particular, the smoothness of the noise-free output Ce At x for x ∈ X. As was noted in [START_REF] Aalto | Convergence of discrete time Kalman filter estimate to continuous time estimate[END_REF] and as seen later in the proof of Theorem 4.1, the effects of the input noise process u and the initial state x can be treated separately. Therefore we shall first derive several results with different assumptions on the system concerning just the effect of the initial state in Section 3. Finally, in Section 4, we shall consider the effect of the input noise. The input noise effect is shown with the assumption of admissibility of the observation operator C.

Notation and standing assumptions. We denote by {e k } ∞ k=1 ⊂ D(A) an orthonormal basis for the state space X. The operator A generates a strongly continuous semigroup that is bounded by e At L(X) ≤ µ for t ∈ [0, T ].

Background

The idea of the proofs is exactly the same as in [START_REF] Aalto | Convergence of discrete time Kalman filter estimate to continuous time estimate[END_REF], but here we need to deal with many more technical issues. That is, we define a dense, numerable subset of the time interval [0, T ] and apply the discrete-time Kalman filter in this subset. Then we compute an upper bound for each increment in the state estimate and finally sum up these bounds. So let us define the time points t j for j = 1, 2, ... through the dyadic division (2.1) t j = T n j, j = 1, ..., n T 2 K n (2j -1), j = 2 K-1 n + 1, ..., 2 K n, for some K = 1, 2, ... The time point definition is illustrated in Figure 1. Then define T j := {t i } j i=1 and the X-valued martingale zj = E(z(T )|{y(t), t ∈ T j }). Define also the shorthand notation y(T j ) = {y(t), t ∈ T j }. Now it holds that zn = ẑT,n , and as discussed in [1, Section 2.1], as j → ∞, the martingale zj converges almost surely strongly to ẑ(T ). The idea in the proofs in this paper is to find upper bounds for the increments zj+1 -zj , for j ≥ n. As the martingale zj is square integrable, we have the telescope identity shown in [1, Lemma 1] for L, N ∈ N with L ≥ N :

(2.2) E ||z L -zN || 2 X = L-1 j=N E ||z j+1 -zj || 2 X .
We remark that setting N = n and letting

L → ∞ gives E ||ẑ(T ) -ẑT,n || 2 X = ∞ j=n E ||z j+1 -zj || 2 X .
Let us then establish an expression for one increment E ||z j+1 -zj || 2 X . Say [ξ, ξ 1 ] is a jointly Gaussian random variable in some product space. Denote ξ1 := E(ξ|ξ 1 ) and P 1 := Cov ξ1 -ξ, ξ1 -ξ . Then say

(2.3) ξ 2 = Hξ + w
where H is a bounded operator and w is a Gaussian random variable with mean zero and covariance R > 0 and w is independent of [ξ,

ξ 1 ]. Then it holds that ξ2 := E(ξ|{ξ 1 , ξ 2 }) = ξ1 + P 1 H * (HP 1 H * + R) -1 (ξ 2 -H ξ1 ). Then ξ2 -ξ1 = P 1 H * (HP 1 H * + R) -1 (H(ξ -ξ1 ) + w) from which it directly follows that Cov ξ2 -ξ1 , ξ2 -ξ1 = P 1 H * (HP 1 H * + R) -1 HP 1 ,
and further,

(2.4) E ξ2 -ξ1 2 = tr P 1 H * (HP 1 H * + R) -1 HP 1 .
The basic ingredients for the proofs in this paper are now presented, namely the martingale zj defined as E(z(T )|y(T j )) with T j defined in (2.1), the telescope identity (2.2), and equation (2.4) for the increment norm. Later we sometimes need the assumption that x ∈ D(A) almost surely. With Gaussian random variables this means that x is actually a D(A)-valued random variable. Proposition 2.1. Let ξ be an X-valued Gaussian random variable s.t. ξ ∈ X 1 almost surely where X 1 ⊂ X is another Hilbert space with continuous and dense embedding. Then ξ is an X 1 -valued Gaussian random variable.

Proof. Pick h ∈ X 1 . We intend to show that ξ, h X 1 is a real-valued Gaussian random variable. For h ∈ X 1 there exists

h ′ ∈ X ′ 1 , the dual space of X 1 , s.t. ξ, h X 1 = ξ, h ′ (X 1 ,X ′ 1 )
and further, there exists a sequence

{h i } ∞ i=1 ⊂ X such that ξ, h ′ (X 1 ,X ′ 1 ) = lim i→∞ ξ, h i X . Now ξ, h i X
is a pointwise converging sequence of Gaussian random variables and so the limit is also Gaussian.

Fernique's theorem [START_REF] Da Prato | Stochastic equations in infinite dimensions[END_REF]Theorem 2.6] can be applied to note that if ξ is an X 1 -valued Gaussian random variable then ξ ∈ L p (Ω; X 1 ) for any p > 0. In particular, E ||ξ|| 2 X 1 < ∞ and if A ∈ L(X 1 , X) then Aξ is an X-valued Gaussian random variable.

Convergence results without input noise

Assume now that there is no input noise in the system (1.1), that is, u = 0. The output is then given by

y(t) = C t 0 e As x ds + w(t)
where w is a Brownian motion. In a sense, the output is parameterized by the initial state x, and therefore we also define the martingale xj := E(x|y(T j ). Then in the absence of input noise it holds that zj = e AT xj .

In order to use (2.4) to compute one increment E ||z j+1 -zj || 2 X , we need to consider how to take into account an intermediary observation y(t j+1 ) in the state estimate. Obviously the noise process value w(t j+1 ) is not independent of the measurements y(t i ) with i = 1, ..., j as required in order to use (2.4). As in [START_REF] Aalto | Convergence of discrete time Kalman filter estimate to continuous time estimate[END_REF], the dependence of the noise term w(t j+1 ) on y(t i ) with i = 1, ..., j is removed by subtracting from y(t j+1 ) the linear interpolant 1 2 y(t j+1 -h) + y(t j+1 + h) where h = T 2 K n with the corresponding K (see (2.1) and Figure 1). Note that t j+1 ± h ∈ T j . When there is no input noise, this new output is obtained as

(3.1) ỹj+1 = C h (t j+1 )x + wj+1
where C h (t) is defined for x ∈ X and t ≥ h by

(3.2) C h (t)x := C 2 t t-h
e As x ds -t+h t e As x ds , and wj+1 ∼ N 0, h 2 R is independent of y(t i ) and hence of ỹi for i = 1, ..., j. Note that we can first define C h (t)x only for x ∈ D(A). However, it holds that t+h t e As x ds

2 D(A) = t+h t e As x ds 2 X + A t+h t e As x ds 2 X ≤ h 2 µ 2 ||x|| 2 X + e A(t+h) -e At x 2 X ≤ h 2 + 4 µ 2 ||x|| 2 X .
Therefore, we can uniquely extend C h (t) to a continuous operator from X to Y, and

||C h (t)|| L(X,Y) ≤ µ √ h 2 + 4 ||C|| L(D(A),Y) . In addition, using Ce As ξ = Ce At ξ + s t CAe Ar ξdr in (3.2) yields a useful bound (3.3) ||C h (t)ξ|| Y ≤ 1 2 t t-h t s CAe Ar ξ Y dr ds + t+h t s t CAe Ar ξ Y dr ds provided that CAe Ar ξ ∈ L 1 (t -h, t + h; Y). Now the form (3.1) is exactly as (2.
3) and so (recalling zj = e AT xj ) we can use (2.4) to obtain

E ||z j+1 -zj || 2 X = tr e AT P j C h (t j+1 ) * × (3.4) × C h (t j+1 )P j C h (t j+1 ) * + h 2 R -1 C h (t j+1 )P j e A * T .
Finally we are able to show the main lemma, which links the convergence of the state estimate to properties of the operator C h (t) and the smoothness of the output y. This lemma serves as the basis for all the proofs of our main theorems when there is no input noise.

Lemma 3.1. Let X 1 be a dense Hilbert subspace of X with a continuous embedding. Assume that there exist M > 0 and k > 1 such that for any K ∈ N and ξ ∈ X 1 , it holds that

(3.5) 2 K n j=2 K-1 n+1 ||C h (t j )ξ|| 2 Y ≤ M h k ||ξ|| 2 X 1
where h = T 2 K n . Assume also that x ∈ X 1 almost surely and u = 0. Then for ẑT,n and ẑ(T ) defined in (1.2), it holds that

E ||ẑ T,n -ẑ(T )|| 2 X ≤ 2M E ||x|| 2 X 1 E ||ẑ T,n -z(T )|| 2 X (2 k-1 -1) min(eig(R)) T k-1 n k-1 .
Note that E ||x|| 2 X 1 is well defined and finite by Proposition 2.1. A strict a priori result is obtained by replacing

E ||ẑ T,n -z(T )|| 2 X ≤ µ 2 E ||x|| 2 X .
Proof. We will use a shorter notation C h = C h (t j+1 ) in the proof. Note that

C h P j C * h + h 2 R -1 ≤ 2 h min(eig(R)) =: C R h .
One increment can be bounded by (3.4):

tr e AT P j C * h C h P j C * h + h 2 R -1 C h P j e A * T = ∞ k=1 C h P j e A * T e k , C h P j C * h + h 2 R -1 C h P j e A * T e k Y ≤ C R h ∞ k=1 C h P j e A * T e k 2 Y = C R h ∞ k=1 E C h (x j -x) e AT (x j -x), e k X 2 Y ≤ C R h E ||C h (x j -x)|| 2 Y ∞ k=1 E e AT (x j -x), e k 2 X ≤ C R h tr(C h P j C * h )E ||z j -z(T )|| 2 X .
Now we use this bound, P j ≤ P n ≤ P 0 , and assumption (3.5) to sum up the increments corresponding to h = T 2 K n , that is, j = 2 K-1 n, ..., 2 K n -1:

2 K n-1 j=2 K-1 n E ||z j+1 -zj || 2 X ≤ M C R E ||x|| 2 X 1 E ||z n -z(T )|| 2 X T 2 K n k-1
and finally, summing up over K = 1, 2, ... gives the result.

3.1. Diagonalizable main operator. We proceed to prove a convergence result for systems with unbounded observation operator C -provided that A is (unitarily) diagonalizable. The proof is based on Lemma 3.1. To get a useful bound for ||C h ξ|| 2 Y , some assumptions on the degree of unboundedness of C and the spectral asymptotics of A are required. Theorem 3.1. Let ẑT,n and ẑ(T ) be as defined above in (1.2). Denote by {λ k } ∞ k=1 ⊂ C the spectrum of A ordered so that |λ k | is non-decreasing and let {e k } ∞ k=1 ⊂ D(A) be the corresponding set of eigenvectors that give an orthonormal basis for X. Make the following assumptions on x, A, and C:

(i) x ∈ D(A) almost surely; (ii) There exists δ > 1/2 such that lim k→∞ |λ k | k β = 0 when β > δ, ∞ when β < δ;
(iii) There exists γ ∈ [0, 1) such that 2γ + 1/δ < 2 and

sup k ||Ce k || Y |λ k | γ < ∞.
Then the following holds:

• If lim k→∞ |λ k | k δ = Γ ∈ (0, ∞), then E ||ẑ T,n -ẑ(T )|| 2 X ≤ M T 3-2γ-1/δ n 2-2γ-1/δ
where the constant M is given below in (3.7).

• If either this limit does not exist, or it is 0 or ∞, then for all ǫ ∈ 0, δ -

1 2-γ E ||ẑ T,n -ẑ(T )|| 2 X ≤ M ǫ T 3-2γ-1/(δ+ǫ) n 2-2γ-1/(δ-ǫ)
where the ǫ-dependent constant M ǫ is given below also in (3.7) but with different, ǫ-dependent parameters (see the last paragraph of the proof ).

For example, 1D wave equation on interval [0, L] with Dirichlet boundary conditions in the natural state space where some pointwise value of the state is observed, satisfies the assumptions of the above theorem with δ = 1 and γ = 0. The limit of

|λ k | k as k → ∞ exists and it is Γ = π 2L . This would imply convergence rate E ||ẑ T,n -ẑ(T )|| 2 X ≤ M T 2 n .
Proof. Assume first that lim k→∞

|λ k | k δ = Γ ∈ (0, ∞). Denote ξ = ∞ k=1 α k e k ∈ D(A) which is equivalent to ∞ k=1 |λ k | 2 α 2 k ≤ ∞. Now (3.6) C h ξ = ∞ k=1 α k 2 t t-h e λ k s ds - t+h t e λ k s ds Ce k .
For the term inside parentheses, we have

||C h ξ|| Y ≤ n(h) k=1 |α k | 2 µ ||Ce k || Y |λ k |h 2 + ∞ k=n(h)+1 |α k | ||Ce k || Y 2µ |λ k | =: (I) + (II).
We then proceed to find upper bounds for the two parts. Using Cauchy-Schwartz inequality and denoting Γ := sup k

|λ k | k δ gives (I) ≤ µ h 2 2   n(h) k=1 α 2 k ||Ce k || 2 Y |λ k | 2-2γ   1/2   n(h) k=1 |λ k | 2γ   1/2 ≤ µh 2 Γγ 2 M I   n(h) k=1 k 2γδ   1/2
where

M I = n(h) k=1 α 2 k ||Ce k || 2 Y |λ k | 2-2γ 1/2
. The sum inside the parentheses can be bounded from above by the integral

n(h)+1 0 x 2γδ dx to get (I) ≤ µh 2 Γγ 2 √ 2γδ + 1 M I (n(h) + 1) 2γδ+1 ≤ 3 δ µ Γγ 2 √ 2γδ + 1 M I h 2-γ-1 2δ ≤ 3 δ µ Γγ 2 M I h 2-γ-1 2δ
where the last row follows from the facts that

(n(h) + 1) 2γδ+1 ≤ (h -1/δ + 2) 2γδ+1 = (1+2h 1/δ ) γδ+ 1 2 h -γ-1 2δ ≤ 3 δ h -γ-1 2δ
if h ≤ 1, and that 2γδ + 1 > 1.

For the second part, assume |λ k | ≥ Γk δ for k ≥ n(h) + 1 where Γ = 0.9Γ for example. Again, using Cauchy-Schwartz inequality yields

(II) ≤ 2µ   ∞ k=n(h)+1 α 2 k ||Ce k || 2 Y |λ k | 2-2γ   1/2   ∞ k=n(h)+1 1 |λ k | 4-2γ   1/2 ≤ 2µ Γ2-γ M II   ∞ k=n(h)+1 1 k (4-2γ)δ   1/2 where M II = ∞ k=n(h)+1 α 2 k ||Ce k || 2 Y |λ k | 2-2γ 1/2
. Now the sum inside the parentheses can be bounded from above by the integral

∞ n(h) 1
x (4-2γ)δ dx. Note that our assumptions on γ and δ imply (4 -2γ)δ > 2. So we get

(II) ≤ 2µM II Γ2-γ (4 -2γ)δ -1 1 n(h) (4-2γ)δ-1 1/2 ≤ 2µ Γ2-γ M II h 2-γ-1 2δ
where in the last row we have used n(h) ≥ h -1/δ .

Combining the bounds gives

||C h ξ|| 2 Y ≤ 2 (I) 2 + (II) 2 ≤ 2µ 2 M 2 I + M 2 II max 9 δ Γ2γ 4 , 4 Γ4-2γ h 4-2γ-1/δ ≤ 2µ 2 ||Aξ|| 2 X sup k ||Ce k || 2 Y |λ k | 2γ max 9 δ Γ2γ 4 , 4 Γ4-2γ h 4-2γ-1/δ
where we have used

M 2 I + M 2 II = ∞ k=1 α 2 k ||Ce k || 2 Y |λ k | 2-2γ ≤ ∞ k=1 |λ k | 2 α 2 k sup j ||Ce j || 2 Y |λ j | 2γ .
Note that we assumed that we could choose for example Γ = 0.9Γ. In some sense this is not our choice but we need to make sure that the "original" h = T 2n is small enough so that n(T /(2n)) = 2n T 1/δ is such that there exists

Γ > 0 for which |λ k | k δ ≥ Γ for k ≥ n(T /(2n)).
To get a bound for the sum in (3.5), we simply multiply the bound obtained for

||C h ξ|| 2 Y by 2 K-1 n = T 2h to get 2 K n j=2 K-1 n+1 ||C h (t j )ξ|| 2 Y ≤ sup k ||Ce k || 2 Y |λ k | 2γ max 9 δ Γ2γ 4 , 4 Γ4-2γ µ 2 T ||ξ|| 2 D(A) h 3-2γ-1/δ
and so the result follows by Lemma 3.1 with (3.7)

M = 2µE ||ẑ T,n -z(T )|| 2 X E ||x|| 2 D(A) (2 2-2γ-1/δ -1) min(eig(R)) sup k ||Ce k || 2 Y |λ k | 2γ max 9 δ Γ2γ 4 , 4 Γ4-2γ .
In the case that lim k→∞ |λ k | k δ is 0, ∞, or it does not exist, some modifications are required to the bounds of (I) and (II). In the bound for (I), δ needs to be replaced by δ +ǫ and then Γǫ = sup k

|λ k | k δ+ǫ < ∞.
In the bound for (II), δ needs to be replaced by δ -ǫ and then Γǫ = inf k≥n(h)+1

|λ k | k δ-ǫ > 0.
The assumption (iii) in the theorem differs from our minimal assumption

C ∈ L(D(A), Y) which is equivalent to ||Ce k || Y |λ k |
∈ l 2 for unitarily diagonalizable A. It is possible to construct a system for which C ∈ L(D(A), Y) but (iii) does not hold. Remark 3.1. Theorem 3.1 can be extended to γ < 0. In that case, when determining the bounds for (I) and (II), the computations are carried out as if γ were zero. This eventually leads to a bound

E ||ẑ T,n -ẑ(T )|| 2 X ≤ M T 3-1/δ
n 2-1/δ . Note that if assumption (iii) holds for γ < -1 2δ then C is actually bounded.

Admissible observation operator.

In the next result we assume that the observation operator C is admissible in the sense of Weiss [START_REF] Weiss | Admissible observation operators for linear semigroups[END_REF]. One good example of systems that satisfy assumption (iii) in the following theorem is provided by scattering passive boundary control systems, see the article [START_REF] Malinen | Conservative boundary control systems[END_REF] by Malinen and Staffans. For a more extensive background, we refer to [START_REF] Staffans | Well-posed Linear Systems[END_REF] by Staffans. Theorem 3.2. Let ẑT,n and ẑ(T ) be as defined above in (1.2) and u = 0. Make the following assumptions:

(i) x ∈ D(A) almost surely;
(ii) The orthonormal basis {e k } ⊂ X is such that e k ∈ D(A 2 ) for every k ∈ N and there exists δ > 1/2 such that for ξ = ∞ k=1 α k e k the norm given by

∞ k=1 k 2δ α 2 k is equivalent to the D(A)-norm and ∞ k=1 k 4δ α 2 k is equivalent to the D(A 2 )-norm; (iii)
The observation operator is admissible, that is, for any T ≥ 0 there exists

H T ≥ 0, such that Ce A(•) x L 2 ((0,T );Y) ≤ H T ||x|| X .
Then

E ||ẑ T,n -ẑ(T )|| 2 X ≤ M (T ) T 2-1/2δ n 1-1/2δ with M (T ) = 2E ||ẑT,n-z(T )|| 2 X E(||x|| 2 D(A) ) (2 1-1/2δ -1) min(eig(R)) max 3 2δ+1 T µ||C|| 2 L(D(A),Y) 8δ+4 , H T 2δ-1 .
Proof. In this proof, the aforementioned norms are used in D(A) and D(A 2 ). We need to utilize the global output bound Ce

A(•) x L 2 ((0,T );Y) ≤ H T ||x|| X .
To this end, define a stacked operator

C h := [C h (h), C h (3h), . .

. , C h (T -h)]

T for h = T 2 K n mapping to a product space Y 2 K-1 n . Then the sum on the left hand side of (3.5) is obtained as C h ξ 2 Y 2 K-1 n . In this proof, a i N i=1 is used to denote an augmented vector with N components a i .

The proof proceeds similarly as the proof of Theorem 3.1 but the sum in (3.6) is split using the index n(h

) = ⌈h -1/(2δ) ⌉ to get C h ξ Y 2 K-1 n ≤ (I) + (II) where (I) = n(h) k=1 α k 2 C (2j-1)h (2j-2)h
e As e k ds -C 2jh (2j-1)h

e As e k ds

2 K-1 n i=1 Y 2 K-1 n ≤ n(h) k=1 |α k | 2 T 2h h 2 µ ||C| | L(D(A),Y) k 2δ ≤ T h 3 8 µ ||C| | L(D(A),Y)   n(h) k=1 k 2δ α 2 k   1/2   n(h) k=1 k 2δ   1/2
where the first inequality is obtained using (3.3) with t = (2j -1)h and bounding the derivative of Ce Ar e k by

CAe Ar e k Y ≤ µ ||C| | L(D(A),Y) ||e k || D(A 2 ) = µ ||C|| L(D(A),Y) k 2δ
and noting that 2 K-1 n = T 2h .

For the remaining part it holds that

(II) = ∞ k=n(h)+1 α k 2 C (2i-1)h (2i-2)h
e As e k ds -C 2ih (2i-1)h

e As e k ds

2 K-1 n i=1 Y 2 K-1 n ≤ ∞ k=n(h)+1 |α k | 2 √ 2hH T ≤ h 2 H T   ∞ k=n(h)+1 k 2δ α 2 k   1/2   ∞ k=n(h)+1 k -2δ   1/2
since it holds that

C (2i-1)h (2i-2)h
e As e k ds -C 2ih (2i-1)h

e As e k ds

2 K-1 n i=1 Y 2 K-1 n ≤   2 K-1 n i=1 2ih (2i-2)h Ce As e k Y ds 2   1/2 ≤ √ 2h Ce As e k L 2 ((0,T );Y)
where the last inequality follows from Cauchy-Schwartz inequality. Finally, the sums n(h) k=1 k 2δ and ∞ k=n(h)+1 k -2δ are bounded by integrals of x 2δ and x -2δ , respectively, as in the proof of Theorem 3.1 and then the result is obtained by Lemma 3.1 and

C h ξ 2 Y 2 K-1 n ≤ 2 (I) 2 + (II) 2 . 3.3. Analytic semigroup.
In this section we show the convergence estimate when A is the generator of an analytic semigroup. One result is first shown without additional assumptions for bounded and unbounded observation operator C. Then we assume further that -A is a sectorial operator in X which enables us to treat non-integer powers (-A) η for η ≥ 0. An example of such case is provided by heat equation treated below in Example 3.1.

An important tool here is that for analytic semigroups it holds that [START_REF] Tanabe | Equations of Evolution[END_REF]Theorem 3.3.1]). Using this to bound the derivative of Ce At ξ, that is CAe At ξ, gives by (3.3),

(3.8) A κ e At L(X) ≤ c(κ) t κ , t > 0, κ ∈ N (see
(3.9) ||C h (t)ξ|| Y ≤ c(1) ||C|| L(X/D(A),Y) ||ξ|| X/D(A) 2(t -h) h 2 , t > h.
Theorem 3.3. Let ẑT,n and ẑ(T ) be as defined above in (1.2). Assume A is the generator of an analytic C 0 -semigroup and assume either

(i) C ∈ L(X, Y), or (ii) C ∈ L(D(A), Y) and x ∈ D(A) almost surely. Then E ||ẑ T,n -ẑ(T )|| 2 X ≤ M T n where M = 2||C|| 2 L(X/D(A),Y)
min(eig(R))

µ 2 + c(1) 2 π 2 96 E ||z(T ) -ẑT,n || 2 X E ||x|| 2 X/D(A) .
Proof. The proofs for the two cases are identical so only the case (i) is presented. In the second case just replace X by D(A) in ||ξ|| X and ||C|| L(X,Y) .

The proof is based on Lemma 3.1, so we need to find a bound for (3.10)

2 K n j=2 K-1 n+1 ||C h (t j )ξ|| 2 Y = 2 K-1 n l=1 C h (2l -1)h ξ 2 Y For l = 1, we use ||C h (h)ξ|| Y ≤ hµ ||C|| L(X,Y) ||ξ|| X which is clear from the definition of C h (t) in (3.
2). For l > 1, we use (3.9) where the denominator becomes 2 (2l -1)h -h = 4h(l -1) and so,

2 K n j=2 K-1 n+1 ||C h (t j )ξ|| 2 Y ≤ ||C|| 2 L(X,Y) ||ξ|| 2 X   µ 2 + c(1) 2 16 2 K-1 n l=2 1 (l -1) 2   h 2 ≤ ||C|| 2 L(X,Y) ||ξ|| 2 X µ 2 + c(1) 2 π 2 96 h 2 .
The result follows by Lemma 3.1.

One more case is treated where A is as before and, in addition, -A is a sectorial operator, see [START_REF] Arendt | Vector-valued Laplace Transforms and Cauchy Problems[END_REF]Section 3.8] for definitions. Then it is possible to define non-integer powers (-A) η where η ∈ R and spaces D((-A) η ) equipped with the corresponding graph norm. Also (3.8) holds then for non-integer κ ≥ 0 if A is replaced by -A, see [START_REF] Tanabe | Equations of Evolution[END_REF]Thm. 3.3.3]. In particular, if A is strictly negative definite, then it is sectorial. This type of systems are also studied in [START_REF] Da Prato | Riccati equations with unbounded coefficients[END_REF] and [START_REF] Flandoli | Direct solution of a Riccati equation arising in a stochastic control problem with control and observation on the boundary[END_REF]. Theorem 3.4. Let ẑT,n and ẑ(T ) be as defined above in (1.2). Assume that A is the generator of an analytic C 0 -semigroup and, in addition, -A is a sectorial operator. Then assume C ∈ L(D((-A) ν ), Y) and x ∈ D((-A) η ) almost surely where ν ∈ R and η ∈ R are such that |η -ν| < 1/2. Then1

E ||ẑ T,n -ẑ(T )|| 2 X ≤ M T 1+2(η-ν) n 1+2(η-ν)
where M is given below in (3.12).

Proof. This is done exactly as the proof of Theorem 3.3 above. Just the bounds for ||C h ((2l -1)h)ξ|| Y in the summation (3.10) are computed differently. To begin with, we note that (3.8) with non-integer κ = 1 -η + ν yields,

(3.11) CAe At ξ Y ≤ ||C|| L(D((-A) ν ),Y) ||ξ|| D((-A) η ) c(1 -η + ν) t 1-η+ν .
When treating the term with l = 1 in (3.10), the cases ν ≥ η and ν < η have to be considered separately. First for ν ≥ η,

Ce At ξ Y = C(-A) -ν (-A) ν-η e At (-A) η ξ Y ≤ ||C|| L(D((-A) ν ),Y) ||ξ|| D((-A) η ) c(ν -η) t ν-η .
Then directly by the definition of C h in (3.2) (recalling 1 + η -ν > 0),

||C h (h)ξ|| Y ≤ ||C|| L(D((-A) ν ),Y) ||ξ|| D((-A) η ) c(ν -η) 2h 0 1 s ν-η ds ≤ ||C|| L(D((-A) ν ),Y) ||ξ|| D((-A) η ) c(ν -η) 1 + η -ν (2h) 1+η-ν .
For ν < η < 1+ν, use (3.3) with t = h. Using (3.11) to bound the derivative norm CAe At ξ Y and computing the integrals yields

||C h (h)ξ|| Y ≤ ||C|| L(D((-A) ν ),Y) ||ξ|| D((-A) η ) c(1 -η + ν) 2(η -ν) 2 1+η-ν -2 1 + η -ν h 1+η-ν ≤ ||C|| L(D((-A) ν ),Y) ||ξ|| D((-A) η ) 2 ln 2 c(1 -η + ν) 1 + η -ν h 1+η-ν .
For the terms with l > 1, it holds by (3.3) and (3.11), that

||C h ((2l -1)h)ξ|| Y ≤ h 2 ||C| | L(D((-A) ν ),Y) c(1 -η + ν) 2(2h) 1-η+ν (l -1) 1-η+ν ||ξ|| D((-A) η ) .

Now summing up the bounds for ||C

h ((2l -1)h)ξ|| 2 Y for l = 1, ..., 2 K-1 n and using ∞ l=2 1 (l-1) 2+2(ν-η) ≤ 2-2(η-ν)
1-2(η-ν) yields a bound for the sum in (3.5), and finally Lemma 3.1 can be used to get the result with

M = 2 ||C|| 2 L(D((-A) ν ),Y) E ||z(T ) -ẑT,n || 2 X E ||x|| 2 D((-A) η ) (2 1+2(η-ν) -1) min(eig(R)) × (3.12) × M ν,η + c(1 -η + ν) 2 2 4-2(ν-η) 2 -2(η -ν) 1 -2(η -ν)
where the term with l = 1 gives

M ν,η =    4(ln 2) 2 c(1-η+ν) 2 (1+η-ν) 2 if η > ν, 2 2+2(η-ν) c(ν-η) 2 (1+η-ν) 2 if η ≤ ν. Example 3.1. Consider the 1D heat equation          ∂ ∂t z(x, t) = ∂ 2 ∂x 2 z(x, t), x ∈ [0, 1], z(0, t) = z(1, t) = 0, z(x, 0) = z 0 , dy(t) = ∂ ∂x z(0, t) dt + dw(t) with state space X = L 2 (0, 1) and D(A) = H 2 0 [0, 1]. Assume z 0 ∈ D(A) almost surely. Now the spectrum of A is {-π 2 k 2 }
and the corresponding eigenvectors are e k = sin(πkx). Then it is easy to see that the assumptions of Theorem 3.1 are satisfied with δ = 2 and γ = 1/2 and thus the theorem implies convergence rate O(n -1/2 ) for E ||ẑ T,n -ẑ(T )|| 2 X . The assumptions of Theorem 3.2 are satisfied with δ = 2 implying convergence rate O(n -3/4 ), and Theorem 3.3 clearly implies convergence rate O(n -1 ) but we can do even better.

Denoting z = ∞ k=1 α k e k we have ||z|| 2 D((-A) ν ) = ∞ k=1 k 4ν α 2 k .
For the output it holds that

|Cz| 2 = ∞ k=1 πkα k 2 ≤ π ∞ k=1 1 k 1+ǫ ∞ k=1 k 3+ǫ α 2 k from which it can be deduced that C ∈ L(D((-A) ν ), Y)) for ν > 3/4. Now Theorem 3.4 implies convergence rate O(n -3/2+ǫ ) for E ||ẑ T,n -ẑ(T )|| 2 X
with ǫ > 0 -of course, with a multiplicative constant that tends to infinity as ǫ → 0.

Input noise

In this section, we shall study the effect of the input noise u in (1.1). This is done only under the assumption of an admissible observation operator C. The main theorem is a generalization of [1, Theorem 2] where C was assumed to be bounded, and the proof follows the same outline.

In the cases without input noise, the state z(t) was parameterized by the initial state through z(t) = e At x. Now we define the solution operator S(t) through

S(t) : [x, u] → e At x + t 0 e A(t-s) Bdu(s).
Formally we first define S(t) for u ∈ H 1 (0, T ; U) and then extend it for Brownian motion as a Wiener integral. Then the solution to (1.1) is given by z(t) = S(t)[x, u] and the conditional expectation over a given sigma algebra σ by E(z(t)|σ) = S(t)E([x, u]|σ). In the following theorem, we consider virtually estimating [x, u] using the outputs y(t j ) and then the state estimate is obtained by S(T ). The well-posedness of S(t)[0, ũj ] for ũj := E(u|{y(t i ), i = 1, ..., j}) is established next. Lemma 4.1. For fixed j, it holds that ũj (•) ∈ H 1 (0, T ; U).

Proof. For fixed t ∈ [0, T ] and j ∈ N, it holds that

E(u(t)|y(T j )) = Cov [u(t), y(T j )] Cov [y(T j ), y(T j )] -1 y(T j ).
Here Cov [y(T j ), y(T j )] -1 y(T j ) is a well defined finite-dimensional vector, and so we only need to differentiate Cov [u(t), y(T j )] ∈ L(Y j , U) which is equivalent to the differentiability of its adjoint, Cov [y(T j ), u(t)] which is a block operator with components Cov [y(t i ), u(t)], for i = 1, ..., j. We show that each block has a bounded derivative.

Since We are now ready to proceed to the main result of the paper.

y(t i ) = C t i 0 e As x
Theorem 4.1. Assume that C is an admissible observation operator, that is, for all T ≥ 0 there exists H T ≥ 0 such that

T 0 Ce At x 2 Y dt ≤ H 2 T ||x|| 2 X .
In addition, assume that one set of assumptions in Theorems 3.1-3.4 are satisfied. Then for ẑT,n and ẑ(T ) defined in (1.2) it holds that

E ||ẑ T,n -ẑ(T )|| 2 X ≤ M 1 (T ) n + M 2 (T )
n 2 + err x where M (T ) is given below in (4.5) and err x is given by the respective Theorem 3.1-3.4.

Proof. The first part of proof follows essentially the same outline as the proof of Lemma 3.1. Say we are estimating [x, u] and we have [x j , ũj ] := E([x, u]|F j ) and the corresponding error covariance P j . Then the state estimate and the corresponding error covariance are given by zj = S(T )[x j , ũj ] and S(T )P j S(T ) * -although in the last equation the formal adjoint S(T ) * is only defined in connection with the covariance P j , namely P j S(T

) * h = E [x j -x, ũj -u] S(T )[x j -x, ũj -u], h X .
The output in (1.1) is given by As with the noiseless case, when we are including the measurement y(t j+1 ), in order to get rid of the output noise correlation, we shall subtract the linear interpolant from y(t j+1 ), namely define

y(t) = C
ỹj+1 = y(t j+1 ) - 1 2 y(t j+1 + h) - 1 2 y(t j+1 -h)
where h = T 2 K n for the corresponding K, given in (2.1). Now this output can be written as ỹj+1 = C (t j+1 )[x, u] T + wj+1 where wj+1 ∼ N (0, h/2 R) is independent of w(t i ) and hence, wi with i = 1, ..., j, and for t ≥ h. Now the error covariance increment is obtained just as in Section 3, and by the same computation as in the proof of Lemma 3.1 (just replacing e AT by S(T ), P j by P j , and C h by C ), we get

C (t) := [C h (t), C h,u (t)] with C h (t) defined in (3.
(4.2) E ||z j+1 -zj || 2 X ≤ 2 h min(eig(R)) E ||C [x, u]|| 2 Y E ||ẑ T,n -z(T )|| 2 X .
In order to get a suitable bound for the increment, we must find a bound for the term In (I), r ∈ [0, t j+1 ] and in (II), r ∈ [t j+1 , t j+1 + h], and thus they are independent. Therefore,

E ||C [x, u]|| 2 Y = E ||C h x|| 2 Y + E ||C h,
(t j+1 + h -s)(t j+1 + h -r) × Ce A(s-t) ABQB * A * e A * (
E ||(I) + (II)|| 2 Y ≤ h 4 6 tr(Q) Ce As AB 2 L 2 (0,T ;L(U,Y)) .

It is well known that

Cov h 0 Bu(t) dt, h 0 Bu(t) dt = h 3 3 BQB *
and so

E   C 2 t j+1 +h t j+1 B(u(t) -u(t j+1 ))dt 2 Y   ≤ h 3 12 ||C|| 2 L(D(A),Y) ||B| | 2 L(U,D(A)) tr(Q).
In (4.3), the two B(u(t) -u(t j+1 ))-terms are independent (because in the first one, t ≤ t j+1 and in the second, t ≥ t j+1 ) and by utilizing this and gathering the above bounds, we get

E ||C h,u (t j+1 )u|| 2 Y ≤ 6E ||(I) + (II)|| 2 Y + 6 h 3 12 ||C|| 2 L(D(A),Y) ||B|| 2 L(U,D(A)) tr(Q) ≤ h 4 Ce As AB 2 L 2 (0,T ;L(U,Y)) tr(Q) + h 3 2 ||C|| 2 L(D(A),Y) ||B|| 2 L(U,D(A)) tr(Q).
Combining this with (4.2) gives

E ||z j+1 -zj || 2 X ≤ 2h 3 min(eig(R)) Ce As AB 2 L 2 (0,T ;L(U,Y)) tr(Q)E ||ẑ T,n -z(T )|| 2 X + h 2 min(eig(R))
||C|| where incr [x,j+1] is the contribution of the initial state x obtained from one of the Theorems 3.1-3.4.

As was done in the proof of Theorem 3.1, the bound obtained for increment E ||z j+1 -zj || 2 X is multiplied by T 2h to get an upper bound for increments corresponding to one value h, that is, one K (see Figure 1): (4.4)

2 K n-1 j=2 K-1 n E ||z j+1 -zj || 2 X ≤ M 1 (T ) 2 K n + 3M 2 (T ) (2 K n) 2 + incr [x,

Discussion

In this paper we extended the convergence results presented by the author in [START_REF] Aalto | Convergence of discrete time Kalman filter estimate to continuous time estimate[END_REF]. There the convergence rate estimates for E ||ẑ T,n -ẑ(T )|| 2 X were shown for finite-dimensional systems and infinite-dimensional systems with bounded observation operator C. Now convergence rate estimates were found for systems with unbounded observation operators with some additional assumptions on the system operators. Firstly, a result was shown for systems with diagonalizable main operators A. In this case, some additional assumptions were needed, including a slightly nonstandard assumption on the output operator (assumption (iii) in Theorem 3.1). In the problems arising from PDEs on one-dimensional spatial domains, this is not a big problem but unfortunately with more complicated systems, finding a suitable γ might be close to a mission impossible. The spectral asymptotics, on the other hand, is an extensively studied field -so much so that it has even been a subject of a few books, such as [START_REF] Levendorskiì | Asymptotic Distribution of Eigenvalues of Differential Operators[END_REF] by Levendorskiì and [18] Safarov and Vassiliev. Theorem 3.2 treats the case with assuming essentially just the admissibility of the observation operator C. Two results were shown for systems with analytic semigroups, in which case the other technical assumptions were not needed. The effect of the input noise was studied in Theorem 4.1, which extends the corresponding earlier result [START_REF] Aalto | Convergence of discrete time Kalman filter estimate to continuous time estimate[END_REF]Theorem 2] to admissible observation operators.

In all results of the paper -except for the analytic semigroup case without input noise -the convergence rate estimates are of the form M T k+1 n k = M T ∆t k , meaning that the estimates deteriorate as T grows. In general, this cannot be completely avoided because there is no guarantee that the output sampling does not cause essential loss of information. A result where the bound would not deteriorate as T grows could be possible under some additional assumptions (like exponential stability of the system), but the long time behaviour should be anyway studied by comparing the solutions of the corresponding discrete-and continuous-time algebraic Riccati equations.

We remark that in this paper, as well as in [START_REF] Aalto | Convergence of discrete time Kalman filter estimate to continuous time estimate[END_REF], it has been assumed that the time discretization can be done perfectly and the only error source is the time-sampling of the output signal y that is defined in continuous time. Further research would be needed to estimate the error caused by approximate discretization schemes.

( 1 . 2 )

 12 ẑT,n := E z(T ) y iT n n i=1

r r r r r r r r r r j = 1 ,

 1 ..., n j = n + 1, ..., 2n j = 2n + 1, ..., 4n j = 4n + 1, ..., 8n

Figure 1 .

 1 Figure1. Illustration of the time point addition scheme in the construction of the martingales xj and zj (from[START_REF] Aalto | Convergence of discrete time Kalman filter estimate to continuous time estimate[END_REF]).

t 0 es 0 e

 00 As x ds + C t 0 A(s-r) B du(r) ds + w(t).

  s-r) Bdu(r) ds

2 L 2

 22 r-t) C * dr ds dt and, as with the first term, E ||(II)|| 2 Y = tr Cov [(II), (II)] ≤ h 4 12 tr(Q) Ce As AB (0,T ;L(U,Y)) .

2 L

 2 (D(A),Y) ||B|| 2 L(U,D(A)) tr(Q)E ||ẑ T,n -z(T )|| 2 X + incr [x,j+1]

M 1 (

 1 T ) = T 2 tr(Q) 2 min(eig(R)) ||C| | 2 L(D(A),Y) ||B|| 2 L(U,D(A)) E ||ẑ T,n -z(T )|| 2 X . M 2 (T ) = T 3 tr(Q) 3 min(eig(R)) H 2 T ||B|| 2 L(U,D(A)) E ||ẑ T,n -z(T )|| 2X , Now summing up (4.4) for K = 1, 2, ... yields the result.

  At L(X) ≤ µ. On the other hand, computing the integrals yields Now the idea is to bound the sum in (3.6) by using the first bound for small k and the latter for large k. Define the index n(h) := ⌈h -1/δ ⌉ for splitting the sum to get

	t t-h	e λ k s ds -	t	t+h	e λ k s ds ≤	4µ |λ k |	.

t t-h e λ k s ds -t+h t e λ k s ds ≤ h 2 sup s≥0 d ds e λ k s ≤ µh 2 |λ k |, since e

  ds + C

	and therefore
	d dt which is bounded for t, t i ∈ [0, T ], concluding the result. Cov [y(t i ), u(t)] = t i t Ce A(s-t) BQds, if t < t i , 0, if t ≥ t i
				t i	s	e A(s-r) B du(r) ds + w(t i )
			0	0
	= C	t i	e As x ds + C	t i	t i	e A(s-r) Bds du(r) + w(t i )
	0		0	r
	it holds that
	Cov [y(t i ), u(t)] = C	t∧t i		t i	e A(s-r) BQds dr
			0			r

  u u|| 2 and u are independent). The first term is now bounded as in one of the Theorems 3.1-3.4, so then remains the input noise induced term. To evaluate C h,u (t j+1 )u, note that for t ≥ t j+1 . In order to get a similar expression for t < t j+1 , just change t j+1 ↔ t in the bounds of the last two integrals and put minus signs in front them. Of course the last term is just Schwartz inequality. Now make a change of variables τ = s-t and note that for a fixed s ∈ [t j+1 , t j+1 + h], it holds that τ ∈ [s -t j+1 , s] ⊂ [0, T ], yielding finally

	and from this,
	E ||(I)|| 2 Y = tr Cov [(I), (I)]
	t ≤ 0	e A(t-s) Bdu(s) = tr(Q) 4 t j+1 0 t j+1 +h t j+1 0 t j+1 (t j+1 + h -s) Ce A(s-t) AB L(U,Y) ds e A(t-s) Bdu(s)	2	dt
	≤	h 3 12	tr(Q)		0	t j+1	+ t j+1 +h t t j+1 t j+1	A Ce A(s-t) AB s 0 e A(s-r) Bdu(r) ds + 2 L(U,Y) ds dt	t t j+1	Bdu(s)
	t t j+1 Bdu(s) = B(u(t) -u(t j+1 )). where the second inequality follows by Cauchy-E ||(I)|| 2 Inserting this to (4.1) gives C h,u (t j+1 )u Y ≤ h 4 2 tr(Q) Ce As AB L 2 (0,T ;L(U,Y)) 12 (4.3) where Ce As AB 2 L 2 (0,T ;L(U,Y)) ≤ H 2 L(U,D(A)) due to the admissibility T ||B|| 2 assumption.
	= -For the second term we have t j+1 t j+1 C 2 t j+1 -h t j+1 +h t t j+1 +h t j+1 +h 1 Cov [(II), (II)] = 4 t j+1 t	t	t j+1 +h
										A
		t j+1					t j+1
	These two terms are very similar by nature so it suffices to find a bound for
	one of them and use the same bound for both terms. Thus, let us consider
	the first part of the latter term, namely
		C 2	t j+1 t j+1 +h	t j+1 t A	0	s	e A(s-r) Bdu(r) ds dt
		=	C 2		t j+1 t j+1 +h	t j+1 t A	0	t j+1	e A(s-r) Bdu(r) ds dt
				+	C 2	t j+1 t j+1 +h	t j+1 t A	t j+1 s	e A(s-r) Bdu(r) ds dt
		=	C 2		0	t j+1	t j+1 +h t j+1 (t j+1 + h -s)Ae A(s-r) Bds du(r)
				+	C 2	t j+1 +h t j+1	r	t j+1 +h (t j+1 + h -s)Ae A(s-r) Bds du(r)
		= (I) + (II).
	Then
	Cov [(I), (I)] =	1 4		0	t j+1	t j+1 +h t j+1	t j+1 +h t j+1

Y (recall that

x t A s 0 e A(s-r) Bdu(r) ds -B u(t) -u(t j+1 ) dt + s 0 e A(s-r) Bdu(r) ds + B u(t) -u(t j+1 ) dt . (t j+1 + h -s)(t j+1 + h -r)

× Ce A(s-t) ABQB * A * e A * (r-t) C * dr ds dt,

This result extends to η -ν = 1/2 in which case the convergence rate is O(T

n -2 ln n).


Atef Masmoudi
email: atef.masmoudi@lirmm.fr

William Puech

Lossless chaos-based crypto-compression scheme for image protection

In this study, the authors proposed a new scheme which performs both lossless compression and encryption of images. Lossless compression is done by arithmetic coding (AC) while encryption is based on a chaos-based pseudorandom bit generator. Hence, they proposed to incorporate recent results of chaos theory into AC in order to shuffle the cumulative frequency vector of input symbols chaotically to make AC secure and the decoding process completely key-dependent. Many other techniques based on varying the statistical model used by AC have been proposed in literature, however, these techniques suffer from losses in compression efficiency that result from changes in entropy model statistics and are weak against known attacks. The proposed compression-encryption techniques were developed and discussed. The numerical simulation analysis indicates that the proposed scheme is highly satisfactory for image encryption without any AC compression efficiency loss. In addition, it can be incorporated into any image compression standard or algorithm employing AC as entropy coding stage, including static, adaptive and context-based adaptive models, and at any level, including bit, pixel and predictive error pixel levels.

Introduction

In this paper, the combination of compression and encryption by using arithmetic coding (AC) and chaotic systems for image protection was particularly given much interest. Several chaotic systems have been recently used for image encryption [START_REF] Li | Improving security of a chaotic encryption approach[END_REF][START_REF] Yang | A survey of chaotic secure communication systems[END_REF][START_REF] Wu | Analyzing and improving a chaotic encryption method[END_REF][START_REF] Zhang | An image encryption approach based on chaotic maps[END_REF][START_REF] Wong | An efficient diffusion approach for chaos-based image encryption[END_REF][START_REF] Masmoudi | A new image cryptosystem based on chaotic map and continued fractions[END_REF] and some of these novel chaotic systems have designed pseudorandom bit generators (PRBGs) for stream cipher applications [START_REF] Kanso | Logistic chaotic maps for binary numbers generations[END_REF][START_REF] Patidar | A novel pseudo random bit generator based on chaotic standard map and its testing[END_REF][START_REF] Masmoudi | An efficient prbg based on chaotic map and engel continued fractions[END_REF] since they have interesting cryptographic properties such as ergodicity, sensitivity to initial values and control parameters.

The AC has been widely used as an efficient entropy coding technique by most interesting image and video standards such as JBIG2, JPEG2000 and H.264/AVC, and by a variety of lossless data compression methods operating in bit- [START_REF] Marpe | Context-based adaptive binary arithmetic coding in the h.264/avc video compression standard[END_REF][START_REF] Alcaraz-Corona | Bi-level image compression estimating the Markov order of dependencies[END_REF], pixel- [START_REF] Carpentieri | Lossless compression of continuous-tone images[END_REF][START_REF] Masmoudi | Efficient adaptive arithmetic coding based on updated probability distribution for lossless image compression[END_REF][START_REF] Zhang | Segmentation of source symbols for adaptive arithmetic coding[END_REF][START_REF] Masmoudi | A new arithmetic coding model for a block-based lossless image compression based on exploiting inter-block correlation[END_REF] or predictive error pixel-level [START_REF] Motta | Lossless image coding via adaptive linear prediction and classification[END_REF][START_REF] Kuroki | Adaptive arithmetic coding for image prediction errors[END_REF][START_REF] Matsuda | A lossless coding scheme using adaptive predictors and arithmetic code optimized for each image[END_REF].

The AC is also considered for data encryption. Grangetto et al. [START_REF] Grangetto | Multimedia selective encryption by means of randomized arithmetic coding[END_REF] proposed a novel multimedia security framework by means of the randomised AC (RAC). This framework is based on a random organisation of the encoding intervals using a secret key. A chaos-based zeroth-order adaptive AC (AAC-0) [START_REF] Bose | A novel compression and encryption scheme using variable model arithmetic coding and coupled chaotic system[END_REF] technique was proposed. The statistical model is made varying in nature according to a coupled chaotic systems-based PRBG which makes AC compression efficiency loss of about 6%. Wen et al. [START_REF] Wen | Binary arithmetic coding using key-based interval splitting[END_REF] designed the binary AC (BAC) with key-based interval splitting. They proposed to use a key to split the interval associated with the symbol to be encoded. Thus, the traditional assumption in AC that a single continuous interval is used for each symbol is not preserved. The repeated splitting at each encoding/ decoding step allows both encryption and compression. Kim et al. [START_REF] Kim | Secure arithmetic coding[END_REF] demonstrated the insecurity of the interval splitting AC against a known plain-text attack and a chosen plain-text attack. They also provided an improved version called the secure AC applying a series of permutations at the input symbol sequence and output codeword. It should be noticed that because of the permutations process, it is difficult to extend the secure AC to the adaptive and context-based adaptive models which exploit the input symbol redundancy to encode messages. Mi et al. [START_REF] Mi | A novel chaotic encryption scheme based on arithmetic coding[END_REF] proposed a chaotic encryption scheme based on RAC using the logistic map for PRBG. Zhou et al. [START_REF] Zhou | Joint security and performance enhancement for secure arithmetic coding[END_REF] presented a scheme for joint security and performance enhancement of secure AC. They proposed to incorporate the interval splitting AC [START_REF] Wen | Binary arithmetic coding using key-based interval splitting[END_REF] scheme with the bit-wise XOR operation. This scheme can be extended to any adaptive model because of the elimination of the input symbol permutation step, and its implementation is lower in complexity than the original secure AC. A secure BAC scheme based on non-linear dynamic filter with changeable coefficients was presented and discussed in [START_REF] Li | A secure and efficient entropy coding based on arithmetic coding[END_REF]. In a previous work, Masmoudi et al. [START_REF] Masmoudi | A new joint lossless compression and encryption scheme combining a binary arithmetic coding with a pseudo random bit generator[END_REF] proposed a joint lossless compression and encryption scheme combining BAC with a PRBG. They proposed to randomly exchange the two intervals of the cumulative distribution vector which conserves the BAC compression efficiency while using the same probabilities with and without the exchanging intervals process. Duan et al. [27] introduced a RAC based on order-1 Markov model with a higher compression efficiency and good security. Recently, the proof of insecurity of the RAC based on Markov model has been proved and an improvement has been suggested in [28].

In this paper, we suggested a new scheme to achieve encryption by controlling certain operations in the AC compression algorithm using a secret key. The proposed scheme can be applied with any statistical model, including static, adaptive and context-based adaptive models, and to a variety of lossless data compression methods.

The rest of this paper is organised as follows. In Section 2, we briefly overview AC. Section 3 details the proposed AC-encryption scheme. In Section 4, we analyse the compression efficiency and the security of the proposed scheme and discuss the experimental results. Finally, our conclusion is discussed in Section 5.

Overview of AC

The AC [29, 30] is a very efficient statistical coding technique for data compression. The AC principle consists to assign one codeword to the entire input data symbols. To calculate the appropriate codeword for input data symbols, the AC works with a statistical model that estimates the probability of each symbol at the encoding/decoding process. The model used by AC can be static (SAC) [31,32], adaptive (AAC) [31,32] or context-based adaptive [33,[START_REF] Marpe | Context-based adaptive binary arithmetic coding in the h.264/avc video compression standard[END_REF]. In the static model, AC uses a fixed frequency count vector to encode all the symbols, however, in adaptive and context-based adaptive models, the AC works with a variable frequency count vector which is estimated/updated before/after encoding each symbol.

Let X be a sequence of m events X = x 1 , …, x m taking values from an alphabet of n symbols {α 1 , …, α n }. We denoted by Freq = [ f 1 , …, f n ], in Algorithm 1 (see Fig. 1) line 11, the frequency count vector where f k represents the frequency count of the corresponding symbol α k from the total length of the sequence X.

During the encoding process, AC firstly estimates the probability (frequency count) of each symbol and then calculates the cumulative frequency vector cum_freq, in Algorithm 1 (Fig. 1) line 12, by assigning a subinterval for each symbol α k with a size proportional to its frequency count in the interval [0, R) (R in Algorithm 1 Fig. 1) line 13). The cum_freq vector is calculated as

cum freq[a k ] = a k-1 k=1 f k (1)
Next, for any new symbol α k from the input sequence, AC selects the subinterval for α k and defines it as the new current interval. The AC iterates this step until all input sequence would be processed and finally generates the codeword that uniquely identifies the input data, which corresponds to any number within the final coding interval.

Let range be the size of the current subinterval, range in Algorithm 2 (see Fig. 2) line 12 and Algorithm 3 (see Fig. 3) line 15. The range is always equal to the product of the probabilities of the individual symbols encoded so far, and can be mathematically defined by [START_REF] Yang | A survey of chaotic secure communication systems[END_REF]. The 'range' value monotonically decreases at each iteration and therefore the number of bits needed to specify the founded interval of that length increase. It corresponds to the likelihood function given the observed data X in a specific model

M range = P M (X ) = m i=1 P M (x i ) (2)
If we let L(X, M ) denotes the likelihood function given the observed data X in the model M, then L(X, M ) = range. Therefore we can bound the number of bits required to represent the message X by the following equation [29]

-log 2 (L(X , M )) + 2 bits ( 3
)
where ⌊•⌋ is the floor function. The AC is therefore most useful when there are large probabilities in the corresponding model. The techniques, which secretly modifying the statistical model, suffer from losses in compression efficiency. In fact, our work which consists in randomly permuting the symbols' positions in the cumulative frequency count vector used by AC to make it secure and the decoding process completely key-dependent is within this framework. The proposed scheme can be seen as a coder-based encryption method.

In the encryption and decryption algorithms, we suggest to use new variable S called the symbols vector, in Algorithm 1 (Fig. 1) line 10, which is initially rearranged into the ascii order. For example, for grey-scale image data, where each pixel takes a value from {0, …, 255}, the symbols vector S is initially equal to S = [0, . . . , 255] Therefore, the cum_freq vector, used by AC during the encoding/decoding process, is initially equal to cum Freq = [0, Freq[S [START_REF] Li | Improving security of a chaotic encryption approach[END_REF]], Freq[S [START_REF] Li | Improving security of a chaotic encryption approach[END_REF]]

+ Freq[S[2]], . . . , n i=1

Freq[S[i]]

If we suggested swapping randomly the position of only two symbols in S, and the two symbols happen to be of the same probability, then nothing would change in the cum_freq vector and consequently all the other symbols would keep the same lower and upper bounds. In this case, the new symbol would be encoded without encryption and the generated bitstream remains the same with and without the swapping step which leads to a non-secured scheme. Similarly, if we wish to encrypt a data source and a part of its distribution It is also possible to randomly choose a position j and place its corresponding symbol at this position S[ j] in the beginning of the vector S and simultaneously shift all the other symbols S [START_REF] Li | Improving security of a chaotic encryption approach[END_REF], …, S[ j -1] to the right. This allows the modification of the cum_freq vector. However, all the symbols that were initially situated at on the right of S[ j], which are S[ j + 1], …, S[n], would keep the same low and high limits. Therefore, if the next symbol to encoded is within S[ j + 1], …, S[n], it will be encoded without encryption even after modifying the cum_freq vector. Thus, the adequate solution would be to apply a permutation of all the symbols positions at each encoding/decoding step, The idea presented in this section consists in chaotically reordering the symbols' positions in S, in order to randomly modify the lower and upper bound of all symbols prior to performing the AC encoding/decoding process. In fact, we propose to calculate a random position j to partition the symbols vector S in two parts to be interchanged. Thus, S becomes after permutaion step

S = [S[j], S[j + 1], . . . , S[n], S[1], . . . , S[j -1]
and the corresponding cum_freq vector becomes cum freq

= [0, Freq[S[j]], Freq[S[j]] + Freq[S[j + 1]], ..., × n i=j Freq[S[i]] + Freq[S[1]], ..., n i=1

Freq[S[i]]]

Thanks to this solution, we chaotically modify the lower and upper bounds of each symbol which introduces randomisation in the AC procedure. The position j is calculated from a chaos-based PRBG, which leads to perform maximum randomisation in the proposed encryption approach. We also propose to use two other chaos-based PRBGs in the proposed compression-encryption scheme. The first PRBG is used to changing dually the S and cum_freq vectors in order to reduce the additional processing time needed by updating the cum_freq vector after the random organisation of the symbols vector S, and to provide well protected results. The second PRBG is used to avoid XOR operation of output bits to enhance security and to make the proposed AC-encryption scheme more resistant against known attacks.

In this work, we have used, with no loss of generality, the logistic map (LM) as chaos-based PRBG, since it is fast, easy to implement and widely used to design chaotic-based cryptographic system. The known one-dimensional LM is defined as

y n+1 = ly n * (1 -y n )
where λ ∈ [0, 4], the control parameter, and y n ∈ (0, 1). The true chaotic behaviour of LM appears for λ = 4, this is why we suggest to use this value when designing LM-based PRBGs. It should be noted that multiple chaotic systems, high-dimensional chaotic systems, multiple iterations of chaotic systems and many other techniques can be used to design PRBG.

In this work, there are three PRBGs used to carry out the dual random organisation and the XOR operation: † PRBG 1 : randomly generates a controlling bit to control whether the symbols vector S and the cum_freq vector are modified or not. † PRBG 2 : randomly generates a value j that partition the symbols vector S in two parts to be interchanged. † PRBG 3 : randomly generates a bit to do XOR operation with the AC output bits.

Three initial values (or seeds) y 0 1 , y 0 2 and y 0 3 [Algorithm 1 (Fig. 1) lines 1, 2 and 3, repectively] are needed to initialise the three PRBGs. Therefore the proposed AC-encryption scheme works with a key K = (y 0 1 , y 0 2 , y 0 3 ) composed of three floating value numbers. Given the same K, both the encoder and the decoder generate the same controlling bit sequence. However, as long as one of the three key components is unknown or incorrectly given, the decoder cannot decode the compressed and encrypted data properly and the decoded data will be completely useless.

The encryption and decryption procedures are illustrated in Algorithms 2 and 3 (Figs. 2 and3), respectively. Algorithm 1 presents the Initialise_AC() function used both by the encoder and the decoder to initialise all necessary variables. For more details about the AC implementation used in this work, please refer to [31]. In addition, AC is very sensitive to errors in the compressed data, and this undesired property ameliorates the security of the encryption-based AC method.

The proposed method for random interchanging the two parts of the symbols vector is very important as it guarantees the modification of the lower and upper bounds of all the symbols in the cum_freq vector.

Our compression-encryption scheme can be applied at bit-, pixel-or predictive error pixel-level. Even in the case of a binary source, we do not really need to omit the PRBG 2 generator as the permutation position for a binary source noted as j will usually be equal to 2 which guarantees the modification of the lower and upper bounds of the two binary symbols. Therefore our encryption algorithm remains the same whether the source is binary or not.

Fig. 4a presents a short example to illustrate the swapping position procedure of the input symbols and their corresponding frequencies. The conventional SAC encoding of the sequence 'CACBCC', obtained from an alphabet {A, B, C} with known probabilities {0.2, 0.3, 0.5}, is depicted in Fig. 4b. The encoding of the same sequence using our approach is shown in Fig. 4c. The dual swapping of the cum_freq vector is determined according to two random and secret keystreams '010110' and ' -2 -23 -', which are generated by PRBG 1 and PRBG 2 , respectively. It should be noticed that, when randomly interchanging the symbols' positions in S, we need to calculate the cum_freq vector according to these new positions, prior to performing the AC encoding process.

It should be noted that the proposed scheme conserves the compression efficiency of the AC because we use the same probabilities when encoding the input sequence with and without the dual swapping process described above. In addition, Figs 4b andc show that AC encodes the input sequence with and without encryption in two different intervals, but with the same size.

Experimental results

Analysis of compression efficiency

We propose to analyse the compression efficiency of our method on some well-known 8-bit grey-scale images, which are available from the University of Waterloo Greyset2 collection (http://links.uwaterloo.ca/repository.html), both on static and adaptive order-0 models. These 12 images range in size from 464 × 352 to 672 × 498 pixels.

To verify the effect of the interchanging symbols' positions procedure in compression efficiency, Table 1 gives the compression results of all test images, both on SAC and AAC-0.

When encoding an input sequence with SAC, the exact frequency count vector of all the existing symbols will be used. Consequently, the likelihood function given an image X and the static model based on true symbols probabilities www.ietdl.org

(called P) is calculated as L(X , P) = m i=1 P(x i ) = 255 k=0 f (k) m f (k) (4)
However, AAC-0 assumes that all pixels have initially an equal frequency count, and after encoding/decoding each new pixel, it updates this pixel occurrence by adding one to its frequency count. Let U positive integers, so that n = ba + 1 and each of the n values has an equal probability (1/n). The probability of observing an image X whith respect to the model used by AAC-0, which corresponds to the likelihood function given that image, is calculated as

L(X , U [0,255] ) = m i=1 P(x i |x 1 , . . . , x i-1 ) = 255 k=0 f (k) j=1 (j) m i=1 (n + i -1) ( 5
)
where m is the total image size and f is the vector of true frequency counts. The number of bits required to represent the message X, either in SAC or AAC-0, can be bounded according to (3). The proposed encryption scheme has the effect of finding randomly a position j between all possible positions in the symbols vector, which will be subdivided in two parts, and then swapping the order of these two parts both in symbols vector and cum freq vector. The swapping procedure does not make any change to the subinterval width associate with each symbol, thereby allowing encryption without any loss in coding efficiency. It should be noticed that during experiments, we used an integer AC implementation with 32-bit precision. In this case, changing the frequency counts positions in cum_freq vector may cause a slight difference in AC encoding and decoding procedures. From Table 1, the obtained compression results, both on SAC and AAC-0, with and without encryption were slightly different but there is no loss in compression efficiency.

In addition, Figs. 5-7 show three grey-scale images and their corresponding compressed-encrypted images in grey-scale form. The compressed-encrypted images are obtained by using the key K, see Table 2, both on SAC and AAC-0 models.

Security analysis

The security analysis of the proposed AC-encryption scheme includes the key sensitivity test, the statistical analysis, the entropy analysis and the key space analysis. Table 2 lists four different secret keys used for the security analysis.

Statistical analysis:

An ideal encryption algorithm should resist statistical attacks. This has been shown using the histogram analysis, correlation analysis and information entropy analysis of the proposed scheme for both SAC and AAC-0.

(a) The histogram analysis: We have analysed the histograms of all the test images and their corresponding compressed-encrypted images using the key K. Fig. 8 presents the histograms of three original images. We presented the histogram analysis of the Lena, France and Frog images in Figs. 9, 10 and 11, respectively. It should be noticed that the unit of output acts as byte and therefore the range of output is 0 -255.

From Figs. 9a andb, we can see that for Lena image, the histogram of the AC output, both on SAC and AAC-0, has uniform distribution. However, for some images, such as France and Frog images, the output value of SAC has a non-uniform distribution, as illustrated in Figs. 10a and11a, respectively. Similar results were obtained for AAC-0. The histograms of the compressed-encrypted outputs were also calculated. From Figs. 9c, 10c and11c, we can see that the histograms of the compressed-encrypted outputs for SAC have always uniform distribution. We have found similar results for the histogram analysis of the all test images. (b) The χ 2 test: We have also computed the χ 2 values of both AC output and ciphertexts, on SAC and AAC-0. The χ 2 value is described as follows where i is the number of grey levels (0, …, 255), f i is the observed frequency count of each grey level and n c is the length, in bytes, of the obtained bitsream when applying AC. Assuming a significant level of 0.05, then χ 2 (255, 0.05) = 293. The χ 2 value for all test images and their corresponding bitstream with and without encryption both on static and zero-order adaptive models are presented in Table 3. The results shown in boldface are the worst. We clearly show that all the observed values of the ciphertext, both on static and adaptive models, are less than 293. This implies that the distribution of the compressed-encrypted histogram is uniform which points out good quality of the AC-encryption scheme. Therefore the encrypted image does not provide any clue to employ any statistical attack on the proposed image encryption scheme, which makes any statistical attacks difficult.

x 2 = f i - n c 256 2 / n c
(c) The correlation test: For real images, each pixel is usually highly correlated with its adjacent pixels either in horizontal, vertical or diagonal directions. However, an efficient encryption scheme should generate a ciphertext with low correlation between adjacent values [START_REF] Wu | Analyzing and improving a chaotic encryption method[END_REF]. Horizontal, vertical and diagonal correlation coefficients γ xy of two adjacent pixels can be calculated using the following equation where N is the total number of duplets (x i , y i ) obtained from the image, and E(x) and E(y) are the mean values of x and y, respectively. Table 4 shows the correlations of the original Lena image and its corresponding ciphertext found using the key K. Similar results for all test images were obtained.

E(x) = 1 N N i=1 x i D(x) = 1 N N i=1 x i -E(x) 2 cov(x, y) = 1 N N i=1 x i -E(x) y i -E(y) g xy = cov(x, y) D(x) √ D(y) √ (6)
In addition, Figs. 12 and 13 show the correlation distribution of 1000 pairs of adjacent pixels randomly selected from the original images of France and Frog, and that in their compressed-encrypted images produced using the key K.

Fig. 14 reports the correlation achieved by a decoder which attempts to decode the Lena image using slightly different initial values, of which only one is correct. It can be seen that a very poor correlation is obtained among the decoded images obtained using slightly different initial values.

Both the correlation coefficients and the figures show that the proposed AC-encryption scheme decorrelates effectively the neighbouring pixels of the plain images, either in SAC or AAC-0.

When an image is encrypted, its entropy, expressed in bits/ pixel, should ideally be equal to 8 bits/pixel. If the output of an image encryption algorithm [35] generates symbols with entropy less than 8 bits/pixel, there exists certain degree of predictability and the encryption algorithm will be threatened in its security. Table 5 presents the average entropy obtained for the ciphertexts of all test images, both on SAC and ACC-0, which is equal to 7.9990 and 7.9991 bits/pixel, respectively. The average entropy found is very close to the theoretical value ≃ 8, which confirms that the proposed AC-encryption scheme is secure against the entropy attack. In Table 5, the entropies shown in boldface represent the worst results, which confirms that the AC cannot be considered as an encryption scheme.

Key space:

From a strong cryptographic point of view, the data security should be based on an algorithm where a brute force attack has a minimum of 2 128 different keys that can be used in the encryption [START_REF] Schneier | Applied cryptography: protocols, algorithms, and source code in C[END_REF]. In our scheme, the secret key K includes three floating point numbers with precision of 10 14 , and so the key space of the key is 10 42 ≃ 2 140 which is far enough to resist the brute-force attack.

Key sensitivity test:

According to the basic principle of cryptology, a cryptosystem should be sensitive to the key. Thus, we propose the following tests [START_REF] Maniccam | Lossless image compression and encryption using SCAN[END_REF]:

(a) Assume that the encryption key used is K presented in Table 2. First, the original image of Lena is encrypted using The unit of output acts as byte, so the range of output is 0-255 the test key. Next, the same original image is encrypted with three slightly different keys described in Table 2.

We propose for each of the used secret keys, to conserve two parts of K and to change the third one by 10 -14 . Suppose C 1 and C 2 are two ciphertexts having the same length N, then the number of pixel change rate (NPCR) and the unified average changing intensity (UACI) [START_REF] Chen | A symmetric image encryption scheme based on 3D chaotic cat maps[END_REF][START_REF] Masmoudi | A new image cryptosystem based on chaotic map and continued fractions[END_REF] can be mathematically defined by ( 9) and [START_REF] Marpe | Context-based adaptive binary arithmetic coding in the h.264/avc video compression standard[END_REF], respectively, where |.| denotes the absolute value function. The NPCR and UACI between various ciphertexts produced using slightly different keys on SAC, are calculated and the results are given Table 6. Similar results were found when using AAC-0. Table 6 shows that ciphertexts obtained with keys K and K 1 have 99.63% pixels which are not identical.

It should be noticed that the NPCR value should be as close to 100% as possible, and the UACI value should be as close to 33% as possible

D i = 0, if C 1 i = C 2 i 1, if C 1 i = C 2 i (8) NPCR = i D i N × 100% (9) UACI = 1 N i |C 1 i -C 2 i | 255 × 100% (10)
(b) In addition, to test the key sensitivity of our encryption algorithm, we propose to encrypt plain Lena image with the key K and to decrypt it using a key with a difference of 10 -14 . Fig. 15 shows the experimental results. As observed from Fig. 15, the decrypted images are totally different from the plain one even where there is only 10 -14 difference in the decryption key. Thus, having a perfect approximation of the encryption secret key makes decryption impossible.

4.2.4

Comparison: Tables 7 and8 present the comparison between the proposed method and the other methods in terms of encryption and compression efficiency, respectively. From Table 7, its is shown that most of the other researchers have not reported the histogram, χ 2 , correlation, entropy, NPCR and UACI tests, and suffer from vulnerability to some attacks. In addition, most of the other techniques are based on BAC (each data symbol is coded bit by bit) by swapping the two mapping intervals randomly using different methods. It should be noted that the other methods are restricted to be applied only to some types of statistical models.

Conclusion

In this paper, we proposed an efficient and secure encryption scheme which combines AC, both on static and adaptive models, with chaos-based PRBGs to perform both lossless compression and encryption of images. In our scheme, we exploit both the efficiency of the AC in lossless image compression and the advantages of the chaos theory in data encryption to provide a scheme which can be very useful in many applications, especially for image protection. The proposed scheme presents several interesting features like: extendable to any context-based AAC approach, lossless encryption of images, the same sizes for the output bitstream and the ciphertext, for both static and adaptive models, key space is large enough, extendable to any PRBG, fast and easily implemented at a low cost. expected value uniform <293

≃ 0 ≃ 8 >99.58 [START_REF] Wu | Npcr and uaci randomness tests for image encryption[END_REF] 33.37 < . < 33.55 [START_REF] Wu | Npcr and uaci randomness tests for image encryption[END_REF] The missing values are represented by the symbol NA (not available).

Fig. 1

 1 Fig. 1 Initialise_AC()

Fig. 2

 2 Fig. 2 AC-encryption algorithm

Fig. 3

 3 Fig. 3 AC-decryption algorithm

Fig. 4

 4 Fig. 4 Swapping position procedure of the input symbols and their corresponding frequencies a Symbols and their corresponding probability positions are exchanged in the same order and are key-dependent (y 1 = 1 and j = 2) Interval partitions for b Conventional SAC c AC with the proposed dual permutation scheme of symbol positions

Fig. 6 Fig. 5 Fig. 7

 657 Fig. 6 Three grey-scale images of France a Original image of France and its corresponding compressed-encrypted images using the key K both on b SAC c AAC-0

Fig. 9

 9 Fig. 9 Histogram analysis of the AC outputs with and without encryption both on SAC and AAC-0 for original image of Lena Histogram of a Output of SAC b Output of AAC-0, compressed-encrypted output obtained by using the key K on c SAC d AAC-0

256Fig. 10

 10 Fig. 10 Histogram analysis of the AC outputs with and without encryption both on SAC and AAC-0 for original image of France Histogram of a Output of SAC b Output of AAC-0, compressed-encrypted output obtained by using the key K on c SAC d AAC-0

Fig. 11

 11 Fig. 11 Histogram analysis of the AC outputs with and without encryption both on SAC and AAC-0 for original image of Frog Histogram of a Output of SAC b Output of AAC-0, compressed-encrypted output obtained by using the key K on c SAC d AAC-0

  (d) The information entropy test: The information entropy is initially defined by Shannon [34]. To calculate the entropy H of a source of length n, we should use the

Fig. 12

 12 Fig. 12 Original images of France a Correlations of 1000 pairs of adjacent pixels randomly selected from the original image of France and their corresponding compressed-encrypted images obtained by b SAC c AAC-0 The unit of output acts as byte, so the range of output is 0-255

Fig. 13

 13 Fig. 13 Original images of Frog a Correlations of 1000 pairs of adjacent pixels randomly selected from the original image of Frog and their corresponding compressed-encrypted images obtained by b SAC c AAC-0 The unit of output acts as byte, so the range of output is 0-255

Fig. 14

 14 Fig. 14 Correlation test obtained by decoding the compressed-encrypted Lena image with different combinations of the keys Generated using different initial values a y 0 1 b y 0 2 c y 0 3 on AAC-0 The plots are measured against initial values obtained from the key K

Fig. 15

 15 Fig. 15 Decrypted image with a K b K 1 c K 2 d K 3

Table 1

 1 Compression results (in byte) of different images both on static and adaptive order-0 models

	Image	SAC	SAC-encryption	AAC-0	AAC0-encryption
	Lena	278	244 278	231 622	231 621
	Goldhill	269	245 269	237 882	237 882
	Barbara	944	244 945	238 748	238 749
	Boat	705	233 705	227 049	227 050
	Zelda	565	238 564	232 816	232 816
	Mandrill	387	241 387	237 049	237 050
	Peppers	285	248 284	238 376	238 376
	France	626	261 623	258 217	258 216
	Mountain	836	239 837	232 293	232 294
	Library	564	119 564	112 108	112 109
	Washsat	95	95 286	93 071	93 071
	Frog	203	193 201	180 982	180 981

Table 2

 2 Keys used for security analysis

		K	K 1	K 2	K 3
	y 0 1 y 0 2 y 0 3	0.23951648742195 0.54397486939831 0.83215648972136	0.23951648742194 0.54397486939831 0.83215648972136	0.23951648742195 0.54397486939832 0.83215648972136	0.23951648742195 0.54397486939831 0.83215648972135

Table 3 χ

 3 2 values for test images of the AC output and

	ciphertext both on static and adaptive models
		SAC	SAC-encryption AAC0 AAC0-encryption
	Lena	272	250	260	260
	Goldhill	293	290	250	251
	Barbara	264	268	261	232
	Boat	260	238	270	257
	Zelda	223	216	271	249
	Mandrill	242	266	264	244
	Peppers	256	284	309	241
	France	18101	259	24734	245
	Mountain	3859	283	447	261
	Library	51292	269	510	290
	Washsat	240	222	259	278
	Frog	1348	272	1038	258

Table 4

 4 Correlation of the original Lena image and its corresponding ciphertext found using the key K

		Direction	Correlation coefficients
	original	horizontal	0.9714
	images	vertical	0.9725
		diagonal	0.9531
	output of SAC	horizontal	0.0156
	ciphertext of SAC	horizontal	0.0003
	output of AAC-0	horizontal	0.0054
	ciphertext of AAC-0	horizontal	0.0009

Table 5

 5 Entropy of the AC outputs with and without

	encryption for both static and adaptive models
		SAC	SAC-encryption AAC0 AAC0-encryption
	Lena	7.9992	7.9991	7.9991	7.9992
	Goldhill	7.9992	7.9993	7.9992	7.9992
	Barbara	7.9992	7.9992	7.9993	7.9993
	Boat	7.9993	7.9991	7.9993	7.9993
	Zelda	7.9991	7.9991	7.9992	7.9993
	Mandrill	7.9992	7.9993	7.9993	7.9993
	Peppers	7.9992	7.9993	7.9992	7.9993
	France	7.9648	7.9994	7.9573	7.9994
	Mountain 7.9962	7.9993	7.9960	7.9992
	Library	7.8524	7.9983	7.9858	7.9985
	Washsat	7.9982	7.9981	7.9981	7.9980
	Frog	7.9921	7.9991	7.9932	7.9990
	average	7.9832	7.9990	7.9937	7.9991

Table 6

 6 Difference between ciphertexts obtained by using

	keys with slight differences
	Test item	Test results between ciphertexts with tiny
			change in the key
		k 1	k 2	k 3
	NPCR, %	99.63	99.52	99.61
	UACI, %	33.42	33.51	33.48

Table 7

 7 Comparison of encryption efficiency of the proposed method and the other methods

	Methods	Histogram	χ 2 Correlation	Entropy	Key sensitivity	Key	Proof of
								space	insecurity
						NPCR, %	UACI, %
	Grangetto et al. [19] Bose et al. [20]	NA NA	NA NA	NA NA	NA NA	NA NA	NA NA	NA 2 30	[39]
	Wen et al. [21]	NA	NA	NA	NA	NA	NA	NA	[39, 40]
	Kim et al. [22]	NA	NA	NA	>0.99 (binary	NA	NA	NA	[41, 42]
	Mi et al. [23] Zhou et al. [24] Li et al. [25]	NA uniform uniform	NA NA NA	NA NA NA	data) NA NA >0.99 (binary	NA NA 50	NA NA NA	2 32 2 128 2 199
	Masmoudi et al. [26]	passed all NIST	2 157		data)
	Duan et al. [27]	statistical tests [43] NA	NA	NA	>0.99 (binary	99.60	NA	2 256	[28]
	proposed method	uniform		<0.0009	data) > 7.999	99.59	33.47	2 140
			<272

Table 8

 8 Comparison of compression efficiency of the proposed method and the other methodsThe symbol 'any' in the Data type column corresponds to binary and non binary data. However, in the 'statistical model' column it corresponds to static, adaptive and context-based adaptive models.

	Methods	Loss in, %	Data	Statistical model
			type
	Grangetto et al.	0	binary	MQ-coder for
	[19]			JPEG 2000
	Bose et al. [20]	6	any	zeroth-order
	Wen et al. [21]	<1	binary	static
	Kim et al. [22]	<1	binary	static
	Mi et al. [23]	0	binary	static
	Zhou et al. [24]	<1	binary	any
	Li et al. [25]	0	binary	any
	Masmoudi	0	binary	any
	et al. [26]
	Duan et al. [27]	improve RAC	binary	order-1 Markov
		[19] by 33%		model
	proposed	0	any	any
	method
	expected value	0	any	any

& The Institution of Engineering and Technology 2014

& The Institution of Engineering and Technology 2014 IET Image Process., pp. 1-16 doi: 10.1049/iet-ipr.2013.0598