Almost-cyclostationary signal processing

3 Almost-cyclostationary signal processing Contents 3.1

Introduction.......................................................................................................

82

3.2

Linear filtering ...................................................................................................

82

3.2.1

Linear almost-periodically time-variant systems .............................................

82

3.2.2

Input/output relations in terms of cyclic statistics...........................................

83

3.3

Products of ACS signals .......................................................................................

88

3.4

Supports of cyclic spectra of band limited signals ...................................................

90

3.5

Rice’s representation ...........................................................................................

93

3.6

Sampling and aliasing ......................................................................................... 105

3.7

3.6.1

Cyclic statistics of the sampled signal ........................................................... 105

3.6.2

(Almost-)cyclostationarity of a sampled cyclostationary signal ......................... 108

3.6.3

Bandpass sampling.................................................................................... 109

Multirate processing of discrete-time ACS signals ................................................... 111 3.7.1

Expansion (upsampling) ............................................................................. 113

3.7.2

Sampling ................................................................................................. 116

3.7.3

Decimation (downsampling)....................................................................... 117

3.8

Special topics ..................................................................................................... 119

3.9

Summary........................................................................................................... 120

3.10 Proofs ............................................................................................................... 121 Abstract The processing of almost-cyclostationary signals is addressed. The linear almost-periodically timevariant filtering is considered and input/output relations are derived in terms of cyclic crosscorrelations and cyclic cross-spectra. The cyclic statistics of the signals involved in the Rice’s representation are obtained as special cases. The discrete-time signal obtained by uniformly sampling a continuous-time almost-cyclostationary signal is considered. The cyclic statistical functions of the discrete-time signal are linked to those of the continuous-time signal and aliasing issues are addressed. Multirate processing of discrete-time almost-cyclostationary signals are also considered. Keywords Linear almost-periodically time variant (LAPTV) filtering, Frequency shift (FRESH) filtering, Uniform sampling, Bandpass sampling, Rice’s representation, Multirate systems Cyclostationary Processes and Time Series. https://doi.org/10.1016/B978-0-08-102708-0.00014-5 Copyright © 2020 Elsevier Ltd. All rights reserved.

81

82

Cyclostationary Processes and Time Series

3.1 Introduction Cyclostationarity properties of almost-cyclostationary (ACS) signals are transformed by signal processing. In particular, the linear almost-periodically time-variant (LAPTV) transformations are of great interest. In fact, many man-made signals are obtained by LAPTV transformations of random (possibly wide-sense stationary) data. In addition, the class of the ACS processes is closed under LAPTV transformations. The linear time-invariant (LTI) filtering is as a special case of LAPTV filtering. The LAPTV filtering of ACS signals is introduced in (Dragan, 1969b), (Gardner, 1972a), (Pourahmadi and Salehi, 1983), (Gardner, 1985), (Gardner, 1987f), and (Brown, 1987). It is also considered in follow-on work in (McLernon, 1998), (Franks, 1994), (Spooner, 1992), (Spooner and Gardner, 1994a), and (Napolitano, 1995). Properties of linear periodically time-varying systems are analyzed in (Voychishin and Dragan, 1972b), (Meyer and Burrus, 1975), (Ding and Fahmy, 1990), (Chen and Qiu, 1997), (Mehr and Chen, 2001), and (Mehr and Chen, 2002). Applications are considered in (Ericson, 1981). More general kinds of linear time-variant (LTV) filtering are addressed in (Izzo and Napolitano, 2002a), (Izzo and Napolitano, 2002b), where linear systems that are deterministic in the fraction-of-time probability sense are considered. In general, LTV transformations of ACS signals produce oscillatory ACS signals (Section 14.2.1), (Napolitano, 2016a, Sec. 6.2.1). The ACS signal obtained by the product of ACS signals are characterized in (Gardner, 1987f), (Gardner et al., 2006). The Rice’s representation models real signals as linear periodically time-variant filtering of the complex envelope or, equivalently, of the in-phase and in-quadrature components. See (Gardner, 1987c) for the second-order statistics and (Izzo and Napolitano, 1997) for the higher-order statistics. The problem of sampling a continuous-time ACS signal is treated in (Napolitano, 1995), (Izzo and Napolitano, 1996), where sufficient conditions on the sampling frequency are derived to reconstruct cyclic statistics of the continuous-time signal starting from the discrete-time sequence obtained by uniformly sampling the continuous-time signal. Multirate processing of discrete-time jointly ACS signals is addressed in (Sathe and Vaidyanathan, 1993), (Izzo and Napolitano, 1995), (Izzo and Napolitano, 1998b), (Akkarakaran and Vaidyanathan, 2000), (Aach, 2007), (Napolitano, 2012, Sec. 4.10). In particular, expansion (upsampling) and decimation (downsampling) operations are considered. Note that such operations are liner time variant but not periodically or almostperiodically time variant.

3.2 Linear filtering 3.2.1 Linear almost-periodically time-variant systems Definition 3.1. A linear time-variant (LTV) system with input x(t), output y(t), impulseresponse function h(t, u), and input-output relation

Chapter 3 • Almost-cyclostationary signal processing 83  y(t) =

R

h(t, u) x(u) du

(3.1)

is said to be linear almost-periodically time-variant (LAPTV) if the impulse-response function admits the Fourier series expansion  hσ (t − u) ej 2πσ u (3.2) h(t, u) = σ ∈J

where J is a countable set of possibly incommensurate frequency shifts. By substituting (3.2) into (3.1) we see that the output y(t) can be expressed in the two equivalent forms (Gardner, 1987f), (Franks, 1994)    y(t) = hσ (t) ⊗ x(t) ej 2πσ t (3.3a) σ ∈J

=



 gσ (t) ⊗ x(t) ej 2πσ t

(3.3b)

σ ∈J

where gσ (t)  hσ (t) e−j 2πσ t . Eqs. (3.3a) and (3.3b) can be re-expressed in the frequency domain:  Y (f ) = Hσ (f ) X(f − σ )

(3.4)

(3.5a)

σ ∈J

=



Gσ (f − σ ) X(f − σ )

(3.5b)

σ ∈J

where Hσ (f ) and Gσ (f ) are the Fourier transforms of hσ (t) and gσ (t), respectively. From (3.3a) it follows that a LAPTV system combines LTI filtered versions of frequencyshifted versions of the input signal. For this reason LAPTV filtering is also referred to as frequency-shift (FRESH) filtering (Gardner, 1993). Equivalently, from (3.3b) it follows that a LAPTV system combines frequency shifted versions of LTI filtered versions of the input (Fig. 3.1). In the special case where J ≡ {k/T0 }k∈Z for some period T0 , the system is said to be linear periodically time-variant (LPTV). If J contains only the element σ = 0, then the system is linear time-invariant.

3.2.2 Input/output relations in terms of cyclic statistics Let us consider two LAPTV systems (Fig. 3.2) with input/output relationships  yi (t) = hi (t, u) xi (u) du i = 1, 2 R

(3.6)

84

Cyclostationary Processes and Time Series

FIGURE 3.1 LAPTV systems: Two equivalent realizations and corresponding input/output relations in time and frequency domains (see (3.3a)–(3.5b)).

FIGURE 3.2 LAPTV systems.

whose impulse-response functions admit the Fourier series expansions  hi (t, u) = hσi (t − u) ej 2πσi u i = 1, 2 .

(3.7)

σi ∈Ji

Assumption 3.2. The LAPTV systems (3.6) are such that   hσ (τ ) ∈ L1 (R) i = 1, 2 a. i

(3.8a)

σi ∈Ji

b.

 hσ i

σi ∈Ji

∞

<∞

i = 1, 2 .

(3.8b)

Assumption 3.3. The input signals x1 (t) and x2 (t) have almost-periodic expected values  E {xi (t)} = x ηi ej 2πηi t i = 1, 2 (3.9) ηi ∈E1i

with absolutely convergent Fourier series. Assumption 3.4. The input signals x1 (t) and x2 (t) are jointly ACS with (conjugate) cyclic cross-correlation functions such that

Chapter 3 • Almost-cyclostationary signal processing 85  R α

a.

α∈A12

(∗) x1 x2

 R α

b.

α∈A12

 (τ ) ∈ L1 (R)

(3.10a)



(3.10b)

(∗) x1 x2

∞

<∞

Assumptions 3.2, 3.3, and 3.4 are used to allow the interchange of integral and sum operations in the derivations of the following input/output relationships. Theorem 3.5. Under Assumptions 3.2a and 3.3, the expected value of yi (t) is given by   Hσi (ηi + σi ) x ηi ej 2π(ηi +σi )t i = 1, 2 (3.11) E {yi (t)} = σi ∈Ji ηi ∈E1i

where

 Hσi (f ) 

R

hσi (τ ) e−j 2πf τ dτ .

(3.12)

Proof. See Section 3.10. Theorem 3.6. Under Assumptions 3.2a and 3.4b or under Assumptions 3.2b and 3.4a, the (conjugate) cross-correlation of the outputs yi (t) is given by

 Ry y (∗) (t, τ )  E y1 (t + τ ) y2(∗) (t) 1 2

   α j 2πσ1 τ R (∗) (τ ) e = x1 x2 α∈A12 σ1 ∈J1 σ2 ∈J2 1 +(−)σ2 ⊗ rσα+σ (τ ) ej 2π(α+σ1 +(−)σ2 )t 1 σ2 (∗) τ

(3.13)

denotes convolution with respect to τ , (−) is an optional minus sign that is linked where ⊗ τ to (∗), and  γ −j 2πγ s ds . (3.14) rσ1 σ2 (∗) (τ )  hσ1 (s + τ ) h(∗) σ2 (s) e R

Proof. See Section 3.10. From Theorem 3.5 it follows that E {yi (t)} is an almost-periodic function of t. From Theorem 3.6 it follows that Ry y (∗) (t, τ ) is an almost-periodic function of t whose 1 2

(conjugate) cycle frequencies are obtained by sum and differences of those of Rx

(∗) 1 x2

(t, τ )

and the frequencies of the Fourier series of the two impulse-response functions h1 (t, u) and h2 (t, u). Therefore, from Theorems 3.5 and 3.6 considered for x1 (t) ≡ x2 (t), we have that LAPTV systems transform wide-sense ACS processes into wide-sense ACS processes (with different cyclic features). That is, wide-sense ACS processes are closed under LAPTV transformations.

86

Cyclostationary Processes and Time Series

Corollary 3.7. Under Assumptions 3.2 and 3.4a, the (conjugate) cyclic cross-correlation and the (conjugate) cyclic cross-spectrum of the outputs yi (t) are given by 

  (∗) R α (∗) (τ )  E y1 (t + τ ) y2 (t) e−j 2παt y1 y2 t

  α−σ −(−)σ j 2πσ1 τ 1 2 = (τ ) e ⊗ rσα1 σ2 (∗) (τ ) (3.15) R (∗) τ x1 x2

σ1 ∈J1 σ2 ∈J2



Sα

(∗)

y1 y2

R α (∗) (τ ) e−j 2πf τ dτ y1 y2   α−σ −(−)σ 2 S (∗)1 (f − σ1 ) Hσ1 (f ) Hσ(∗) ((−)(α − f )) = 2

(f ) 

R

x1 x2

σ1 ∈J1 σ2 ∈J2

(3.16)

where in the sums in (3.15) and (3.16), only those σ1 ∈ G1 and σ2 ∈ G2 such that α − σ1 − ∈ A12 give nonzero contribution.

(−)σ2

Proof. See Section 3.10. The results of Corollary 3.7 can be specialized to several cases of interest. a. By specializing (3.15) and (3.16) to the case of two LTI systems, we have   (∗) R α (∗) (τ ) = R α (∗) (τ ) ⊗ h1 (τ ) ⊗ h2 (−τ ) ej 2πατ y1 y2

F ←→ S α

(∗)

y1 y2

x1 x2

(f ) = S α

(∗)

(∗)

x1 x2

(f ) H1 (f ) H2 ((−)(α − f )) .

(3.17) (3.18)

b. If in (3.15) and (3.16) x1 = x2 = x, h1 = h2 = h, y1 = y2 = y, and (∗) is conjugation, then we obtain the input-output relations for LAPTV systems in terms of cyclic autocorrelation functions and cyclic spectra:

  α−σ +σ α j 2πσ1 τ 1 2 (τ ) = (τ ) e ⊗ rσα1 σ2 ∗ (τ ) (3.19) R Ryy ∗ xx ∗ τ σ1 ∈J σ2 ∈J

α Syy ∗ (f ) =

 

σ1 ∈J σ2 ∈J

α−σ1 +σ2 Sxx (f − σ1 ) Hσ1 (f ) Hσ∗2 (f − α) . ∗

(3.20)

c. By specializing (3.16) to the case x1 = x2 = x, h1 = h2 = h LTI system, and y1 = y2 = y, we obtain y(t) = x(t) ⊗ h(t)

(3.21)

α α (∗) ((−)(α − f )) . Syy (∗) (f ) = Sxx (∗) (f ) H (f ) H

(3.22)

and

Analogously, by specializing (3.15), we obtain the inverse Fourier transform of both sides of (3.22) α α α Ryy (∗) (τ ) = Rxx (∗) (τ ) ⊗ rhh(∗) (τ )

(3.23)

Chapter 3 • Almost-cyclostationary signal processing 87

α where the asymmetric (conjugate) ambiguity function rhh (∗) (τ ) can be written in one of the following equivalent forms (see Section 3.10)  α rhh h(s + τ ) h(∗) (s) e−j 2παs ds (∗) (τ )  R   = h(τ ) ⊗ h(∗) (−τ ) ej 2πατ (3.24a)  (3.24b) = H (λ) H (∗) ((−)(α − λ)) ej 2πλτ dλ R

For α = 0 and (∗) = ∗, (3.22) and (3.23) reduce to the input/output relations for LTI systems in terms of power spectral densities and averaged autocorrelations: 0 0 2 Syy ∗ (f ) = Sxx ∗ (f ) |H (f )|

(3.25)

0 0 0 Ryy ∗ (τ ) = Rxx ∗ (τ ) ⊗ rhh∗ (τ ) .

(3.26)

d. By specializing (3.15) and (3.16) to the case x1 = x2 = x, h2 (t + τ, t) = δ(τ ), y1 = y, y2 = x, so that J2 = {0} and  α rσ1 σ2 (∗) (τ ) = hσ1 (s + τ ) δ(s) e−j 2παs ds = hσ1 (τ ) (3.27) R

