Stationary Signal Models Versus Cyclostationary Signal Models

What happens when a cyclostationary time-series is treated as if it were stationary?

In this post let’s consider the difference between modeling a communication signal as stationary or as cyclostationary.

There are two contexts for this kind of issue. The first is when someone recognizes that a particular signal model is cyclostationary, and then takes some action to render it stationary (sometimes called ‘stationarizing the signal’). They then proceed with their analysis or algorithm development using the stationary signal model. The second context is when someone applies stationary-signal processing to a cyclostationary signal model, either without knowing that the signal is cyclostationary, or perhaps knowing but not caring.

At the center of this topic is the difference between the mathematical object known as a random process (or stochastic process) and the mathematical object that is a single infinite-time function (or signal or time-series).

A related paper is The Literature [R68], which discusses the pitfalls of applying tools meant for stationary signals to the samples of cyclostationary signals.

Stationarizing a Random Process Using Phase Randomization

Suppose we have our old friend the textbook pulse-amplitude modulated (PAM) signal given by

\displaystyle s(t) = \sum_{k=-\infty}^\infty a_k p(t - kT_0), \hfill (1)

where p(t) is the pulse function (for example, rectangular or square-root raised-cosine), k is the symbol index, and \{a_k\} is the symbol sequence, typically drawn from a finite alphabet such as \{+1, -1\}. Textbook signals employ independent and identically distributed symbols. Figure 1 shows some examples of pulse functions that we’ll encounter throughout this post.

Figure 1. Five pulse-shaping functions that could be used in digital QAM signals (1). Each pulse here corresponds to an intended symbol rate of 1/20 (normalized frequency units).

The signal (time-series) s(t) in (1) is cyclostationary, by our lights, as we’ve seen in many posts at the CSP Blog. For example, for the binary alphabet a_k \in \{+1, -1\},  s(t) is the complex envelope of an RF BPSK signal. When the alphabet is \{+1, -1, +i, -i\}, the signal is QPSK.

In conventional stochastic process theory, s(t) is conceived of as an ensemble of signals in some probabilistic sample space. For each possible sequence \{a_k\}_{k=-\infty}^\infty, there is a corresponding s_j(t), called a sample path (or just ‘signal’). When you calculate probabilistic parameters, such as the mean, moments, and cumulants, you do so by averaging the sample paths–which is called averaging over the ensemble. This is done by applying the expectation operator E[\cdot]. For s(t) in (1) above, averaging over the ensemble produces periodically time-varying results (cyclostationarity). Let’s look at the autocorrelation for the BPSK signal version, for example,

\displaystyle R_s(t+t_1, t+t_2) = E[s(t+t_1)s^*(t+t_2)], \hfill (2)

\displaystyle = \sum_{j=-\infty}^\infty \sum_{k=-\infty}^\infty E[a_j a^*_k] p(t-jT_0+t_1)p^*(t-kT_0+t_2) \hfill (3)

\displaystyle = \sum_{k=-\infty}^\infty (1) p(t-kT_0+t_1)p^*(t-kT_0 + t_2) \hfill (4)

which follows because E[a_j a^*_k] = 1 if j=k and is zero otherwise for BPSK, which itself is a consequence of our assumption that the symbols are independent and identically distributed random variables. Is R_s(t+t_1, t+t_2) periodic? Yes, with period T_0. Replace t with t + T_0 and rename the summation index, and you have the same quantity. If the autocorrelation was instead independent of time, and the mean was too, the process is referred to as stationary.

People have lots of tools to apply when they have a stationary process, so in the past many looked for ways to render a non-stationary process stationary (‘stationarize’) to enable the fruitful use of those tools. We can stationarize our PAM signal by adding a random variable to the pulse train. Let’s call it \tau_0 (because it relates to time). The new signal is given by

\displaystyle s(t) = \sum_{k=-\infty}^\infty a_k p(t - kT_0 + \tau_0), \hfill (5)

where the random variable \tau_0 has a uniform distribution on the interval [-T_0/2, T_0/2] for example. Is this signal stationary? To check, let’s compute the first- and second-order moments both with and without the phase-randomizing variable \tau_0.

Stationarizing the Mean

We start with the mean, or first-order moment, of s(t). Using the original signal, we have

\displaystyle E[s(t)] = E\left[ \sum_{k=-\infty}^\infty a_k p(t-k T_0) \right], \hfill (6)

 \displaystyle = \sum_{k=-\infty}^\infty E[a_k] p(t - kT_0) \hfill (7)

\displaystyle = \sum_{k=-\infty}^\infty M_a p(t-kT_0). \hfill (8)

So this mean value depends on time t and is also periodic for all p(t) and nontrivially periodic for some choices for p(t), assuming that M_a \neq 0. For most pulse-amplitude-modulated signals, it is indeed true that M_a = 0, but not for some, such as on-off keying (OOK), where the symbols are drawn from the set \{0, 1\}.  In any case, what do I mean by nontrivially periodic?

Let’s think about several different pulse functions p(t), as shown in Figure 1.

First, consider the highly textbook-ish case of rectangular p(t). Here p(t) is 1 for t\in [-T_0/2, T_0/2], and is zero otherwise. In such a case E[s(t)] = M_a for all t, and therefore the mean value is independent of time. But for a more realistic pulse like the half-cosine pulse, the sum of shifted pulses in (8) leads to a periodic function that is not simply a constant: it is nontrivially periodic. For T_0 = 20 and M_a = 1, we have the mean values (8) for our considered pulses in Figures 2-6.

Figure 2. Sum of shifted pulses from (8) with M_a = 1 and rectangular pulse function.

 

Figure 3. Sum of shifted pulses from (8) with M_a = 1 and square-root raised-cosine pulse function.

 

Figure 4. Sum of shifted pulses from (8) with M_a = 1 and half-cosine pulse function.

 

Figure 5. Sum of shifted pulses from (8) with M_a = 1 and Manchester pulse function.

 

Figure 6. Sum of shifted pulses from (8) with M_a = 1 and raised-cosine (Nyquist) pulse function.

For the practical square-root raised-cosine pulse, and the related raised-cosine (‘Nyquist’) pulse, it turns out that the sum of shifted pulses is very close to a constant, but for the half-cosine and Manchester pulses, the sum is clearly nontrivially periodic. They are all periodic, but some are either exactly a constant or very close to a constant, and those are trivially periodic.

