# Adaptive Algorithm for Multichannel Autoregressive Estimation in Spatially Correlated Noise

Authors: Alimorad Mahmoudi
Journal: Journal of Stochastics
 (2014)

---

## Abstract

This paper addresses the problem of multichannel autoregressive (MAR) parameter estimation in the presence of spatially correlated noise by steepest descent (SD) method which combines low-order and high-order Yule-Walker (YW) equations. In addition, to yield an unbiased estimate of the MAR model parameters, we apply inverse filtering for noise covariance matrix estimation. In a simulation study, the performance of the proposed unbiased estimation algorithm is evaluated and compared with existing parameter estimation methods.


---

## Body

## 1. Introduction

The noisy MAR modeling has many applications such as high resolution multichannel spectral estimation [1], parametric multichannel speech enhancement [2], MIMO-AR time varying fading channel estimation [3], and adaptive signal detection [4].When the noise-free observations are available, the Nuttall-Strand method [5], the maximum likelihood (ML) estimator [6], and the extension of some standard schemes in scalar case to multichannel can be used for estimation of MAR model parameters. The relevance of the Nuttall-Strand method is explained in [7] by carrying out a comparative study between these methods. The noise-free MAR estimation methods are sensitive to the presence of additive noise in the MAR process which limits their utility [1].The modified Yule-Walker (MYW) method is a conventional method for noisy MAR parameter estimation. This method uses estimated correlation at lags beyond the AR order [1]. The MYW method is perhaps the simplest one from the computational point of view [8], but it exhibits poor estimation accuracy and relatively low efficiency due to the use of large-lag autocovariance estimates [8]. Moreover, numerical instability issues may occur when it is used in online parameter estimation [9].The least-squares (LS) method is another method for noisy MAR parameter estimation. The additive noise causes the least-squares (LS) estimates of MAR parameters to be biased. In [10] an improved LS (ILS) based method has been developed for estimation of noisy MAR signals. In this method, bias correction is performed using observation noise covariance estimation.The method proposed in [10], denoted by vector ILS based (ILSV) method, is an extension of Zheng’s method [11] to the multichannel case. In the ILSV method, the channel noises can be correlated and no constraint is imposed on the covariance matrix of channel noises. Nevertheless this method has poor convergence when the SNR is low.In [12], the ILSV algorithm is modified using symmetry property of the observation covariance matrix which is named an advanced least square vector (ALSV) algorithm. One step called symmetrisation is added to the ILSV algorithm to estimate the observation noise covariance matrix.In [13, 14], methods based on errors-in-variables and cross coupled, Kalman and H
∞, filters are suggested, respectively. These methods are also an extension of the scalar methods presented in [15] and [16, 17], respectively. In these methods, the observation noise covariance matrix is assumed to be a diagonal matrix. This assumption means that the channel noises are uncorrelated. These methods are not suitable for spatially correlated noise case.Note that cross coupled filter idea consists of using two mutually interactive filters. Each time a new observation is available, the signal is estimated using the latest estimate value of the parameters, and conversely the parameters are estimated using the latest a posteriori signal estimate.In this paper, we combine low-order and high-order YW equations in a new scheme and use adaptive steepest descent (SD) algorithm to estimate the parameters of the MAR process. We also apply the inverse filtering idea for observation noise covariance matrix estimation. Finally we can theoretically perform the convergence analysis of the proposed method. Using computer simulations, the performance of the proposed algorithm is evaluated and compared with that of the LS, MYW, ILSV, and ALSV methods in spatially correlated noise case. Note that we do not do comparisons with [13, 14] because they assume that the channel noises are uncorrelated.The proposed estimation algorithm provides more accurate parameter estimation. The convergence behaviour of the proposed method is better than those of the ILSV and the ALSV methods.


## 2. Data Model

