# A Gallery of Cyclic Cumulants

The third in a series of posts on visualizing the multidimensional functions characterizing the fundamental statistics of communication signals.

Let’s continue our progression of galleries showing plots of the statistics of communication signals. So far we have provided a gallery of spectral correlation surfaces and a gallery of cyclic autocorrelation surfaces. Here we introduce a gallery of cyclic-cumulant matrices.

When we look at the spectral correlation or cyclic autocorrelation surfaces for a variety of communication signal types, we learn that the cycle-frequency patterns exhibited by modulated signals are many and varied, and we get a feeling for how those variations look (see also the Desultory CSP posts). Nevertheless, there are large equivalence classes in terms of spectral correlation. That simply means that a large number of distinct modulation types map to the exact same second-order statistics, and therefore to the exact same spectral correlation and cyclic autocorrelation surfaces. The gallery of cyclic cumulants will reveal, in an easy-to-view way, that many of these equivalence classes are removed once we consider, jointly, both second- and higher-order statistics.

### Review of Cyclic Temporal Cumulants

We develop cyclic cumulants from random-process theory in the cyclic temporal cumulants post. Let’s review the basic formula for an $n$th-order cyclic cumulant before turning to the main topic of this post, which is the visualization of cyclic cumulants.

A cyclic cumulant is expressible in terms of cyclic moments, and both cyclic cumulants and cyclic moments are intimately related to the various $n$th-order joint probability density functions of a cyclostationary random process or time-series. We can sidestep the density functions by just focusing on the moments and cumulants.

The most familiar (to signal processors) moment is the second-order moment $\displaystyle R_x(t,\tau) = E[x(t+\tau/2)x^*(t-\tau/2)], \hfill (1)$

for some complex-valued signal $x(t)$. For a stationary signal, this moment is not a function of $t$, and it is usually renamed to the autocorrelation function. For a cyclostationary signal, this moment is periodic (or almost periodic) in $t$, and is expressible as a generalized Fourier series, the Fourier frequencies of which are the cycle frequencies, and the Fourier amplitudes are the cyclic autocorrelation functions.

To generalize the second-order moment to arbitrary orders $n$, we consider a slightly more general version of (1), $\displaystyle R_x(t, [\tau_1\ \tau_2]) = E[x(t+\tau_1)x^*(t+\tau_2)], \hfill (2)$

then we also realize that the moment changes if we delete the conjugation on the second factor or introduce one to the first factor. So we consider $m$ conjugations by using an optional conjugation notation $(*)_j$ as in $\displaystyle R_x(t, [\tau_1, \tau_2]; m) = E[x^{(*)_1} (t+\tau_1) x^{(*)_2}(t+\tau_2)]. \hfill (3)$

Now we can more easily generalize this second-order moment to orders $n$ by introducing an $n$-vector of delays $\boldsymbol{\tau} = [\tau_1\ \tau_2\ \ldots\ \tau_n]$ as in $\displaystyle R_x(t, \boldsymbol{\tau};n,m) = E\left[ \prod_{j=1}^n x^{(*)_j}(t+\tau_j) \right]. \hfill (4)$

Here the variable $m$ denotes the number of optional conjugations that are selected–it must range from zero to $n$. However, as I’ve noted elsewhere, this notation is not satisfactory because it does not specify exactly which factors in the product are conjugated. To remedy this, in the post on symmetries of higher-order cyclic functions, I introduce a binary vector $\boldsymbol{b} = [b_1\ b_2\ \ldots\ b_n]$ that has elements equal to one for the indices of $x^{(*)_j}(t)$ that are conjugated and elements equal to zero for those indices corresponding to unconjugated factors. So for $n=2$ we get the non-conjugate cyclic autocorrelation with $\boldsymbol{b} = [0\ 1]$, the conjugated non-conjugate cyclic autocorrelation with $\boldsymbol{b} = [1\ 0]$, the conjugate cyclic autocorrelation with $\boldsymbol{b} = [0\ 0]$, and the conjugated conjugate cyclic autocorrelation with $\boldsymbol{b} = [1\ 1]$.

Now for cyclostationary communication signals $x(t)$, the moments (4) are periodic for some (usually all) choices of even $n$, and so we can express those moments in Fourier series $\displaystyle R_x(t, \boldsymbol{\tau};n,m) = \sum_\alpha R_x^\alpha(\boldsymbol{\tau}; n,m) e^{i 2 \pi \alpha t}, \hfill (5)$

or if you like with the added specification of $\boldsymbol{b}$, $\displaystyle R_x(t, \boldsymbol{\tau};n,m, \boldsymbol{b}) = \sum_\alpha R_x^\alpha(\boldsymbol{\tau}; n,m,\boldsymbol{b}) e^{i 2 \pi \alpha t}. \hfill (6)$