we obtain the cyclic cross-statistics of the output and input signals

 α−σ α j 2πσ1 τ 1 ⊗ hσ1 (τ ) Rxx (∗) (τ ) e Ryx (∗) (τ ) = τ

(3.28)

σ1 ∈J1

α Syx (∗) (f ) =



σ1 ∈J1

α−σ1 Sxx (∗) (f − σ1 ) Hσ1 (f ) .

(3.29)

e. By specializing (3.16) to the case h1 and h2 LTI and (∗) absent, for α = 0 we obtain Sy01 y2 (f ) = Sx01 x2 (f ) H1 (f ) H2 (−f ) F ←→ Ry01 y2 (τ ) = Rx01 x2 (τ ) ⊗ h1 (τ ) ⊗ h2 (−τ )

(3.30) (3.31)

Example 3.8. Time Delay. Let us consider the LTI system that introduces a time delay d: y(t) = x(t − d) = x(t) ⊗ δ(t − d)

(3.32)

The (conjugate) cyclic statistics of the output signal can be expressed in terms of those of the input signal: α α −j 2παd Ryy (∗) (τ ) = Rxx (∗) (τ ) e

(3.33)

α α −j 2παd Syy (∗) (f ) = Sxx (∗) (f ) e

(3.34)

88

Cyclostationary Processes and Time Series

α α Ryx (∗) (τ ) = Rxx (∗) (τ − d)

(3.35)

α α −j 2πf d Syx (∗) (f ) = Sxx (∗) (f ) e

(3.36)

From (3.33) it follows that the (conjugate) cyclic autocorrelation of y(t) for α = 0 contains a phase term depending on d. In contrast, for α = 0 the time averaged (conjugate) 0 0 autocorrelation Ryy (∗) (τ ) is equal to Rxx (∗) (τ ) and is independent of the time shift d. Relationships (3.33)–(3.36) are exploited in (Gardner and Chen, 1988a), (Gardner and Chen, 1992), (Gardner and Chen, 1992) to design signal-selective time-difference-of arrival (TDOA) estimation algorithms. Example 3.9. Spectral Coherence and LTI Filtering. The magnitude of the spectral coherence (1.60) is invariant to LTI filtering, provided that the spectral components are not annihilated by the filter harmonic response. Let y(t) be given by (3.21). Accounting for (3.22), we have α γyy (∗) (f )  

=

α Syy (∗) (f )

1/2 0 (f ) S 0 ((−)(α − f )) Syy ∗ yy ∗ α (∗) ((−)(α − f )) Sxx (∗) (f ) H (f ) H 0 (f ) |H (f )|2 S 0 ((−)(α − f )) |H ((−)(α − f ))|2 Sxx ∗ xx ∗

α = γxx (∗) (f )

Thus,

H (f ) H (∗) ((−)(α − f )) |H (f )| |H ((−)(α − f ))|

1/2

(3.37)

      α   α γyy (∗) (f ) = γxx (∗) (f )

(3.38)

H (f ) H (∗) ((−)(α − f )) = 0 .

(3.39)

for all (α, f ) such that

3.3 Products of ACS signals Let x1 (t), x2 (t), c1 (t), and c2 (t) be jointly ACS complex-valued signals with (conjugate) cross-correlation functions  R αx (∗) (τ ) ej 2παx t (3.40) Rx x (∗) (t, τ ) = 1 2

αx ∈A

(∗) x1 x2

Rc

(∗) 1 c2

(t, τ ) =

 αc ∈A

(∗) c1 c2

x1 x2

R αc (∗) (τ ) ej 2παc t . c1 c2

(3.41)

If x1 (t) and x2 (t) are statistically independent of c1 (t) and c2 (t), then the eighth-order joint probability density function of their real and imaginary parts factors into the prod-

Chapter 3 • Almost-cyclostationary signal processing 89

uct of the two fourth-order joint probability densities of real and imaginary parts of x1 , x2 and of c1 , c2 , (Gardner, 1987f), (Brown, 1987), (Gardner, 1994a). Moreover, the (conjugate) cross-correlation function of the product waveforms yi (t) = ci (t) xi (t)

i = 1, 2

(3.42)

also factors: Ry

(∗) 1 y2

(t, τ ) = Rc

(∗) 1 c2

(t, τ ) Rx

(∗) 1 x2

(t, τ ) .

(3.43)

Theorem 3.10. Under Assumption 3.4 (and a similar assumption on the (conjugate) cyclic cross-correlations of c1 and c2 ) the (conjugate) cyclic cross-correlation function and the (conjugate) cyclic cross-spectrum of y1 (t) and y2 (t) are given by  R α (∗) (τ ) = R αx (∗) (τ ) R α−α(∗)x (τ ) (3.44) y1 y2

Sα

x1 x2

αx ∈A

(∗) x1 x2

(∗) (f ) =

y1 y2





c1 c2

S αx

(∗)

R x1 x2

αx ∈A

(∗) x1 x2

x (λ) S α−α (∗) (f − λ)dλ

c1 c2

where, in the sums, only those (conjugate) cycle frequencies αx such that α − αx ∈ Ac nonzero contribution.

(3.45)

(∗) 1 c2

give

Proof. See Section 1.9. From (3.44) and (3.45), it follows that the set of (conjugate) cycle frequencies of the (conjugate) cross-correlation of y1 (t) and y2 (t) is

 {α} = α = αc + αx , αc ∈ Ac c(∗) , αx ∈ Ax x (∗) . (3.46) 1 2

1 2

In the special case where c1 (t) and c2 (t) are almost-periodic functions with (generalized) Fourier series  ci,γi ej 2πγi t i = 1, 2 (3.47) ci (t) = γi ∈Gi

then Rc

(∗) 1 c2

(∗)

(t, τ ) = c1 (t + τ ) c2 (t)   (∗) = c1,γ1 c2,γ2 ej 2πγ1 τ ej 2π(γ1 +(−)γ2 )t

(3.48)

γ1 ∈G1 γ2 ∈G2

and Rα

(∗) c1 c 2

(τ ) =

 γ1 ∈G1

(∗)

c1,γ1 c2,(−)(α−γ1 ) ej 2πγ1 τ

(3.49)

90

Cyclostationary Processes and Time Series

Sα

(∗) c1 c 2

(f ) =



(∗)

γ1 ∈G1

c1,γ1 c2,(−)(α−γ1 ) δ(f − γ1 )

(3.50)

where in the sum give nonzero contribution only those frequencies γ1 such that (−)(α − γ1 ) ∈ G2 . From (3.49) and (3.50) it follows that the set of the (conjugate) cycle-frequencies of the (conjugate) cross-correlation of c1 (t) and c2 (t) is {α} = {α = γ1 + (−)γ2 , γ1 ∈ G1 , γ2 ∈ G2 } .

(3.51)

In the special case x1 ≡ x2 and c1 ≡ c2 (3.48)–(3.50) reduce to (corrected versions of) (3.93)–(3.95) in (Gardner et al., 2006). By substituting (3.49) and (3.50) into (3.44) and (3.45), respectively, and making the variable change γ2 = (−)(α − αx − γ1 ), we get (see (Napolitano, 1995, Eqs. (46)–(48)))   α−(γ +(−)γ2 ) (∗) c1,γ1 c2,γ2 R (∗)1 (τ ) ej 2πγ1 τ (3.52) R α (∗) (τ ) = y1 y2

Sα

(∗)

y1 y2

x1 x2

γ1 ∈G1 γ2 ∈G2

(f ) =

 

γ1 ∈G1 γ2 ∈G2

(∗)

c1,γ1 c2,γ2 S

α−(γ1 +(−)γ2 ) (f (∗) x1 x2

− γ1 ) ,

(3.53)

where, in the sums, give nonzero contribution only those values of γ1 and γ2 such that α − (γ1 + (−)γ2 ) ∈ Ax x (∗) . From (3.52) and (3.53), it follows that the set of the (conjugate) 1 2

cycle-frequencies of the (conjugate) cross-correlation of y1 (t) and y2 (t) is

 {α} = α = γ1 + (−)γ2 + αx , αx ∈ Ax x (∗) , γ1 ∈ G1 , γ2 ∈ G2 . 1 2

(3.54)

Example 3.11. Modulation. Let us consider the system that performs the modulation of a complex carrier: y(t) = x(t) ej 2πf0 t .

(3.55)

The (conjugate) cyclic statistics of the output signal can be expressed in terms of those of the input signal: α−κf

α j 2πf0 τ 0 Ryy (∗) (τ ) = Rxx (∗) (τ ) e α Syy (∗) (f ) = α Ryx (∗) (τ ) = α Syx (∗) (f ) =

α−κf Sxx (∗) 0 (f − f0 ) α−f Rxx (∗)0 (τ ) ej 2πf0 τ α−f Sxx (∗)0 (f − f0 )

(3.56) (3.57) (3.58) (3.59)

where κ  1 + (−)1, that is, κ = 0 if (∗) = ∗ and κ = 2 if (∗) is absent.

3.4 Supports of cyclic spectra of band limited signals Definition 3.12. A signal x(t) is said to be a strictly band-limited low-pass signal with monolateral bandwidth B if and only if x(t) passes unaltered through the ideal low-pass

Chapter 3 • Almost-cyclostationary signal processing 91

filter with monolateral bandwidth B. That is, x(t) ≡ x(t) ⊗ hLPF (t)

(3.60)

where hLPF (t) is the impulse-response function of the ideal low-pass filter with harmonicresponse function    1 |f |  B f HLPF (f ) = rect  (3.61) 2B 0 |f | > B .

Fact 3.13. Accounting for (3.25), for a strictly band-limited low-pass signal we have 0 0 2 Sxx ∗ (f ) = Sxx ∗ (f ) |H (f )|

(3.62)

0 (f ) ≡ 0 for f ∈ (−B, B), which is commonly used as definition of strictly bandthat is, Sxx ∗ limited low-pass process in the stationary case.

From (3.22) it follows that the support in the (α, f ) plane of the (conjugate) cyclic spectrum of y(t) is such that 

  α α (f )  (α, f ) ∈ R × R : S (f ) = 0 supp Syy (∗) yy (∗) 

(3.63) ⊆ (α, f ) ∈ R × R : H (f ) H (∗) ((−)(α − f )) = 0 . For the LTI filtering (3.60), accounting for the fact that rect(f ) is real and even, the input/output relation (3.22) specializes to α α Sxx (∗) (f ) = HLPF (f ) HLPF (f − α) Sxx (∗) (f ) .

(3.64)

Hence, (3.63) specializes into   α supp Sxx (∗) (f ) ⊆ {(α, f ) ∈ R × R : HLPF (f ) HLPF (f − α) = 0} ≡ {(α, f ) ∈ R × R : |f |  B, |α − f |  B}

(3.65a)

⊆ {(α, f ) ∈ R × R : |f |  B, |α|  2B} .

(3.65b)

In Fig. 3.3, set (3.65a) is represented. By common convention, f and α variables are on the abscissa and ordinate axis, respectively of the (α, f )-plane. From the conditions characterizing the bigger set (3.65b) it results that α Sxx (∗) (f ) = 0

|f | > B, |α| > 2B .

(3.66)

Let x(t) be a strictly band-limited high-pass signal, that is, x(t) ≡ x(t) ⊗ hHPF (t)

(3.67)

92

Cyclostationary Processes and Time Series

FIGURE 3.3 Cyclic-spectrum support of a low-pass signal. Source: (Gardner et al., 2006) Copyright of Elsevier.

FIGURE 3.4 Cyclic-spectrum support of a high-pass signal. Source: (Gardner et al., 2006) Copyright of Elsevier.

where hHPF (t) is the impulse-response function of the ideal high-pass filter with harmonicresponse function   f . (3.68) HHPF (f ) = 1 − rect 2b 0 (f ) ≡ 0 for f ∈ (−b, b). In such a case, Sxx ∗ From (3.63), we have (see Fig. 3.4)

  α supp Sxx (∗) (f ) ⊆ {(α, f ) ∈ R × R : HHPF (f ) HHPF (f − α) = 0} .

(3.69)

Chapter 3 • Almost-cyclostationary signal processing 93

FIGURE 3.5 Cyclic-spectrum support of a band-pass signal. Source: (Gardner et al., 2006) Copyright of Elsevier.

Finally, let x(t) be a strictly band-limited band-pass signal, that is, x(t) ≡ x(t) ⊗ hBPF (t)

(3.70)

where hBPF (t) is the impulse-response function of the ideal band-pass filter with harmonicresponse function HBPF (f ) = HLPF (f ) HHPF (f )

(3.71)

where HLPF (f ) and HHPF (f ) are given by (3.61) and (3.68), respectively. In such a case, 0 (f ) ≡ 0 for f ∈ (−B, −b) ∪ (b, B), where 0 < b < B. Sxx ∗ Accounting for (3.63), we have (see Fig. 3.5)   α supp Sxx (∗) (f ) ⊆ {(α, f ) ∈ R × R : HBPF (f ) HBPF (f − α) = 0} = {(α, f ) ∈ R × R : HLPF (f ) HLPF (f − α) = 0} ∩ {(α, f ) ∈ R × R : HHPF (f ) HHPF (f − α) = 0} .

(3.72)

It is noted that supports for symmetric definitions of cyclic spectra (1.35) are “diamond shaped” and are given in (Gardner, 1986c), (Gardner, 1987f, Chap. 11), (Gardner, 1991a).

3.5 Rice’s representation In this section, cyclic statistical functions of real signal, analytic signal, complex envelope, and in-phase and in-quadrature components are linked among them (Gardner, 1987c), (Izzo and Napolitano, 1997). For related results see also (Rybakov, 1981), (Rybakov, 1986), (Gurevich and Korshunov, 1987). Expressions specialized for α = 0 reveal that cyclosta-

94

Cyclostationary Processes and Time Series

tionarity cannot be neglected even if only the time-averaged autocorrelation or the power spectrum is computed. Let x(t) be a real-valued signal. Assumption 3.14. xdc  E{x(t)}t = 0, that is, the power spectrum Sxx (f ) does not contain the additive term xdc δ(f ). Assumption 3.15. The signal x(t) does not contain any finite-strength additive sine-wave α (f ) do not contain Dirac deltas (Section 1.4). component. That is, cyclic spectra Sxx Assumption 3.15 avoids in the following expressions the presence of undetermined products of sign(·), u(·), and rect(·) functions with Dirac deltas centered in the discontinuity points of these functions. In the special case α = 0, only Assumption 3.14 is needed.

Hilbert transform The Hilbert transform of a signal x(t) is defined by the LTI transformation (Section D.1)  x (t)  x(t) ⊗ hH (t)