A real noisym-channel AR signal can be modeled as
(1)
z
(
n
)
=
A
1
z
(
n
-
1
)
+
⋯
+
A
p
z
(
n
-
p
)
+
e
(
n
)
,
y
(
n
)
=
z
(
n
)
+
w
(
n
)
,
where e
(
n
) is the m
×
1 input noise vector and z
(
n
) is the m
×
1 MAR output noise vector. A
i’s (
i
=
1,2
,
…
,
p
) are the m
×
m MAR parameter matrices, and y
(
n
) is the m
×
1 observed noise vector. The noises e
(
n
) and w
(
n
) are assumed to be zero-mean stationary white noises with the m
×
m covariance matrices given by
(2)
R
e
(
l
)
=
E
{
e
(
n
)
e
(
n
-
l
)
T
}
=
Σ
e
δ
l
,
R
w
(
l
)
=
E
{
w
(
n
)
w
(
n
-
l
)
T
}
=
Σ
w
δ
l
,
where δ
l is the Kronecker delta and Σ
e and Σ
w are the unknown input noise covariance matrix and the observation noise covariance matrix which are generally nondiagonal. This means that the observation noise is spatially correlated. E
{
·
} is the expectation operator and T denotes the transpose operation.In this data model, we have the following assumptions.(i)
The AR orderp and the number of channels m are known.
(ii)
The MAR parameter matricesA
i are constrained such that the roots of
(3)
det
⁡
{
A
p
(
z
)
}
=
0lie inside the unit circle, where A
p
(
z
)
=
I
m
-
∑
i
=
1
p
A
i
z
-
i, I
m is the m
×
m identity matrix, and det denotes determinant of a matrix. This condition on the matrices A
i guarantees the stability of the MAR model [1]. Note that the roots of (3) are poles of the MAR model.(iii)
The observation noisew
(
n
) is uncorrelated with the input noise e
(
n
); that is, E
{
w
(
n
)
e
(
k
)
T
}
=
0 for all n and k.Using (1), we can write the following linear regression for y
(
n
):
(4)
y
(
n
)
=
A
T
g
(
n
)
+
ε
(
n
)
,
where A is the m
p
×
m parameter matrix defined as
(5)
A
=
[
A
1
⋯
A
p
]
T
and g
(
n
) is the m
p
×
1 regression vector defined as
(6)
g
(
n
)
=
[
y
(
n
-
1
)
T
⋯
y
(
n
-
p
)
T
]
T
.
The m
×
1 equation error vector ε
(
n
) is
(7)
ε
(
n
)
=
-
A
T
n
(
n
)
+
w
(
n
)
+
e
(
n
)
,
where
(8)
n
(
n
)
=
[
w
(
n
-
1
)
T
⋯
w
(
n
-
p
)
T
]
T
.The covariance matrices ofε
(
n
) and y
(
n
) are given as
(9)
R
ε
(
l
)
=
E
{
ε
(
n
)
ε
(
n
-
l
)
T
}
=
∑
j
=
l
p
A
j
Σ
w
A
j
-
l
T
+
Σ
e
δ
l
;
0
≤
l
≤
p
,
(10)
R
y
(
l
)
=
E
{
y
(
n
)
y
(
n
-
l
)
T
}
=
R
z
(
l
)
+
Σ
w
δ
l
,
where A
0
=
-
I
m.The main objective in this paper is to estimateA, Σ
w, and Σ
e using the given samples {
y
(
1
)
,
y
(
2
)
,
…
,
y
(
N
)
}.


## 3. Steepest Descent Algorithm for Noisy MAR Models

The proposed estimation algorithm uses steepest descent method. It iterates between the following two steps until convergence: (I) given noise covariance matrix, estimate MAR matrices using steepest descent method; (II) given MAR matrices, estimate noise covariance matrix based on inverse filtering idea.


### 3.1. MAR Parameter Estimation Given Observation Noise Covariance Matrix

