We assume a signal model as follows:
xn=∑
k
=1K
A
k
ei
ω
k
n+zn
x
n
k
K
1
A
k
ω
k
n
z
n
(1)
where
A
k
∈C
A
k
is a complex number representing the magnitude and phase of
the
kk-th frequency component. We
will examine the
subspace structure of signals
composed of several frequency components, starting with
examining its autocorrelation matrix.
Recall the autocorrelation of a signal
xn
x
n
is defined as:
r
x
k=Exnxn−k¯
r
x
k
E
x
n
x
n
k
(2)
and the autocorrelation matrix of
xn
x
n
is defined as:
R
x
=ExxH=(
r
x
0…
r
x
M−1
⋮⋱⋮
r
x
M−1…
r
x
0
)
R
x
E
x
x
r
x
0
…
r
x
M
1
⋮
⋱
⋮
r
x
M
1
…
r
x
0
(3)
Now we can eigen-decompose the correlation matrix
R
x
R
x
:
R
x
=UH
Λ
x
U
R
x
U
Λ
x
U
(4)
But what does this give us? Let us examine a simple case
first: a single frequency component in noise.
In this case we have:
xn=
A
1
ei
ω
1
k+zn
x
n
A
1
ω
1
k
z
n
(5)
where we assume as usual that
zn
z
n
is white noise. It can be shown that the
autocorrelation of
Equation 2
will be:
r
x
k=|
A
1
|2ei
ω
1
k+
σ
ω
2δk
r
x
k
A
1
2
ω
1
k
σ
ω
2
δ
k
(6)
Where the first term in the sum in the above equation,
containing the exponential, is the signal component and the
second term is the noise component. This equation gives a
decomposition of the autocorrelation matrix of
Equation 3:
R
x
=
R
s
+
R
n
R
x
R
s
R
n
(7)
where
R
s
R
s
and
R
n
R
n
are the signal and noise contributions
respectively. Writing out
Equation 6 in matrix form, we see that
R
x
=(
1e(−i)
ω
1
e(−i)
ω
1
2…e(−i)
ω
1
(M−1)
e(−i)
ω
1
⋮⋱⋱⋮
⋮⋱⋱⋮e(−i)
ω
1
e(−i)
ω
1
(M−1)…e(−i)
ω
1
2e(−i)
ω
1
1
)
R
x
1
ω
1
ω
1
2
…
ω
1
M
1
ω
1
⋮
⋱
⋱
⋮
⋮
⋱
⋱
⋮
ω
1
ω
1
M
1
…
ω
1
2
ω
1
1
and
R
n
=
σ
ω
2I
R
n
σ
ω
2
I
Now define
e
1
=(
1ei
ω
1
ei
ω
1
2…e(−i)
ω
1
(M−1)
)
e
1
1
ω
1
ω
1
2
…
ω
1
M
1
. Then:
R
s
=|
A
1
|2
e
1
e
1
H
R
s
A
1
2
e
1
e
1
By inspection,
rank
R
s
=1
rank
R
s
1
, so there is only one nonzero eigenvalue
λ
1
≠0
λ
1
0
. Moreover,
R
s
e
1
=|
A
1
|2
e
1
e
1
H
e
1
=M|
A
1
|2
e
1
R
s
e
1
A
1
2
e
1
e
1
e
1
M
A
1
2
e
1
(8)
In other words,
e
1
e
1
is an eigenvector of
R
s
R
s
with eigenvalue
λ
1
=M|
A
1
|2
λ
1
M
A
1
2
.
Here we extend the treatment of the previous section to include
multiple frequency components. First recall Equation 1:
xn=∑
k
=1K
A
k
ei
ω
k
n+zn
x
n
k
K
1
A
k
ω
k
n
z
n
(9)
Using the decomposition idea of
Equation 7, it can be shown (for example, see
Hayes:96) that similar to the decomposition of a single
tone in
Equation 8:
R
x
=
R
s
+
R
n
R
x
R
s
R
n
(10)
R
x
=∑
k
=1K|
A
1
|2
e
k
e
k
H+
σ
k
2I
R
x
k
K
1
A
1
2
e
k
e
k
σ
k
2
I
where
e
k
=(
1ei
ω
k
ei
ω
k
2…e(−i)
ω
k
(M−1)
)
e
k
1
ω
k
ω
k
2
…
ω
k
M
1
. Again, the first term of the sum is the signal
component and the second is the noise component. We can
rewrite
Equation 10 compactly as:
R
x
=EΛEH+ω2I
R
x
E
Λ
E
ω
2
I
(11)
E=(
e
1
…
e
k
)
E
e
1
…
e
k
where
EE is a M
× K matrix.
Λ=(
|
A
1
|2 0
|
A
2
|2 0
⋱ ⋮
|
A
k
|20
00…00
)
Λ
A
1
2
0
A
2
2
0
⋱
⋮
A
k
2
0
0
0
…
0
0
where
ΛΛ is
a M × M matrix. Note how the autocorrelation matrix
nicely decomposes, into the
signal subspace and
noise subspace. We can write
Equation 11 in a more algebraic form:
R
x
=∑
i
=1K(
λ
i
+
σ
ω
2)
u
i
u
i
H+∑
i
=K+1M
σ
ω
2
u
i
u
i
H
R
x
i
K
1
λ
i
σ
ω
2
u
i
u
i
i
M
K
1
σ
ω
2
u
i
u
i
from which we can similarly collect the terms as:
U
s
=
u
1
…
u
K
, signal eigenvectors M × K
U
s
u
1
…
u
K
, signal eigenvectors M × K
U
n
=
u
K
+
1
…
u
M
, signal eigenvectors M × K
U
n
u
K
+
1
…
u
M
, signal eigenvectors M × K
Note that what we have done so far is identify the
subspace where the signal belongs to. Recall
from linear algebra that there can be infinitely many
different sets of basis which can span the same subspace
(for example, think of rotation in the subspace, which we
will take advantage of later in ESPRIT). So we still need
a little more work before we can get the tones out of the
autocorrelation matrix.
Think of how we can take advantage of this decomposition
idea by using subspace projection. For example, define
P
s
=
U
s
U
s
H
P
s
U
s
U
s
(12)
P
n
=
U
n
U
n
H
P
n
U
n
U
n
(13)
But first let's explore a few techniques that have been
well-developed in the literature.
This was developed in 1973 by Pisarenko (Pisarenko:73),
published in a geophysics journal in the UK. The main idea
uses Caratheodory's Theorem (this theorem gives a bound on how
big the set we need to be able to capture the dynamics of
parameters of interest), which is of fundamental importance in
estimation theory (Billingsley:86). This algorithm is
intuitive, but is very sensitive to noise. First set
M=K+1
M
K
1
such that
dim
U
s
=K
dim
U
s
K
and
dim
U
n
=1
dim
U
n
1
.
Suppose we have the eigenvalue decomposition of the
autocorrelation matrix as in Equation 4. There is only one noise eigenvalue/vector,
denoted by
λ
min
=
σ
ω
2
λ
min
σ
ω
2
and
u
min
u
min
. This eigenvector is orthogonal to subspace of the
signal, or in other words, to every basis function of the
signal subspace:
u
min
⊥
u
signal
⇒
u
min
annihilates the signal components
u
min
⊥
u
signal
u
min
annihilates the signal components
and vice versa. Recall the shorthand that we have used
e
k
=(
1ei
ω
k
ei
ω
k
2…e(−i)
ω
k
(M−1)
)
e
k
1
ω
k
ω
k
2
…
ω
k
M
1
. Then the above statement is equivalent to:
e
k
H
u
min
=∑
k
=0K
u
min
ke−(i
ω
i
k)=0
e
k
u
min
k
K
0
u
min
k
ω
i
k
0
(14)
or,
U
min
eiω=∑
k
=0K
u
min
ke−(i
ω
i
k)
U
min
ω
k
K
0
u
min
k
ω
i
k
(15)
has zeros at frequencies
ω
k
ω
k
, for
1…K−1
1
…
K
1
. In other words:
U
min
z=∑
k
=0K
u
min
kz−k=∏
k
=0K1−eiωkz-1
U
min
z
k
K
0
u
min
k
z
k
k
K
0
1
ω
k
z
-1
(16)
which is usually called the
annihilating filter.
The annihilating filter of Equation 16 has zeros on the unit circle at angles
corresponding to the frequencies of the signal. We can
find the frequencies
ω
k
ω
k
, for
1…K−1
1
…
K
1
, by finding the zeros of
U
min
eiω=∑
k
=0K
u
min
ke−(i
ω
i
k)
U
min
ω
k
K
0
u
min
k
ω
i
k
.
Suppose that the eigenvectors
u
i
u
i
are unit norm.
R
x
u
i
=
λ
i
u
i
R
x
u
i
λ
i
u
i
u
i
H
R
x
u
i
=
λ
i
u
i
H
u
i
=
λ
i
u
i
R
x
u
i
λ
i
u
i
u
i
λ
i
u
i
H
R
x
u
i
=
u
i
(∑
k
|
A
k
|2
e
k
e
k
H+
ω
w
2I)=
λ
i
u
i
R
x
u
i
u
i
k
A
k
2
e
k
e
k
ω
w
2
I
λ
i
∑
k
=1K|
A
k
|2|
e
k
H
u
k
|=
λ
i
−
σ
w
2
k
K
1
A
k
2
e
k
u
k
λ
i
σ
w
2
We can estimate the magnitude
|
A
i
|
A
i
from the above equation.
Let
xn=Asin
ω
1
n+φ+wn
x
n
A
ω
1
n
φ
w
n
where the noise is white. First estimate the
autocorrelation sequence
r
⁁
x
k
r
⁁
x
k
from the available data for
k=012
k
0
1
2
. Now form the autocorrelation matrix and
eigen-decomposition:
R
x
=(
r
x
0
r
x
1
r
x
2
r
x
1
r
x
0
r
x
1
r
x
2
r
x
1
r
x
0
)
R
x
r
x
0
r
x
1
r
x
2
r
x
1
r
x
0
r
x
1
r
x
2
r
x
1
r
x
0
(17)
R
x
=UΛUH
R
x
U
Λ
U
U=
u
1
u
2
u
3
U
u
1
u
2
u
3
Suppose
λ
1
≥
λ
2
>
λ
3
λ
1
λ
2
λ
3
. From the annihilating filter from the eigenvector
of the minimum eigenvalue
u
min
=
u
3
u
min
u
3
:
U
min
z=∑
k
=02
u
min
kz−k
U
min
z
k
2
0
u
min
k
z
k
find the roots which will be at (or hopefully near)
z=e±i
ω
1
z
±
ω
1
From Prop. 3.1, we
can create an artificial spectra by evaluating the
annihilating filter at different frequencies. So we construct
e
i
ω
0
e
ω
0
for different frequencies
ω
0
ω
0
and evaluate:
Pei
ω
0
=1|
e
ω
0
H
u
min
|2
P
ω
0
1
e
ω
0
u
min
2
(18)
If we take different values of
ω
0
ω
0
, then we can plot what is called a
pseudo
spectra, which has peaks where the zeros are. This is
a good way to see how the spectra might look, instead of using
the nonparametric approach such as the periodogram.
This was developed by Schmidt (Schmidt:81), originally
developed around 1979. The main idea is really just using
averaging to improve the performance of the estimator of
Pisarenko. Now we pick
M>K+1
M
K
1
. Let
(
λ
1
≥
λ
2
≥…≥
λ
K
)≥(
λ
K
+
1
≥…≥
λ
M
)
λ
1
λ
2
…
λ
K
λ
K
+
1
…
λ
M
along with their associated eigenvectors. In the above
statement, the left-hand side is the K
K signal eigenvalues and the right-hand side is the
M−K
M K
noise eigenvalues. Now instead of just one
annihilating filter, we construct M−K
M K such
noise eigenfilters:
U
i
z=∀
i
,i=K+1…M:∑
m
=0M−1
u
i
mz−m
U
i
z
i
i
K
1
…
M
m
M
1
0
u
i
m
z
m
(19)
Observe that every single one of these filters will have
M−1 M
1 roots by the fundamental
theorem of calculus,
KK of which
are common, and are at the angles corresponding to the
frequencies of the signal in consideration. The other
M−K−1
M
K
1
zeros are spurious, and some
of them can be close to the unit circle.
This is more appropriately a "trick" than a
proposition. All the
u
i
u
i
's are orthogonal to the signal subspace, so we
find the KK common zeros by
averaging!
After averaging, we simply do the same computation as we did
for Pisarenko.
Similarly we can also generate the pseudo spectra based on
MUSIC. So we construct
e
i
ω
0
e
ω
0
for different frequencies
ω
0
ω
0
and evaluate:
Pei
ω
0
=1∑
k
=K+1M|
e
ω
0
H
u
k
|2
P
ω
0
1
k
M
K
1
e
ω
0
u
k
2
(20)
or using our projection idea in
Equation 12, we get a nice compact formula:
Pei
ω
0
1
e
ω
0
H
P
n
e
ω
0
P
ω
0
1
e
ω
0
P
n
e
ω
0
(21)
If we take different values of
ω
0
ω
0
we can plot what is called a pseudo spectra, which
has peaks where the zeros are. This is a good way to see
how the spectra might look, instead of using the
nonparametric approach such as the periodogram.
We give an example of the pseudo spectra
when there are 3 exponential components in the deisred
signal in Figure 1. Note that
they all have common zeros at the frequency of the desired
signal.
There is a long story about this one. It is not clear whether
Paulraj or Roy (Roy:89) was the first inventor of ESPRIT. We
start our development by considering a cute "trick" that is
the basic idea behind ESPRIT.
Consider the case when
K=1
K
1
, and no noise. Consider a vector of
observations of
xn=
A
1
ei
ω
1
n
x
n
A
1
ω
1
n
, called x
x:
x=
x
0
x
1
…
x
N
−
1
x
x
0
x
1
…
x
N
−
1
x=
A
1
1ei
ω
k
ei
ω
k
2…ei
ω
k
M−1
x
A
1
1
ω
k
ω
k
2
…
ω
k
M
1
Now we partition this vector as follows:
x=
x
0
x
1
…
x
N
−
2
x
N
−
1
x
x
0
x
1
…
x
N
−
2
x
N
−
1
where we will set
s
1
s
1
eqaul to all the terms but the last,
s
1
=
x
0
x
1
…
x
N
−
2
s
1
x
0
x
1
…
x
N
−
2
and we will also set
s
2
s
2
eqaul to all the terms but the first in the
above expresion,
s
2
=
x
1
…
x
N
−
2
x
N
−
1
s
2
x
1
…
x
N
−
2
x
N
−
1
Now note that:
s
1
=ei
ω
1
s
2
s
1
ω
1
s
2
So in the noisy case, we can use some sort of averaging
to get a better estimate of what
ω
1
ω
1
is. How can we extend this to the multiple
signal case?
To extend it to the multiple signal case, we start again with
the autocorrelation matrix:
R
x
=UH
Λ
x
U
R
x
U
Λ
x
U
(22)
Now let us define two matrices, each having
(M−1)
×
M
M
1
×M dimensions:
Γ
1
=
I
M
−
1
0
(
M
−
1
)
×
1
Γ
1
I
M
−
1
0
(
M
−
1
)
×
1
(23)
Γ
2
=
0
(
M
−
1
)
×
1
I
M
−
1
Γ
2
0
(
M
−
1
)
×
1
I
M
−
1
(24)
which we will use to select the first and last
M−1
M
1
columns of an (M × M) matrix respectively. We
call these the
selector matrices.
We use them as follows:
S
1
=
Γ
1
U
S
1
Γ
1
U
(25)
S
2
=
Γ
2
U
S
2
Γ
2
U
(26)
Now we use a very very very important property.
For the matrices defined in Equation 25, for every k k
denoting the different frequency components, we have:
(
Γ
1
e
ω
k
)ei
ω
k
=
Γ
2
e
ω
k
Γ
1
e
ω
k
ω
k
Γ
2
e
ω
k
which we can compactly write as:
(
Γ
1
U)Φ=
Γ
2
U
Γ
1
U
Φ
Γ
2
U
(27)
Where
Φ=diagei
ω
1
ei
ω
2
…ei
ω
k
Φ
diag
ω
1
ω
2
…
ω
k
.
From here we can see the light at the end of the tunnel.
Recall Remark 2.1,
where we argued that we have so far only identified the signal
up to the subspace and that it is invariant to transformation
by a unitary matrix on the subspace. Let us do a little more
work to tie Theorem 5.2
(Rational Invariance) with the favorite engineering
tool in linear algebra, namely the eigenvalue decomposition,
using some unitary matrix
T
T.
Γ
1
(UT)Φ=
Γ
2
UT
Γ
1
U
T
Φ
Γ
2
U
T
(28)
Γ
1
U(TΦTH)=
Γ
2
U
Γ
1
U
T
Φ
T
Γ
2
U
(29)
where we can label
(TΦTH)
T
Φ
T
as eig. In the first line we apply a unitary
transformation on
U
U, and in the second line we use the property that
TTH=I
T
T
I
.
Now we are ready to unconfuse ourselves by summarizing the
steps in the algorithm:
-
Estimate
r
x
k
r
x
k
and form
R
xx
R
xx
of size at least K × K.
-
Calculate eigendecomposition
R
x
=UH
Λ
x
U
R
x
U
Λ
x
U
. Pick the eigenvectors corresponding to the
KK largest eigenvalues and
form the matrix
U
s
=
u
1
…
u
K
U
s
u
1
…
u
K
. The signal subspace is spanned by this matrix.
-
Solve the equation
(
Γ
1
U
s
)Φ=
Γ
2
U
s
Γ
1
U
s
Φ
Γ
2
U
s
, get Φ
Φ.
-
Compute the KK eigenvalues
of Φ
Φ. In the ideal case the eigenvalues are:
∀
k
,k=1…K:ei
ω
k
k
k
1
…
K
ω
k
.
In solving the rotational invariance formula of Theorem 5.2, we can also
use the least-squares solution. The original problem is:
(
Γ
1
U
s
)Φ=
Γ
2
U
s
Γ
1
U
s
Φ
Γ
2
U
s
which least-squares solution is given by:
Φ=
U
s
H
T
1
H
T
1
U
s
-1
U
s
T
1
T
2
U
s
Φ
U
s
T
1
T
1
U
s
-1
U
s
T
1
T
2
U
s
(30)
And we can apply what we know about least-squares solutions
to analytically derive the mean-square error of this estimate.
We can use the frequency esitmation algorithms which we have
developed in the context of smart antennas, especially
applicable to narrowband wireless
communication systmes.
Consider Figure 2. Suppose we have
MM equi-spaced antennas on a line
configuration of distance dd from
each other, and an incoming waveform at angle
θθ. We take a spatial
snapshot of the incoming signal, by taking
MM simultaneous samples at the
MM antennas in the array. The
dashed lines illustrate the incident wavefront, in which the
angles are equal. It can then be shown that the spatial samples
will be a function of the center frequency, the speed of light
and the incident angle θθ.
We can write the samples at antenna
mm at time
nn when there are
KK incoming signals as:
∀
m
,m=0…M−1:xmn=∑
k
=1K
s
k
nei
ω
k
m+wmn
m
m
0
…
M
1
x
m
n
k
K
1
s
k
n
ω
k
m
w
m
n
(31)
where
ω
k
=2πλdcos
θ
k
ω
k
2
λ
d
θ
k
and
s
k
n
s
k
n
is the message of the
kk-th signal.
Note that by using Equation 31 we have
related the problem of finding the direction of arrival (DoA)
of the different signals. This is especially useful for
beamforming and interference rejection, where we would like to
be able to separate out interfering users by the DoAs.
When the signals of interest are narrowband, that means that
the message component
s
k
n
s
k
n
changes very slowly with time, and thus spatial translation.
Therefore we can consider the message component to be a
constant:
∀
m
,m=0…M−1:xm=∑
k
=1K
s
k
ei
ω
k
m+wm
m
m
0
…
M
1
x
m
k
K
1
s
k
ω
k
m
w
m
(32)
which is a very familiar problem statement for frequency
estimation.
Note that from the expression
ω
k
=2πλdcos
θ
k
ω
k
2
λ
d
θ
k
we can see that we cannot disambiguate signals
coming in from the "front" of the array and from the "back"
of the array. Thus in real-world implementation, other
configurations are also considered (such as circular
{Zoltowski:94} or triangular {Rappaport:99}).
Now we are ready to go back to our task of separating the
different users. WLOG consider the
k=1
k
1
user to be the "good" user, and all the other
signals to be the interfering users. What we would like to do
is to be able to filter reject the signals of the bad users
while passing the signal of the good user.
Before we continue further, consider the solution if the form
of linear combining of the input signals from the different
receive antennas:
yn=∑
m
=0M−1hmn¯xmn
y
n
m
M
1
0
h
m
n
x
m
n
(33)
or in vector form:
yn=
h
M
H
x
n
y
n
h
M
x
n
(34)
We can recognize this as the ever-familiar FIR filter,
except that now it operates in space and time. Therefore this
is called a
space-time FIR filter.
Consider the effect on the samples of signals:
yn=
h
M
H
x
n
y
n
h
M
x
n
(35)
yn=
s
k
n(
h
M
H
e
ω
k
)+
h
M
H
v
n
y
n
s
k
n
h
M
e
ω
k
h
M
v
n
(36)
where
e
ω
k
e
ω
k
is defined as usual, and contains the complex
exponential components.
We first write our problem as a constrained optimization
problem:
min
h
M
h
M
E|yn|2
h
M
E
y
n
2
(37)
subject to
h
M
H
e
ω
1
=1
h
M
e
ω
1
1
where the first user is the desired user. First we rewrite
the optimization problem as follows:
min
h
M
h
M
E
h
M
H
R
x
h
M
h
M
E
h
M
R
x
h
M
(38)
subject to
h
M
H
e
ω
1
=1
h
M
e
ω
1
1
From here using the Lagrange method we obtain the solution,
which is given by
h
M
opt
=
R
x
-1
e
ω
1
e
ω
1
H
R
x
e
ω
1
h
M
opt
R
x
-1
e
ω
1
e
ω
1
R
x
e
ω
1
(39)
We can also show the pseudo-spectra given by the
minimum-variance beamforming algorithm:
S
⁁
ei
ω
0
=1
e
ω
0
H
R
x
-1
e
ω
0
S
⁁
ω
0
1
e
ω
0
R
x
-1
e
ω
0
(40)
We can use our previously developed frequency estimators to
find the angles of arrival, and then designing filters as we
like. We give an example of the performance in the below
figure.
As with many other fields, there is a high degree of
factionalism in frequency estimation. European researchers
tend to be biased towards MUSIC, and Americans towards
ESPRIT. For evidence, see any frequency estimation textbook!
A good paper which compares both in the context of direction
estimation is by Kangas, Stoica, and Soderstrom (Kangas:94).
It is shown that in the intermediate range of modeling error,
ESPRIT outperforms MUSIC. In general MUSIC is more accurate
than ESPRIT.
-
P.Billingsley, Probability and Measure, second edition,
John Wiley and Sons, New York, 1986.
-
M. Hayes, Digital Signal Processing and Modeling,
Wiley,New Yourk, NY, 1996.
-
V.F. Pisarenko, The Retrieval of Harmonics from a
Covariance Function, Geophysics J. Roy. Astron. Soc. 33
(1973), 347-366.
-
R. Roy and T. Kailath, ESPRIT - Estimation of Signal
Parameters via Rotational Invariance Techniques, IEEE
Transactions on Acoustics, Speech, and Signal Processing.
ASSP-37 (1989), 984-995.
-
R.O Schmidt, A Signal Subspace Approach to Multiple
Emitter Location and Spectral Estimation, Ph.D. thesis,
Stanford University, Stanford, CA, 1981.
