Giannakis, G.B. “Cyclostationary Signal Analysis”
Digital Signal Processing Handbook
Ed. Vijay K. Madisetti and Douglas B. Williams
Boca Raton: CRC Press LLC, 1999
c

1999byCRCPressLLC
17
Cyclostationary Signal Analysis
Georgios B. Giannakis
University of Virginia
17.1 Introduction
17.2 Definitions, Properties, Representations
17.3 Estimation, Time-Frequency Links, Testing
Estimating Cyclic Statistics
•
Links with Time-Frequency Rep-
resentations
•
Testing for Cyclostationarity
17.4 CS Signals and CS-Inducing Operations
Amplitude Modulation
•
Time Index Modulation
•
Fractional
Sampling andMultivariate/Multirate Processing
•
Periodically
Varying Systems
17.5 Application Areas

× 1 vector θ
s
(with d
θ
s
<N), the problem is ill-posed because
c

1999 by CRC Press LLC
adding a new datum, say x(n
0
), adds a new unknown, s(n
0
), tobe determined. Thus, only structured
nonstationarities can be handled when rapid variations are present; and only for classes of finitely
parameterized nonstationary processes can reliable statistical descriptors be computed using a single
time series. One such class is that of (wide-sense) cyclostationary processes which are characterized
by the periodicity they exhibit in their mean, correlation, or spectral descriptors.
An overview of cyclostationary signal analysis and applications are the main goals of this sec-
tion. Periodicity is omnipresent in physical as well as manmade processes, and cyclostationary
signals occur in various real life problems entailing phenomena and operations of repetitive nature:
communications [15], geophysical and atmospheric sciences (hydrology [66], oceanography [14],
meteorology [35], and climatology [4]), rotating machinery [43], econometrics [50], and biological
systems [48].
In 1961 Gladysev [34] introduced key representations of cyclostationary time series, while in 1969
Hurd’s thesis [38] offered an excellent introduction to continuous time cyclostationary processes.
Since 1975 [22], Gardner and co-workers have contributed to the theory of continuous-time cyclo-
stationary signals, and especially their applications to communications engineering. Gardner [15]
adopts a “non-probabilistic” viewpoint of cyclostationarity (see [19] for an overview and also [36]
and [18] for comments on this approach). Responding to a recent interest in digital periodically

x
(n) = µ
x
(n + lP), c
xx
(n; τ) = c
xx
(n + lP; τ),or,¯c
xx
(n; τ) =¯c
xx
(n + lP; τ),
∀n, l ∈ Z. The smallest of all such P s is called the period. Being periodic, they all accept Fourier
Series expansions over complex harmonic cycles with the set of cycles defined as: A
c
xx
:= {α
k
=
2πk/P , k= 0,...,P− 1}; e.g., c
xx
(n; τ)and its Fourier coefficients called cyclic correlations are
related by:
c
xx
(n; τ)=
P−1

k=0
C

P
kn
.
(17.1)
Strict sense cyclostationarity, or, periodic (non-) stationarity, can also be defined in terms of
probability distributions or density functions when these functions vary periodically (in n). But
c

1999 by CRC Press LLC
the focus in engineering is on periodically and almost periodically correlated
1
time series, since real
data are often zero-mean, correlated, and with unknown distributions. Almost periodicity is very
common in discrete-time because sampling a continuous-time periodic process will rarely yield a
discrete-time periodic signal; e.g., sampling cos(ω
c
t + θ)every T
s
seconds results in cos(ω
c
nT
s
+ θ)
for which an integer period exists only if ω
c
T
s
= 2π/P. Because 2π/(ω
c
T

k
; τ) = lim
N→∞
1
N
N−1

n=0
c
xx
(n; τ)e
−jα
k
n
.
(17.2)
The set of cycles, A
c
xx
(τ ) := {α
k
: C
xx
(α
k
; τ) = 0 , −π<α
k
≤ π}, must be countable and the
limit is assumed to exist at least in the mean-square sense [9, Thm. 1.15].
Definition 17.2 and Eq. (17.2) for ACS, subsume CS Definition 17.1 and Eq. (17.1). Note that

N
N−1