The Yule-Walker equations for the noisy MAR process given by (4) are
(11)
R
y
(
l
)
-
∑
i
=
1
p
A
i
R
y
(
l
-
i
)
+
∑
i
=
1
p
A
i
Σ
w
δ
l
-
i
=
0
,
l
≥
1
.
Evaluating (11) for l
=
1
,
…
,
p and using R
y
(
l
)
=
R
y
(
-
l
)
T give the following system of linear equations:
(12)
(
R
y
-
R
w
)
A
-
r
y
=
0
,
where
(13)
R
y
=
[
R
y
(
0
)
R
y
(
1
)
⋯
R
y
(
p
-
1
)
R
y
(
-
1
)
R
y
(
0
)
⋯
R
y
(
p
-
2
)
⋮
⋮
⋱
⋮
R
y
(
1
-
p
)
R
y
(
2
-
p
)
⋯
R
y
(
0
)
]
,
r
y
=
[
R
y
(
1
)
R
y
(
2
)
⋯
R
y
(
p
)
]
T
,
and R
w is a block-diagonal m
p
×
m
p matrix given by
(14)
R
w
=
[
Σ
w
0
⋯
0
0
Σ
w
⋮
⋮
⋮
0
⋱
0
0
⋯
⋯
Σ
w
]
=
Σ
w
⊗
I
p
.Evaluating (11) for l
=
p
+
1
,
…
,
p
+
q and using R
y
(
l
)
=
R
y
(
-
l
)
T give the following modified Yule-Walker (MYW) equations:
(15)
R
y
q
A
-
r
y
q
=
0
,
where
(16)
R
y
q
=
[
R
y
(
-
p
)
R
y
(
-
p
+
1
)
⋯
R
y
(
-
1
)
R
y
(
-
p
-
1
)
R
y
(
-
p
)
⋯
R
y
(
-
2
)
⋮
⋮
⋱
⋮
R
y
(
-
p
-
q
)
R
y
(
-
p
-
q
+
1
)
⋯
R
y
(
-
q
)
]
,
r
y
q
=
[
R
y
(
p
+
1
)
R
y
(
p
+
2
)
⋯
R
y
(
p
+
q
)
]
T
,
and q is an arbitrary nonnegative integer.Combining the low-order Yule-Walker equations given by (12) and the high-order Yule-Walker equations given by (15), we obtain
(17)
H
y
A
-
H
w
A
-
h
y
=
0
,
where
(18)
H
y
=
[
R
y
R
y
q
]
,
h
y
=
[
r
y
r
y
q
]
,
(19)
H
w
=
[
R
w
0
]
.Premultiplying (17) by (
H
y
T
H
y
)
-
1
H
y
T, we obtain
(20)
A
-
A
H
-
R
H
H
w
A
=
0
,
where
(21)
A
H
=
(
H
y
T
H
y
)
-
1
H
y
T
h
y
,
R
H
=
(
H
y
T
H
y
)
-
1
H
y
T
.
And R
H of dimension m
p
×
m
(
p
+
q
) can be partitioned as
(22)
R
H
=
[
R
H
1
R
H
2
]
,
where the dimensions of R
H
1 and R
H
2 are m
p
×
m
p and m
p
×
m
q, respectively.Substituting (19) and (22) into (20), we get
(23)
A
-
A
H
-
R
H
1
R
w
A
=
0
.
If R
w is assumed to be known, we can iteratively estimate A using steepest descent algorithm as follows:
(24)
A
(
i
)
=
A
(
i
-
1
)
+
μ
(
A
H
+
R
H
1
R
w
A
(
i
-
1
)
-
A
(
i
-
1
)
)
,
where μ is the step size parameter and A
H, R
H
1 can be estimated from the observations. We can also control the stability and rate of convergence of the algorithm by changing μ. The above equation is converged if μ is selected between zero and 2
/
(
1
-
λ
min
⁡
) where λ
min
⁡ is the maximum eigenvalue of the R
H
1
R
w matrix. Furthermore if we set μ
=
1 and q
=
0 in (24), it is changed to
(25)
A
(
i
)
=
A
l
s
+
R
y
-
1
R
w
A
(
i
-
1
)
,
where A
l
s
=
R
y
-
1
r
y, which is previously derived in [8].


### 3.2. Application of Inverse Filtering for Estimation of Observation Noise Covariance Matrix Given MAR Matrices

If the observed process is filtered by MAR matrices, the output process of the filter is obtained as(26)
ε
(
n
)
=
y
(
n
)
-
A
T
g
(
n
)
.
By evaluating (9) for l
=
1, we can obtain
(27)
R
ε
(
1
)
=
∑
j
=
1
p
A
j
Σ
w
A
j
-
1
T
.
We can rewrite (27) as
(28)
(
∑
j
=
1
p
A
j
⊗
A
j
-
1
)
vec
⁡
(
Σ
w
T
)
=
vec
⁡
(
R
ε
(
1
)
T
)
.
We solve (28) as follows:
(29)
vec
⁡
(
Σ
w
T
)
=
(
∑
j
=
1
p
A
j
⊗
A
j
-
1
)
-
1
r
ε
,
where r
ε
=
vec
⁡
(
R
ε
(
1
)
T
) and vec is a column vectorization operator; that is,
(30)
vec
⁡
(
[
a
1
a
2
a
3
a
4
]
)
=
[
a
1
a
3
a
2
a
4
]
.
By assuming Σ
w
=
Σ
w
T, symmetry property of Σ
w, we can improve estimate as [12]
(31)
Σ
w
=
unvec
(
Σ
w
T
)
+
unvec
(
Σ
w
T
)
T
2
,
where unvec is the inverse operation of vec.Using (26) R
ε
(
1
) is evaluated from the observed signal as
(32)
R
ε
(
1
)
=
E
{
ε
(
n
)
ε
(
n
-
1
)
T
}
=
E
{
(
y
(
n
)
-
A
T
g
(
n
)
)
(
y
(
n
-
1
)
-
A
T
g
(
n
-
1
)
T
)
}
.
In (32), A is unknown. So, we will use a recursive algorithm for estimating Σ
w and A. This algorithm is discussed in Section 3.3.Now, we assume thatΣ
w and A are known and derive an estimate for Σ
e.Evaluating (9), for l
=
0, we obtain
(33)
Σ
e
=
R
ε
(
0
)
-
∑
j
=
0
p
A
j
Σ
w
A
j
T
,
where
(34)
R
ε
(
0
)
=
E
{
ε
(
n
)
ε
(
n
)
T
}
=
E
{
(
y
(
n
)
-
A
T
g
(
n
)
)
(
y
(
n
)
-
A
T
g
(
n
)
T
)
}
.


### 3.3. The Algorithm