The cyclic cumulant is a specific nonlinear function of cyclic moments. This follows directly from the well-known relationship between the $n$th-order cumulant of a set of random variables and all the lower-order joint moments of subsets of those random variables (see the cyclic cumulant post for details). The cyclic cumulant is given by the following expression involving a sum of products of cyclic moments, where each product of cyclic moments has a set of cycle frequencies that sum to the cyclic-cumulant cycle frequency, $\displaystyle C_x^\alpha (\boldsymbol{\tau}; n,m, \boldsymbol{b}) = \sum_{P} \left[ (-1)^{p-1}(p-1)! \sum_{\boldsymbol{\beta}^\dagger \boldsymbol{1}=\alpha} \prod_{j=1}^p R_x^{\beta_j} (\boldsymbol{\tau}_j; n_j, m_j, \boldsymbol{b}_j) \right]. \hfill (7)$

In this post, then, what we want to do is compute and display, somehow, the multidimensional function $\displaystyle C_x^\alpha(\boldsymbol{\tau};n,m,\boldsymbol{b})$ for a wide variety of communication signals for which these cumulants are not zero. And that is potentially helpful to us visually oriented humans because these cyclic cumulant functions are highly useful as features in signal-processing-based and machine-learning-based modulation-recognition algorithms and systems (My Papers [17,25,26,28,30,38,43,50,51,52,54,55]).

### Visualizing Cyclic Cumulants

So how should we visualize these complicated multidimensional functions? One way I’ve tried in the past is to plot two-dimensional slices of the function and arrange them in a sequence indexed by the values of a third dimension, then show the sequence in a video. I do this in the cyclic-cumulants and higher-order symmetries posts. That works OK for fourth-order cumulants, but the dimensionality of sixth- and higher-order functions would lead to some very long videos indeed. In the fourth-order case, I can set $\tau_4 = 0$, leading to the reduced-dimension cyclic temporal cumulant, fix $\tau_3$, the cycle frequency, and the conjugation configuration $\boldsymbol{b}$, and then make plots of the cyclic-cumulant magnitude as a function of $\tau_1$ and $\tau_2$. But for sixth-order cumulants, I’ve got too many delays to deal with for that approach to be viable.

Another possibility is to simply arrange all the cyclic cumulants in a giant matrix and plot that matrix using, say MATLAB’s imagesc.m. Let’s look at that for rectangular-pulse BPSK. Here the sampling rate is unity (as usual), the bit rate is 1/8, and the carrier offset is zero. I use 32,768 samples to estimate each cyclic cumulant (32,768/8 = 4,096 bits), delays ranging from -8 to 8, and cycle frequencies $\{k/8\}_{k=-7}^7$. Note that this set of cycle frequencies is sufficient to cover all BPSK cycle frequencies $(n-2m)f_c + kf_{bit}$ because $f_c = 0$.

I’ll consider cyclic-cumulant orders of 2 and 4 separately. The basic reason I don’t plot them all on one set of axes is that the dimensionalities of the two cyclic cumulant functions differ due to the delay vector. For $N=2$, it is simply $\boldsymbol{\tau} = [\tau_1]$ and for $N=4$ it is $\boldsymbol{\tau} = [\tau_1\ \tau_2\ \tau_3]$.

#### All N=2 BPSK Cyclic Cumulants

Let’s first look at the cyclic cumulants for the ideal BPSK signal, and then we’ll take a quick look at the same functions for that BPSK signal after it experiences a simple multipath-channel impulse response.

All of the non-zero second-order ( $N=2$) cyclic cumulant magnitudes and phases are shown in Figures 1 and 2, respectively. To check my work, I plot a column of the matrix in Figure 1 in Figure 3. I chose the column corresponding to $\boldsymbol{b} = [0\ 1]$ and $\alpha = 0$, which is column 23. This column corresponds to the conventional autocorrelation, and we know that the autocorrelation for our rectangular-pulse BPSK signal is a triangle with width equal to twice the symbol interval, or 16. The height should be just greater than one, since the BPSK signal has unit power and there is a small amount of noise present. And from Figure 3 this all checks out. Figure 1. All non-zero second-order cyclic-cumulant magnitudes for the ideal rectangular-pulse BPSK signal with rate 1/8 and carrier 0.0. The y-axis shows the value of the single delay $\latex \tau_1$ and the x-axis is an index into the cycle-frequency and conjugation-configuration vector . For example, column eight corresponds to and . Kinda looks like a stained-glass window. Figure 2. The phases of the second-order cyclic cumulants for all non-zero second-order cyclic-cumulant magnitudes for the ideal rectangular-pulse BPSK signal with rate 1/8 and carrier 0.0. Figure 3. A sanity check: This is a plot of the 23rd column of the magnitude matrix in Figure 1. It should be equal to the conventional autocorrelation function for a rectangular-pulse BPSK signal with symbol interval equal to 8, unit power, and high SNR. Which it is.

#### All N=4 BPSK Cyclic Cumulants