The periodicities give rise to spectral lines, which we can visualize by taking the Fourier transform of the mean (8) and plotting its magnitude. For the rectangular and half-cosine pulses, we obtain the plots in Figures 7 and 8.

Figure 7. FFT magnitude for the sum of shifted pulses in (8) with M_a = 1 and rectangular pulse function. Note the lack of any peaks other than for frequency zero.

 

Figure 8. FFT magnitude for the sum of shifted pulses in (8) with M_a = 1 and half-cosine pulse function. Note the peaks at harmonics of the symbol rate 1/20.

We can force the probabilistic mean to be independent of time for any pulse using the phase-randomizing variable \tau_0. Let’s go through the math.

Our random process that includes a phase-randomizing random variable is now given by

\displaystyle s(t) = \sum_{k=-\infty}^\infty a_k p(t - kT_0 + \tau_0) \hfill (9)

where the probability density function for \tau_0 is p_{\tau_0}(x), which we will make specific later. The goal is to find the expected value of s(t), and then see how the distribution of \tau_0 (that is, p_{\tau_0}(x)) affects that expected value–for which distributions is the expected value of s(t) a constant function of time? Applying the expectation operator to the signal starts us off,

\displaystyle E[s(t)] = \sum_{k=-\infty}^\infty E[a_k]E[p(t-kT_0+\tau_0)] \hfill (10)

which follows by statistical independence between \{a_k\} and \tau_0 and the linearity of the expected value operator.

\displaystyle E[s(t)] = \sum_{k=-\infty}^\infty M_a E[p(t-kT_0+\tau_0)] \hfill (11)

\displaystyle = M_a E\left[ \sum_{k=-\infty}^\infty p(t - kT_0 + \tau_0) \right] \hfill (12)

For any value of \tau_0, the function within the backets in (12) is periodic in t with period T_0, so it can be exactly represented as a Fourier series

\displaystyle r(t) = \sum_{k=-\infty}^\infty p(t - kT_0 + \tau_0) = \sum_{j=-\infty}^\infty v_j e^{i 2 \pi j t /T_0} \hfill (13)

with Fourier coefficients

\displaystyle v_j = \frac{1}{T_0} \int_{\cal{I}} r(t) e^{-i 2 \pi j t/T_0} \, dt \hfill (14)

where \cal{I} represents some interval of length T_0,

\displaystyle {\cal{I}} = [t_1, t_2], \ \ \ t_2 - t_1 = T_0 \hfill (15)

Continuing with the Fourier coefficient,

\displaystyle v_j = \frac{1}{T_0} \int_{\cal{I}} \sum_{k=-\infty}^\infty p(t - kT_0 + \tau_0) e^{-i 2 \pi j t/T_0} \, dt \hfill (16)

\displaystyle = \frac{1}{T_0} \sum_{k=-\infty}^\infty \int_{\cal{I}} p(t - kT_0 + \tau_0) e^{-i 2 \pi j t/T_0} \, dt \hfill (17)

\displaystyle = \frac{1}{T_0} \sum_{k=-\infty}^\infty \int_{t_1}^{t_2} p(t - kT_0 + \tau_0) e^{-i 2 \pi j t/T_0} \, dt \hfill (18)

Using the substitution s = t + \tau_0 and then renaming s as t we obtain the following expression for the Fourier coefficient v_j

\displaystyle v_j = \frac{e^{i 2 \pi j \tau_0/T_0}}{T_0} \sum_k \int_{t_1+\tau_0}^{t_2 + \tau_0} p(t - kT_0) e^{-i 2 \pi j t /T_0} \, dt \hfill (19)

Next, make the substitution for t given by s = t - kT_0, then rename s by t again, yielding

\displaystyle v_j = \frac{e^{i 2 \pi j \tau_0/T_0}}{T_0} \sum_{k=-\infty}^\infty \int_{t_1+\tau_0 -kT_0}^{t_2+\tau_0-kT_0} p(t) e^{-i 2 \pi jt/T_0}\, dt \hfill(20)

\displaystyle = \frac{e^{i 2 \pi j \tau_0/T_0}}{T_0} \int_{-\infty}^\infty p(t) e^{-i 2 \pi j t /T_0}\, dt \hfill (21)

Let’s define a new symbol for the integral of the sine-wave-weighted pulse function,

\displaystyle P(j) = \int_{-\infty}^\infty p(t) e^{-i 2 \pi j t /T_0}\, dt \hfill (22)

which means we have an expression for the Fourier series representation of r(t),

\displaystyle r(t) = \sum_{j=-\infty}^\infty v_j e^{i 2 \pi j t/T_0} = \sum_{j=-\infty}^\infty \frac{e^{i 2 \pi j \tau_0/T_0}}{T_o} P(j) e^{i 2 \pi j t /T_0} \hfill (23)

With this expression we can return to the expected value of the signal,

\displaystyle E[s(t)] = M_a E[r(t)] \hfill (24)

\displaystyle = M_a E\left[ \sum_{j=-\infty}^\infty v_j e^{i 2 \pi j t/T_0} = \sum_{j=-\infty}^\infty \frac{e^{i 2 \pi j \tau_0/T_0}}{T_o} P(j) e^{i 2 \pi j t /T_0} \right] \hfill (25)

\displaystyle = \frac{M_a}{T_0} \sum_j E[e^{i 2 \pi j \tau_0/T_0}] P(j) e^{i 2 \pi j t /T_0} \hfill (26)

So, what is the expected value E[e^{i 2 \pi j \tau_0/T_0}]? Let the random variable \tau_0 have uniform distribution on the interval [0, T_0). Then we compute the expected value easily as

\displaystyle E[e^{i 2 \pi j \tau_0/T_0}] = \int_{0}^{T_0} \frac{1}{T_0} e^{i 2 \pi j x /T_0}\, dx \hfill (27)

\displaystyle = \frac{1}{T_0} \left. \frac{e^{i 2 \pi j \tau_0/T_0}}{i2\pi j/T_0}\right|_{\tau_0=0}^{\tau_0 = T_0} = 0 \hfill (j \neq 0) \hfill (28)

For j=0,

\displaystyle E[e^{i 2 \pi j \tau_0/T_0}] = E[e^0] = E[1] = 1. \hfill (29)