The results derived in previous sections and subsections are summarized to propose a recursive algorithm for estimatingΣ
w and A
i’s  (
i
=
1,2
,
…
,
p
). The improved least-squares algorithm for multichannel processes based on inverse filtering (IFILSM) is as follows.Step 1.
Compute autocovariance estimates by means of given samples{
y
(
1
)
,
y
(
2
)
,
…
,
y
(
N
)
}:
(35)
R
^
y
(
k
)
=
1
N
∑
n
=
p
+
q
+
1
N
y
(
n
)
y
(
n
-
k
)
T
,
k
=
0
,
1
,
2
,
…
,
p
+
q
,
and use them to construct the covariance matrix estimates R
^
y, r
^
y, R
^
y
q, and r
^
y
q.Step 2.
Initializei
=
0, where i denotes the iteration number. Then, compute the estimate
(36)
A
^
(
i
)
=
A
^
H
=
(
H
^
y
T
H
^
y
)
-
1
H
^
y
T
h
^
y
.Step 3.
Seti
=
i
+
1 and calculate the estimate of observation noise covariance matrix Σ
w:
(37)
R
^
ε
(
1
)
=
1
N
∑
n
=
p
+
2
N
{
(
y
(
n
)
-
A
^
IFILSM
(
i
-
1
)

T
g
(
n
)
)
m
m
m
m
m
×
(
y
(
n
-
1
)
-
A
^
IFILSM
(
i
-
1
)

T
g
(
n
-
1
)
T
)
A
^
IFILSM
(
i
-
1
)

T
g
}
vec
⁡
(
Σ
^
w
T
)
=
(
∑
j
=
0
p
A
^
j
⊗
A
^
j
-
1
)
-
1
r
^
ε
Σ
^
w
(
i
)
=
unvec
(
Σ
^
w
T
)
+
unvec
(
Σ
^
w
T
)
T
2
R
^
w
(
i
)
=
Σ
^
w
⊗
I
p
.Step 4.
Perform the bias correction(38)
A
^
IFILSM
(
i
)
=
A
^
IFILSM
(
i
-
1
)
+
μ
(
A
^
H
+
R
^
H
1
R
^
w
(
i
)
A
^
IFILSM
(
i
-
1
)
-
A
^
IFILSM
(
i
-
1
)
)
.Step 5.
If(39)
∥
A
^
IFILSM
(
i
)
-
A
^
IFILSM
(
i
-
1
)
∥
∥
A
^
IFILSM
(
i
-
1
)
∥
≤
δ
,
where ∥
·
∥ is the Euclidean norm and δ is an appropriate small positive number, the convergence is achieved and the iteration process must be terminated; otherwise, go to Step 3.Step 6.
EstimateΣ
e via (33):
(40)
Σ
^
e
=
R
^
ε
(
0
)
-
∑
j
=
0
p
A
^
j
Σ
^
w
A
^
j
T
,
R
^
ε
(
0
)
=
1
N
∑
n
=
p
+
1
N
{
(
y
(
n
)
-
A
^
T
g
(
n
)
)
(
y
(
n
)
-
A
^
T
g
(
n
)
T
)
}
.


## 3.1. MAR Parameter Estimation Given Observation Noise Covariance Matrix

