Adaptive Filtering

Recall optimal ltering: Given

x(n) = d(n) +v(n),
estimate and extract d(n) from the current and past values of x(n).
EE 524, # 11 1
Let the lter coecients be
w =
_

_
w
0
w
1
.
.
.
w
N1
_

_
.
Filter output:
y(n) =
N1

k=0
w

k
x(n k) = w
H
x(n) =

d(n),
where
x(n) =
_

_
x(n)
x(n 1)
.
.
.
x(n N + 1)
_

_
.
EE 524, # 11 2
Wiener-Hopf equation:
R(n)w(n) = r(n) w
opt
(n) = R(n)
1
r(n),
where
R(n) = E{x(n)x(n)
H
},
r(n) = E{x(n)d(n)

}.
EE 524, # 11 3
Adaptive Filtering (cont.)
Example 1: Unknown system identication.
EE 524, # 11 4
Adaptive Filtering (cont.)
Example 2: Unknown system equalization.
EE 524, # 11 5
Adaptive Filtering (cont.)
Example 3: Noise cancellation.
EE 524, # 11 6
Adaptive Filtering (cont.)
Example 4: Signal linear prediction.
EE 524, # 11 7
Adaptive Filtering (cont.)
Example 5: Interference cancellation without reference input.
EE 524, # 11 8
Adaptive Filtering (cont.)
Idea of the Least-Mean-Square (LMS) algorithm:
w
k+1
= w
k
(
w
E{|e
k
|
2
})

, ()
where the indices are given as subscripts [e.g. d(k) = d
k
], and
E{|e
k
|
2
} = E{|d
k
w
H
k
x
k
|
2
}
= E{|d
k
|
2
} w
H
k
r r
H
w
k
+w
H
k
Rw
k
,
(
w
E{|e
k
|
2
})

= Rw r.
Use single-sample estimates of R and r:

R = x
k
x
H
k
, r = x
k
d

k
,
EE 524, # 11 9
and insert them into ():
w
k+1
= w
k
+x
k
e

k
, e
k
= d
k
w
H
k
x
k
LMS alg.
EE 524, # 11 10
Adaptive Filtering: Convergence Analysis
Convergence analysis: Subtract w
opt
from both sides of the previous
equation:
w
k+1
w
opt
. .
v
k+1
= w
k
w
opt
. .
v
k
+x
k
(d

k
x
H
k
w
k
) ()
and note that
x
k
(d

k
x
H
k
w
k
) = x
k
d

k
x
k
x
H
k
w
k
= x
k
d

k
x
k
x
H
k
w
k
+x
k
x
H
k
w
opt
x
k
x
H
k
w
opt
= (x
k
d

k
x
k
x
H
k
w
opt
) x
k
x
H
k
v
k
.
EE 524, # 11 11
Observe that
E
_
x
k
(d

k
x
H
k
w
k
)
_
= r Rw
opt
. .
0
RE{v
k
} = RE{v
k
}.
Let c
k
= E{v
k
}. Then
c
k+1
= [I R]c
k
( )
Sucient condition for convergence:
c
k+1
< c
k
k.
EE 524, # 11 12
Adaptive Filtering: Convergence Analysis
Let us premultiply both parts of the equation ( ) by the matrix U
H
of
the eigenvectors of R, where
R = UU
H
.
Then, we have
U
H
c
k+1
. .

c
k+1
= U
H
[I R] UU
H
. .
I
c
k
,
and, hence
c
k+1
= [I ]c
k
.
Since
c
k

2
= c
H
k
c
k
= c
H
k
UU
H
. .
I
c
k
= c
H
k
c
k
= c
k

2
,
EE 524, # 11 13
the sucient condition for convergence can be rewritten as
c
k+1

2
< c
k

2
k.
Let us then require that the absolute value of each component of the vector
c
k+1
is less than that of c
k
:
|1
i
| < 1, i = 1, 2, . . . , N.
The condition
|1
i
| < 1, i = 1, 2, . . . , N,
is equivalent to
0 < <
2