At long last, then, we have

\displaystyle E[s(t)] = \frac{M_a}{T_0} P(0) e^0 = \frac{M_a}{T_0} \int_{-\infty}^\infty p(t) \, dt. \hfill (30)

So, if the distribution of the phase-randomizing random variable \tau_0 is uniform on an interval with length equal to the symbol interval T_0, then the resulting random process has a mean value that is independent of time for all pulses p(t).

Let’s check this result before proceeding (and there is a lot more to proceed with). Keep it simple: Does this result make sense for a rectangular pulse p(t)? Such a pulse is equal to one on the interval [-T_0/2, T_0/2], so we can evaluate (30) by inspection: E[s(t)] = (M_a/T_0)\times T_0 = M_a. So the mean of the signal is the average value of the symbols, which makes good sense.

Let’s now look at the case in which the distribution of the phase-randomizing random variable \tau_0 is uniform on the interval [0, x) for x \leq T_0 and x > 0. So we go back to the general formula (26) and evaluate the expectation with this new distribution for \tau_0.

We want to evaluate \displaystyle E[e^{i 2 \pi j \tau_0/T_0}] for the more general uniform distribution. For j=0, the expected value is one, as before, because E[C] = C for any non-random constant C. For j\neq 0, we have

\displaystyle E[e^{i 2 \pi j \tau_0/T_0}] = \int_0^x (1/x) e^{i2 \pi j s/T_0} \, ds \hfill (31)

\displaystyle = \frac{1}{x} \left. \frac{e^{i 2 \pi j s/T_0}}{i2\pi j/T_0} \right|_{s=0}^{s=x} \hfill (32)

\displaystyle = \frac{T_0 e^{i \pi j x/T_0}}{i 2 \pi j x} (2i)\sin(\pi j x/T_0) \hfill (33)

(which checks because it is zero for x=T_0)

\displaystyle = e^{i \pi jx/T_0} \mbox{\rm sinc}(jx/T_0). \hfill (34)

Putting it all together, we have the mean value

\displaystyle E[s(t)] = \frac{M_a}{T_0}P(0) + \frac{M_a}{T_0} \sum_{j\neq 0} e^{i 2 \pi j t /T_0} P(j) e^{i \pi jx/T_0} \mbox{\rm sinc}(jx/T_0). \hfill (35)

A remaining check on the validity of (35) is that it matches our previous result when x = 0 (no phase randomization). I’ll leave that to you.

In conclusion, the mean is a non-constant function of time, generally, when the phase-randomizing random variable is uniform but on a smaller interval than the symbol interval T_0, and is a constant function of time if the phase-randomizing random variable is uniform on [0, T_0). Figure 9 shows an example featuring OOK (M_a = 0.5) and half-cosine pulses. As x increases to T_0 = 20, the time-variation of the mean vanishes. The final constant value (for x = T_0) is equal to M_a P(0) / T_0 = (0.5)(12.73)/20 = 0.32. The lack of symmetry is due to the lack of symmetry in the density function p_{\tau_0}(x).

Figure 9. Illustration of the effect of a phase-randomizing delay variable on the first-order moment (mean value) of the pulse-amplitude-modulated signal (1) for the case of half-cosine pulses, M_a = 0.5, and a symbol rate of 1/20. The mean becomes time-invariant when the distribution of the random delay is uniform on an interval of width 20.

Stationarizing the Second-Order Moment

Since most of the communication signals we observe, capture, and process have zero mean values, what we really want to understand is how phase randomization is used to stationarize the second-order moment–the autocorrelation function. That’s where the cyclostationarity property that we make such good use of here at the CSP Blog meets its end. So let’s do the math for stationarizing the second-order moment here. The form is similar to what we went through for the first-order moment (the mean value) above, so I’ll be a little less detailed.

Our signal is (9), which is repeated here for convenience

\displaystyle s(t) = \sum_{k=-\infty}^\infty a_k p(t - kT_0 + \tau_0). \hfill (36)

The non-conjugate (conventional) autocorrelation function, defined within the random-process framework, is given by

\displaystyle R_s(t, \tau) = E\left[s(t+\tau/2)s^*(t-\tau/2)\right] \hfill (37)

Assuming the pulse function itself is non-random, and that the symbols \{a_k\} are statistically independent of the phase-randomizing random variable \tau_0 (assumptions we also made in the case of the mean value analysis above), we have

\displaystyle R_s(t,\tau) = \sum_{k_1,k_2=-\infty}^\infty E[a_{k_1} a_{k_2}^*] E\left[p(t+\tau/2 -k_1T_0 + \tau_0)p^*(t-\tau/2-k_2 T_0 + \tau_0)\right]. \hfill (38)

Assuming, as before, independent and identically distributed symbols \{a_k\} and E[a_k] = 0 (which rules out OOK here, which is OK, OK?), the symbol expectation simplifies to