(3.73)

where 1 F ←→ HH (f ) = −j sign(f ) . πt Under Assumption 3.15, and using (3.22) and (3.23), we have hH (t) 

(3.74)

   α α S x x (f ) = Sxx (f ) −j sign(f ) −j sign(α − f ) α = Sxx (f ) sign(f ) sign(f − α)   f − α/2 α (f ) 1 − 2 rect = Sxx |α|

1 1 F α α j 2πατ ←→ R (τ ) = R (τ ) ⊗ ⊗ − e xx x x πτ πτ   α = Rxx (τ ) ⊗ δ(τ ) − 2|α| sinc(|α|τ ) ej πατ

  α α S x x (f ) = Sxx (f ) −j sign(f ) 1 F α α α xx (τ ) R ←→ R x x (τ ) = Rxx (τ ) ⊗ πτ   α α Sx x (f ) = Sxx (f ) −j sign(α − f )

1 j 2πατ F α α (τ ) = R (τ ) ⊗ − e ←→ Rx xx x πτ

(3.75)

(3.76) (3.77) (3.78) (3.79) (3.80)

In the special case α = 0, under Assumption 3.14 it results    S x x (f ) = Sxx (f ) −j sign(f ) −j sign(−f ) = Sxx (f )

(3.81)

Chapter 3 • Almost-cyclostationary signal processing 95

F ←→ R x x (τ ) = Rxx (τ )

(3.82)

  S x x (f ) = Sxx (f ) −j sign(f )

(3.83)

1 F xx (τ ) R ←→ R x x (τ ) = Rxx (τ ) ⊗ πτ   Sx x (f ) = Sxx (f ) −j sign(−f ) = j sign(f ) Sxx (f ) 1 F xx (τ ) ←→ Rx  −R x (τ ) = − Rxx (τ ) ⊗ πτ

(3.84) (3.85) (3.86)

Analytic signal Let us consider the LTI system with impulse response function

1 2 1 F hA (t)  δ(t) = δ(t) + j ←→ HA (f ) = u(f ) k πt k ◦

(3.87)

√ where k is a normalizing factor which is chosen differently from different authors (k = 2 in (Van Trees, 1971b), k = 1 in (Papoulis, 1983), (Gardner, 1987c), (Gardner, 1987f, Chap. 11, Sec. D-2), (Papoulis, 1991), and k = 2 in (Izzo and Napolitano, 1997)). The analytic signal of a signal x(t) is defined by the LTI transformation (Section D.1) ◦

x(t)  x(t) ⊗ hA (t) =

 1 x(t) + j x (t) . k ◦

(3.88) ◦

Under Assumption 3.15, and using (3.22) and (3.23) and δ ∗ (−t) = δ(t), we have     ◦ ◦ R α◦ ◦ ∗ (τ )  E x(t + τ ) x ∗ (t) e−j 2παt xx t ◦  ◦ α j 2πατ = Rxx (τ ) ⊗ δ(τ ) ⊗ δ(τ ) e 4 F α (f ) 2 u(f ) u(f − α) ←→ S α◦ ◦ ∗ (f ) = Sxx xx k  u(f − α) α0 4 α (f ) 2 = Sxx k u(f ) α<0  ◦ j 2πατ δ(τ ) e α0 2 α F (τ ) ⊗ ←→ R α◦ ◦ ∗ (τ ) = Rxx ◦ xx k δ(τ ) α<0     ◦ ◦ R α◦ ◦ (τ )  E x(t + τ ) x(t) e−j 2παt xx t ◦  ◦ α = Rxx (τ ) ⊗ δ(τ ) ⊗ δ(−τ ) ej 2πατ 4 F α (f ) 2 u(f ) u(α − f ) ←→ S α◦ ◦ (f ) = Sxx xx k

(3.89) (3.90) (3.91) (3.92)

(3.93) (3.94)

96

Cyclostationary Processes and Time Series ⎧   ⎨ rect f −α/2 α0 4 α α (f ) 2 = Sxx k ⎩ 0 α<0  α sinc(ατ ) ej πατ α0 4 α F ←→ R α◦ ◦ (τ ) = 2 Rxx (τ ) ⊗ xx k 0 α<0

(3.95)

(3.96)

Using the symmetry relationships (Section 1.5.1) and the fact that x(t) is real, we have S α◦ ∗ ◦ (f ) = S α◦ ◦ ∗ (α − f ) x

x

xx

4 u(α − f ) u(−f ) k2

α (α − f ) = Sxx α (f ) = Sxx

S α◦ ∗ x

◦∗

x

4 u(α − f ) u(−f ) k2

(3.97)

∗ (f ) = S −α ◦ ◦ (−f ) xx

−α = Sxx (−f )∗

4 u(−f ) u(−α + f ) k2

4 u(−f ) u(−α + f ) k2

(3.98)

◦  k ◦  ◦ x(t) = k Re x(t) = x(t) + x ∗ (t) . 2

(3.99)

α (f ) = Sxx

From (3.88) we have

Thus, k2 4 2 k F α ←→ Sxx (f ) = 4 α Rxx (τ ) =

 R α◦ ◦ ∗ (τ ) + R α◦ ∗ ◦ (τ ) + R α◦ ◦ (τ ) + R α◦ ∗ xx

x

x

xx

x

 S α◦ ◦ ∗ (f ) + S α◦ ∗ ◦ (f ) + S α◦ ◦ (f ) + S α◦ ∗ xx

x

x

xx

x

 (τ )

(3.100)

 (f )

(3.101)

◦ x∗ ◦ x∗

In (3.101), the cyclic spectra in the sum have non-overlapping supports in the (α, f ) plane (Gardner, 1986c), (Gardner, 1987f, Chap. 11, Sec. D-2), (Izzo and Napolitano, 1997). In fact, S α◦ ◦ ∗ (f ) is nonzero for f  0 and f  α (see (3.90)), S α◦ ◦ (f ) is nonzero for f  0 xx xx and f  α (see (3.94)), S α◦ ∗ ◦ (f ) is nonzero for f  0 and f  α (see (3.97)), and S α◦ ∗ ◦ ∗ (f ) is x x x x nonzero for f  0 and f  α (see (3.98)). In the special case α = 0, under Assumption 3.14 it results  ◦  ◦ Rx◦ x◦ ∗ (τ )  E x(t + τ ) x ∗ (t) t 4 4 F ←→ Sx◦ x◦ ∗ (f ) = Sxx (f ) 2 u(f ) u(−(−f )) = Sxx (f ) 2 u(f ) k k

2 1 1 2 ◦ F ←→ Rx◦ x◦ ∗ (τ ) = Rxx (τ ) ⊗ δ(t) + j = Rxx (τ ) k k πt k

(3.102) (3.103) (3.104)

Chapter 3 • Almost-cyclostationary signal processing 97  ◦  ◦ Rx◦ x◦ (τ )  E x(t + τ ) x(t) t

(3.105)

4 F ←→ Sx◦ x◦ (f ) = Sxx (f ) 2 u(f ) u(−f ) = 0 k F ←→ Rx◦ x◦ (τ ) = 0 Sx◦ ∗ x◦ (f ) = Sxx (f )

(3.106) (3.107)

4 u(−f ) k2

(3.108)

Sx◦ ∗ x◦ ∗ (f ) = 0

(3.109)

Complex envelope The complex envelope of a signal x(t) is defined as (Section D.1) ◦

 x (t)  x(t) e−j 2πf0 t .

(3.110)

Under Assumption 3.15 we have    −j 2παt  α ∗ R (τ )  E  x (t + τ )  x (t) e ∗ x x   t  ◦ ◦ = E x(t + τ ) x ∗ (t) e−j 2παt e−j 2πf0 τ t

= R α◦ ◦ xx

∗

(τ ) e

−j 2πf0 τ

(3.111)

F α α ←→ S x x ∗ (f ) = S ◦ ◦ ∗ (f + f0 ) xx