The Yule-Walker equations for the noisy MAR process given by (4) are
(11)
R
y
(
l
)
-
∑
i
=
1
p
A
i
R
y
(
l
-
i
)
+
∑
i
=
1
p
A
i
Σ
w
δ
l
-
i
=
0
,
l
≥
1
.
Evaluating (11) for l
=
1
,
…
,
p and using R
y
(
l
)
=
R
y
(
-
l
)
T give the following system of linear equations:
(12)
(
R
y
-
R
w
)
A
-
r
y
=
0
,
where
(13)
R
y
=
[
R
y
(
0
)
R
y
(
1
)
⋯
R
y
(
p
-
1
)
R
y
(
-
1
)
R
y
(
0
)
⋯
R
y
(
p
-
2
)
⋮
⋮
⋱
⋮
R
y
(
1
-
p
)
R
y
(
2
-
p
)
⋯
R
y
(
0
)
]
,
r
y
=
[
R
y
(
1
)
R
y
(
2
)
⋯
R
y
(
p
)
]
T
,
and R
w is a block-diagonal m
p
×
m
p matrix given by
(14)
R
w
=
[
Σ
w
0
⋯
0
0
Σ
w
⋮
⋮
⋮
0
⋱
0
0
⋯
⋯
Σ
w
]
=
Σ
w
⊗
I
p
.Evaluating (11) for l
=
p
+
1
,
…
,
p
+
q and using R
y
(
l
)
=
R
y
(
-
l
)
T give the following modified Yule-Walker (MYW) equations:
(15)
R
y
q
A
-
r
y
q
=
0
,
where
(16)
R
y
q
=
[
R
y
(
-
p
)
R
y
(
-
p
+
1
)
⋯
R
y
(
-
1
)
R
y
(
-
p
-
1
)
R
y
(
-
p
)
⋯
R
y
(
-
2
)
⋮
⋮
⋱
⋮
R
y
(
-
p
-
q
)
R
y
(
-
p
-
q
+
1
)
⋯
R
y
(
-
q
)
]
,
r
y
q
=
[
R
y
(
p
+
1
)
R
y
(
p
+
2
)
⋯
R
y
(
p
+
q
)
]
T
,
and q is an arbitrary nonnegative integer.Combining the low-order Yule-Walker equations given by (12) and the high-order Yule-Walker equations given by (15), we obtain
(17)
H
y
A
-
H
w
A
-
h
y
=
0
,
where
(18)
H
y
=
[
R
y
R
y
q
]
,
h
y
=
[
r
y
r
y
q
]
,
(19)
H
w
=
[
R
w
0
]
.Premultiplying (17) by (
H
y
T
H
y
)
-
1
H
y
T, we obtain
(20)
A
-
A
H
-
R
H
H
w
A
=
0
,
where
(21)
A
H
=
(
H
y
T
H
y
)
-
1
H
y
T
h
y
,
R
H
=
(
H
y
T
H
y
)
-
1
H
y
T
.
And R
H of dimension m
p
×
m
(
p
+
q
) can be partitioned as
(22)
R
H
=
[
R
H
1
R
H
2
]
,
where the dimensions of R
H
1 and R
H
2 are m
p
×
m
p and m
p
×
m
q, respectively.Substituting (19) and (22) into (20), we get
(23)
A
-
A
H
-
R
H
1
R
w
A
=
0
.
If R
w is assumed to be known, we can iteratively estimate A using steepest descent algorithm as follows:
(24)
A
(
i
)
=
A
(
i
-
1
)
+
μ
(
A
H
+
R
H
1
R
w
A
(
i
-
1
)
-
A
(
i
-
1
)
)
,
where μ is the step size parameter and A
H, R
H
1 can be estimated from the observations. We can also control the stability and rate of convergence of the algorithm by changing μ. The above equation is converged if μ is selected between zero and 2
/
(
1
-
λ
min
⁡
) where λ
min
⁡ is the maximum eigenvalue of the R
H
1
R
w matrix. Furthermore if we set μ
=
1 and q
=
0 in (24), it is changed to
(25)
A
(
i
)
=
A
l
s
+
R
y
-
1
R
w
A
(
i
-
1
)
,
where A
l
s
=
R
y
-
1
r
y, which is previously derived in [8].


## 3.2. Application of Inverse Filtering for Estimation of Observation Noise Covariance Matrix Given MAR Matrices

If the observed process is filtered by MAR matrices, the output process of the filter is obtained as(26)
ε
(
n
)
=
y
(
n
)
-
A
T
g
(
n
)
.
By evaluating (9) for l
=
1, we can obtain
(27)
R
ε
(
1
)
=
∑
j
=
1
p
A
j
Σ
w
A
j
-
1
T
.
We can rewrite (27) as
(28)
(
∑
j
=
1
p
A
j
⊗
A
j
-
1
)
vec
⁡
(
Σ
w
T
)
=
vec
⁡
(
R
ε
(
1
)
T
)
.
We solve (28) as follows:
(29)
vec
⁡
(
Σ
w
T
)
=
(
∑
j
=
1
p
A
j
⊗
A
j
-
1
)
-
1
r
ε
,
where r
ε
=
vec
⁡
(
R
ε
(
1
)
T
) and vec is a column vectorization operator; that is,
(30)
vec
⁡
(
[
a
1
a
2
a
3
a
4
]
)
=
[
a
1
a
3
a
2
a
4
]
.
By assuming Σ
w
=
Σ
w
T, symmetry property of Σ
w, we can improve estimate as [12]
(31)
Σ
w
=
unvec
(
Σ
w
T
)
+
unvec
(
Σ
w
T
)
T
2
,
where unvec is the inverse operation of vec.Using (26) R
ε
(
1
) is evaluated from the observed signal as
(32)
R
ε
(
1
)
=
E
{
ε
(
n
)
ε
(
n
-
1
)
T
}
=
E
{
(
y
(
n
)
-
A
T
g
(
n
)
)
(
y
(
n
-
1
)
-
A
T
g
(
n
-
1
)
T
)
}
.
In (32), A is unknown. So, we will use a recursive algorithm for estimating Σ
w and A. This algorithm is discussed in Section 3.3.Now, we assume thatΣ
w and A are known and derive an estimate for Σ
e.Evaluating (9), for l
=
0, we obtain
(33)
Σ
e
=
R
ε
(
0
)
-
∑
j
=
0
p
A
j
Σ
w
A
j
T
,
where
(34)
R
ε
(
0
)
=
E
{
ε
(
n
)
ε
(
n
)
T
}
=
E
{
(
y
(
n
)
-
A
T
g
(
n
)
)
(
y
(
n
)
-
A
T
g
(
n
)
T
)
}
.