max
where
max
is the maximum eigenvalue of R. In practice, even a stronger
EE 524, # 11 14
condition is (often) used:
0 < <
2
tr{R}
,
where tr{R} >
max
.
EE 524, # 11 15
Normalized LMS
A promising variant of LMS is the so-called Normalized LMS (NLMS)
algorithm:
w
k+1
=w
k
+

x
k

2
x
k
e

k
, e
k
= d
k
w
H
k
x
k
NLMS alg.
The sucient condition for convergence:
0 < < 2.
In practice, at some time points x
k
can be very small. To make the
NLMS algorithm more robust, we can modify it as follows:
w
k+1
= w
k
+

x
k

2
+
x
k
e

k
,
so that the gain constant cannot go to innity.
EE 524, # 11 16
Recursive Least Squares
Idea of the Recursive Least Squares (RLS) algorithm: use sample estimate

R
k
(instead of true covariance matrix R) in the equation for the weight
vector and nd w
k+1
as an update to w
k
. Let

R
k+1
=

R
k
+x
k+1
x
H
k+1
r
k+1
= r
k
+x
k+1
d

k+1
,
where 1 is the (so-called) forgetting factor. Using the matrix inversion
lemma, we obtain

R
1
k+1
= (

R
k
+x
k+1
x
H
k+1
)
1
=
1

R
1
k

R
1
k
x
k+1
x
H
k+1

R
1
k
+x
H
k+1

R
1
k
x
k+1
_
.
EE 524, # 11 17
Therefore,
w
k+1
=

R
1
k+1
r
k+1
=
_

R
1
k

R
1
k
x
k+1
x
H
k+1

R
1
k
+x
H
k+1

R
1
k
x
k+1
_
r
k
+
1

R
1
k

R
1
k
x
k+1
x
H
k+1

R
1
k
+x
H
k+1

R
1
k
x
k+1
_
x
k+1
d

k+1
= w
k
g
k+1
x
H
k+1
w
k
+g
k+1
d

k+1
,
where
g
k+1
=

R
1
k
x
k+1
+x
H
k+1

R
1
k
x
k+1
.
EE 524, # 11 18
Hence, the updating equation for the weight vector is
w
k+1
= w
k
g
k+1
x
H
k+1
w
k
+g
k+1
d

k+1
= w
k
+g
k+1
(d

k+1
x
H
k+1
w
k
)
. .
e

k,k+1
= w
k
+g
k+1
e

k,k+1
.
EE 524, # 11 19
RLS algorithm:
Initialization: w
0
= 0, P
0
=
1
I
For each k = 1, 2, . . ., compute:
h
k
= P
k1
x
k
,

k
= 1/( +h
H
k
x
k
),
g
k
= h
k

k
,
P
k
=
1
_
P
k1
g
k
h
H
k

,
e
k1,k
= d
k
w
H
k1
x
k
,
w
k
= w
k1
+g
k
e

k1,k
,
e
k
= d
k
w
H
k
x
k
.
EE 524, # 11 20
Example
LMS linear predictor of the signal
x(n) = 10e
j2fn
+e(n)
where f = 0.1 and
N = 8,
e(n) is circular unit-variance white noise,

1
= 1/[10 tr(R)],
2
= 1/[3 tr(R)],
3
= 1/[tr(R)].
EE 524, # 11 21
EE 524, # 11 22
Adaptive Beamforming
The above scheme describes narrowband beamforming, i.e.
conventional beamforming if w
1
, . . . , w
N
do not depend on the
EE 524, # 11 23
input/output array signals,
adaptive beamforming if w
1
, . . . , w
N
are determined and optimized based
on the input/output array signals.
Input array signal vector:
x(i) =
_

_
x
1
(i)
x
2
(i)
.
.
.
x
N
(i)
_

_
.
Complex beamformer output:
y(i) = w
H
x(i).
EE 524, # 11 24
Adaptive Beamforming (cont.)
Input array signal vector:
x(k) =
_

_
x
1
(k)
x
2
(k)
.
.
.
x
N
(k)
_