4 u(f + f0 ) u(f + f0 − α) k2  u(f + f0 − α) α0 4 α (f + f0 ) 2 = Sxx k α<0 u(f + f0 )   α x (t + τ )  x (t)} e−j 2παt R x x (τ )  E { t    −j 2π[α+2f ]t  −j 2πf τ ◦ ◦ 0 0 = E x(t + τ ) x(t) e e α (f + f0 ) = Sxx

(3.112) (3.113)

t

α+2f = R ◦ ◦ 0 (τ ) e−j 2πf0 τ xx

(3.114)

F α+2f0 α ←→ S (f + f0 ) x x (f ) = S ◦ ◦ xx

α+2f0

(f + f0 )

α+2f0

(f + f0 )

= Sxx = Sxx

4 u(f + f0 ) u(α + 2f0 − (f + f0 )) k2 ⎧   f −α/2 α  −2f 4 ⎨ rect k2 ⎩ 0

α+2f0

0

α < −2f0

(3.115) (3.116)

Using the symmetry relationships (Section 1.5.1) and the fact that x(t) is real, we have α α S x ∗ x x ∗ (α − f ) x (f ) = S

98

Cyclostationary Processes and Time Series

α = Sxx (α − f + f0 ) α (f − f0 ) = Sxx

4 u(α − f + f0 ) u(α − f + f0 − α) k2

4 u(α − f + f0 ) u(−f + f0 ) k2

(3.117)

−α α ∗ S x ∗ x ∗ (f ) = S x x (−f ) −α+2f0

= Sxx

α−2f0

= Sxx

(−f + f0 )∗

(f − f0 )

4 u(−f + f0 ) u(−α + f0 + f ) k2

4 u(−f + f0 ) u(−α + f0 + f ) k2

(3.118)

From (3.88) and (3.110) we have   k   x(t) = k Re  x (t) ej 2πf0 t = x ∗ (t) e−j 2πf0 t .  x (t) ej 2πf0 t +  2

(3.119)

Thus (see Example 3.11), α (τ ) = Rxx

k2  α j 2πf0 τ α −j 2πf0 τ + R R x (τ ) e x x ∗ (τ ) e x∗  4 α−2f0

+ R x x

−j 2πf0 τ 0 (τ ) ej 2πf0 τ + R x∗  x ∗ (τ ) e α+2f

 (3.120)

k2  α F α α (f ) = Sx  ←→ Sxx x (f + f0 ) x ∗ (f − f0 ) + S x∗  4  α−2f0

+ S x x

 α+2f0 (f − f0 ) + S x∗  x ∗ (f + f0 )

(3.121)

In the special case α = 0, under Assumption 3.15 it results    x (t + τ ) x ∗ (t) t R x x ∗ (τ )  E   ◦  ◦ = E x(t + τ ) x ∗ (t) t e−j 2πf0 τ = Rx◦ x◦ ∗ (τ ) e−j 2πf0 τ

(3.122)

F ←→ S x x ∗ (f ) = Sx◦ x◦ ∗ (f + f0 ) 4 = Sxx (f + f0 ) 2 u(f + f0 ) k

(3.123)

x (t + τ ) x (t)}t R x x (τ )  E {     ◦ ◦ = E x(t + τ ) x(t) e−j 2π2f0 t e−j 2πf0 τ t

2f = R ◦ ◦0 (τ ) e−j 2πf0 τ xx

F 2f0 ←→ S x x (f ) = S ◦ ◦ (f + f0 ) xx 2f

= Sxx 0 (f + f0 )

4 u(f + f0 ) u(f0 − f ) k2

(3.124)

Chapter 3 • Almost-cyclostationary signal processing 99

2f

= Sxx 0 (f + f0 )

  4 f rect 2f0 k2

4 u(−f + f0 ) k2   4 f −2f0 S x ∗ x ∗ (f ) = Sxx (f − f0 ) 2 rect 2f0 k S x ∗ x (f ) = Sxx (f − f0 )

(3.125) (3.126) (3.127)

Supports of cyclic spectra of band pass signals Definition 3.16. The signal x(t) is bandpass with central frequency f0 and strictly bandlimited with monolateral bandwidth W , with f0 > W , if and only if it passes unaltered through the ideal bandpass filter with harmonic response  f − f  f + f  1 ||f | − f0 |  W 0 0 HBPF (f )  rect + rect = (3.128) 2W 2W 0 otherwise. That is, x(t) = x(t) ⊗ hBPF (t)

(3.129)

where the impulse-response function hBPF (t) is the inverse Fourier transform of HBPF (f ). If x(t) is ACS, then accounting for (3.128) and (3.129), and the input/output relation in terms of cyclic spectra for linear time-invariant (LTI) filters (3.22), the support in the (α, f )-plane of its cyclic spectra is given by  α  supp Sxx (f ) ⊆ {(α, f ) ∈ R × R : HBPF (f ) HBPF (α − f ) = 0} = {(α, f ) ∈ R × R : ||f | − f0 |  W, ||α − f | − f0 |  W } .

(3.130)

Remember that f and α variables are on the abscissa and ordinate axis, respectively of the (α, f )-plane. From (3.112) and (3.130) it follows that   α   α supp S x x ∗ (f ) = supp Sxx (f + f0 ) ∩ {(α, f ) ∈ R × R : f  −f0 , α  f + f0 } ⊆ D1 (α, f )

(3.131)

where D1 (α, f ) is the domain defined as (Fig. 3.3, Eq. (3.65a) with B = W ) D1 (α, f )  {(α, f ) ∈ R × R : |f |  W, |α − f |  W }

(3.132)

which is obtained when the asymmetric definition for the cyclic spectrum (1.30) is used. Note that a diamond-shaped domain is obtained for D1 (α, f ) (Gardner, 1991a, Fig. 2b) if the symmetric definition for the cyclic spectrum (1.35) (Gardner, 1991a, Eqs. (44)–(48)) is adopted.

100

Cyclostationary Processes and Time Series

FIGURE 3.6 Set D2 (α, f ) defined in (3.137) containing the support of the cyclic spectrum of a bandpass ACS signal. ˘ Source: (Napolitano and Dogançay, 2018) Copyright of Elsevier.

Analogously, it can be shown that  α  supp S x (f ) ⊆ D1 (α, f ) x ∗  α  supp S x x (f ) ⊆ D1 (α, f )  α  supp S x ∗ (f ) ⊆ D1 (α, f ) . x ∗

(3.133) (3.134) (3.135)

Therefore, from (3.121) and (3.131)–(3.135) we have (Fig. 3.6)   α (f ) ⊆ D2 (α, f ) supp Sxx

(3.136)

where D2 (α, f )  D1 (α, f − f0 ) ∪ D1 (α, f + f0 ) ∪ D1 (α − 2f0 , f − f0 ) ∪ D1 (α + 2f0 , f + f0 ) .

(3.137)

In Fig. 3.6, the set D2 (α, f ) defined in (3.137) is represented with axes α and f normalized to a reference frequency fr . It is coincident with the set represented in Fig. 3.5 (Eq. (3.72) with b = f0 − W and B = f0 + W ). In (3.121), the cyclic spectrum of x(t) is expressed as the superposition of the contributions of cyclic spectra of the complex envelope  x (t) with different conjugation choices. In Fig. 3.6 the location of such contributions in the (α, f ) plane is evidenced.

Chapter 3 • Almost-cyclostationary signal processing

101

In-phase and in-quadrature components The in-phase and in-quadrature components are defined as (Section D.1)  k   x (t) +  x ∗ (t) 2  k  x (t)] = −j  x (t) −  x ∗ (t) xs (t)  k Im [ 2

xc (t)  k Re [ x (t)] =

(3.138) (3.139)

respectively. Under Assumption 3.15 we have   Rxαc xc (τ )  E {xc (t + τ ) xc (t)} e−j 2παt

t

= = F ←→ Sxαc xc (f ) = Rxαs xs (τ )  = = F ←→ Sxαs xs (f ) = Rxαc xs (τ ) 

   k 2   x (t) +  x ∗ (t) e−j 2παt E  x (t + τ ) +  x ∗ (t + τ )  t 4 2   k α α α α R x x ∗ (τ ) + R x ∗ x x (τ ) + R x ∗ x (τ ) + R x ∗ (τ ) 4  k2  α α α α S x (f ) + S x ∗ (f ) x x ∗ (f ) + S x ∗ x x (f ) + S x ∗ 4   E {xs (t + τ ) xs (t)} e−j 2παt

(3.140) (3.141)

t

   k 2   − x (t) −  x ∗ (t) e−j 2παt E  x (t + τ ) −  x ∗ (t + τ )  t 4  k2  α α α α R x (τ ) − R x ∗ (τ ) x x ∗ (τ ) + R x ∗ x x (τ ) − R x ∗ 4  k2  α α α α S x x ∗ (f ) + S x ∗ x x (f ) − S x ∗ x (f ) − S x ∗ (f ) 4   E {xc (t + τ ) xs (t)} e−j 2παt

(3.142) (3.143)

t

   k 2   = −j x (t) −  x ∗ (t) e−j 2παt E  x (t + τ ) +  x ∗ (t + τ )  t 4  k2  α α α α =j R x x ∗ (τ ) − R x ∗ x x (τ ) + R x ∗ x (τ ) − R x ∗ (τ ) 4  k2  α F α α α ←→ Sxαc xs (f ) = j S x (f ) − S x ∗ (f ) x x ∗ (f ) − S x ∗ x x (f ) + S x ∗ 4   Rxαs xc (τ )  E {xs (t + τ ) xc (t)} e−j 2παt

(3.144) (3.145)

t

= −j

k2

4 k2 = −j 4 2 k F ←→ Sxαs xc (f ) = − j 4

     x (t) +  x ∗ (t) e−j 2παt E  x (t + τ ) −  x ∗ (t + τ ) 

t

 α α α α R x (τ ) + R x ∗ (τ ) x x ∗ (τ ) − R x ∗ x x (τ ) − R x ∗

(3.146)

 α α α α S x x ∗ (f ) − S x ∗ x x (f ) − S x ∗ x (f ) + S x ∗ (f )

(3.147)

 

102

Cyclostationary Processes and Time Series  1  −j 2παt {(x (t + τ ) + j x (t + τ )) (x (t) + (−)j xs (t))} e E c s c t k2  1 α α α α = 2 Rxc xc (τ ) + j Rxs xc (τ ) + (−)j Rxc xs (τ ) − (−)Rxs xs (τ ) k  1 α F α α α α ←→ S (f ) = (f ) + j S (f ) + (−)j Sx x (f ) − (−)Sx x (f ) S (∗) x x x x c c s c c s s s x x k2 α R (τ ) = x x (∗)

(3.148) (3.149)

In the special case α = 0, under Assumption 3.15 and accounting for (D.17) and (D.18), it results Sxc xc (f ) = Sxx (f + f0 ) u(f + f0 ) + Sxx (f − f0 ) u(−f + f0 )  f   f  2f −2f + Sxx 0 (f − f0 ) rect + Sxx 0 (f + f0 ) rect 2f0 2f0     F 2f0 −2f ←→ Rxc xc (τ ) = Rxx (τ ) c + Rxx (τ ) e−j 2πf0 τ + Rxx 0 (τ ) ej 2πf0 τ ⊗ 2f0 sinc(2f0 τ ) Sxs xs (f ) = Sxx (f + f0 ) u(f + f0 ) + Sxx (f − f0 ) u(−f + f0 )  f   f  2f −2f − Sxx 0 (f − f0 ) rect − Sxx 0 (f + f0 ) rect 2f0 2f0     F 2f −2f ←→ Rxs xs (τ ) = Rxx (τ ) c − Rxx0 (τ ) e−j 2πf0 τ + Rxx 0 (τ ) ej 2πf0 τ ⊗ 2f0 sinc(2f0 τ )  Sxc xs (f ) = j Sxx (f + f0 ) u(f + f0 ) − Sxx (f − f0 ) u(−f + f0 )  f   f  2f −2f − Sxx 0 (f + f0 ) rect + Sxx 0 (f − f0 ) rect 2f0 2f0     F 2f −2f ←→ Rxc xs (τ ) = − Rxx (τ ) s − j Rxx0 (τ ) e−j 2πf0 τ − Rxx 0 (τ ) ej 2πf0 τ ⊗ 2f0 sinc(2f0 τ )  Sxs xc (f ) = − j Sxx (f + f0 ) u(f + f0 ) − Sxx (f − f0 ) u(−f + f0 )  f   f  2f −2f + Sxx 0 (f + f0 ) rect − Sxx 0 (f − f0 ) rect 2f0 2f0     F 2f −2f ←→ Rxs xc (τ ) = Rxx (τ ) s − j Rxx0 (τ ) e−j 2πf0 τ − Rxx 0 (τ ) ej 2πf0 τ ⊗ 2f0 sinc(2f0 τ )

(3.150)

(3.151)

(3.152)

(3.153)

(3.154)

(3.155)

(3.156)

(3.157)

where [·]c and [·]s denote the operators that provide the in-phase and in-quadrature components of their argument according to (D.15) and (D.16), respectively. From (3.150)–(3.156) it follows that the power spectra and cross-power spectra (α = 0) of in-phase and in-quadrature components of a cyclostationary signal depend not only

Chapter 3 • Almost-cyclostationary signal processing

103

on the power spectrum of the real signal, but also from cyclic spectra at cycle frequencies α = ±2f0 of the real signal. Thus, neglecting these cyclic spectra leads to wrong results.

Wide-sense stationary signals Let x(t) be a zero mean wide-sense stationary signal. Under the zero-mean assumption we have   ◦        E x(t) = 0 ⇒ E x(t) = E  x (t) = E xc (t) = E xc (t) = 0 . (3.158) From the wide-sense stationarity assumption it follows that ±2f0

Sxx

F ±2f (f ) ≡ 0 ←→ Rxx 0 (τ ) ≡ 0 .

(3.159)

In addition, under Assumption 3.14, it results that F S x x (f ) = S x ∗ x ∗ (f ) ≡ 0 ←→ R x x (τ ) = R x ∗ x ∗ (τ ) ≡ 0

(3.160)

  Rxc xc (τ ) = Rxs xs (τ ) = Rxx (τ ) c

xx (τ ) sin(2πf0 τ ) = Rxx (τ ) cos(2πf0 τ ) + R

(3.161)

F ←→ Sxc xc (f ) = Sxs xs (f ) = Sxx (f + f0 ) u(f + f0 ) + Sxx (f − f0 ) u(−f + f0 )  Rxs xc (τ ) = −Rxc xs (τ ) = Rxx (τ ) s xx (τ ) cos(2πf0 τ ) = − Rxx (τ ) sin(2πf0 τ ) + R

(3.162)



(3.163)



F ←→ Sxs xc (f ) = −Sxc xs (f ) = − j Sxx (f + f0 ) u(f + f0 )

 − Sxx (f − f0 ) u(−f + f0 )

(3.164)

xx (τ ) xx (−τ ) = R R

(3.165)

Rxs xc (τ ) = Rxc xs (−τ )     Rxx (τ ) = Rxx (τ ) c cos(2πf0 τ ) − Rxx (τ ) s sin(2πf0 τ )

(3.166) (3.167)

Stationary white noise In this section, the complex envelope signal w (t) =

  1 1 w(t) + j w (t) e−j 2πf0 t = wc (t) + j ws (t) k k

(3.168)

is analyzed, where w(t) is wide-sense stationary white noise (Kailath, 1966). For related results see also (Viswanathan, 2006). It results that Rww (τ ) = E{w(t + τ ) w(t)} =

N0 δ(τ ) 2

(3.169)

104

Cyclostationary Processes and Time Series

N0 F ←→ Sww (f ) = 2

(3.170)

Rw w  (τ ) = Rww (τ ) =

N0 δ(τ ) 2

(3.171)

N0 F ←→ Sw w  (f ) = 2

(3.172)

N0 1 N0 1 δ(τ ) ⊗ = 2 πτ 2 πτ N0 F ←→ Sw w (f ) = − Sw w (f ) = −j sign(f ) 2   ww (τ ) sin(2πf0 τ ) Rwc wc (τ ) = Rww (τ ) c = Rww (τ ) cos(2πf0 τ ) + R N0 1 N0 δ(τ ) cos(2πf0 τ ) + sin(2πf0 τ ) = 2 2 πτ N0 N0 1 = δ(τ ) + sin(2πf0 τ ) 2 2 πτ   Rws wc (τ ) = − Rwc ws (τ ) = Rww (τ ) s ww (τ ) cos(2πf0 τ ) = − Rww (τ ) sin(2πf0 τ ) + R ww (τ ) = Rw w (τ ) = − Rw w (τ ) = R

N0 1 N0 δ(τ ) sin(2πf0 τ ) + cos(2πf0 τ ) = − 2 2 πτ N0 1 = cos(2πf0 τ ) 2 πτ    2  Rw w(t + τ ) w ∗ (t)} = 2 Rww (τ ) c + j Rww (τ ) s w ∗ (τ ) = E{ k   N0 1 N0 1 2 N0 δ(τ ) + sin(2πf0 τ ) + j cos(2πf0 τ ) = 2 2 2 πτ 2 πτ k     2 N0 2 1 = 2 δ(τ ) + j e−j 2πf0 τ  Rww (τ ) lp πτ k k 2   2 N0 F 1 + j (−j sign(f )) ⊗ δ(f + f0 ) ←→ Sw w ∗ (f ) = 2 k 2 4 N0 u(f + f0 ) = 2 k 2

(3.173) (3.174) (3.175a)

(3.175b)

(3.176a)

(3.176b)

(3.177)

(3.178)

where [·]lp denotes the operator that provides the complex envelope of its argument according to (D.11).   Sws wc (f ) = − j Sww (f + f0 ) u(f + f0 ) − Sww (f − f0 ) u(−f + f0 )  N0  u(f + f0 ) − u(−f + f0 ) = −j 2  +∞ Sws wc (f ) df = 0 . lim Rws wc (τ ) = p.v. τ →0

−∞

(3.179) (3.180)

Chapter 3 • Almost-cyclostationary signal processing

105

From the wide-sense stationarity of w(t) it follows that E{ w(t + τ ) w (t)} = 0 .

(3.181)

3.6 Sampling and aliasing Let x(n)  xa (t)|t=nTs

n∈Z

(3.182)

be the discrete-time process obtained by uniformly sampling with period Ts = 1/fs the continuous-time ACS signal xa (t) with (conjugate) autocorrelation 

 Rxαa (τ ) ej 2παt (3.183) E xa (t + τ ) xa(∗) (t) = α∈Ax a

and Loève bifrequency spectrum (see (1.56))

  E Xa (f1 ) Xa(∗) (f2 ) = Sxαa (f1 ) δ(f2 − (−)(α − f1 )) ,

(3.184)

α∈Ax a (∗)

where subscript x a ≡ xa xa . The relationship between the Fourier transform Xa (f ) of the continuous-time signal xa (t) and the Fourier transform X(ν) of the discrete-time sampled signal x(n) is (Section C.1) (Lathi, 2002, Sec. 9.5) X(ν) =

+∞ 1  Xa ((ν − k)fs ) . Ts

(3.185)

k=−∞

If the signal xa (t) is strictly band limited, that is Xa (f ) = 0 for |f | > B, the following two conditions are equivalent, provided that fs  2B: 1. Replicas in (3.185) do not overlap; 2. In the principal frequency domain ν ∈ [−1/2, 1/2), the Fourier transform X(ν) of the discrete-time sampled signal is an amplitude and frequency scaled version of only one replica (the one with k = 0) of the Fourier transform Xa (f ) of the continuous time signal. That is 1 Xa (f )|f =νfs ν ∈ [−1/2, 1/2) . (3.186) X(ν) = Ts In the following, it is shown that in aliasing formulas for cyclic spectra, conditions for not overlap of replicas and conditions for the presence of a unique replica in the principal cycle- and spectral-frequency domain ( α , ν) ∈ [−1/2, 1/2) × [−1/2, 1/2) are not equivalent (see the proofs of Theorem 3.22 and 3.25 in Section 3.10).

3.6.1 Cyclic statistics of the sampled signal Lemma 3.17. Aliasing Formula for the Loève Bifrequency Spectrum of a Sampled ACS Signal. Let xa (t) be a continuous-time ACS signal with Loève bifrequency spectrum (3.184). The

106

Cyclostationary Processes and Time Series

Loève bifrequency spectrum of the sampled signal x(n) (3.182) is given by

 1   α E X(ν1 ) X (∗) (ν2 ) = Sx a ((ν1 − k1 )fs ) Ts α∈Ax a k1 ∈Z  δ ((ν2 − k2 ) − (−)(α/fs − (ν1 − k1 )))

(3.187)

k2 ∈Z

Proof. See Section 3.10. Corollary 3.18. From (3.187), it follows that the discrete-time signal x(n)  xa (t)|t=nTs obtained by uniformly sampling a continuous-time ACS signal xa (t) is a discrete-time ACS signal. In fact, its Loève bifrequency spectrum has support concentrated on lines with slope ±1 in the bifrequency plane (ν1 , ν2 ) and can be written in the form (1.42).  in (1.42) is linked to Ax a by the relationships The countable set A 

= (3.188)  α ∈ (−1/2, 1/2] :  α = (α/fs ) mod 1 A α∈Ax a



Ax a ⊆

α ∈ R : α = α fs − pfs

 (3.189)

 p∈Z  α ∈A

with mod 1 denoting the modulo 1 operation with values in [−1/2, 1/2). The (conjugate) cyclic spectrum and the (conjugate) cyclic autocorrelation function of x(n) are given in the following result. Lemma 3.19. Aliasing Formulas for Cyclic Statistics of ACS Signals (Napolitano, 1995), (Izzo and Napolitano, 1996). Let xa (t) be a continuous-time ACS signal with (conjugate) autocorrelation (3.183) and Loève bifrequency spectrum (3.184). The (conjugate) autocorrelation, cyclic autocorrelation, and cyclic spectrum of the sampled signal x(n) (3.182) are given by E{x(n + m) x (∗) (n)} = E{xa (t + τ ) xa (t)}|t=nTs ,τ =mTs ,   α−pf  α  Rx a s (τ ) R x (m) = p∈Z

τ =mTs ,α= α fs

 1   α−pfs  α  (ν) = Sx a (f − qfs ) , Sx f =νfs ,α= α fs Ts

(3.190) (3.191) (3.192)

p∈Z q∈Z

where subscript x ≡ x x (∗) . Proof. See Section 3.10. Assumption 3.20. The continuous-time process xa (t) is strictly bandlimited, i.e., Sx0a x ∗ (f ) = 0 for |f | > B (Fact 3.13). a

Chapter 3 • Almost-cyclostationary signal processing

107

According to (3.65a)–(3.66), we have the following result. Lemma 3.21. Support of the (Conjugate) Cyclic Spectrum. Under Assumption 3.20 (xa (t) strictly band-limited), it results   supp Sxαa (f ) ⊆ {(α, f ) ∈ R × R : |f |  B, |α − f |  B} ⊆ {(α, f ) ∈ R × R : |f |  B, |α|  2B} .

(3.193a) (3.193b)

Therefore, Sxαa (f ) = 0

|f | > B, |α| > 2B .

(3.194)

Moreover, in (3.184) it results   supp Sxαa (f )δ(f2 − (−)(α − f1 )) ⊆ {(f1 , f2 ) ∈ R × R : |f1 |  B, |α − f1 |  B, f2 = (−)(α − f1 )}

(3.195a)

= {(f1 , f2 ) ∈ R × R : |f1 |  B, |f2 |  B, f2 = (−)(α − f1 )} .

(3.195b)

Theorem 3.22. Let xa (t) be a continuous-time ACS signal with Loève bifrequency spectrum (3.184). Under Assumption 3.20 and for fs  4B, it results 

1  α Sx a (ν1 fs ) δ (ν2 − (−)(α/fs − ν1 )) E X(ν1 ) X (∗) (ν2 ) = Ts α∈Ax a

1 1 |ν1 |  , |ν2 |  . 2 2

(3.196)

Moreover, 

E X(ν1 ) X (∗) (ν2 ) = 0 for 1/4  |ν1 |  1/2 and 1/4  |ν2 |  1/2 .

(3.197)

Proof. See Section 3.10. Theorem 3.22 provides a lower bound to the sampling frequency fs such that the Loève bifrequency spectrum of the discrete-time signal x(n) in the principal frequency domain (ν1 , ν2 ) ∈ [−1/2, 1/2]×[−1/2, 1/2] is identical, but for frequency scaling factors, to the Loève bifrequency spectrum of the continuous-time signal xa (t). Note that the condition fs  4B is more restrictive with respect to fs  2B that is sufficient for having non-overlapping replicas in (3.187) or for the mean-square reconstruction of the process by the Shannon interpolation formula (Gardner, 1972a) (see also the proof of Theorem 3.22 in Section 3.10). Theorem 3.23. Sampling Theorem for (Conjugate) Cyclic Autocorrelation of ACS Signals (Napolitano, 1995), (Izzo and Napolitano, 1996). Let xa (t) be a continuous-time ACS signal

108

Cyclostationary Processes and Time Series

with (conjugate) autocorrelation (3.183) and Loève bifrequency spectrum (3.184). Under Assumption 3.20 and for fs  4B, it results:  1 α α   | α|  R x (m) = Rx a (τ ) τ =mTs ,α= α fs 2 ⎧  fs ⎨  α    |α|  Rx (m)  α =α/fs 2 Rxαa (τ )τ =mT = s ⎩ 0 otherwise    τ α  α  −m ∀τ ∈ R, α =  α fs Rx a (τ ) = Rx (m) sinc Ts

(3.198a) (3.198b)

(3.199)

m∈Z

(∗)

where subscripts x a ≡ xa xa and x ≡ x x (∗) . Proof. The results are consequence of (3.191) and (3.193b). Theorem 3.24. Sampling Theorem for (Conjugate) Cyclic Spectrum of ACS Signals (Napolitano, 1995), (Izzo and Napolitano, 1996). Let xa (t) be a continuous-time ACS signal with (conjugate) autocorrelation (3.183) and Loève bifrequency spectrum (3.184). Under Assumption 3.20 and for fs  4B, it results:  1 α α  (ν) = Sx a (f )f =νf ,α= Sx α fs s T ⎧s  ⎨ T   α  s Sx (ν) ν=f/fs , α =α/fs Sxαa (f ) = ⎩ 0

1 1 | α |  , |ν|  2 2 fs fs |α|  , |f |  2 2 otherwise

(3.200a) (3.200b)

By Fourier transforming both sides in (3.199) we get the following expression for the cyclic spectra of a strictly band-limited signal  α −j 2πf mTs  Sxαa (f ) = Ts rect(f Ts ) , α = α fs (3.201) R x (m) e m∈Z

which is equivalent to (3.200b). Proof. The results are consequence of (3.192) and (3.193b).

3.6.2 (Almost-)cyclostationarity of a sampled cyclostationary signal Let xa (t) be cyclostationary with period T0 . It results that (Theorem 1.4) Sxαa (f ) ≡ 0 for α = m/T0 , m ∈ Z .

(3.202)

In such a case, in the double sum in (3.192), only terms such that ( α − p)fs is an integer multiple of 1/T0 give nonzero contribution. Thus, (3.192) specializes into 1   ( α −p)fs α  (ν) = Sx a ((ν − q)fs ) δ( Sx α −p) mod (Ts /T0 ) Ts p∈Z q∈Z

Chapter 3 • Almost-cyclostationary signal processing

=

1   m/T0 Sx a ((ν − q)fs ) δ( α −mTs /T0 ) mod 1 Ts

109

(3.203)

m∈Z q∈Z

where δα is the Kronecker delta, that is, δα = 1 for α = 0 and δα = 0 for α = 0. From (3.203) it follows that, for a given conjugation configuration (xx ∗ or xx), the cyclostationarity properties of the discrete-time signal x(n) depend on the ratio Ts /T0 . Let Ts /T0 be a rational number, say q0  n(q0 )/d(q0 ) with numerator n(q0 ) and denominator d(q0 ) relatively prime integers. We have that α (ν) ≡ 0 only for Sx if d(q0 ) = 1 (i.e., Ts is an integer multiple of T0 ), then x(n) is WSS, i.e.,   α ∈ Z; • if d(q0 ) > 1, then x(n) is cyclostationary with period d(q0 ) i.e., all cycle frequencies are multiple of the rational number 1/d(q0 ). That is,  α = k/d(q0 ), k = 0, 1, . . . , d(q0 ) − 1.

•

If Ts and T0 are incommensurate, the discrete-time signal x(n) exhibits cyclostationarity, for the considered conjugation configuration, with cycle frequencies  α = (mTs /T0 ) mod 1, m ∈ Z. Therefore the cycle frequencies are not all integer multiple of a fundamental rational cycle frequency, that is, x(n) is not cyclostationary but, rather, almost-cyclostationary. If the cyclostationary continuous-time signal xa (t) is strictly band limited to [−B, B], then it exhibits cyclostationarity only in correspondence of the finite number of cycle frequencies α = m/T0 , m ∈ Z, such that |α|  2B (see Fig. 3.3). Therefore, if x(n) is obtained by sampling xa (t) at rate fs  4B, then from (3.200a) we have that x(n) exhibits cyclostationarity with cycle frequencies (considered in the interval [0, 1[)  α=

Ts α =m , fs T0

0 < m  int(4BT0 ) 

T0 Ts

(3.204)

where int(·) denotes integer part. As already observed, x(n) is cyclostationary if Ts /T0 is rational and is almost-cyclostationary otherwise. In general, no information on the cyclostationarity properties of a continuous-time signal xa (t) can be achieved by the knowledge of the cyclostationarity properties of the sampled signal x(n). However, if x(n) is obtained by sampling a strictly band limited signal at a rate fs  4B, then (3.200b) holds. Consequently, if x(n) exhibits cyclostationarity with cycle frequencies given by (3.204), then the continuous-time signal xa (t) is wide-sense cyclostationary with period T0 , whether T0 /Ts is a rational number or not.

3.6.3 Bandpass sampling From Theorems 3.23 and 3.24 it follows that the continuous-time signal should be sampled at twice the Nyquist frequency in order to avoid aliasing in both cycle and spectral frequency domains. Thus, for bandpass signals in communications and optics applications, very high sampling rates are required. The problem of sub-Nyquist sampling for cyclic spectrum reconstruction can be addressed resorting to compressive sampling (Section 10.2). An alternative way to reduce the

110

Cyclostationary Processes and Time Series

sampling rate without demodulating the bandpass signal is to resort to bandpass sampling (Napolitano and Do˘gançay, 2018). It consists in sampling a bandpass signal with a sampling frequency significantly lower than the Nyquist frequency of the bandpass signal and comparable with the Nyquist frequency of the complex envelope. Sufficient conditions involving sampling frequency and carrier frequency and the bandwidth of the bandpass signal must be satisfied in order to avoid the overlap of replicas in the Fourier transform of the sampled signal and to allow one to reconstruct the continuous-time complex envelope. Conditions for the reconstruction of the cyclic spectrum can be derived under assumptions less restrictive than those made for compressive sampling based techniques (Napolitano and Do˘gançay, 2018). By substituting the expression (3.121) (with x replaced by xa ) of the cyclic spectrum of xa (t) into the expression (3.192) xa (t) expressed in terms of the (conjugate) cyclic spectra of  of the cyclic spectrum of x(n) expressed in terms of replicas of cyclic spectra of xa (t) leads to k 2   α−pfs α−pfs  α  Sxx (ν) = S xa  xa∗ (f − f0 − qfs ) + S xa∗ xa (f + f0 − qfs ) 4 p∈Z q∈Z

α−2f0 −pfs

+ S xa  xa

α+2f0 −pfs

(f − f0 − qfs ) + S x ∗ x∗ a a

 (f + f0 − qfs )

f =νfs ,α= α fs

.

(3.205)

In the following, sufficient conditions on the sampling frequency fs , the frequency f0 , and the signal bandwidth W are obtained such that the double sum in (3.205) involves terms with non-overlapping supports in the ( α , ν) plane. In addition, more restrictive conditions are derived such that, for each  α ∈ [−1/2, 1/2), only one replica in (3.205) gives nonzero contribution for ν ∈ [−1/2, 1/2). Theorem 3.25. Band Pass Sampling (Napolitano and Do˘gançay, 2018). Let xa (t) be a band pass ACS signal with central frequency f0 and strictly band limited with monolateral bandwidth W , with f0 > W (Definition 3.16). 1. For every integer m and sampling frequency fs such that  m  f0 /(2W ) − 1/2 fs = 2f0 /(m + 1/2) is verified, replicas in (3.205) do not overlap. 2. For every integer m and sampling frequency fs such that  m  f0 /(4W ) − 1/2 fs = 2f0 /(m + 1/2)

(3.206)

(3.207)

is verified, only one replica in (3.205), for each  α ∈ [−1/2, 1/2), gives nonzero contribution for ν ∈ [−1/2, 1/2). Proof. See Section 3.10.

Chapter 3 • Almost-cyclostationary signal processing

111

The smallest sampling frequency satisfying (3.206) is obtained by taking m = mmax  floor(f0 /(2W ) − 1/2)

(3.208)

and the admissible values for fs that allow bandpass sampling are fs = 2f0 /(m + 1/2)

m = 1, 2, . . . , mmax .

(3.209)

The smallest sampling frequency satisfying (3.207) is obtained by m = mmax  floor(f0 /(4W ) − 1/2) .

(3.210)

The set of all admissible sampling frequencies is obtained by replacing mmax with mmax in (3.209). If condition (3.207) is satisfied, then for each fixed  α ∈ [−1/2, 1/2), the cyclic spectrum  α (ν), with ν ∈ [−1/2, 1/2), is proportional to only one of the four (frequency scaled) cyclic  Sxx spectra of the complex envelope  xa (t) in the right-hand side of (3.205). In Fig. 3.7, an example is illustrated corresponding to two different values of sampling frequency. Let fr be a reference frequency. In the first case (Fig. 3.7 (top)), fs = fr /4 and condition (3.207) is satisfied. That is, in the principal domain ( α , ν) = (α/fs , f/fs ) ∈ 2 α only one cyclic spectrum in (3.205) gives nonzero contribution. In [−1/2, 1/2) , for each  the second case (Fig. 3.7 (bottom)), fs = fr /6 and condition (3.206) (but not (3.207)) is satisfied. In such a case, replicas in (3.205) do not overlap, but for each  α ∈ [−1/2, 1/2) either one or two cyclic spectra in (3.205) can give nonzero contribution for ν ∈ [−1/2, 1/2). Note that the reciprocal positions of the replicas of the four parallelograms depend on the value of fs and are not necessarily coincident with the reciprocal positions as in Fig. 3.6. In practice, since f0  W , the above results on bandpass sampling allow one to choose a sampling frequency fs for the continuous-time signal xa (t) which is significantly smaller than the sampling frequency fN = 4(f0 + W ) (i.e., twice the Nyquist frequency) which is necessary to have no aliasing in both f and α domains (Theorem 3.24) when the bandpass nature of xa (t) is not exploited. In particular, for m = mmax we have fs =

2f0 4W floor(f0 /(2W ) − 1/2) + 1/2

(3.211)

which is the sampling frequency to avoid aliasing in both α and f for a low-pass ACS signal with monolateral bandwidth W (Theorem 3.24).

3.7 Multirate processing of discrete-time ACS signals In multirate digital signal processing, the main building blocks are constituted by expansors (upsamplers) and decimators (downsamplers). These devices are linear time-variant systems that can be classified as deterministic in the FOT probability framework (Izzo

112

Cyclostationary Processes and Time Series

FIGURE 3.7 Set D3 (α, f ) defined in (3.274) containing the support of the cyclic spectrum of the bandpass sampled bandpass signal. (Top) sampling frequency fs = fr /4 (condition (3.207) satisfied). (Bottom) fs = fr /6 (condition (3.206) ˘ satisfied and (3.207) not satisfied). Source: (Napolitano and Dogançay, 2018) Copyright of Elsevier.

and Napolitano, 2002a), (Napolitano, 2012, Sec. 6.3.8). In particular, these systems are not LAPTV. The effects of multirate systems on second-order cyclostationary processes are analyzed in (Sathe and Vaidyanathan, 1993) and (Akkarakaran and Vaidyanathan, 2000). A cyclostationary analysis of multirate filter banks is made in (Aach, 2007). Higher-order statistics of expanded or decimated ACS time series are considered in (Izzo and Napolitano, 1998b).

Chapter 3 • Almost-cyclostationary signal processing

113

In (Sathe and Vaidyanathan, 1993), (Izzo and Napolitano, 1998b), and (Akkarakaran and Vaidyanathan, 2000), multirate operations on a single signal are addressed and no cross-statistics of a signal and its expanded or decimated version are considered. Moreover, no cross-statistics are considered between signals obtained by multirate elaborations with different rates of the same signal or of jointly ACS signals. In these cases, since the elaborations are not LAPTV, the output signals are not jointly ACS but, rather, jointly SC (Section 13.4). In (Napolitano, 2012, Sec. 4.10), multirate processing of jointly SC processes is considered. Results for jointly ACS processes reported here are obtained as special case of the results in (Napolitano, 2012, Sec. 4.10). Let x1 (n) and x2 (n) be discrete-time jointly ACS processes with (conjugate) crosscorrelation function  

(∗) α j 2π αn  (3.212) R E x1 (n + m) x2 (n) = x (m) e   α ∈A (∗)

where subscript x ≡ x1 x2 , the Fourier coefficients N

  1 (∗) αn E x1 (n + m) x2 (n) e−j 2π N→∞ 2N + 1

α  R x (m)  lim

(3.213)

n=−N

are the (conjugate) cyclic cross-correlation functions, and   α    A α ∈ [−1/2, 1/2) : R x (m) ≡ 0

(3.214)

is the countable set of (conjugate) cycle frequencies  α in the principal domain [−1/2, 1/2). By double Fourier transforming both sides of (3.212) with n1 = n + m and n2 = n, one obtains the Loève bifrequency cross-spectrum  

(∗) α  (ν1 )  δ (ν2 − (−)( α − ν1 )) (3.215) Sx E X1 (ν1 ) X2 (ν2 ) =   α ∈A

where  δ (·) denotes the periodic delta train with period 1 (C.18) and  α α −j 2πνm   (ν) = R Sx x (m) e

(3.216)

m∈Z

are the (conjugate) cyclic cross-spectra.

3.7.1 Expansion (upsampling) Let xI i (n) be the Li -fold expanded (upsampled) version of xi (n), (i = 1, 2) (see Fig. 3.8). That is ⎧ $ % ! " ⎪ ⎨ x n n = kLi , k ∈ Z n i Li (3.217) xI i (n) = xi  ⎪ Li ⎩ 0 otherwise

114

Cyclostationary Processes and Time Series

FIGURE 3.8 Expansion of two processes with different expansion factors.

whose Fourier transform is XI i (ν) = Xi (νLi ) .

(3.218)

Theorem 3.26. (Napolitano, 2012, Sec. 4.10.1.2). Let x1 (n) and x2 (n) be jointly ACS with Loève bifrequency cross-spectrum (3.215). The Loève bifrequency cross-spectrum of the expanded processes xI 1 (n) and xI 2 (n) is given by 

(∗) E XI 1 (ν1 ) XI 2 (ν2 ) 

(∗) = E X1 (ν1 L1 ) X2 (ν2 L2 )  α  = (ν1 L1 )  δ (ν2 L2 − (−)( α − ν1 L1 )) (3.219a) Sx   α ∈A

$ $ %% 1  h L1  α  α  = δ ν2 − − (−) − ν1 Sx (ν1 L1 ) L2 L2 L2 L2  h∈Z  α ∈A $ $ %% L 2 −1  1 p L  α 1 α   (ν1 L1 ) − (−) − ν1 δν2 ν2 − Sx . = L2 L2 L2 L2 

  α ∈A

(3.219b)

(3.219c)

p=0

Proof. It is a consequence of the scaling property of the Dirac delta (Zemanian, 1987, Sec. 1.7). See also (Napolitano, 2012, Sec. 5.11). α (ν L ) are periodic in ν with period 1/L , and In (3.219a)–(3.219c), the cyclic spectra  Sx 1 1 1 1 the delta train is periodic in ν1 with period 1/L1 and is periodic in ν2 with period 1/L2 . Notation  δν2 (·) clarifies that the periodicity with period 1 must be considered with respect to the variable ν2 . (∗) Therefore, E{XI 1 (ν1 ) XI 2 (ν2 )} is periodic in ν1 with period 1/L1 and in ν2 with period 1/L2 . In addition, from (3.219c) it follows that the Loève bifrequency cross-spectrum of the expanded processes xI 1 (n) and xI 2 (n) in the principal domain (ν1 , ν2 ) ∈ [−1/2, 1/2)2 has support (for the case (−) present)

ν2 = (L1 /L2 )ν1 −  α /L2 + p/L2 mod 1

 p ∈ {0, 1, . . . , L2 − 1}  α ∈ A,

(3.220)

that is, lines with slope L1 /L2 . Consequently, xI 1 (n) and xI 2 (n) are jointly SC (Section 13.4). In particular, for a single process x1 (n) ≡ x2 (n) ≡ x(n), by taking L1 = 1, that is xI 1 (n) ≡ x(n), it follows that the ACS process x(n) and its expanded version xI 2 (n) are jointly SC with

Chapter 3 • Almost-cyclostationary signal processing

115

support lines with slope 1/L2 . In the special case of x(n) WSS, x(n) and xI 2 (n) are jointly SC with support constituted by a unique line. Theorem 3.27. (Napolitano, 2012, Sec. 4.10.1.2). Let x1 (n) and x2 (n) be jointly ACS with (conjugate) cross-correlation (3.212). The (conjugate) cross-correlation function of the expanded processes xI 1 (n) and xI 2 (n) is given by

 (∗) E xI 1 (n + m) xI 2 (n) 

 (∗) = E XI 1 (ν1 ) XI 2 (ν2 ) ej 2π[ν1 (n+m)+(−)ν2 n] dν1 dν2 [−1/2,1/2]2 ! "  /L )n + m (1 − L 1 2 α /L2 )n  = ej 2π( (3.221) Rxα δn mod L2 L1   α ∈A

α  where the Kronecker delta δn mod L2 is defined in (C.17) and R x [ /L2 ] is defined according to (3.217), that is, ! "  α  = kL, k ∈ Z R x ( /L) α   R (3.222) x L 0 otherwise

Proof. See (Napolitano, 2012, Sec. 5.11.2). α  When L1 = L2 , the argument of the cyclic cross-correlation function R x (·) contains the α  time variable n. Consequently, under the mild assumption that R x (·) ∈ 1 (Z), it follows that, for every fixed m, the cross-correlation function in (3.221) as a function of n does not con(∗) tain any finite-strength additive sinewave component. That is, xI 1 (n) and xI 2 (n) do not exhibit joint cyclostationarity at any cycle frequency: N

  1 (∗) αn E xI 1 (n + m) xI 2 (n) e−j 2π =0 (∗) (m)  lim xI 1 xI 2 N→∞ 2N + 1

α  R

∀ α ∈ [−1/2, 1/2) .

n=−N

(3.223) 0 In particular, their time-averaged (conjugate) cross-correlation R

(∗)

xI 1 xI 2

(m) is identically

zero. Moreover, from (3.221) it follows that xI 1 (n) and xI 2 (n) are jointly ACS if and only if L1 = L2 = L. In such a case, their cyclic cross-correlation function is (Napolitano, 2012, Sec. 5.11.3) ! " 1 m  α  α L  .  (∗) (m) = R  α L mod 1 ∈ A (3.224) R xI 1 xI 2 L x L Furthermore, since from (3.224) it follows α /L α   R (∗) (mL) x (m) = LR xI 1 xI 2

(3.225)

116

Cyclostationary Processes and Time Series

in accordance with the Fourier-transform pair (3.217) and (3.218), the cyclic cross-spectra of xI 1 (n) and xI 2 (n) are α  S

(∗)

xI 1 xI 2

(ν) =

1  S α L (νL) L x

.  α L mod 1 ∈ A

(3.226)

In the special case of x1 ≡ x2 , (3.224) and (3.226) reduce to (Izzo and Napolitano, 1998b, eqs. (32) and (34)) specialized to second order. The problem of image suppression in (3.226) is addressed in (Napolitano, 2012, Sec. 4.10.6).

3.7.2 Sampling Let xδi (n) be the sampled version of xi (n), with sampling period Mi (i = 1, 2). That is xδi (n)  xi (n)

+∞ 

δn− Mi =

=−∞

+∞ 

xi ( Mi ) δn− Mi .

(3.227)

=−∞

Accounting for (C.17), the Fourier transforms of the sampled processes are Mi −1  Mi −1  1  1  p  p   Xi ν − = δ ν− Xδi (ν) = Xi (ν) ⊗ Mi Mi Mi Mi p=0

(3.228)

p=0

where ⊗ denotes periodic convolution with period 1. Theorem 3.28. (Napolitano, 2012, Sec. 4.10.2). Let x1 (n) and x2 (n) be jointly ACS, with Loève bifrequency cross-spectrum (3.215). The Loève bifrequency cross-spectrum of the sampled processes xδ1 (n) and xδ2 (n) is given by

 (∗) E Xδ1 (ν1 ) Xδ2 (ν2 ) =

M 1 −1 M 2 −1     1 p1  α  ν1 − Sx M1 M2 M1  p1 =0 p2 =0  α ∈A    p2 p1  − (−)  α − ν1 + . δ ν2 − M2 M1

(3.229)

Proof. It is an immediate consequence of (3.215) and (3.228). From (3.229) it follows that the sampled processes are in turn jointly ACS for every value of M1 and M2 . Such a result is not surprising since sampling is a linear periodically timevariant transformation (Napolitano, 1995). Theorem 3.29. (Napolitano, 2012, Sec. 4.10.2). For M1 = M2 = M, the (conjugate) cyclic cross-correlation function of xδ1 (n) and xδ2 (n) is given by α  R

(∗) (m) =

xδ1 xδ2

M−1 1   α −q/M (m) δm mod M . Rx M q=0

(3.230)

Chapter 3 • Almost-cyclostationary signal processing

117

FIGURE 3.9 Decimation of two processes with different decimation factors.

Proof. See (Napolitano, 2012, Sec. 5.11.4). Theorem 3.30. (Napolitano, 2012, Sec. 4.10.2). For M1 = M2 = M, the (conjugate) cyclic cross-spectrum of the sampled processes xδ1 (n) and xδ2 (n) is given by α  S

(∗) xδ1 xδ2

(ν) =

M−1 M−1 1   1   p α −q/M  (ν) ⊗ δ ν− Sx M M M q=0

=

1 M2

M−1  M−1  p=0 q=0

p=0

 p  α −q/M  . ν− Sx M

(3.231)

Proof. In (3.231), ⊗ denotes periodic convolution with period 1. By Fourier transforming both sides of (3.230), using (C.17) and the fact that  δ (ν) is the unit element for the periodic convolution with period 1, we obtain (3.231).

3.7.3 Decimation (downsampling) Let xDi (n) be the decimated (downsampled) version of xi (n), with decimation factor Mi (i = 1, 2) (Fig. 3.9). That is, xDi (n) = xi (nMi ) = xδi (nMi ) .

(3.232)

The (conjugate) cross-correlation functions of decimated, sampled, and original signals are linked by the relationship 



(∗) (∗) E xD1 (n + m) xD2 (n) = E x1 ((n + m)M1 ) x2 (nM2 ) 

(∗) (3.233) = E xδ1 ((n + m)M1 ) xδ2 (nM2 ) . Since the sampled process is an expanded version of the decimated process, that is xδi (n) = xDi [n/Mi ], accounting for (3.218), the Fourier transform of the sampled and decimated processes are linked by Xδi (ν) = XDi (νMi ). Thus, using (3.228), the Fourier transforms of the decimated processes are XDi (ν) = Xδi

Mi −1  ν  ν − p 1  Xi = . Mi Mi Mi p=0

(3.234)

118

Cyclostationary Processes and Time Series

Theorem 3.31. (Napolitano, 2012, Sec. 4.10.3.2). Let x1 (n) and x2 (n) be jointly ACS, with Loève bifrequency cross-spectrum (3.215). The Loève bifrequency cross-spectrum of the decimated processes xD1 (n) and xD2 (n) is given by 

(∗) E XD1 (ν1 ) XD2 (ν2 ) =

M1 −1      1  M2  α ν1 − p1   δν2 ν2 − (−)  α M2 − (ν1 − p1 ) . Sx M1 M1 M1

(3.235)

 p1 =0  α ∈A

Proof. It is consequence of (3.234), the scaling property of Dirac delta (Zemanian, 1987, Sec. 1.7), and identity (Napolitano, 2012, Eq. (5.192)). α ((ν − p )/M ) are periodic in ν with period 1 and the delta train is The cyclic spectra  Sx 1 1 1 1 periodic in both ν1 and ν2 with period 1. In addition, from (3.235) it follows that the Loève bifrequency cross-spectrum of the decimated processes xD1 (n) and xD2 (n) in the principal domain (ν1 , ν2 ) ∈ [−1/2, 1/2)2 has support (in the case (−) present)

ν2 = (M2 /M1 )ν1 −  α M2 − (M2 /M1 )p1 mod 1   α ∈ A, p1 ∈ {0, 1, . . . , M1 − 1}

(3.236)

that is, lines with slope M2 /M1 . Consequently, unlike the case of sampled processes, the decimated processes xD1 (n) and xD2 (n) are jointly SC. In particular, for a single process x1 (n) ≡ x2 (n) = x(n), by taking M1 = 1, that is xD1 (n) ≡ x(n), it follows that the ACS process x(n) and its decimated version xD2 (n) are jointly SC with support lines with slope M2 . In the special case of x(n) WSS, the support line is unique. Theorem 3.32. (Napolitano, 2012, Sec. 4.10.3.2). Let x1 (n) and x2 (n) be jointly ACS with (conjugate) cross-correlation (3.212). The (conjugate) cross-correlation function of the decimated processes xD1 (n) and xD2 (n) is given by 

(∗) E xD1 (n + m) xD2 (n) 

 (∗) = E XD1 (ν1 ) XD2 (ν2 ) ej 2π[ν1 (n+m)+(−)ν2 n] dν1 dν2 [−1/2,1/2]2    α M2 n  = ej 2π (3.237) Rxα (M1 − M2 )n + M1 m .   α ∈A

Proof. See (Napolitano, 2012, Sec. 5.11.6). α 1  By reasoning as for (3.221), we have that if R x (·) ∈ (Z) then for M1 = M2 the right-hand side of (3.237), as a function of n, does not contain any finite-strength additive sinewave (∗) component. That is, for M1 = M2 the processes xD1 (n) and xD2 (n) do not exhibit joint cyclostationarity at any cycle frequency. From (3.235) and (3.237) it follows that xD1 (n) and xD2 (n) are jointly ACS if and only if M1 = M2 = M.

Chapter 3 • Almost-cyclostationary signal processing

119

Theorem 3.33. (Napolitano, 2012, Sec. 4.10.3.2). For M1 = M2 = M, the (conjugate) cyclic cross-correlation function of the decimated processes xD1 (n) and xD2 (n) is given by α  R

(∗) (m) =

xD1 xD2

M−1 

α −q)/M x( (mM) . R

(3.238)

q=0

Proof. See (Napolitano, 2012, Sec. 5.11.7). Theorem 3.34. (Napolitano, 2012, Sec. 4.10.3.2). For M1 = M2 = M, the (conjugate) cyclic cross-spectrum of the decimated processes xD1 (n) and xD2 (n) is given by α  S

(∗) (ν) =

xD1 xD2

M−1 M−1   1   ( α −q)/M ν − p . Sx M M

(3.239)

p=0 q=0

Proof. Eq. (3.239) is obtained by Fourier transforming both sides of (3.238) and accounting for the Fourier-transform pair (3.232) and (3.234). In the special case x1 ≡ x2 , (3.238) and (3.239) reduce to (Izzo and Napolitano, 1998b, eqs. (22) and (24)) specialized to second order. In addition, by comparing (3.230) with (3.238) the following relationship between (conjugate) cyclic cross-correlations of decimated and sampled signals can be established α  R

(∗)

xD1 xD2

α /M  (m) = M R (∗) (mM)

or, equivalently α  R (∗) (m) = xδ1 xδ2

xδ1 xδ2

! " 1  m αM R . (∗) M xD1 xD2 M

Moreover, by comparing (3.231) with (3.239) it results $ % ν  α /M  α   . S (∗) (ν) = M S (∗) xD1 xD2 xδ1 xδ2 M

(3.240)

(3.241)

(3.242)

3.8 Special topics Special topics related to filtering and sampling of cyclostationary and almost-cyclostationary signals are: •

Sampling with jitter (Dehay and Monsan, 2007); • Models for time variant systems (Duverdier et al., 1999), (Duverdier and Lacaze, 2000), (Huang et al., 2002b), (Guidorzi and Diversi, 2003), (Sadeghi and Yu, 2016);

120

• •

Cyclostationary Processes and Time Series

Digital-to-analog conversion (Mandyam, 2002); Cyclic spectral analysis with missing observations (Dunsmuir, 1984), (Dandawaté and Giannakis, 1993), (Giannakis and Zhou, 1994), (Drake et al., 2014), (Drake et al., 2015).

3.9 Summary LAPTV systems, also referred to as FRESH filters, transform ACS processes into ACS processes with different cyclic features. That is, the class of the ACS processes is closed under LAPTV transformations. Input/output relations in terms of (conjugate) cross-correlation (Theorem 3.6), (conjugate) cyclic cross-correlations and (conjugate) cyclic cross-spectra (Corollary 3.7) are derived for two LAPTV systems with two jointly ACS inputs. Input/output relations for a single system with a single input are obtained a special case. Cyclic statistics of the product of ACS processes are also derived (Theorem 3.10). Cyclic statistics of signals involved in the Rice’s representation of a real-valued signal are derived (Section 3.5). Specifically, cyclic statistics of real-valued signal, analytic signal, complex envelope, and in-phase and in-quadrature components are linked among them. Cyclic statistics of the discrete-time signal obtained by uniformly sampling with sampling frequency fs a continuous-time ACS signal are derived. It is shown that the discretetime signal is ACS (Corollary 3.18) and its cyclic statistics are obtained by the superposition of replicas of those of the continuous-time ACS signal (Lemmas 3.17 and 3.19). For a continuous-time signal with bandwidth B, the sufficient condition fs  4B is derived such that the Loève bifrequency spectrum and the (conjugate) cyclic spectrum of the discrete-time sampled signal are amplitude- and frequency-scaled versions of the corresponding functions of the continuous-time signal in the whole main cycle and spectral frequency domains (Theorems 3.22 and 3.24). A similar sampling theorem is proved for the (conjugate) cyclic autocorrelation function (Theorem 3.23). Condition fs  4B is more restrictive than condition fs  2B which is sufficient to have non-overlapping replicas in the aliasing formulas or to reconstruct the process by the Shannon interpolation formula. For bandpass signals adopted in communications, the anti-aliasing condition requires very high sampling rates. This problem can be circumvented by bandpass sampling. It consists in sampling the continuous-time bandpass signal by a sampling frequency significantly lower than the Nyquist frequency 2B and comparable with the Nyquist frequency of the complex envelope. Admissible values of the sampling frequency for bandpass sampling are derived such that (conjugate) cyclic spectra of the complex envelope can be reconstructed by samples of the continuous-time signal (Theorem 3.25). The effects of expansors (upsamplers) (Theorems 3.26 and 3.27) and decimators (downsamplers) (Theorems 3.31, 3.32, 3.33, 3.34) on ACS processes are analyzed since these devices are the main building blocks of multirate systems. These devices are linear timevariant systems that are not LAPTV. The expanded (upsampled) version of an ACS process is ACS. Similarly, the decimated (downsampled) version of an ACS process is ACS. How-

Chapter 3 • Almost-cyclostationary signal processing

121

ever, expanded or decimated versions of the same ACS signal with different expansion or decimation factors are not jointly ACS but, rather, jointly spectrally correlated.

3.10 Proofs Proof of Theorem 3.5 Here the proof of the FOT counterpart of Theorem 3.5 is reported, since the proof in the FOT approach is more challenging than the analogous proof in the stochastic process framework. {α}

E

{α}

{yi (t)} = E

R

{α}

=E =E







R σ ∈J i i

 

R σ ∈J i i {α}

=E =



hi (t, u) xi (u) du

R σ ∈J i i

{α}

=

&





& hσi (t − u) e

j 2πσi u

xi (u) du &

hσi (s) e

j 2πσi (t−s)

xi (t − s) ds

hσi (s) ej 2πσi (t−s) E{α} {xi (t − s)} ds 

R σ ∈J i i

 

hσi (s) e 

x ηi

σi ∈Ji ηi ∈E1i

R

j 2πσi (t−s)



& x ηi e

j 2πηi (t−s)

ds

ηi ∈E1i

hσi (s) e−j 2π(ηi +σi )s ds ej 2π(ηi +σi )t .

(3.243)

In the third equality the variable change t − u = s is made. This allows to interchange the order of integral and almost-periodic component extraction operator (Theorem 2.32 and Remark 2.30). Eq. (3.5) follows accounting for (3.12). Assumption 3.2a assures that 

hσi (s) e−j 2πσi s ∈ L1 (R)

(3.244)

σi ∈Ji

and Assumption 3.3 assures that    x η  < ∞ . i ηi ∈E1i

Thus the assumptions for Theorem 2.32 are satisfied.

(3.245)

122

Cyclostationary Processes and Time Series

Proof of Theorem 3.6 Here the proof of the FOT counterpart of Theorem 3.6 is reported, since the proof in the FOT approach is more challenging than the analogous proof in the stochastic process framework. By substituting (3.7) into (3.6) we have E{α} {y1 (t + τ ) y2 (t)}    (∗) (∗) {α} =E h1 (t + τ, u1 ) x1 (u1 ) du1 h2 (t, u2 ) x2 (u2 ) du2 R ⎧R ⎨  = E{α} hσ1 (t + τ − u1 ) ej 2πσ1 u1 x1 (u1 ) du1 ⎩ R σ1 ∈J1 ⎫   ⎬ (−)j 2πσ2 u2 (∗) h(∗) (t − u ) e x (u ) du 2 2 2 σ2 2 ⎭ R σ ∈J 2 2 ⎧ ⎨  hσ1 (s1 + τ ) ej 2πσ1 (t−s1 ) x1 (t − s1 ) ds1 = E{α} ⎩ R σ1 ∈J1 ⎫   ⎬ (−)j 2πσ2 (t−s2 ) (∗) h(∗) (s ) e x (t − s ) ds 2 2 σ2 2 2 ⎭ R (∗)

(3.246)

σ2 ∈J2

where, in the last equality, the variable changes s1 = t − u1 and s2 = t − u2 are made. Under Assumptions 3.2a and 3.4b, or alternatively under Assumptions 3.2b and 3.4a, the order of integration and summation operations can be interchanged. In fact, the Fubini and Tonelli theorem (Champeney, 1990, Chapter 3) and the dominated convergence theorem (Champeney, 1990, Chapter 4) can be applied. In addition, under the above assumptions also the almost-periodic component extraction operator can be interchanged with the integration and summation operations (Theorem 2.32). Finally, by observing that almost-periodic functions can be led outside the almost-periodic component extraction operator (see also Remark 2.30), we get E{α} {y1 (t + τ ) y2 (t)}     −j 2πσ1 s1 −(−)j 2πσ2 s2 = hσ1 (s1 + τ ) h(∗) e σ2 (s2 ) e (∗)

σ1 ∈J1 σ2 ∈J2 R R



(∗) E{α} x1 (t − s1 ) x2 (t − s2 ) ds1 ds2 ej 2π(σ1 +(−)σ2 )t      = R α (∗) (s2 − s1 ) hσ1 (s1 + τ ) h(∗) σ2 (s2 ) α∈A12 σ1 ∈J1 σ2 ∈J2 R R

x1 x2

e−j 2πσ1 s1 e−(−)j 2πσ2 s2 e−j 2παs2 ds1 ds2 ej 2π(α+σ1 +(−)σ2 )t

Chapter 3 • Almost-cyclostationary signal processing

=

    α∈A12 σ1 ∈J1 σ2 ∈J2



R

R

Rα

(∗)

x1 x2

123

(u) ej 2πσ1 u

−j 2π(α+σ1 +(−)σ2 )s2 hσ1 (s2 − u + τ ) h(∗) ds2 du ej 2π(α+σ1 +(−)σ2 )t σ2 (s2 ) e

(3.247)

where in the second equality, according to (1.74) (with the replacements t − s1  t + τ and t − s2  t), we put 

 (∗) E{α} x1 (t − s1 ) x2 (t − s2 ) = Rα α∈A12

(∗)

x1 x2

(s2 − s1 ) ej 2πα(t−s2 )

(3.248)

and in the last equality the variable change s2 − s1 = u is made. By substituting (3.14) into (3.247), we get the FOT counterpart of (3.13) which is formally the same as (3.13). By replacing (3.14) into (3.247), for the generic term of the multiple sum in (3.247) it results that  α R

  1 +(−)σ2 (τ ) ej 2πσ1 τ ⊗ rσα+σ (τ ) 1 σ2 (∗) τ     1 +(−)σ2 (τ − s) ds =  R α (∗) (s) ej 2πσ1 s rσα+σ  1 σ2 (∗) x x 1 2 R      1 +(−)σ2  ds (τ − s)  R α (∗) (s) rσα+σ σ (∗)  1 2 x1 x2 R      1 +(−)σ2  R α (∗) (s) ds rσα+σ 1 σ2 (∗) (∗)

x1 x2

R

∞

x1 x2

(3.249)

Using Assumption 3.2 we have      (s)  ds hσ1 (τ + s) h(∗) σ2 ∞ R     hσ2 (s) ds hσ1 ∞ < ∞ .

α+σ1 +(−)σ2 rσ1 σ2 (∗)



R

(3.250)

Therefore, accounting for Assumptions 3.2 and 3.4, the series in the right-hand side of (3.13) is convergent. In addition, the assumptions for (a bi-dimensional version of) Theorem 2.32 hold. In fact, accounting for Assumptions 3.2 and 3.4 we have   1. hσ1 (s1 ) hσ2 (s2 ) ∈ L1 (R2 ). σ1

σ2

2. The function t → x1 (t − s1 ) x2 (t − s2 ) ∈ L1loc (R) ∩ L∞ (R) and R α (∗) (τ ) exists ∀α ∈ A. x1 x2  R α (∗) ∞ < ∞. 3. (∗)

α

x1 x2

124

Cyclostationary Processes and Time Series

Proof of Corollary 3.7 By taking the Fourier coefficient at frequency α of the almost periodic-function in the FOT counterpart of (3.13), we have Rα

(∗) y1 y2

  (∗) (τ )  y1 (t + τ ) y2 (t) e−j 2παt t

   β β+σ1 +(−)σ2 = rσ1 σ2 (∗) (τ ) R (∗) (τ ) ej 2πσ1 τ ⊗ τ β∈A12 σ1 ∈J1 σ2 ∈J2

x1 x2

  ej 2π(β+σ1 +(−)σ2 )t e−j 2παt .

(3.251)

  ej 2π(β+σ1 +(−)σ2 )t e−j 2παt = δα−(β+σ1 +(−)σ2 )

(3.252)

t

Since t

the FOT counterpart of (3.15), which is formally identical to (3.15), easily follows. Under Assumptions 3.2 (a and b) and 3.4a, accounting for inequalities (3.249), (3.250), we have that R α (∗) (τ ) ∈ L1 (R). Thus, its Fourier transform S α (∗) (f ) exists in ordinary y1 y2 y1 y2 sense.

Proof of (3.24a) and (3.24b) Accounting for the Fourier pairs F h(−τ ) ←→ H (−f ) F h(∗) (−τ ) ←→ H (∗) ((−)(−f )) F h(∗) (−τ ) ej 2πατ ←→ H (∗) ((−)(−(f − α))) = H (∗) ((−)(α − f )) we have  α rhh (∗) (τ ) 



R

h(s + τ ) h(∗) (s) e−j 2παs ds 

h(τ − s  ) h(∗) (−s  ) ej 2παs ds  R   = h(τ ) ⊗ h(∗) (−τ ) ej 2πατ

=



h(s + τ ) h(∗) (s) e−j 2παs ds R   F h(τ + s) h(∗) (s) =s→α   H (∗) ((−)α) = H (α) ej 2πατ ⊗ α

α rhh (∗) (τ ) 

(3.253a)

Chapter 3 • Almost-cyclostationary signal processing

125

 =

H (λ) H (∗) ((−)(α − λ)) ej 2πλτ dλ

R

(3.253b)

Therefore, we have the Fourier transform pair  α rhh (∗) (τ ) 

R

F h(t + τ ) h(∗) (t) e−j 2παt dt ←→ H (f ) H (∗) ((−)(α − f ))

(3.254)

Proof of Theorem 3.10 By using the FOT counterparts of (3.40) and (3.41), which are formally identical to (3.40) and (3.41), we get Rα

(∗) y1 y2

  (∗) (τ )  y1 (t + τ ) y2 (t) e−j 2παt t  

 (∗) {α} = E y1 (t + τ ) y2 (t) e−j 2παt   

t  (∗) (∗) {α} {α} = E c1 (t + τ ) c2 (t) E x1 (t + τ ) x2 (t) e−j 2παt t *  +  αc αx j 2παc t j 2παx t −j 2παt = R (∗) (τ ) e R (∗) (τ ) e e αc ∈A

(∗) c1 c2



=

αx ∈A

(∗) x1 x2

c1 c2

x1 x2

αx ∈A

(∗) x1 x2



R αx

αc ∈A

(∗) c1 c2

(∗)

x1 x2

t

(τ ) R αc (∗) (τ ) c1 c2

  ej 2παx t ej 2παc t e−j 2παt

(3.255)

t

from which the FOT counterpart of (3.44), formally identical to (3.44), follows observing that 

ej 2παx t ej 2παc t e−j 2παt

Under Assumption 3.4a for R α

(∗)

x1 x2

  α R

(∗) y1 y2

  (τ ) 

αx ∈A

(∗) x1 x2

Thus, the Fourier transform of R α

(∗)

y1 y2

t

= δα−αx −αc .

(3.256)

(τ ) and Assumption 3.4b for R α

(∗)

c1 c 2

  αx R





(∗) x1 x2

  (τ )

sup

αc ∈A

(∗) c1 c2

αc R

(∗) c1 c2



∞

(τ ), we have

∈ L1 (R) .

(3.257)

(τ ) exists in the ordinary sense and is given by (3.45).

The FOT counterpart of (3.45) easily follows by Fourier transforming both sides of the FOT counterpart of (3.44).

126

Cyclostationary Processes and Time Series

Proof of Lemma 3.17 Accounting for (3.185) and (3.184) we have    

1 1  (∗) (∗) Xa ((ν1 − k1 )fs ) Xa ((ν2 − k2 )fs ) E X(ν1 ) X (ν2 ) = E Ts Ts k1 ∈Z k2 ∈Z  1  

= 2 E Xa ((ν1 − k1 )fs ) Xa(∗) ((ν2 − k2 )fs ) Ts k1 ∈Z k2 ∈Z

1    α = 2 Sx a ((ν1 − k1 )fs ) Ts k1 ∈Z k2 ∈Z α∈Ax a

δ ((ν2 − k2 )fs − (−)(α − (ν1 − k1 )fs )) 1   α = Sx a ((ν1 − k1 )fs ) Ts α∈Ax a k1 ∈Z  δ ((ν2 − k2 ) − (−)(α/fs − (ν1 − k1 )))

(3.258)

k2 ∈Z

where, in the last equality, the scaling property of the Dirac delta fs δ(νfs ) = δ(ν) is used (Zemanian, 1987, Sec. 1.7).

Proof of Lemma 3.19 The proof of (3.190) in the stochastic approach is straightforward. In contrast, the proof of the FOT Counterpart of (3.190) is more elaborated. In the following, a formal derivation of the FOT Counterpart of (3.190) is reported. Moreover, the general case of continuous-time GACS signals (Chapter 12) is considered which include the ACS signals as special case. For continuous-time GACS signals, the set of cycle frequencies depends on the lag parameter. The FOT counterpart of (3.190) is   α} E{ {x(n + m) x (∗) (n)} = E{α} {xa (t + τ ) xa (t)} (3.259) t=nTs ,τ =mTs

where E{α} {·} is the continuous-time almost-periodic component extraction operator (see α } {·} is the discrete-time almost-periodic component extraction operator (2.22)) and E{ which is defined similarly to its continuous-time counterpart. (∗) The second-order lag product xa (t + τ ) xa (t) can be decomposed into the sum of its almost-periodic component and a residual term x a (t, τ ) not containing any finitestrength additive sine-wave component (Section 2.3.1.6). For a GACS signal it results 

xa (t + τ ) xa(∗) (t) = E{α} xa (t + τ ) xa(∗) (t) + x a (t, τ )  = Rxαa (τ ) ej 2παt + x a (t, τ ) (3.260) α∈Aτ

Chapter 3 • Almost-cyclostationary signal processing

127

(∗)

where subscript x a ≡ xa xa , Aτ  Ax a ,τ , and   x a (t, τ ) e−j 2παt ≡ 0 t

∀α ∈ R .

(3.261)

For any 0 <  < Ts we formally have N   1 αn x a (t, τ )t=nT ,τ =mT e−j 2π s s 2N + 1 n=−N



1 = 2N + 1 =

N Ts +

−N Ts −

x a (t, mTs )

+∞ 

1 (2N + 1)Ts p=−∞



+∞ 

α t/Ts δ(t − nTs ) e−j 2π dt

n=−∞ N Ts + −N Ts −

α −p)t/Ts x a (t, mTs ) e−j 2π( dt

(3.262)

where, in the second equality, the Poisson formula (C.8a) +∞ 1  j 2πpt/Ts δ(t − nTs ) = e Ts p=−∞ n=−∞ +∞ 

(3.263)

is accounted for. Taking the limit for N → ∞ in (3.262) and accounting for (3.261) one has N   1 αn x a (t, τ )t=nT ,τ =mT e−j 2π =0. s s N→∞ 2N + 1

lim

(3.264)

n=−N

That is, if the continuous-time residual term x a (t, τ ) does not contain any (continuoustime) finite-strength additive sine-wave component (see (3.261)), then also the sampled residual term x a (nTs , mTs ) does not contain any discrete-time finite-strength additive sine-wave component. Thus, α} {x(n + m) x (∗) (n)} E{     α} = E{ xa (t + τ ) xa(∗) (t) t=nTs ,τ =mTs     α} = E{ E{α} {xa (t + τ ) xa(∗) (t)} t=nTs ,τ =mTs

  { α} x a (t, τ )t=nT ,τ =mT +E s s     { α} α j 2παt  =E Rx a (τ ) e 

=

 α∈Aτ

α∈Aτ

  α j 2παt  Rx a (τ ) e 

t=nTs ,τ =mTs

t=nTs ,τ =mTs

(3.265)

128

Cyclostationary Processes and Time Series

where in the third equality (3.264) is accounted for and in the last equality the fact that the right-hand side of (3.265) is a discrete-time almost periodic function is exploited. Eq. (3.259) is then proved observing that the right-hand side of (3.265) is coincident with the right-hand side of (3.259). Note that uniformly sampling a continuous-time GACS signal leads to a discrete-time ACS signal (Section 12.4).

Proof of the FOT counterpart of (3.191) Accounting for (3.265), we have   α −j 2π αn  R x (m) = x(n + m) x(n) e n ,  αn = Rxαa (mTs ) ej 2παnTs e−j 2π α∈AmTs

=



Rxαa (mTs )

α∈AmTs

=





j 2π(αTs − α )n



n

e n . /0 1  δαTs − α +p p∈Z

( α −p)fs Rx a (mTs )

(3.266)

p∈Z

where in the last equality we use the fact that there is nonzero contribution only for those cycle frequencies α such that α = ( α − p)/Ts , p ∈ Z.

Proof of the FOT counterpart of (3.192) Using (3.266) we have α  (ν) = Sx



α −j 2πνm  R x (m) e

m∈Z

=

!  

" ( α −p)fs Rx a (mTs )

e−j 2πνm

m∈Z p∈Z

=

!  

" ( α −p)fs Rx a (mTs ) e−j 2πνfs mTs

p∈Z m∈Z

=

1   ( α −p)fs Sx a ((ν − q)fs ) Ts p∈Z q∈Z

where, in the last equality, the second Poisson’s summation formula (C.12) is used.

(3.267)

Chapter 3 • Almost-cyclostationary signal processing

129

Proof of Theorem 3.22 From (3.195a) and (3.195b), under Assumption 3.20 (xa (t) strictly band-limited), and accounting for (3.193a), the support of the replica with k1 = k2 = 0 in the aliasing formula (3.187) is given by   supp Sxαa (ν1 fs ) δ(ν2 − (−)(α/fs − ν1 )) ⊆ {(ν1 , ν2 ) ∈ R × R : |ν1 fs |  B, |α − ν1 fs |  B, ν2 − (−)(α/fs − ν1 ) = 0}

(3.268a)

= {(ν1 , ν2 ) ∈ R × R : |ν1 |  B/fs , |ν2 |  B/fs , α = (ν1 + (−)ν2 )fs }

(3.268b)

Replicas in (3.187) are separated by 1 in both ν1 and ν2 variables. Thus, from (3.268b) it follows that B/fs  1/2, that is, fs  2B

(3.269)

is a sufficient condition such that replicas do not overlap. Note that, however, (3.269) does α fs in not assure that the mappings f1 = ν1 fs and α =   1 α α  (ν1 ) = Sx a (f1 )f =ν f ,α= Sx α fs 1 1 s Ts

(3.270)

holds ∀ α ∈ [−1/2, 1/2) and ∀ν1 ∈ [−1/2, 1/2). For example, for (∗) present and α = −B, equality (3.270) holds only for ν1 ∈ [−1/2, 0] and not for ν1 ∈ [0, 1/2]. In fact, for ν1 ∈ [0, 1/2] the density of the replica with n1 = 0, n2 = 1 is present in the right-hand side of (3.270). In contrast, condition B/fs  1/4, that is, fs  4B

(3.271)

assures in (3.268b) |ν1 |  1/4,

|ν2 |  1/4,

|α/fs | = |ν1 + (−)ν2 |  |ν1 | + |ν2 |  1/2

(3.272)

α fs in (3.270) hold ∀ α ∈ [−1/2, 1/2) and and, consequently, the mappings f1 = ν1 fs and α =  ∀ν1 ∈ [−1/2, 1/2). Moreover, (3.196) holds. Finally, note that the effect of sampling at two times the Nyquist rate (see (3.271)) leads also to (3.197).

Proof of Theorem 3.25 From (3.205), it follows that the support of the rescaled cyclic spectrum of the sampled signal is such that    α (ν)|ν=f/fs , (3.273) supp  Sxx α =α/fs ⊆ D3 (α, f ) where D3 (α, f ) 

D2 (α − pfs , f − qfs ) p∈Z q∈Z

130

Cyclostationary Processes and Time Series

FIGURE 3.10 Reciprocal positions of replicas of D1 (α, f ) in D3 (α, f ) defined in (3.274). Only the replicas closest to D2 (α, f ) defined in (3.137) are represented. In this example, 2f0 − fs /2 = mfs , fs = 4W . Source: (Napolitano and ˘ Dogançay, 2018) Copyright of Elsevier.



=

D1 (α − pfs , f − f0 − qfs ) ∪ D1 (α − pfs , f + f0 − qfs )

p∈Z q∈Z

 ∪ D1 (α − 2f0 − pfs , f − f0 − qfs ) ∪ D1 (α + 2f0 − pfs , f + f0 − qfs ) (3.274)

with D1 (α, f ) and D2 (α, f ) defined in (3.132) and (3.137), respectively. In the following, sufficient conditions are derived such that the set D3 (α, f ) in (3.274) is constituted by the union of non-overlapping translations of the set D1 (α, f ). In (3.274), in order to avoid overlap along the α direction, it is sufficient that there exists an integer m  1 such that 2f0 − fs /2 = mfs and, moreover (Fig. 3.10), D1 (α − mfs , f − f0 )  D1 (α − 2f0 , f − f0 )  D1 (α − (m + 1)fs , f − f0 ) α

α

(3.275)

or, equivalently, D1 (α − mfs , f + f0 )  D1 (α + 2f0 , f + f0 )  D1 (α − (m + 1)fs , f + f0 ) . α

α

(3.276)

In (3.275) and (3.276), the notation D  D means that the domain D is contained in a α

horizontal strip of the (α, f ) plane with cycle frequencies less than those of the domain D . Consequently, accounting for (3.132), we must have fs  4W . Analogously, in order to avoid overlap along the f direction, it is again sufficient that there exists an integer m  1 such that 2f0 − fs /2 = mfs and, moreover (Fig. 3.10), D1 (α, f − mfs − f0 )  D1 (α, f − f0 )  D1 (α, f − (m + 1)fs − f0 ) f

f

(3.277)

Chapter 3 • Almost-cyclostationary signal processing

131

or, equivalently, D1 (α, f − mfs + f0 )  D1 (α, f + f0 )  D1 (α, f − (m + 1)fs + f0 ) . f

(3.278)

f

In (3.277) and (3.278), the notation D  D means that the domain D is contained in a f

vertical strip of the (α, f ) plane with spectral frequencies less than those of the domain D . Thus we have no aliasing in α and f provided that (Fig. 3.10)  2f0 − fs /2 = mfs (3.279) fs  4W which is equivalent to (3.206). If, for the purpose of cyclic spectral analysis, it is required that for each  α ∈ [−1/2, 1/2) only one replica (i.e., only one cyclic spectrum) in (3.205) gives nonzero contribution if ν ∈ [−1/2, 1/2), then, in all the above relations involving supports, the domain D1 (α, f ) must be replaced by (see (3.65a) and (3.65b) with B = W ) D1 (α, f )  {(α, f ) ∈ R × R : |f |  W, |α|  2W } .

(3.280)

In such a case, the condition to avoid the overlap of domains along the f direction remains unaltered. In contrast, in order to avoid the overlap of domains along the α direction, the inequality in (3.279) must be replaced by the more stringent condition fs  8W . Consequently, (3.206) must be replaced by (3.207).