## 3.3. The Algorithm

The results derived in previous sections and subsections are summarized to propose a recursive algorithm for estimatingΣ
w and A
i’s  (
i
=
1,2
,
…
,
p
). The improved least-squares algorithm for multichannel processes based on inverse filtering (IFILSM) is as follows.Step 1.
Compute autocovariance estimates by means of given samples{
y
(
1
)
,
y
(
2
)
,
…
,
y
(
N
)
}:
(35)
R
^
y
(
k
)
=
1
N
∑
n
=
p
+
q
+
1
N
y
(
n
)
y
(
n
-
k
)
T
,
k
=
0
,
1
,
2
,
…
,
p
+
q
,
and use them to construct the covariance matrix estimates R
^
y, r
^
y, R
^
y
q, and r
^
y
q.Step 2.
Initializei
=
0, where i denotes the iteration number. Then, compute the estimate
(36)
A
^
(
i
)
=
A
^
H
=
(
H
^
y
T
H
^
y
)
-
1
H
^
y
T
h
^
y
.Step 3.
Seti
=
i
+
1 and calculate the estimate of observation noise covariance matrix Σ
w:
(37)
R
^
ε
(
1
)
=
1
N
∑
n
=
p
+
2
N
{
(
y
(
n
)
-
A
^
IFILSM
(
i
-
1
)

T
g
(
n
)
)
m
m
m
m
m
×
(
y
(
n
-
1
)
-
A
^
IFILSM
(
i
-
1
)

T
g
(
n
-
1
)
T
)
A
^
IFILSM
(
i
-
1
)

T
g
}
vec
⁡
(
Σ
^
w
T
)
=
(
∑
j
=
0
p
A
^
j
⊗
A
^
j
-
1
)
-
1
r
^
ε
Σ
^
w
(
i
)
=
unvec
(
Σ
^
w
T
)
+
unvec
(
Σ
^
w
T
)
T
2
R
^
w
(
i
)
=
Σ
^
w
⊗
I
p
.Step 4.
Perform the bias correction(38)
A
^
IFILSM
(
i
)
=
A
^
IFILSM
(
i
-
1
)
+
μ
(
A
^
H
+
R
^
H
1
R
^
w
(
i
)
A
^
IFILSM
(
i
-
1
)
-
A
^
IFILSM
(
i
-
1
)
)
.Step 5.
If(39)
∥
A
^
IFILSM
(
i
)
-
A
^
IFILSM
(
i
-
1
)
∥
∥
A
^
IFILSM
(
i
-
1
)
∥
≤
δ
,
where ∥
·
∥ is the Euclidean norm and δ is an appropriate small positive number, the convergence is achieved and the iteration process must be terminated; otherwise, go to Step 3.Step 6.
EstimateΣ
e via (33):
(40)
Σ
^
e
=
R
^
ε
(
0
)
-
∑
j
=
0
p
A
^
j
Σ
^
w
A
^
j
T
,
R
^
ε
(
0
)
=
1
N
∑
n
=
p
+
1
N
{
(
y
(
n
)
-
A
^
T
g
(
n
)
)
(
y
(
n
)
-
A
^
T
g
(
n
)
T
)
}
.


## 4. Simulation Result

In simulation study, we have compared the IFILSM method with the other existing methods. Two examples have been considered: first frequency estimation and second synthetic noisy MAR model estimation. In all simulations, we setμ
=
1 for the proposed IFILSM algorithm.


### 4.1. Estimating the Frequency of Multiple Sinusoids Embedded in Spatially Colored Noise