_
.
Complex beamformer output:
y(k) = w
H
x(k),
x(k) = x
s
(k)
. .
signal
+x
N
(k)
. .
noise
+ x
I
(k)
. .
interference
.
The goal is to lter out x
I
and x
N
as much as possible and, therefore,
EE 524, # 11 25
to obtain an approximation x
S
of x
S
. Most popular criteria of adaptive
beamforming:
MSE minimum
min
w
MSE, MSE = E{|d(i) w
H
x(i)|
2
}.
Signal-to-Interference-plus-Noise-Ratio (SINR)
max
w
SINR, SINR =
E{|w
H
x
s
|
2
}
E{|w
H
(x
I
+x
N
)|
2
}
.
EE 524, # 11 26
Adaptive Beamforming (cont.)
EE 524, # 11 27
Adaptive Beamforming (cont.)
In the sequel, we consider the max SINR criterion. Rewrite the snapshot
model as
x(k) = s(k)a
s
+x
I
(k) +x
N
(k),
where a
S
is the known steering vector of the desired signal. Then
SINR =

2
s
|w
H
a
s
|
2
w
H
E{(x
I
+x
N
)(x
I
+x
N
)
H
}w
=

2
s
|w
H
a
s
|
2
w
H
Rw
where
R = E{(x
I
+x
N
)(x
I
+x
N
)
H
}
is the interference-plus-noise covariance matrix.
Obviously, SINR does not depend on rescaling of w, i.e. if w
opt
is an
optimal weight, then w
opt
is such a vector too. Therefore, max SINR is
EE 524, # 11 28
equivalent to
min
w
w
H
Rw subject to w
H
a
S
= const.
Let const = 1. Then
H(w) = w
H
Rw +(1 w
H
a
s
) +

(1 a
H
s
w)

w
H(w) = (Rw a
s
)

= 0 =
Rw = a
s
=w
opt
= R
1
a
s
.
This is a spatial version of the Wiener-Hopf equation!
From the constraint equation, we obtain
=
1
a
H
s
R
1
a
s
EE 524, # 11 29
and therefore
w
opt
=
1
a
H
s
R
1
a
s
R
1
a
s
MVDR beamformer.
Substituting w
opt
into the SINR expression, we obtain
max SINR = SINR
opt
=

2
s
(a
H
s
R
1
a
s
)
2
a
H
s
R
1
RR
1
a
s
=
2
s
a
H
s
R
1
a
s
.
If there are no interference sources (only white noise with variance
2
):
SINR
opt
=

2
s

2
a
H
s
a
s
=
N
2
s

2
.
EE 524, # 11 30
Adaptive Beamforming (cont.)
Let us study what happens with the optimal SINR if the covariance matrix
includes the signal component:
R
x
= E{xx
H
} = R +
2
s
a
s
a
H
s
.
Using the matrix inversion lemma, we have
R
1
x
a
s
= (R +
2
s
a
s
a
H
S
)
1
a
s
=
_
R
1

R
1
a
s
a
H
s
R
1
1/
2
s
+a
H
s
R
1
a
s
_
a
s
=
_
1
a
H
s
R
1
a
s
1/
2
s
+a
H
s
R
1
a
s
_
R
1
a
s
= R
1
a
s
.
EE 524, # 11 31
Optimal SINR is not aected!
However, the above result holds only if
there is an innite number of snapshots and
a
S
is known exactly.
EE 524, # 11 32
Adaptive Beamforming (cont.)
Gradient algorithm maximizing SNR (very similar to LMS):
w
k+1
= w
k
+(a
s
x
k
x
H
k
w
k
),
where, again, we use the simple notation w
k
= w(k) and x
k
= x(k). The
vector w
k
converges to w
opt
R
1
a
s
if
0 < <
2

max
= 0 < <
2
tr{R}
.
The disadvantage of the gradient algorithms is that the convergence may
be very slow, i.e. it depends on the eigenvalue spread of R.
EE 524, # 11 33
Example
N = 8,
single signal from
s
= 0

and SNR = 0 dB,

single interference from
I
= 30

and INR= 40 dB,

1
= 1/[50 tr(R)],
2
= 1/[15 tr(R)],
3
= 1/[5 tr(R)].
EE 524, # 11 34
EE 524, # 11 35
Adaptive Beamforming (cont.)
Sample Matrix Inversion (SMI) Algorithm:
w
SMI
=

R
1
a
S
,
where

R is the sample covariance matrix

R =
1
K
K