Turning to fourth-order, the same three kinds of plots that we see in Figures 1-3 are shown for $N=4$ in Figures 4-6. Now you’ve literally ‘seen’ a higher-order cyclic cumulant! (I’ve posted the jpg images and the corresponding MATLAB .fig files on the Downloads page.) To check my work, I plot the 53rd column of the magnitude matrix in Figure 6 (see caption). Figure 4. All non-zero fourth-order cyclic-cumulant magnitudes for the ideal rectangular-pulse BPSK signal with rate 1/8 and carrier 0.0. The y-axis shows the index of the three-vector latex \boldsymbol{b}\$. Figure 5. The phases of the fourth-order cyclic cumulants for all non-zero fourth-order cyclic-cumulant magnitudes for the ideal rectangular-pulse BPSK signal with rate 1/8 and carrier 0.0. Figure 6. A sanity check: This is a plot of the 53rd column of the magnitude matrix in Figure 4. The x-axis is an index into the three-vector . Near the middle of the plot, there is an index that corresponds to , and column 53 corresponds to a conjugation configuration of . Therefore, the peak of this plot should be , which we know should be equal to 2.

#### All Filtered-BPSK Cyclic Cumulants

The second- and fourth-order cyclic-cumulant functions for the ideal rectangular-pulse BPSK signal are estimated for a filtered version of the signal and shown in Figures 7-12. The filter is a four-ray multipath channel with delays 0.0, 2.0, 3.3, and 5.5 samples and associated complex-valued gains 1.0, -0.4+0.25i, 0.2-0.2i, and -0.1-0.07i.

So we can see that this method of visualization has a drawback–too much information is plotted. We can’t make much sense out of the fourth-order plots, so you can imagine that we’d have lots of trouble with the sixth- and eighth-order plots.

Fortunately there is another way.

### The Gallery

The important, at least to me, aspect of the multidimensional cyclic-cumulant function is what I call the cycle-frequency pattern. We know that cycle frequencies depend on the order $n$, the number of conjugated factors $m$, the signal’s probability structure (for example, the moments and cumulants of the symbol random variable), the pulse-shaping function, and even the particular values of the delay vector $\boldsymbol{\tau}$ (for example, consider rectangular-pulse signals and delays all equal to zero). For the broad class of digital PSK and QAM signals (and SQPSK and CPM/CPFSK), the cycle-frequency pattern can be discerned from a low-dimensional cut or slice of the full cyclic temporal cumulant function. In particular, we can set $\boldsymbol{\tau} = [0\ 0\ \ldots 0]$ and not render any cycle frequencies invisible (unlike for rectangular-pulse signals). So the cycle-frequency pattern is then some function $\displaystyle \alpha(x(t), n, \boldsymbol{b}). \hfill (8)$

And for almost all of these types of signals, that function is $\displaystyle \alpha(x(t), n, \boldsymbol{b}) = (n-2m)f_c + kf_{sym}. \hfill (9)$

However, the cyclic-cumulant magnitude might be zero for some or even many of these cycle frequencies–the set of the cycle-frequencies corresponding to all of the non-zero cyclic-cumulant magnitudes is the cycle-frequency pattern.

So an alternative visualization approach is to set the delay vector equal to the origin and find the cyclic-cumulant magnitudes and phases for all the cycle frequencies in (9), up to some maximum desired order, and plot those values. This gives rise to the style of cyclic-cumulant visualizations that form the banner of the CSP Blog–they are what I consider the fingerprint of a communication signal (I was using that term before the Machine Learners got their hands on it, sigh).

Of course there are complications even with this drastically dimension-reduced view of the cyclic cumulants. The main one is that typically the cyclic-cumulant magnitudes grow rapidly with $n$, and so the dynamic range of the plot, so to speak, is high. The larger values (for larger $n$) will tend to dominate the plot. So the first thing we can do is raise each cyclic cumulant magnitude to the power $2/n$. This puts all the cyclic cumulant magnitudes on a equal footing with the power ( $(n,m,k) = (2, 1, 0)$). The second thing we can do is plot the base-10 logarithm of that warped cyclic-cumulant magnitude. These steps render the cycle-frequency pattern easy to visually grasp–and that’s why we’re here after all.

Putting it all together results in the gallery of cyclic-cumulant features shown in Video 1.

The constellation definitions are shown in Figure 13. Figure 13. Constellation names and plots. Use these names to identify which cyclic cumulant plot might be of interest to you.

Regarding the CPM and CPFSK feature plots, CPFSK is a CPM signal with a rectangular pulse. All other CPM signals are called CPM, and the pulse type is indicated by a component of the name string. CPM signals can be partial-response signals or full-response signals, depending on whether the frequency pulse extends over more than one symbol interval (partial) or not (full). The underlying pulse-amplitude-modulated signal (that forms part of the argument of the carrier cosine) can have binary or larger alphabets. The alphabet size is also embedded in the filename. So, for example, the feature with title CPMLRC-0.3-3 is a CPM signal with raised-cosine pulses (‘LRC’), modulation index of 0.3, and the pulse spans three symbol intervals. 