Consider several sinusoidal signals in spatially colored noise as follows:(41)
z
(
n
)
=
[
z
1
(
n
)
=
∑
k
=
1
L
a
1
k
cos
⁡
(
ω
1
k
n
+
φ
1
k
)
⋮
z
m
(
n
)
=
∑
k
=
1
L
a
m
k
cos
⁡
(
ω
m
k
n
+
φ
m
k
)
]
,
y
(
n
)
=
z
(
n
)
+
w
(
n
)
.
It is well known that z
(
n
) can be modeled as an MAR (2
L) process with zero driving noise vector (e
(
n
)
=
0). So, we can apply the method proposed in the previous subsection to estimate the frequency of sinusoidal signals corrupted with spatially colored observation noise. Note that the frequency of sinusoids can be estimated by locating the peaks in the MAR spectrum.Example 1.
We considerN
=
500 data samples of a signal comprised of two sinusoids corrupted with spatially colored noise. Consider y
(
n
) as follows:
(42)
y
(
n
)
=
[
a
1
cos
⁡
⁡
(
ω
1
n
+
φ
1
)
a
2
cos
⁡
⁡
(
ω
2
n
+
φ
2
)
]
+
w
(
n
)
,
where ω
1
=
2
π
0.1, ω
2
=
2
π
0.3, φ
1
,
φ
2 are independent random phases uniformly distributed over [
0
,
2
π
], and w
(
n
) is Gaussian spatially colored noise having the following covariance matrix:
(43)
Σ
w
=
[
σ
w
1
2
1
1
σ
w
2
2
]
.
The signal-to-noise ratio (SNR) of each channel is defined as
(44)
SNR
i
=
a
i
2
2
σ
w
i
2
which is set to 5 and 10 dB. The number of trials M is set to 1000 and the parameter q in the IFILSM and MYW methods is set to 4. The δ in the ILSV, ALSV, and IFILSM methods is set to 0.01. We compare the proposed method with the LS, MYW, ALSV, and ILSV methods in terms of mean square error (MSE) of frequency estimates. Note that other existing methods cannot work in the spatially correlated observation noise. In Table 1, we tabulate the MSEs in terms of dB for the LS, MYW, ALSV, ILSV, and proposed methods. We see that, from the table, the proposed method has the best performance.Table 1
Simulation results (M
=
1000, N
=
500, SNR = 5, 10 dB, ω
1
=
2
π
0.1, and ω
2
=
2
π
0.3).
SNR
5 dB
5 dB
10 dB
10 dB
Methods
MSE1 (dB)
MSE2 (dB)
MSE1 (dB)
MSE2 (dB)
IFILSM
-
42.18
-
47.04
-
51.11
-
60.06
ILSV
-
32.64
-
41.91
-
42.22
-
50.51
ALSV
-
32.35
-
42.76
-
42.33
-
51.36
LS
-
5.75
-
23.02
-
32
-
31.27
MYW
-
37.58
-
45.05
-
46.82
-
51.94


### 4.2. Synthetic MAR Processes

Example 2.
Consider a noisy MAR model withp
=
2, m
=
2, and coefficient matrices
(45)
A
1
=
[
1.1625
0.1432
0.1432
1.5443
]
,
A
2
=
[
-
0.96
-
0.0064
-
0.064
-
0.64
]
.
The driving process e
(
n
) and observation noise w
(
n
) covariance are mutually uncorrelated Gaussian temporally white and spatially correlated processes with covariance matrices
(46)
Σ
e
=
[
1
0.2
0.2
0.7
]
,
Case
1
:
Σ
w
≈
[
0.4149
0.01
0.01
0.4689
]
,
Case
2
:
Σ
w
≈
[
1.29
0.01
0.01
1.44
]
.
We compare the ILSV, the ALSV, and the proposed IFILSM method in terms MSE of MAR coefficient estimate and the convergence probability defined as
(47)
P
=
M
M
t
,
where M is the number of runs that an iterative algorithm is converged and M
t is the total number of runs. We set the parameters of the algorithms M
=
1000 and the parameter q in the IFILSM method is set to 4. The δ in the ILSV, ALSV, and IFILSM methods is set to 0.01.
In Case 1 and Case 2 SNR for each channel is set to 15 and 10 dB, respectively. In this example data samplesN are changed from 1000 to 4000 with step size 1000. The MSE and P values are plotted in Figures 1, 2, and 3 versus N, respectively. The probability of convergence is one for all algorithms in Case 1 and the figure is not plotted here. From the figures, we can see that the accuracy and the convergence of the IFILSM are better than the ILSV and ALSV methods. Note that the ILSV and ALSV methods have poor convergence when SNR is lower or equal to 10 dB. It can be seen from the figures that the performance of the all algorithms is constant over the sample size because autocovariance estimates of the observations are not changed in this rage of data sample size.Figure 1
MSE values of MAR matrices estimate versus data samples;M
=
1000, SNR = 15 dB.Figure 2
MSE values of MAR matrices estimate versus data samples;M
=
1000, SNR = 10 dB.Figure 3
Probability of convergence values of MAR matrices estimate versus data samples;M
=
1000, SNR = 10 dB.


## 4.1. Estimating the Frequency of Multiple Sinusoids Embedded in Spatially Colored Noise