n=0
µ
x
(n)e
−jαn
=
µ
s
2
[δ(α− ω
0
) + δ(α + ω
0
)]+µ
v
δ(α) ,
(17.4)
wherein(17.4) we used the definition of Kronecker’s delta
lim
N→∞
1
N
N−1

n=0
e
jαn

E{X
N
(α)};
thus, the cyclic mean can be interpreted as an averaged DFT and ω
0
can be retrieved by picking the
peak of |X
N
(ω)| for ω = 0.
Case 2: µ
s
= 0.From(17.3) we find the correlation c
xx
(n; τ) = c
ss
(τ )[cos(2ω
0
n + ω
0
τ)+
cos(ω
0
τ)]/2 + c
vv
(τ ). Because c
xx
(n; τ) is periodic in n, x(n) is (second-order) CS with cyclic
correlation [c.f. (17.2) and (17.5)]
C
xx


δ(α) .
(17.6)
The set of cycles is A
c
xx
(τ ) ={±2ω
0
, 0} provided that c
ss
(τ ) = 0 and c
vv
(τ ) = 0. The set A
c
xx
(τ )
is lag-dependent in the sense that some cycles may disappear while others may appear for different
τs. To illustrate the τ -dependence, let s(n) be an MA process of order q. Clearly, c
ss
(τ ) = 0 for
|τ| >q, and thus A
c
xx
(τ ) ={0} for |τ| >q.
The CS process in (17.3) is just one example of signals involving products and sums of stationary
processes such as s(n) with (almost) periodic deterministic sequences d(n), or, CS processes x(n).
For such signals, the following properties are useful:
Property 1 Finite sums and products of ACS signals are ACS. If x
i
(n) is CS with period P

(n) and y
2
(n) equals the least common multiple of
the P
i
s. Similarly, finite sums and products of stationary processes with deterministic (almost) periodic
signals are also ACS processes.
As examples of random–deterministic mixtures, consider
x
1
(n) = s(n) + d(n) and x
2
(n) = s(n)d(n) ,
(17.7)
where s(n) is zero-mean, stationary, and d(n) is deterministic (almost) periodic with Fourier Series
coefficients D(α). Time-varying correlations are, respectively,
c
x
1
x
1
(n; τ)= c
ss
(τ ) + d(n)d(n + τ) and c
x
2
x
2
(n; τ)= c
ss

α-domain. To reiterate the dependence on τ, notice that if d(n) is a periodic ±1 sequence, then
c
x
2
x
2
(n; 0) = c
ss
(0)d
2
(n) = c
ss
(0), and hence periodicity disappears at τ = 0.
ACS signals appear often in nature with the underlying periodicity hidden, unknown, or inacces-
sible. In contrast, CS signals are often man-made and arise as a result of, e.g., oversampling (by a
known integer factor P) digital communication signals, or by sampling a spatial waveform with P
antennas (see also Section 17.4).
Both CS and ACS definitions could also be given in terms of the Fourier Transforms (τ → ω)
of c
xx
(n; τ) and C
xx
(α; τ), namely the time-varying and the cyclic spectra which we denote by
S
xx
(n; ω) and S
xx
(α; ω). Suppose c
xx
(n; τ)and C

(α
k
; ω)e
jα
k
n
(17.10)
S
xx
(α
k
; ω) :=
∞

τ=−∞
C
xx
(α
k
; τ)e
−jωτ
= lim
N→∞
1
N
N−1

n=0
S
xx

−jωn
.
(17.12)
If x(n) is complex ACS, then one also needs
¯
S
xx
(α
k
; ω) := lim
N→∞
N
−1
E{X
∗
N
(−ω) X
N
(α
k
− ω)}.
Both S
xx
and
¯
S
xx
reveal presence of spectral correlation. This must be contrasted to stationary
processeswhose spectral components, X
N

s
xx
.
Before dwelling further on spectral characterization of ACS processes, it is useful to note the diver-
sity of tools available for processing. Stationary signals are analyzed with time-invariant correlations
(lag-domain analysis), or with power spectral densities (frequency-domain analysis). However, CS,
ACS, and generally nonstationary signals entail four variables: (n, τ, α, ω) :=(time, lag, cycle, fre-
quency). Grouping two variables at a time, four domains of analysis become available and their
relationship is summarized in Fig. 17.1. Note that pairs (n; τ)↔ (α; τ),or,(n; ω) ↔ (α; ω),have
τ or ω fixed and are Fourier Series pairs; whereas (n; τ) ↔ (n; ω),or,(α; τ) ↔ (α; ω),haven or
α fixed and are related by Fourier Transforms. Further insight on the links between stationary and
FIGURE 17.1: Four domains for analyzing cyclostationary signals.
cyclostationary processes is gained through the uniform shift (or phase) randomization concept. Let
c