k=1
x
k
x
H
k
.
Reed-Mallet-Brennan (RMB) rule: under mild conditions, the mean losses
(relative to the optimal SINR) due to the SMI approximation of w
opt
do
not exceed 3 dB if
K 2N.
Hence, the SMI provides very fast convergence rate, in general.
EE 524, # 11 36
Adaptive Beamforming (cont.)
Loaded SMI:
w
LSMI
=

R
1
DL
a
S
,

R
DL
=

R +I,
where the optimal weight 2
2
. LSMI allows convergence faster than
N snapshots!
LSMI convergence rule: under mild conditions, the mean losses (relative to
the optimal SINR) due to the LSMI approximation of w
opt
do not exceed
few dBs if
K L
where L is the number of interfering sources. Hence, the LSMI provides
faster convergence rate than SMI (usually, 2N L)!
EE 524, # 11 37
Example
N = 10,
single signal from
s
= 0

and SNR = 0 dB,

single interference from
I
= 30

and INR= 40 dB,

SMI vs. LSMI.
EE 524, # 11 38
EE 524, # 11 39
EE 524, # 11 40
EE 524, # 11 41
Adaptive Beamforming (cont.)
Hung-Turner (Projection) Algorithm:
w
HT
= (I X(X
H
X)
1
X
H
)a
S
,
i.e. data-orthogonal projection is used instead of inverse covariance matrix.
For Hung-Turner method, a satisfactory performance is achieved with
K L.
Optimal value of K
K
opt
=
_
(N + 1)L 1.
Drawback: number of sources should be known a priori.
EE 524, # 11 42
EE 524, # 11 43
This eect is sometimes referred to as the signal cancellation phenomenon.
Additional constraints are required to stabilize the mean beam response
min
w
w
H
Rw subject to C
H
w = f.
1. Point constraints: Matrix of constrained directions:
C = [a
S,1
, a
S,2
a
S,M
],
where a
S,i
are all taken in the neighborhood of a
S
and include a
S
as well.
Vector of constraints:
f =
_

_
1
1
.
.
.
1
_

_
.
EE 524, # 11 44
2. Derivative constraints: Matrix of constrained directions:
C =
_
a
S
,
a()

=
S
, ,

M1
a()

M1

=
S
_
,
where a
S,i
are all taken in the neighborhood of a
S
and include a
S
as well.
Vector of constraints:
f =
_

_
1
0
.
.
.
0
_

_
.
Note that

k
a()

=
S
= D
k
a
S
,
where D is the matrix depending on
s
and on array geometry.
EE 524, # 11 45
Adaptive Beamforming (cont.)
w
opt
= R
1
C(C
H
R
1
C)
1
f
and its SMI version:
w
opt
=

R
1
C(C
H

R
1
C)
1
f.
Additional constraints protect the directions in the neighborhood of
the assumed signal direction.
Additional constraints require enough degrees of freedom (DOFs)
number of sensors must be large enough.
Gradient algorithms exist for the constraint adaptation.
EE 524, # 11 46
EE 524, # 11 47
Adaptive Beamforming (cont.)
Generalized Sidelobe Canceller (GSC): Let us decompose
w
opt
= R
1
C(C
H
R
1
C)
1
f
into two components, one in the constrained subspace, and one orthogonal
to it:
w
opt
= (P
C
+P

C
)
. .
I
w
opt
= C(C
H
C)
1
C
H
R
1
C(C
H
R
1
C)
1
. .
I
f
+P

C
R
1
C(C
H
R
1
C)
1
f.
EE 524, # 11 48
Generalizing this approach, we obtain the following decomposition for w
opt
:
w
opt
= w
q
Bw
a
,
where
w
q
= C(C
H
C)
1
f
is the so-called quiescent weight vector,
B
H
C = 0,
B is the blocking matrix, and w
a
is the new adaptive weight vector.
EE 524, # 11 49
Generalized Sidelobe Canceller (GSC):
Choice of B is not unique. We can take B = P