Consider several sinusoidal signals in spatially colored noise as follows:(41)
z
(
n
)
=
[
z
1
(
n
)
=
∑
k
=
1
L
a
1
k
cos
⁡
(
ω
1
k
n
+
φ
1
k
)
⋮
z
m
(
n
)
=
∑
k
=
1
L
a
m
k
cos
⁡
(
ω
m
k
n
+
φ
m
k
)
]
,
y
(
n
)
=
z
(
n
)
+
w
(
n
)
.
It is well known that z
(
n
) can be modeled as an MAR (2
L) process with zero driving noise vector (e
(
n
)
=
0). So, we can apply the method proposed in the previous subsection to estimate the frequency of sinusoidal signals corrupted with spatially colored observation noise. Note that the frequency of sinusoids can be estimated by locating the peaks in the MAR spectrum.Example 1.
We considerN
=
500 data samples of a signal comprised of two sinusoids corrupted with spatially colored noise. Consider y
(
n
) as follows:
(42)
y
(
n
)
=
[
a
1
cos
⁡
⁡
(
ω
1
n
+
φ
1
)
a
2
cos
⁡
⁡
(
ω
2
n
+
φ
2
)
]
+
w
(
n
)
,
where ω
1
=
2
π
0.1, ω
2
=
2
π
0.3, φ
1
,
φ
2 are independent random phases uniformly distributed over [
0
,
2
π
], and w
(
n
) is Gaussian spatially colored noise having the following covariance matrix:
(43)
Σ
w
=
[
σ
w
1
2
1
1
σ
w
2
2
]
.
The signal-to-noise ratio (SNR) of each channel is defined as
(44)
SNR
i
=
a
i
2
2
σ
w
i
2
which is set to 5 and 10 dB. The number of trials M is set to 1000 and the parameter q in the IFILSM and MYW methods is set to 4. The δ in the ILSV, ALSV, and IFILSM methods is set to 0.01. We compare the proposed method with the LS, MYW, ALSV, and ILSV methods in terms of mean square error (MSE) of frequency estimates. Note that other existing methods cannot work in the spatially correlated observation noise. In Table 1, we tabulate the MSEs in terms of dB for the LS, MYW, ALSV, ILSV, and proposed methods. We see that, from the table, the proposed method has the best performance.Table 1
Simulation results (M
=
1000, N
=
500, SNR = 5, 10 dB, ω
1
=
2
π
0.1, and ω
2
=
2
π
0.3).
SNR
5 dB
5 dB
10 dB
10 dB
Methods
MSE1 (dB)
MSE2 (dB)
MSE1 (dB)
MSE2 (dB)
IFILSM
-
42.18
-
47.04
-
51.11
-
60.06
ILSV
-
32.64
-
41.91
-
42.22
-
50.51
ALSV
-
32.35
-
42.76
-
42.33
-
51.36
LS
-
5.75
-
23.02
-
32
-
31.27
MYW
-
37.58
-
45.05
-
46.82
-
51.94


## 4.2. Synthetic MAR Processes

Example 2.
Consider a noisy MAR model withp
=
2, m
=
2, and coefficient matrices
(45)
A
1
=
[
1.1625
0.1432
0.1432
1.5443
]
,
A
2
=
[
-
0.96
-
0.0064
-
0.064
-
0.64
]
.
The driving process e
(
n
) and observation noise w
(
n
) covariance are mutually uncorrelated Gaussian temporally white and spatially correlated processes with covariance matrices
(46)
Σ
e
=
[
1
0.2
0.2
0.7
]
,
Case
1
:
Σ
w
≈
[
0.4149
0.01
0.01
0.4689
]
,
Case
2
:
Σ
w
≈
[
1.29
0.01
0.01
1.44
]
.
We compare the ILSV, the ALSV, and the proposed IFILSM method in terms MSE of MAR coefficient estimate and the convergence probability defined as
(47)
P
=
M
M
t
,
where M is the number of runs that an iterative algorithm is converged and M
t is the total number of runs. We set the parameters of the algorithms M
=
1000 and the parameter q in the IFILSM method is set to 4. The δ in the ILSV, ALSV, and IFILSM methods is set to 0.01.
In Case 1 and Case 2 SNR for each channel is set to 15 and 10 dB, respectively. In this example data samplesN are changed from 1000 to 4000 with step size 1000. The MSE and P values are plotted in Figures 1, 2, and 3 versus N, respectively. The probability of convergence is one for all algorithms in Case 1 and the figure is not plotted here. From the figures, we can see that the accuracy and the convergence of the IFILSM are better than the ILSV and ALSV methods. Note that the ILSV and ALSV methods have poor convergence when SNR is lower or equal to 10 dB. It can be seen from the figures that the performance of the all algorithms is constant over the sample size because autocovariance estimates of the observations are not changed in this rage of data sample size.Figure 1
MSE values of MAR matrices estimate versus data samples;M
=
1000, SNR = 15 dB.Figure 2
MSE values of MAR matrices estimate versus data samples;M
=
1000, SNR = 10 dB.Figure 3
Probability of convergence values of MAR matrices estimate versus data samples;M
=
1000, SNR = 10 dB.


## 5. Conclusion

Steepest descent (SD) algorithm is used for estimating noisy multichannel AR model parameters. The observation noise covariance matrix is nondiagonal; this means that the channel noises are assumed to be spatially correlated. The inverse filtering idea is used to estimate the observation noise covariance matrix. Computer simulations showed that the proposed method has better performance and convergence than the ILSV and the ALSV methods.



---

*Source: 502406-2014-06-19.xml*