1999 by CRC Press LLC
x(n) be CS with period P , and define y(n) := x(n+ θ),whereθ is uniformly distributed in [0,P)
and independent of x(n). With c
yy
(n; τ):= E
θ
{E
x
[x(n+ θ)x(n+ τ + θ)]}, we find:
c
yy
(n; τ)=
1
P
P−1

. Phase randomization of x(n) in (17.3) leads to a stationary y(n) with correlation found by
substituting α = 0 in (17.6). This leads to c
yy
(τ ) = (1/2)c
ss
(τ ) cos(ω
0
τ)+ c
vv
(τ ), and shows that
if s(n) has multiple spectral peaks, or if s(n) is broadband, then multiple peaks or smearing of the
spectral peak hamper estimation of ω
0
(in fact, it is impossible to estimate ω
0
from the spectrum of
y(n) if s(n) is white). In contrast, picking the peak of C
xx
(α; τ)in (17.6) yields ω
0
, provided that
ω
0
∈ (0,π) so that spectral folding is prevented [33]. Equation (17.13) provides a more general
answer. Phase randomization restricts a CS process only to one cycle, namely α = 0. In other words,
the cyclic correlation C
xx
(α; τ) contains the “stationarized correlation” C
xx
(0; τ) and additional

xx
(τ ). Formally, (17.14) implies:
Property 4 Stationary processes can be viewed as ACS or CS with cyclic correlation C
xx
(α; τ) =
c
xx
(τ )δ(α).
Separation of information bearing ACS signals from stationary ones (e.g., noise) is desired in many
applications and can be achieved based on Property 4 by excluding the cycle α = 0.
Next, it is of interest to view CS signals as special cases of general nonstationary processes
with 2-D correlation r
xx
(n
1
,n
2
) := E{x(n
1
)x(n
2
)}, and 2-D spectral densities S
xx
(ω
1
,ω
2
) :=
FT[r
xx

FIGURE 17.2: Support of 2-D spectrum S
xx
(ω
1
,ω
2
) for CS processes.
Property 5 A CS process with period P is a special case of a nonstationary (harmonizable) process with
2-D spectral density given by
S
xx
(ω
1
,ω
2
) =
P−1

k=−(P−1)
S
xx
(
2π
P
k; ω
1
)δ
D
(ω
2

towards representing CS processes as a superposition of stationary processes. Next we will examine
two such representations introduced by Gladysev [34] (see also [22, 38, 49], and [56]).
We can uniquely write n
0
= nP + i and express x(n
0
) = x(nP + i), where the remainder i takes
values0, 1,...,P−1.Foreachi, define the subprocess x
i
(n) := x(nP +i). In multirate processing,
the P × 1 vector x(n) := [x
0
(n)...x
P−1
(n)]

constitutes the so-called polyphase decomposition of
x(n) [51, Ch. 12]. As shown in Fig. 17.3,eachx
i
(n) is formed by downsampling an advanced copy
of x(n). On the other hand, combining upsampled and delayed x
i
(n)s, we can synthesize the CS
process as:
x(n) =
P−1

i=0

l

xx
becomes
independent of n establishing that x
i
1
(n), x
i
2
(n) are (jointly) stationary with correlation:
c
x
i
1
x
i
2
(τ ) = c
xx
(i
1
; i
2
− i
1
+ τP) , i
1
,i
2
∈[0,P− 1] .
(17.17)

2
=0
S
xx

2π
P
k
1
;
ω − 2πk
2
P

e
j[(
ω−2πk
2
P
)(i
2
−i
1
)+
2π
P
k
1
i
1

−i
1
)
e
−j
2π
P
ki
2
.
(17.19)
Basedon(17.16) through (17.19), we infer that cyclostationary signals with period P can be analyzed
as stationary P × 1 multichannel processes and vice versa. In summary, we have:
Representation 1 (Decimated Components) CS process x(n) can be represented as a P-variate sta-
tionary multichannel process x(n) with components x
i
(n) = x(nP + i), i = 0, 1,...,P − 1. Cyclic
spectra and stationary auto- and cross-spectra are related as in (17.18) and (17.19).
An alternative means of decomposing a CS process into stationary components is by splitting the
(−π, π] spectral support of X
N
(ω) into bands each of width 2π/P [22]. As shown in Fig. 17.4, this
can be accomplished bypassing modulated copies of x(n) through an ideal low-pass filter H
0
(ω) with
spectral support (−π/P, π/P]. The resulting subprocesses ¯x
m
(n) can be shifted up in frequency
and recombined to synthesize the CS process as: x(n) =


(ω) = S
xx

2π
P
(m
1
− m
2
); ω +
2π
P
m
1

, |ω| <
π
P
.
(17.20)
Equation (17.20) suggests that cyclostationary signal analysis is linked with stationary subband pro-
cessing.
Representation 2 (Subband Components) CS process x(n) can be represented as a superposition of P
stationary narrowband subprocesses according to: x(n) =

P−1
m=0
¯x
m
(n) exp(−j2πmn/P ). Auto- and

cyclic quantities, C
xx
(α
k
; τ)and S
xx
(α
k
; ω), which are time-invariant. Indeed, in (17.2) and (17.10)
time-variation is assigned to the Fourier basis.
3
Well-separated samples of such processes are asymptotically independent. Sufficient (so-called mixing) conditions
include absolute summability of cumulants and are satisfied by many real life signals (see [5, 12, Ch. 2]).
c

1999 by CRC Press LLC
17.3.1 Estimating Cyclic Statistics
Firstwewill considerACSprocesseswith known cycles α
k
. Simplerestimatorsfor CSprocessesandcy-
cleestimationmethods will be discussedlaterin the section. If x(n) has nonzeromean, weestimatethe
cyclicmean as in Example17.1 using the normalized DFT:
ˆ
C
xx
(α
k
) = N
−1



n=0
x(n)x(n + τ)e
−jα
k
n
,
ˆc
xx
(n; τ) =

α
k
∈A
c
xx
(τ )
ˆ
C
xx
(α
k
; τ)e
jα
k
n
.
(17.21)
Note that
ˆ

W

ω −
2π
N
n

I
xx

α;
2π
N
n

,
ˆ
S
(ii)
xx
(α; ω) =
M

τ=−M
w(τ)
ˆ
C
xx
(α; τ)e
−jωτ

k
n
.
(17.23)
Estimates (17.21) through (17.23) apply to ACS (and hence CS) processes with a finite number of
known cycles, and rely on the following steps: (1) estimate the time-invariant (or “stationary”)
quantities by dropping limits and expectations from the corresponding cyclic definitions, and (2) use
the cyclic estimates to obtain time-varying estimates relying on the Fourier synthesis Eqs. (17.2)
and (17.10). Selection of the windows in (17.22), variance expressions, consistency, and asymptotic
normality of the estimators in (17.21) through (17.23) under mixing conditions can be found in [11,
12, 24, 39] and references therein.
When x(n) is CS with knowninteger period P, estimationof time-varying correlationsand spectra
becomes easier. Recall that thanks to Representations 1 and 2, not only c
xx
(n; τ)and S
xx
(n; ω), but
the process x(n) itself can be analyzed into P stationary components. Starting with (17.16), it can
be shown that c
xx
(i; τ)= c
x
i
x
i+τ
(0),wherei = 0, 1,...,P− 1 and subscript i + τ is understood
mod(P ). Because the subprocesses x
i
(n) and x
i+τ

−1

P−1
k=0
X
P
(ω)X
P
(2πk/P − ω) exp(−j2πkn/P ), and
then smoothed to obtain a consistent estimate of S
xx
(n; ω).
17.3.2 Links with Time-Frequency Representations
Consistency (and hence reliability) of single record estimates is a notable difference between cyclo-
stationary and time-frequency signal analyses. Short-time Fourier transforms, the Wigner-Ville,
and derivative representations are valuable exploratory (and especially graphical) tools for analyz-
ing nonstationary signals. They promise applicability on general nonstationarities, but unless slow
variations are present and multiple independent data records are available, their usefulness in es-
timation tasks is rather limited. In contrast, ACS analysis deals with a specific type of structured
variation, namely (almost) periodicity, but allows for rapid variations and consistent single record
sample estimates. Intuitively speaking, cyclostationarity provides within a single record, multiple
periods that can be viewed as “multiple realizations.” Interestingly, for ACS processes there is a close
relationship between the normalized asymmetric ambiguity function A(α; τ)[37], and the sample
cyclic correlation in (17.21):
N
ˆ
C
xx
(α; τ)= A(α; τ):=
N−1

˜c
xx∗
(n; τ) := c
xx∗
(n + τ;−2τ):= E{x(n+ τ)x
∗
(n − τ)}
= c
ss
(2τ)exp(j4ω
0
τn)+ c
vv∗
(2τ) ,
(17.26)
exhibits (almost) periodicity and its cycliccorrelation is given by:
˜
C
xx∗
(α; τ)= c
ss
(τ )δ(α−4ω
0
τ)+
c
vv∗
(2τ)δ(α). Assuming c
ss
(τ ) = 0, the latter allows evaluation of ω
0

(17.27)
The
ˆ
˜
C
xx∗
(α; τ)estimate in (17.27) is nothing but the symmetric ambiguity function. Because x(n)
is ACS,
ˆ
˜
C
xx∗
can be shown to be consistent. This provides yet one more reason for the success of
c

1999 by CRC Press LLC
time-frequency representations with chirp signals. Interestingly, (17.27) shows that exploitation of
cyclostationarity allows not only for additive noise tolerance [by avoiding the α = 0 cyclein(17.27)],
but also permits parameter estimation of chirps modulated by stationary multiplicative noise s(n).
17.3.3 Testing for Cyclostationarity
In certain applications involving man-made (e.g., communication) signals, presence of cyclostation-
arity and knowledge of the cycles is assured by design (e.g., baud rates or oversampling factors). In
other cases, however, only a time series {x(n)}
N−1
n=0
is given and two questions arise: How does one
detect cyclostationarity, and if x(n) is confirmed to be CS of a certain order, how does one estimate
the cycles present? The former is addressed by testing hypotheses of nonzero
ˆ
C

ˆ
C
R
xx
(α; τ
1
)...
ˆ
C
R
xx
(α; τ
L
);
ˆ
C
I
xx
(α; τ
1
)...
ˆ
C
I
xx
(α; τ
L
)]

where superscript R(I) denotes real (imaginary) part. Similarly, we define the ensemble vector

xx
(α) will have zero
mean, and
ˆ
D
2c
(α) := ˆc

xx
(α)
ˆ

†
c
(α)ˆc
xx
(α) will be central chi-square. For a given false-alarm rate,
we find from χ
2
tables a threshold  and test [10]
H
0
:
ˆ
D
c
xx
(α) ≥  ⇒ α ∈ A
c
xx

+
2πm
M
)X
∗
N
(ω
2
+
2πm
M
) |
2
1
M

M−1
m=0
| X
N
(ω
1
+
2πm
M
) |
2
1
M


(n), s
2
(n), and v(n) zero-mean,
Gaussian, and mutually independent. To test for cyclostationarity and retrieve the possible periods
present, N = 2, 048 samples were generated; s
1
(n) and s
2
(n) were simulated as AR(1) with variances
σ
2
s
1
= σ
2
s
2
= 2, while v(n) was white with variance σ
2
v
= 0.1. Figure 17.5a shows |
ˆ
C
xx
(α; 0)|
peaking at α =±2(π/8),±2(π/4), 0 as expected, while Fig. 17.5b depicts ρ
xx
(ω
1
,ω


π
−π
S
xx
(ω, ω− α)dω, and for each α, we can view Fig. 17.5a as the (normalized) integral (or
projection) of Fig. 17.5b along each parallel line [40]. Although |
ˆ
C
xx
(α; 0)| is simpler to compute
using the FFT of x
2
(n), ρ
xx
(ω
1
,ω
2
) is generally more informative.
Because cyclostationarity is lag-dependent, as an alternative to ρ
xx
(ω
1
,ω
2
) one can also plot
|
ˆ
C