C
. However, in this case
B is not of full rank. More common choice is to assume N (N M)
full-rank matrix B. Then, the vectors z = B
H
x and w
a
both have
shorter length (N M) 1 relative to the N 1 vectors x and w
q
.
Since the constrained directions are blocked by the matrix B, the signal
cannot be suppressed and, therefore, the weight vector w
a
can adapt
EE 524, # 11 50
freely to suppress interference by minimizing the output GSC power
Q
GSC
= (w
q
Bw
a
)
H
R(w
q
Bw
a
)
= w
H
q
Rw
q
w
H
q
RBw
a
w
H
a
B
H
Rw
q
+w
H
a
B
H
RBw
a
.
The solution is w
a,opt
= (B
H
RB)
1
B
H
Rw
q
.
EE 524, # 11 51
Adaptive Beamforming (cont.)
Generalized Sidelobe Canceller (GSC): Noting that
y(k) = w
H
q
x(k), z(k) = B
H
x(k),
we obtain
R
z
= E{z(k)z(k)
H
}
= B
H
E{x(k)x(k)
H
}B
= B
H
RB,
r
yz
= E{z(k)y

(k)}
= B
H
E{x(k)x(k)
H
}w
q
= B
H
Rw
q
.
EE 524, # 11 52
Hence,
w
a,opt
= R
1
z
r
yz
Wiener-Hopf equation!
EE 524, # 11 53
How to Choose B?
Choose N M linearly independent vectors b
i
:
B = [b
1
b
2
b
NM
]
so that
b
i
c
k
, i = 1, 2, . . . , N M, k = 1, 2, . . . , M,
where c
k
is the kth column of C.
There are many possible choices of B!
EE 524, # 11 54
Example: GSC in the Particular Case of Normal
Direction (Single) Constraint and for a Particular Choice
of Blocking Matrix:
EE 524, # 11 55
In this particular example
C =
_

_
1
1
.
.
.
1
_

_
.
B
H
=
_

_
1 1 0 0 0
0 1 1 0 0
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
0 0 0 1 1
_

_
,
and
x(k) =
_

_
x
1
(k)
x
2
(k)
.
.
.
x
N
(k)
_

_
, z(k) =
_

_
x
1
(k) x
2
(k)
x
2
(k) x
3
(k)
.
.
.
x
N1
(k) x
N
(k)
_

_
.
EE 524, # 11 56
Partially Adaptive Beamforming
In many applications, number of interfering sources is much less than the
number of adaptive weights [adaptive degrees of freedom (DOFs)]. In such
cases, partially adaptive arrays can be used.
Idea: use nonadaptive preprocessor reducing the number of adaptive
channels:
y(i) = T
H
x(i),
where
y has a reduced dimension M1 (M < N) compared with N1 vector
x,
T is an N M full-rank matrix.
EE 524, # 11 57
EE 524, # 11 58
EE 524, # 11 59
Partially Adaptive Beamforming
There are two types of nonadaptive preprocessors:
subarray preprocessor,
beamspace preprocessor.
For arbitrary preprocessor:
R
y
= E{y(i)y(i)
H
} = T
H
E{x(i)x(i)
H
}T = T
H
RT.
Recall the previously-used representation:
R = ASA
H
+
2
I.
EE 524, # 11 60
After the preprocessing, we have
R
y
= T
H
ASA
H
T +
2
T
H
T
=

AS

A
H
+Q

A = T
H
A
Q =
2
T
H
T.
Preprocessing changes array manifold.
Preprocessing may lead to colored noise.
Choosing T with orthonormal columns, we have
T
H
T = I,
and, therefore, the eect of colored noise may be removed.
EE 524, # 11 61
Partially Adaptive Beamforming
EE 524, # 11 62
Preprocessing matrix in this particular case:
T
H
=
1

3
_
_
1 1 1 0 0 0 0 0 0
0 0 0 1 1 1 0 0 0
0 0 0 0 0 0 1 1 1
_
_
,
(note that T
H
T = I here!)
In the general case
T =
_

_
a
S,1
0 0
0 a
S,2
0
.
.
.
.
.
.
.
.
.
.
.
.
0 0 a
S,M
_

