Defining Correlation Functions and Power Spectra

for Multirate Random Processes

Charles W. Therrien
Department of Electrical and Computer Engineering
Naval Postgraduate School
Monterey, California 939435000
Email: [email protected]
[W

Abstract This paper proposes a representation for the timelag cross-correlation function of two random processes sampled at
different rates and its related cross-power density spectrum, using
the theory of lattices. The representation of autocorrelation and
the (auto)power spectral density function then follows directly.

I. I NTRODUCTION
\W


The time-lag cross-correlation function for two (continuoustime) random processes is defined as
Rxy (t; ) = E {x(t)y (t )}

(1)

while what will be called the traditional cross-correlation

function is defined as
Rxy (t, t ) = E {x(t)y (t )}

(2)

The time-lag form has a number of advantages. Among these

are that it is easily related to concepts of the spectrogram, various stochastic time-frequency distributions (Wigner, Rihaczek)
as well as the time-frequency correlation or ambiguity
function and cyclic correlation and power density for cyclostationary (or periodically correlated) random processes [1],
[2]. The time-lag form also leads to a clear distinction between
random processes that are wide-sense stationary (wss), those
that are wide-sense cyclostationary (wscs), and those that are
periodic. Stationary processes are independent of t, cyclostationary processes are periodic in t but not in , and periodic
processes are periodic in both variables.
The analysis of random processes sampled at different rates
from a statistical point of view leads to the concept of joint
cyclostationarity between these processes. The representation
for the time-lag cross-correlation function and cross-spectra of
two such random processes, seems to have eluded researchers
in signal processing, however, and has impeded research for
statistical multirate signal processing. We believe that such
representation is fundamental to the statistical analysis of
multirate processes and is necessary for further progress in
the field.
II. C ROSS - CORRELATION
Consider two continuous-time random processes x(t) and
y(t ) as illustrated in Fig. 1. Beginning at t = t = 0, these two
processes are sampled at rates of Fx and Fy samples/second
with sampling intervals of Tx = 1/Fx and Ty = 1/Fy . It

0-7803-8834-8/05/$20.00 2005 IEEE.

Fig. 1. Continuous random processes sampled at different rates. (Note that

the samples do not line up in time except at points indicated by the dashed
lines.)

is reasonable to assume that the sampling rates are integervalued. We further assume that the samples line up at
multiples of T , where T is called the common period and
is defined by the smallest integers Mx and My satisfying
T = Mx T x = My T y

(3)

With some abuse of notation, we denote the resulting discretetime signals as x[mx ] = x(mx Tx ) and y[my ] = x(my Ty )
where mx and my Z (the set of integers).
The traditional and time-lag forms of the cross-correlation
function for the continuous-time random processes are defined
by (1) and (2) above with the relation
= t t

(4)

The definitions apply for all t, t , R (the real numbers).

It is straightforward and intuitive to define a discrete crosscorrelation function for the two sampled random processes in
the traditional form as
Rxy [mx , my ] = E {x[mx ]y [my ]} = Rxy (mx Tx , my Ty )
This cross-correlation function is defined on a rectangular lattice [3] denoted by R and depicted in Fig. 2. The generating
matrix for the lattice is given by



Tx 0

V R = vx vy
=
0
Ty

6010

where mx , ly Z.
Now, in order to represent the discrete, time-lag correlation
function algebraically, observe that the variables t, t and
are related by


 



t
1
0
t
t
=T
=
(8)

1 1

t

7[

The matrix T is an involutory matrix (T1 = T). Using (5),

(8), and (7) in that order, produces




mx
mx
1
= VR TVL
my
ly

Y[
W

Y\
7\

Fig. 2. Cross-correlation function Rxy [mx , my ] defined on rectangular

lattice R .

where vx and vy are the two basis vectors shown in the figure.
A point on the lattice is thus represented by




mx
t
(5)
= VR
= mx v x + m y v y
my
t
Now consider the representation of the time-lag correlation
function Rxy (t; ) in discrete time. The variable t is of the
form t = mx Tx . Also, can only take on values corresponding
to the actual samples available. To see what values are allowed,
we write (4) as
= t t  = mx T x m y T y
Thus for any fixed value of mx , the allowable values of
are separated by integer multiples of Ty . The possible values
for (t, ) can be seen to fall on the (non-rectangular) lattice
L depicted in Fig. 3. The generating matrix for this lattice

It is easily shown by carrying out the matrix multiplication

that


1
0
1
TVL =
=T
VR
1 1
Thus it is seen from the previous equation that my = mx ly .
We can therefore define the discrete, time-lag cross-correlation
function as
Rxy [mx ; l] = E {x[mx ]y [mx l]}
def

where we have dropped the subscript on the lag variable to

simplify notation.
The algebraic definition of cross-correlation in the multirate
case is seen to be the same as would be used for processes
with a single sampling rate, but the interpretation is different.
As seen above, the variables mx and l specify a position on
a lattice L . This lattice degenerates to a rectangular lattice
only in the case of identical sampling rates (Tx = Ty ), where
the generating matrix (6) becomes diagonal.
Taking a similar approach, one can define the crosscorrelation function Ryx as

Ryx [my ; l ] = E {y[my ]x [my l ]}

def

a'

Observe that this function is defined on its own lattice L
and

that the lag variable l is distinct from the lag variable l in Rxy .
The relation between these two cross-correlation functions and
the traditional forms of cross-correlation is given by

b'
Tx

Rxy [mx ; l]
b

vx

vy

Fig. 3. Time-lag cross-correlation function Rxy [m; l] represented on a lattice

L .

VL = vx

vy =

Tx
Tx Ty

0
Ty

Rxy [mx , mx l]

Ryx [mx l, mx ] = Ryx

[mx l; l]

(10)

All four of the forms are equivalent in the sense that any one
of them can be used to specify all the possible correlation
values that exist between the two sampled random processes.

Ty

is given by

(9)


(6)

Any point on the lattice can be represented by coordinates





mx
t
(7)
= mx v x + l y v y
= VL
ly

III. I NDEXED C ORRELATION F UNCTION

While the representation (9) is quite general, an important
property is not directly apparent. That property has to do with
stationarity (or more generally cyclostationarity, if it exists) of
the underlying continuous random processes. Refer to Fig. 3
again, and assume that the random processes x(t) and y(t ) are
jointly wide-sense stationary (wss). In this case the continuoustime cross-correlation function Rxy (t; ) is a function of
alone. Thus Rxy (t; ) is a surface that resembles a tunnel with
its axis parallel to the t axis. The contours of equal correlation

6011

are straight lines parallel to the t axis. Because the function

is constant in t, the periodic sampling of the lattice induces
a periodicity in correlation. For example, the points a and
b shown in Fig. 3 have the same value of correlation as the
points a and b because they are at the same distance from
the t axis. Further, observe that this periodicity is not twodimensional. The periodicity occurs along only one direction,
namely that of the t axis.
The periodicity of the cross-correlation function
Rxy [mx ; ly ] is not immediately evident from (9) because the
lattice basis vectors do not align with the t axis. It can also
be shown algebraically, however, as follows. If x(t) and y(t )
are jointly wss then (9) can be written as

Although this function is also defined on a lattice, the lattice

is rectangular, since there is no diversity in the sampling rates.
The second-moment statistics then are as listed in Table I. The
TABLE I
L ISTING OF TIME - DOMAIN STATISTICS FOR MULTI - RATE SIGNAL
PROCESSING .

Definition

Rate

x [m] = E {x[m]}, y [m] = E {y[m]}

Fx , Fy

Rx [m; l] = E {x[m]x [m l]}

Fx

Rxy [m; l] = E {x[m]y [m l]}

mixed

Rxy [mx ; ly ] = E {x[mx ]y [mx ly ]}

c
= Rxy
(mx Tx (mx ly )Ty )
c
Rxy

where
is written here with a single argument ( ) which
represents the difference of the sample times. To demonstrate
the periodicity, let mx = mx + Mx and ly = ly + Mx My
where Mx and My are defined by (3). Then
Rxy [mx ; ly ]

c
= Rxy
(mx Tx mx Ty + ly Ty

+Mx Tx Mx Ty + (Mx My )Ty )

covariance can be defined in an analogous way and is related

to the other quantities as
Cx [m; l]
Cxy [m; l]

m = mx P Mx

(11)

where P = mx /Mx . In other words, P is the number of

full periods covered by mx and m mx mod Mx . Then
it follows from the arguments in the previous paragraph that
choosing
(12)
l = ly P (Mx My )
forces Rxy [m; l] and Rxy [mx ; ly ] to have the same values.
The equivalent correlation function defined by
Rxy [m; l]

c
= Rxy
(m(Tx Ty ) + l Ty )

(13)

0 m Mx 1 , < l <
will be called an indexed correlation function because it
can be viewed as a family of correlation functions of a single
variable (l) indexed by m.
IV. AUTOCORRELATION , M EAN , AND C OVARIANCE
For a single random process x[n], it is natural to define the
discrete, time-lag autocorrelation function as
Rx [mx ; l] = E {x[mx ]x [mx l]}
def

(14)

= Rxy [m; l]

x [m]y [m

(a)

l] (b)

(15)

Discrete time processes with constant or periodic mean, and

with correlation functions of the form (13) are said to be jointly
wide-sense cyclostationary (wscs).

added terms

By virtue of (3), the added terms sum to 0 and the periodicity

is demonstrated.
For the case of wss processes, the variables mx and ly can
therefore be confined to the first period of the lattice (i.e., the
period containing the point mx = ly = 0). That is, given
arbitrary values for the variables mx and ly , we can find
equivalent values m and l in the first period of the lattice in
the sense that Rxy [m; l] = Rxy [mx ; ly ]. To confine the indices
to the first period, choose m as

= Rx [m; l] x [m]x [m l]

V. S PECTRAL R EPRESENTATION FOR J OINTLY WSS

P ROCESSES
A. Cross-Power Spectral Density
The two-dimensional (2-D) cross-power spectral density
function of the continuous signals x(t) and y(t) can be defined
as
 


c

Rxy (t; ) e2f t e2f dt d (16)
Sxy (f , f ) =

The following observations can be made:

1) If x(t) and y(t) are jointly wss, then Rxy is constant
c
(f  , f ) has
with respect to t. This implies that Sxy
support only on the f axis, where it is singular.
2) If x(t) and y(t) are jointly wscs with period P , then
c
(f  , f )
Rxy (t; ) is periodic in t. This implies that Sxy

has support on an infinite set of lines f = k 2/P for
k Z.
These are the two cases of primary interest since wss or wscs
processes sampled at different rates become jointly wscs.
It follows from the theory of sampling on a lattice (e.g.,
[4]) that the corresponding sampled version of Rxy (t; ) has
its spectrum repeated at points on a lattice in frequency space
with generating matrix U = (VLT )1 . This lattice is known as
the reciprocal lattice of L . From (6), the generating matrix
of the reciprocal lattice is given by

 

1
Ty Ty Tx
F x Fx Fy
=
(17)
U=
Tx
0
Fy
Tx Ty 0
Let us consider the case where the underlying processes
x(t) and y(t) are jointly wss. In this case the spectrum,
before sampling, has support only on the f axis, where it

6012

is singular (impulsive). The spectrum corresponding to the

sampled correlation function thus has the spectrum repeated
along lines parallel to the f axis with centers at points on the
reciprocal lattice. This is illustrated in Fig. 4(a) where the basis
vectors corresponding to (17) are depicted as ux and uy . The
parallelogram defines one period of the spectrum. Notice that
the aliasing that normally occurs when the correlation function
is sampled in one dimension is less severe because the spectral
components are spread out in the plane. This is to be expected
because the sampling in two dimensions effectively increases
c
( ).
the sampling rate of Rxy
For the sampled signals the cross-power density spectrum
can be defined as
def

(k)
Sxy
() =

M
x 1


Rxy [m; l]el e2km/Mx

(18)

m=0 l=

and used in place of (16). This spectral density function is

depicted in Fig. 4(b).

This can be thought of as a set of one-dimensional spectra

defined by transforming the indexed correlation function, and
is analogous to the spectrogram.
B. Auto-Power Spectral Density
When the random process is wss, the autocorrelation function defined by Table I is independent of m. The power
spectral density function is then defined in the usual way as
the (one-dimensional) Fourier transform of the autocorrelation
function. In our two-dimensional k, representation, this
function has support only on the axis, where it is singular.
If x is wscs, then a two-dimensional cyclic spectrum
(k)
Sx () can be defined in analogy with the definition for
(k)
Sxy (). In this definition Rxy is replaced by Rx , and Mx
becomes the cyclic period. This function is periodic in two
dimensions with centers occuring at points on a rectangular
lattice defined by the sampling frequency Fx . In addition, this
function has a positivity constraint which can be written as
Sx(k) () 0;

k=0

(0)

The term Sx () can be thought of as the stationary

component of the spectrum.

I 

VI. D ISCUSSION

)[

The representation of the correlation function in the timelag form (as opposed to the traditional time-time form) has
important advantages for the analysis of random processes
sampled at different rates. While the latter representation
is straightforward, the former has been elusive. The correct
definition is clear however, if one considers correlation as
represented on a two-dimensional lattice. The corresponding
spectral representation fits well with established descriptors
for processes exhibiting cyclostationarity.
In future publications we will demonstrate how this representation naturally fits in to the extension of linear optimal
filtering to the multirate case both in time and frequency
domains. Some previous publications [5], [6] have addressed
these problems without the benefit of the results discussed
here.

X[
)\

X\
I

)[)\