\displaystyle E[a_{k_1}a_{k_2}^*] = \left\{ \begin{array}{ll} 0, & k_1 \neq k_2, \\ \sigma_a^2, & k_1 = k_2 \end{array} \right. \hfill (39)

\displaystyle R_s(t,\tau) = \sigma_a^2 \sum_{k=-\infty}^\infty E\left[p(t+kT_0+\tau/2+\tau_0)p^*(t+kT_0-\tau/2+\tau_0)\right] \hfill (40)

\displaystyle \sigma_a^2 E\left[\sum_{k=\infty}^\infty p(t+kT_0+\tau/2+\tau_0)p^*(t+kT_0-\tau/2+\tau_0)\right] \hfill (41)

The sum over k in the brackets in (41) is a periodic function of time with period T_0.  Unlike the case of the mean value, the periodicity here depends on the autocorrelation lag variable \tau so that, for example, the rectangular-pulse signal can lead to a nontrivial periodicity. Let’s show some plots of the sum over k in (41) for \tau = T_0/4 and \tau_0 = 0 to fix the idea that this autocorrelation is really periodic (and so therefore we must ruin that periodicity using phase randomization to get back to stationarity). That is, for non-random \tau_0 in (41), the expectation of the sum is just the sum, and so the autocorrelation is just the scaled sum. Examples are shown in Figures 10-12.

Figure 10. The sum of products of shifted pulses from (41) for the rectangular pulse and a delay of T0/4.

 

Figure 11. The sum of products of shifted pulses from (41) with the half-cosine pulse and a delay of T0/4.

 

Figure 12. The sum of products of shifted pulses from (41) with the square-root raised-cosine pulse and a delay of T0/4.

The periodicity of these autocorrelation functions is evident by inspection of the graphs in Figures 10-12, and a simple Fourier analysis shows the exact period to be T_0 = 20, as shown in Figures 13-15. These figures are really showing the cycle frequencies for the non-conjugate cyclic autocorrelations, which are, recall, the Fourier coefficients of the Fourier-series decomposition of the time-varying moment.

Figure 13. The FFT of the sum of products of shifted pulses from (41) with the rectangular pulse and a delay of T0/4. Note the clear presence of harmonics of the symbol rate 1/20.

 

Figure 14. The FFT of the sum of products of shifted pulses from (41) with the half-cosine pulse and a delay of T0/4. Note the clear presence of harmonics of the symbol rate 1/20.

 

Figure 15. The FFT of the sum of products of shifted pulses from (41) with the square-root raised-cosine pulse and a delay of T0/4. Note the clear presence of the symbol rate 1/20.

The analysis approach to evaluating the expectation in (41) is the same as we used for the mean-value analysis: represent the sum, which is periodic with period T_0, as a Fourier series, then apply the expectation.

\displaystyle w(t) = \sum_{k=\infty}^\infty  p(t+k T_0+\tau/2+\tau_0) p^*(t+k T_0-\tau/2+\tau_0) \hfill (42)

\displaystyle = w(t + T_0) \ \ \forall t \hfill (43)

The jthe Fourier coefficient for the Fourier-series representation of w(t) is

\displaystyle w_j = \frac{1}{T_0} \int_{t_1}^{t_2} w(t) e^{-i2 \pi j t/T_0} \, dt \hfill (44)

After following steps similar to those for the mean value above, the Fourier coefficient is found to be

\displaystyle w_j = \frac{1}{T_0} e^{i 2 \pi j \tau_0/T_0} \sum_k \int_{t_1+\tau_0-kT_0}^{t_1 + T_0 + \tau_0 - kT_0} p(s + \tau/2) p^*(t-\tau/2) e^{-i 2 \pi j s /T_0} \, ds \hfill (45)

\displaystyle = \frac{1}{T_0} e^{i 2 \pi j \tau_0/T_0} \int_{-\infty}^\infty p(s+\tau/2)p^*(s-\tau/2) e^{-i 2 \pi j s /T_0} \, ds \hfill (46)

\displaystyle = \frac{1}{T_0} e^{i 2 \pi j \tau_0/T_0} Q(j) \hfill (47)

So the Fourier series is given by

\displaystyle w(t) = \sum_{j=-\infty}^\infty w_j e^{i 2 \pi j t /T_0} \hfill (48)

\displaystyle = \sum_{j=-\infty}^\infty \left( \frac{1}{T_0} e^{i2 \pi j \tau_0/T_0} Q(j) \right) e^{i 2 \pi j t /T_0} \hfill (49)

The autocorrelation is

\displaystyle R_s(t, \tau) = \sigma_a^2 E[w(t)] \hfill (50)

\displaystyle = \frac{\sigma_a^2}{T_0} \sum_{j=-\infty}^\infty E\left[ e^{i 2 \pi j \tau_0/T_0} \right] Q(j) e^{i 2 \pi j t /T_0} \hfill (51)

So we’re back to the question: What is \displaystyle E\left[ e^{i 2 \pi j \tau_0/T_0} \right]? The answer hasn’t changed, fortunately for us. Using a uniform distribution for \tau_0 with width x, as before, we obtain

\displaystyle E\left[ e^{i 2 \pi j \tau_0/T_0} \right] = e^{i \pi j x/T_0} \mbox{\rm sinc}(jx/T_0) \hfill (52)

Finally, then, the autocorrelation of the phase-randomized signal is given by

\displaystyle R_s(t, \tau) = \frac{\sigma_a^2}{T_0} \sum_{j=-\infty}^\infty e^{i \pi j x/T_0} \mbox{\rm sinc}(jx/T_0) Q(j) e^{i 2 \pi j t /T_0} \hfill (53)

In the case of x = T_0, only the j = 0 term in the sum in (53) is non-zero (that term is equal to one), so the autocorrelation is not a function of time:

\displaystyle R_s(t, \tau) = \frac{\sigma_a^2}{T_0} Q(0) \hfill (54)

\displaystyle = \frac{\sigma_a^2}{T_0} \int_{-\infty}^\infty p(t+\tau/2) p^*(t-\tau/2) \, dt = R_s(\tau) \hfill (55)

Otherwise, the autocorrelation is dependent on time t.

Let’s try to check this result. We already know the PSD for the signal is simply

\displaystyle S_s^0(f) = \frac{\sigma_a^2}{T_0} \left| P(f) \right|^2 \hfill (56)

where P(f) \Leftrightarrow p(t). Does our result (55) provide the same answer? Recalling that the PSD is the Fourier transform of the autocorrelation,

\displaystyle S_s^0(f) = {\cal{F}} \left[ R_s^0(\tau) \right] = {\cal{F}} \left[ R_s(\tau) \right] \hfill (57)

\displaystyle = {\cal{F}} \left[ \frac{\sigma_a^2}{T_0} \int_{-\infty}^\infty p(t+\tau/2) p^*(t-\tau/2) \, dt \right] \hfill (58)

\displaystyle = {\cal{F}} \left[ \frac{\sigma_a^2}{T_0} \int_{-\infty}^\infty p(r) p^*(r - \tau) \, dr \right] \hfill (59)

\displaystyle = \frac{\sigma_a^2}{T_0} \int_{-\infty}^\infty p(r) \left[ \int_{-\infty}^\infty p^*(r-\tau) e^{-i 2 \pi f \tau} \, d\tau \right] \, dr \hfill (60)

\displaystyle = \frac{\sigma_a^2}{T_0} \int_{-\infty}^\infty p(r) \left[ \int_{-\infty}^\infty p^*(u) e^{-i 2 \pi f (r-u)} \, du \right] \, dr \hfill (61)

\displaystyle = \frac{\sigma_a^2}{T_0} \int_{-\infty}^\infty p(r) e^{-i2\pi f r} \left[ \int_{-\infty}^\infty p^*(u) e^{i 2 \pi f u} \, du \right] \, dr \hfill (62)

\displaystyle = \frac{\sigma_a^2}{T_0} P(f) \left[ \int_{-\infty}^\infty p(u) e^{-i 2 \pi f u} \, du \right]^*\hfill (63)

\displaystyle = \frac{\sigma_a^2}{T_0} \left| P(f) \right|^2 \hfill (64)

which checks.

All we’ve done so far is show that the mean and non-conjugate autocorrelation for a baseband complex PAM signal are time-invariant if you use a phase-randomizing random variable \tau_0 with a uniform distribution. Let’s get a little more general and connect this basic result to cyclostationary signal processing and to more realistic signal models.

Connection to the non-Conjugate and conjugate autocorrelation functions of CSP

Suppose our signal model includes effects from imperfect downconversion (synchronization):

\displaystyle s(t) = \sum_{k=-\infty}^\infty a_k p(t - kT_0 + \tau_0) e^{-i 2 \pi f_c t + i\phi_0} \hfill (65)

Here \tau_0 is the symbol-clock phase, f_c is the carrier frequency offset, and \phi_0 is the carrier phase. When the only random variable in (65) is the symbol a_k, then we know the non-conjugate and conjugate cyclic autocorrelations, the spectral correlation functions, and the higher-order cyclic cumulants. We know that the cycle frequencies exhibited by the signal (65) follow the form

\displaystyle \alpha(n, m, k) = (n - 2m)f_c + k/T_0 \hfill (66)

Now suppose that \tau_0 is a random variable, as before, but that f_c and \phi_0 are non-random (unknown constants). Have we stationarized the non-conjugate and conjugate autocorrelation functions?

\displaystyle s(t+\tau/2)s^*(t-\tau/2) = \sum_{k_1,k_2} a_{k_1}a_{k_2}^* p(t + \tau/2 - k_1 T_0 + \tau_0)

\displaystyle \times p^*(t -\tau/2 - k_2 T_0 + \tau_0) e^{i2 \pi f_c \tau} \hfill (67)

So this is the same kind of expression we already stationarized with \tau_0 except for the factor e^{i2 \pi f_c \tau}, which is non-random by assumption, so if \tau_0 is a uniform random variable on [0, T_0), the non-conjugate autocorrelation (expectation of (67)) will be stationarized.

The conjugate autocorrelation is the expectation of a different delay product

\displaystyle s(t+\tau/2)s(t-\tau/2) = \sum_{k_1,k_2} a_{k_1}a_{k_2} p(t + \tau/2 - k_1 T_0 + \tau_0)

\displaystyle \times p(t -\tau/2 - k_2 T_0 + \tau_0) e^{i2 \pi 2 f_c t + i2\phi_0} \hfill (68)

Let’s let \phi_0 be random. Assuming the symbol-clock-phase variable \tau_0 and the carrier phase variable \phi_0 are independent, we confront three expectations when we apply the expectation to the delay product

\displaystyle R_{s^*}(t, \tau) = E\left[e^{i2 \pi 2f_c t + i2\phi_0}\right] S_a E\left[\sum_{k} p(t + \tau/2 - kT_0 + \tau_0)p(t -\tau/2 - kT_0 + \tau_0) \right] \hfill (69)

where S_a = E[a_k^2] (which is often, but not always, zero). The expectation over \tau_0 will render the third factor time-invariant, as we’ve painstakingly shown already. But we now have to confront the time-variant term corresponding to the expectation over \phi_0.

\displaystyle E\left[e^{i2 \pi 2 f_c t + i2\phi_0}\right] = E\left[e^{ i2\phi_0}\right] e^{i2 \pi 2f_c t} \hfill (70)

The delay variable \tau_0 is not quite enough to render the entire conjugate autocorrelation time-invariant–we need to consider the carrier phase random variable.

What is a condition on \phi_0 so that E[e^{i2\phi_0}] = 0? Assume \phi_0 is uniformly distributed on the interval [0, x] (sound familiar?) Then

\displaystyle E[e^{i2\phi_0}] = \int_0^x (\frac{1}{x}) e^{i2y} \, dy \hfill (71)

\displaystyle = \frac{e^{ix}}{x} \sin(x) \hfill (72)

Now \sin(x) = 0 for x = b\pi where b is an integer. So we can pick b = 1 to yield the uniform distribution on [0, \pi). And with that choice, the conjugate autocorrelation is stationarized.

Stationarizing Higher-Order Moments and Cumulants

We see that if we use two phase-randomizing random variables, the symbol-clock phase \tau_0 and the carrier phase \phi_0, we can render the time-varying non-conjugate and conjugate autocorrelation function time invariant. This also means that the entire non-conjugate spectral correlation function reduces to the PSD. Staying within that random-process framework, then, we can use any methods or results applicable to stationary random processes on our signal.

But what about the higher-order statistics? What about the cyclic moments and cyclic cumulants? (What’s your guess?) 

The time-varying nth-order moment function is the expectation

\displaystyle R_s(t, \tau; n, m) = E \left[ \prod_{j=1}^n s^{(*)_j} (t + \tau_j) \right] \hfill (73)

\displaystyle = E \left[ \prod_{j=1}^n \left( \sum_{k_j} a_{k_j} p(t-k_jT_0 +\tau_j + \tau_0) e^{i2\pi f_c(t+\tau_j)+i\phi_0} \right)^{(*)_j} \right] \hfill (74)

\displaystyle = \sum_{\mathbf{k}} E[A_{\mathbf{k}}] E \left[ \prod_{j=1}^n p^{(*)_j} (t - k_j T_0 + \tau_j + \tau_0) \right] E \left[ \prod_{j=1}^n \left( e^{i 2 \pi f_c (t+\tau_j) + i \phi_0} \right)^{(*)_j} \right] \hfill (75)

\displaystyle = \sum_{\mathbf{k}} E[A_{\mathbf{k}}] E \left[ \prod_{j=1}^n p^{(*)_j} (t - k_j T_0 + \tau_j + \tau_0) \right] E \left[ e^{i 2 \pi (n-2m)f_c t} e^{i(n-2m)\phi_0} e^{i 2\pi f_c(\sum(-)_j \tau_j)} \right] \hfill (76)

If n=2m, the time-dependency related to f_c disappears, as does any dependency on \phi_0 We’re left with

\displaystyle = \sum_{\mathbf{k}} E[A_{\mathbf{k}}] E \left[ \prod_{j=1}^n   p^{(*)_j}(t-k_j T_0 +\tau_j + \tau_0) \right] e^{i2\pi f_c(\sum (-)_j\tau_j)}. \hfill (77)

We can find all the conditions on \mathbf{k} = [k_1\ k_2\ \ldots k_n] such that the expected value over the symbols is non-zero. For each, we can evaluate the expectation over \tau_0 as before, and this will render each term time-invariant.

If n \neq 2m, we’ll always have the factor e^{i(n-2m)\phi_0}, and we can choose the carrier-phase random variable to force this expectation equal to zero.

In short, the same phase-randomizing random variables we used in the second-order (autocorrelation) case will also force all the higher-order temporal moment functions to be either zero or time-invariant.

Another Way to Stationarize a Cyclostationary Communication Signal: Assume Perfect Knowledge of Synchronization ParameterS and Sample it

Sometimes you’ll see in published papers a signal model that corresponds to application of various signal-processing operations on an actual received signal–the received signal is preprocessed before modulation-recognition begins, for example. If we assume we know, or can estimate, the particular value of the symbol-clock phase parameter \tau_0, the carrier frequency offset f_c, and the carrier phase \phi_0, as well as the pulse-shaping function p(t), then we can process the signal to remove the effect of the carrier and carrier phase (multiply by e^{-i2\pi f_c t - i \phi_0}), apply a perfect matched filter, and sample at the optimal sampling instants. This leads to a signal model such as

\displaystyle s(j) = \sum_{k=-\infty}^\infty (a_k + n_k). \hfill (78)

In this way, we obtain a noisy sequence of symbols, which is typically stationary (it might not be, though, if the symbols contain periodically repeated sequences or some other deviation from independent and identically distributed symbols).

I find this kind of move strange because it assumes high SNR  and very long data records (else the estimates of the synchronization parameters cannot be perfect). It also neglects channel effects.

Consequences of Stationarization

When we adopt a phase-randomized random-process model, what we get is simplicity of mathematical analysis–probabilistic parameters like the autocorrelation are time-invariant and so are that much easier to calculate. What we give up is realism. The key important point is that we have stationarized the process, not the sample paths of the process. In other words, we’ve ensured we don’t have ergodicity, which is a property of a random process that says that an average over the ensemble is equal to an average over infinite time for almost every sample path. If we construct estimators using the stationarized random process, and we have appropriate ergodicity, we can apply the estimator to a sample path (the actual simulated or captured signal) and expect good results. If we don’t have ergodicity, then we may not get the expected or good results, and we might then embark on a program to, for instance, robustify the algorithm against nuisance parameters like carrier frequency offset.

Each of the sample paths of (65) is like a single received example of the process–it has a fixed-but-unknown symbol-clock phase \tau_0, a fixed-but-unknown carrier frequency offset f_c, and a fixed-but-unknown carrier phase \phi_0. It is a cyclostationary signal with properties (for example, cycle frequencies) that depend on the nature of the symbol variables \{a_k\} and the pulse function p(t).

When we do our work in the random-process domain, using a model that contains phase-randomizing variables, we see that there aren’t any cyclic cumulants. The only cumulants are ‘stationary-signal cumulants,’ which involve combining lower-order moments in the usual way, but there aren’t any cycle frequencies except zero. Then we switch to Monte Carlo simulations or processing a captured signal, and each time we process, there are cycle frequencies, there are cyclic cumulants, there is time-variation to moments and cumulants. The question is: What happens to those stationary-signal parameters when we apply them to a cyclostationary signal? I attempt to shed some light on this question in the remainder of this post.

Interlude: The Case of OOK Versus BPSK

Thinking about what the stationary-signal cumulants mean when applied to a cyclostationary input is rather difficult when the cumulant orders n involved are larger than two (at least it is for me). So to ease into the discussion about higher-order stationary-signal cumulants versus higher-order cyclic cumulants, let’s look at an example that only involves first- and second-order cumulants: BPSK vs OOK.

OOK and BPSK are PAM signals like (9) and (65), and if you keep the pulse function p(t) the same between them, the only difference is in the values of the symbols \{a_k\}. The finite set of distinct complex numbers (amplitude and phase of the value multiplying a pulse) is called the constellation. BPSK has a constellation of \{-1, 1\} whereas OOK has a constellation of \{0, 1\}. This means OOK is BPSK plus a constant. The implication is that OOK is BPSK with a finite-strength additive sine-wave component with frequency equal to the carrier frequency. Mathematically,

\displaystyle s_{OOK} (t) = \sum_k a_k p(t - k T_0 + \tau_0)e^{i2 \pi f_c t + i\phi_0}  \hfill (79)

\displaystyle = \sum_k 0.5(b_k + 1) p(t - k T_0 + \tau_0)e^{i2 \pi f_c t + i\phi_0} \hfill (80)

where a_k \in \{0, 1\} and b_k \in \{-1, 1\}. This means that for pulse functions like rectangular, square-root raised-cosine, and raised-cosine, the difference between OOK and BPSK is that OOK has a PSD with an impulse at the carrier f_c and BPSK does not.

We are going to look at cyclic cumulants that use the correct lower-order cycle frequencies and those that don’t to illustrate the kind of thing that can happen when one applies stationary-signal cumulants to cyclostationary signals. When we ignore the first-order cycle frequency that OOK possesses (\alpha = f_c), and compute the cumulant for cycle frequency of zero, we’ll see that the stationary-signal cumulants diverge from the cyclic cumulants–only one of them is a true cumulant. The more common and important case of non-OOK signals and larger values of n is more subtle, but we’ll attempt a similar illustration at the close of the post.

The average power of the OOK signal is made twice that of the BPSK signal so that if the OOK sine-wave component were removed from the signal prior to any kind of estimation, the signal would appear statistically identical to a BPSK signal. This power boost makes the plots a bit easier to interpret–when the OOK and BPSK parameters match, the tone has no effect, or else it was properly taken care of by inclusion of the first-order cycle frequency in the cyclic-cumulant formula.

First, let’s look at the PSDs of the OOK and BPSK signals to fix the ideas here. Figure 16 shows PSD estimates for the two generated signals. The frequency-smoothing method (FSM) is used to estimate these PSDs, so the OOK tone appears as a rectangle with width equal to the FSM smoothing-window width I chose.

Figure 16. Measured power spectra for a BPSK signal and an OOK signal with twice the power of the BPSK signal. Note the spectral line at the carrier frequency of 0.05 for the OOK signal.

Figure 17 shows the corresponding autocorrelation functions (non-conjugate cyclic autocorrelation for \alpha = 0).

Figure 17. Autocorrelation function estimates for the OOK and BPSK signals shown in Figure 16. The spectral line of the OOK signal manifests as the never-ending oscillation in the autocorrelation function.

Apart from the additive tone at the carrier of f_c = 0.05, the PSDs are identical. That additive tone component of the PSD corresponds to the sine-wave component of the autocorrelation function for OOK in Figure 17 (the OOK autocorrelation never decays to zero.)  The symbol rate is 1/10 and the pulse function is square-root raised-cosine with roll-off factor of 0.5.

The non-conjugate spectral correlation function for the symbol-rate cycle frequency and the conjugate spectral correlation function for the doubled-carrier cycle frequency are shown in Figures 18 and 19, respectively. Note here that we might expect that theses functions will differ between OOK and BPSK, because the spectral correlation functions don’t include any notion of lower-order combinations.

Figure 18. Estimates of the non-conjugate spectral correlation function for the OOK and BPSK signals shown in Figure 16 using the symbol-rate cycle frequency. 

 

Figure 19. Estimates of the conjugate spectral correlation function for the OOK and BPSK signals shown in Figure 16 using the doubled-carrier cycle frequency. 

The corresponding cyclic autocorrelation functions are shown in Figures 20 and 21.

Figure 20. Estimates of the non-conjugate cyclic autocorrelation function for the OOK and BPSK signals shown in Figure 16 using the symbol-rate cycle frequency.

 

Figure 21. Estimates of the conjugate cyclic autocorrelation function for the OOK and BPSK signals shown in Figure 16 using the doubled-carrier cycle frequency.

Note that the symbol-rate functions (Figures 18 and 20) for OOK and BPSK match–there is no effect of the first-order cycle frequency \alpha = f_c for this cycle-frequency. That is because there is no combination of first-order cycle frequencies that add up to the symbol rate, so the first-order statistics of the signal do not impact those estimates. On the other hand, the doubled-carrier functions are strongly impacted by the presence of the first-order sine wave. That is because each term in the delay product s(t+\tau/2) s(t-\tau/2) contains a sine-wave component with frequency f_c, and these multiply to contribute to the overall statistic for \alpha = 2_fc. (A real-world example of this phenomenon is the case of ATSC-DTV, which has a pilot tone and two strong conjugate features–the doubled tone frequency ends up contributing to one of those conjugate features.)

Now let’s look at the true cyclic cumulants. These cumulants take into account the first-order cycle frequency of OOK (\alpha = f_c). Figures 22-24 show the true cyclic cumulants for both signals–the BPSK and OOK estimates should match every time.

Figure 22. Estimates of the cyclic cumulants for the OOK and BPSK signals in Figure 16 for n=2, m=1, and a cycle frequency of zero. Note the match between the two signals. This is due to the proper use of the OOK signal’s lower-order cyclostationarity.

 

Figure 23. Estimates of the cyclic cumulants for the OOK and BPSK signals in Figure 16 for n=2, m=1, and a cycle frequency equal to the symbol rate. Note the match between the two signals. The match here is independent of the proper use of the OOK signal’s lower-order cyclostationarity because it does not impact this particular second-order parameter.

 

Figure 24. Estimates of the cyclic cumulants for the OOK and BPSK signals in Figure 16 for n=2, m=0, and a cycle frequency equal to the doubled carrier. Note the match between the two signals. This is due to the proper use of the OOK signal’s lower-order cyclostationarity.

Next let’s look at the cyclic cumulants we obtain when we ignore the first-order cycle frequency. The three sets of cumulants are shown in Figures 25-27.

Figure 25. Estimates of the cyclic cumulants for the OOK and BPSK signals in Figure 16 for n=2, m=1, and a cycle frequency equal to zero. Because we are ignoring the OOK signal’s lower-order cyclostationarity, the parameters for the two signals diverge.

 

Figure 26. Estimates of the cyclic cumulants for the OOK and BPSK signals in Figure 16 for n=2, m=1, and a cycle frequency equal to the symbol rate. The match between the two signals is independent of proper use of the OOK signal’s lower-order cyclostationarity because that lower-order cyclostationarity does not impact this particular parameter. 

 

Figure 27. Estimates of the cyclic cumulants for the OOK and BPSK signals in Figure 16 for n=2, m=, and a cycle frequency equal to the doubled-carrier. Because we are ignoring the OOK signal’s lower-order cyclostationarity, the parameters for the two signals diverge.

As expected, the cumulants diverge for \alpha = 0 and \alpha = 2f_c.

This means that even if you focus exclusively on cycle frequencies equal to zero (‘stationary-signal moments and cumulants’), the ignored-but-still-there lower-order cyclostationarity can affect your result. This is the problem with modeling things as stationary random processes (likely with the aid of phase-randomizing random variables) and then applying stationary-signal moments and cumulants to sample paths of those processes, which are not stationary.

Let’s look at this phenomenon in a little more detail before wrapping this post up.

Applying Stationary-Signal Probabilistic-Parameter Definitions to Cyclostationary Signals

So here is where it all pays off. What happens when you assume you have a stationary signal, or you’ve actively tried to create a stationary process, and you then apply stationary-signal moment and cumulant estimators to the signal, but it is actually cyclostationary?

Let’s consider a BPSK signal with square-root raised-cosine pulses with roll-off factor of one, symbol rate of 1/10, and various carrier frequency offsets ranging from zero to a moderate fraction of the symbol rate such as 1/100. The signal has unit power and noise with power 0.1 is added for realism.

We compute and plot the true cyclic cumulants and the stationary-signal cumulants side-by-side. The true cyclic cumulants are simply cyclic cumulants that employ the correct lower-order cycle frequencies in the required combinations of lower-order cyclic moments. The stationary-signal cumulants are obtained by following the moment-to-cumulant formula, but the only cycle frequency used for any combination of order n and number of conjugations m is zero. The true cyclic cumulants and the stationary-signal cumulants will match for a stationary ergodic random process, otherwise they will diverge.

First, consider (n,m,k) = (2, 0, 0). That is, order n=2, number of conjugated factors m =0, and harmonic number k = 0. The harmonic number k will always be zero in these measurements because we are comparing stationary-signal cumulants (always a cycle frequency of zero) with true cyclic cumulants (can correspond to other cycle frequencies too).

Figure 28. Cyclic cumulant estimates for proper use of lower-order cyclostationarity (upper) and for assuming all moments have a single cycle frequency of zero (lower). The true cyclic cumulants use the doubled-carrier as the cycle frequency, whereas the stationary-signal cyclic cumulants are constrained to use zero for all cycle frequencies. Note the two sets of parameters match only for the case in which the carrier frequency offset is actually zero. The theoretical peak of the true cyclic cumulant is 1.0.

The correct peak of the true cyclic cumulants in Figure 28 is 1.0, and we see that the true cyclic cumulants are correct for all values of the CFO–they are using the doubled CFO as the cycle frequency. So when the actual CFO is zero, the true and stationary cumulants match, as expected. When the actual CFO is not zero, the (2, 0, 0) stationary-signal cumulant is not equal to the cumulant when the CFO is zero. There is no (2, 0) cumulant with cycle frequency zero for those cases. The (2, 0, 0) cumulant is useless unless either (1) you know the CFO, or (2) the CFO is zero.

For the (2, 1, 0) cumulants, we obtain the estimates in Figure 29.

Figure 29. Cyclic cumulant estimates for proper use of lower-order cyclostationarity (upper) and for assuming all moments have a single cycle frequency of zero (lower). The true cyclic cumulants use the zero as the cycle frequency here, whereas the stationary-signal cyclic cumulants are constrained to use zero for all cycle frequencies in all cases. The parameters match because the cycle frequencies match, and there is no lower-order cyclostationarity to deal with for n=2. The theoretical peak of the true cyclic cumulant is 1.1.

In Figure 29, the stationary-signal and true cyclic cumulants match independent of the actual value of the carrier frequency offset. The value of the offset is irrelevant to these cumulants, so lack of knowledge of the offset doesn’t affect the computations.

Next, let’s look at (4, 0, 0) in Figure 30.

Figure 30. Cyclic cumulant estimates for proper use of lower-order cyclostationarity (upper) and for assuming all moments have a single cycle frequency of zero (lower). The true cyclic cumulants use the quadrupled carrier as the cycle frequency here, whereas the stationary-signal cyclic cumulants are constrained to use zero for all cycle frequencies in all cases. The parameters do not match in any case, including when the quadrupled carrier is equal to zero, because the stationary-signal cumulants do not appropriately use the signal’s lower-order cyclostationarity. The theoretical peak of the true cyclic cumulant is 2.28.

For (4, 0, 0), none of the stationary-signal cumulants match the true cyclic cumulants, and the true cyclic-cumulant magnitudes are correct since the theoretical peak for this pulse shape and modulation type is 2.28. Here is where things get interesting. Even when the carrier offset is zero, the stationary-signal cumulant does not match the true cyclic cumulant because the appropriate lower-order cyclic moments are not combined in the moment-to-cumulant formula. This didn’t happen in the second-order case (2, 0, 0) because there isn’t any first-order cyclostationarity to account for (that’s why I did the OOK-BPSK Interlude above, where there was lower-order (first-order) cyclostationarity).

Even worse are the cases in Figure 30 where the offset is non-zero. In those cases the measured stationary-signal cumulants are near zero because there is no zero-valued cycle frequency for (4, 0) for those signals, and the stationary-signal cumulants pretend that there is.

The (n,m) = (4,2) cumulants suffer as well because even though the cycle frequencies for BPSK and (4,2) are just harmonics of the symbol rate–no offset enters–the lower-order cyclic moments that are called for in the moment-to-cumulant formula include cycle frequencies related to the doubled carrier offset. So, we see in Figure 31 that the (4,2) stationary-signal cumulants do not match the true cyclic cumulants.

Figure 31. Cyclic cumulant estimates for proper use of lower-order cyclostationarity (upper) and for assuming all moments have a single cycle frequency of zero (lower). The true cyclic cumulants use zero as the cycle frequency here, whereas the stationary-signal cyclic cumulants are constrained to use zero for all cycle frequencies in all cases. The parameters do not match in any case, including when the carrier offset is zero, because the stationary-signal cumulants do not appropriately use the signal’s lower-order cyclostationarity. The theoretical peak of the true cyclic cumulant is 2.28.

The story is similar for (6, m). I’ll just show the results for m=3 in Figure 32.

Figure 32. Cyclic cumulant estimates for proper use of lower-order cyclostationarity (upper) and for assuming all moments have a single cycle frequency of zero (lower). The true cyclic cumulants use zero as the cycle frequency here, whereas the stationary-signal cyclic cumulants are constrained to use zero for all cycle frequencies in all cases. The parameters do not match in any case, including when the carrier offset is zero, because the stationary-signal cumulants do not appropriately use the signal’s lower-order cyclostationarity. The theoretical peak of the true cyclic cumulant is 23.7.

Discussion

In future posts I’ll be looking at some published papers on modulation recognition (both statistics-based and machine-learning-based) that use stationary-signal moments and cumulants. In this post, the idea I’ve tried to get across is that the kinds of processing you might want to apply to extract valuable information from a signal depend on the statistical nature of the data–the captured or simulated signal that you are actually processing. One can get confused and create poor feature extractors if the features are based on a random-process model that lacks ergodic properties because those properties are needed to forge a strong link between the mathematical model and the processed data record.

Another way of saying it is that if you use phase-randomization to render a random process stationary, then you will get unexpected results when you process the sample paths of that stationary process because they will be cyclostationary signals.

Thanks for getting all the way through this long post! As usual, if there are errors or if you want to make comments, leave a message in the Comments section below.

Author: Chad Spooner

I'm a signal processing researcher specializing in cyclostationary signal processing (CSP) for communication signals. I hope to use this blog to help others with their cyclo-projects and to learn more about how CSP is being used and extended worldwide.

Leave a Comment, Ask a Question, or Point out an Error