_
,
where L = N/M is the size of each subarray, and T
H
T = I holds true if
a
H
S,k
a
S,k
= 1, k = 1, 2, . . . , M.
EE 524, # 11 63
Wideband Space-Time Processing
In the wideband case, we must consider joint space-time processing:
EE 524, # 11 64
Wideband Space-Time Processing (cont.)
Wideband case:
Higher dimension of the problem (NP instead of N),
Steering vector depends on frequency.
EE 524, # 11 65
Constant Modulus Algorithm (CMA)
Application: separation of constant-modulus sources.
Narrowband signals: the received signal is an instantaneous linear
mixture:
x
k
= As
k
.
Objective: nd inverse W, so that
y
k
= W
H
x
k
= s
k
.
Challenge: both A and s
k
are unknown!
However, we have side knowledge: sources are phase modulated, i.e.
s
i
(t) = exp(j
i
(t)).
EE 524, # 11 66
EE 524, # 11 67
Constant Modulus Algorithm (cont.)
Simple example: 2 sources, 2 antennas.
EE 524, # 11 68
Let
w =
_
w
1
w
2
_
be a beamformer. Output of beamforming:
y
k
= w
H
x
k
= [w

1
w

2
]
_
x
1,k
x
2,k
_
.
Constant modulus property: |s
1,k
| = |s
2,k
| = 1 for all k.
Possible optimization problem:
minJ(w) where J(w) = E[(|y
k
|
2
1)
2
].
EE 524, # 11 69
EE 524, # 11 70
The CMA cost function as a function of y (for simplicity, y is taken to be
real here).
No unique minimum! Indeed, if y
k
= w
H
x
k
is CM, then another
beamformer is w, for any scalar that satises || = 1.
EE 524, # 11 71
2 (real-valued) sources, 2 antennas
EE 524, # 11 72
Iterative Optimization
Cost function:
J(w) = E[(|y
k
|
2
1)
2
], y
k
= w
H
x
k
.
Stochastic gradient method: w
k+1
= w
k
[J(w
k
)]

, where is step
size, > 0.
Derivative: Use |y
k
|
2
= y
k
y

k
= w
H
xx
H
w.
J = 2E{(|y
k
|
2
1) (w
H
x
k
x
H
k
w)}
= 2E{(|y
k
|
2
1) (x
k
x
H
k
w)

}
= 2E{(|y
k
|
2
1)x

k
y
k
}
EE 524, # 11 73
Algorithm CMA(2,2):
y
k
= w
H
k
x
k
w
k+1
= w
k
x
k
(|y
k
|
2
1)y

k
.
EE 524, # 11 74
Advantages:
The algorithm is extremely simple to implement
Adaptive tracking of sources
Converges to minima close to the Wiener beamformers (for each source)
Disadvantages:
Noisy and slow
Step size should be small, else unstable
Only one source is recovered (which one?)
Possible convergence to local minimum (with nite data)
EE 524, # 11 75
EE 524, # 11 76
EE 524, # 11 77
Other CMAs
Alternative cost function: CMA(1,2)
J(w) = E[(|y
k
| 1)
2
] = E[(|w
H
x
k
| 1)
2
].
EE 524, # 11 78
Corresponding CMA iteration:
y
k
= w
H
k
x
k

k
=
y
k
|y
k
|
y
k
w
k+1
= w
k
+x
k

k
.
Similar to LMS, with update error
y
k
|y
k
|
y
k
. The desired signal is estimated
by
y
k
|y
k
|
.
EE 524, # 11 79
Other CMAs (cont.)
Normalized CMA (NCMA; becomes scaling independent)
w
k+1
= w
k
+

x
k

2
x
k

k
.
Orthogonal CMA (OCMA): whiten using data covariance R
w
k+1
= w
k
+R
1
k
x
k

k
.
EE 524, # 11 80
Least squares CMA (LSCMA): block update, trying to optimize
iteratively
min
w
s
H
w
H
X
2
where X = [x
1
x
2
x
T
] and s
H
is the best blind estimate at step k of
the complete source vector (at all time points t = 1, 2, . . . , T)
s
H
=
_
y
1
|y
1
|
,
y
2
|y
2
|
, . . . ,
y
T
|y
T
|
_
,
where
y
t
= w
H
k
x
t
, t = 1, 2, . . . , K.
and
w
H
k
= s
H
X
H
(XX
H
)
1
.
EE 524, # 11 81