D
kN
0[

R EFERENCES
2
1

E
Fig. 4.

Two-dimensional cross-spectral density function.

An alternative representation for jointly wss random processes with different rates is the time-dependent cross-power
spectrum defined as
Sxy [m, ) =


l=

Rxy [m; l]e2km/Mx

[1] L. L. Scharf, B. Friedlander, P. Flandrin, and A. Hanssen, The Hilbert

space geometry of the stochastic Rihaczek distribution, in Proceedings
of the 35th Asilomar Conference on Signals, Systems, and Computers,
November 2001, pp. 720725, (Pacific Grove, CA).
[2] C. W. Therrien, Some considerations for statistical characterization of
nonstationary random processes, in Proceedings of the 36th Asilomar
Conference on Signals, Systems, and Computers, November 2002, pp.
15541558, (Pacific Grove, CA).
[3] J. W. S. Cassels, An Introduction to the Geometry of Numbers. New
York: Springer-Verlag, 1959.
[4] D. E. Dudgeon and R. M. Mersereau, Multidimensional Digital Signal
Processing. Upper Saddle River, New Jersey: Prentice Hall, Inc., 1984.
[5] R. Cristi, D. A. Koupatsiaris, and C. W. Therrien, Multirate filtering
and estimation: The multirate Wiener filter, in Proceedings of the 34th
Asilomar Conference on Signals, Systems, and Computers, November
2000, pp. 450454, (Pacific Grove, CA).
[6] R. J. Kuchler and C. W. Therrien, Optimal filtering with multirate observations, in Proceedings of the 37th Asilomar Conference on Signals,
Systems, and Computers, November 2002, pp. 12081212, (Pacific Grove,
CA).

6013