This is a long post. To make it worthwhile, I’ve placed some movies of cyclic-cumulant estimates at the end. Or just skip to the end now if you’re impatient!
In my work on cyclostationary signal processing (CSP), the most useful tools are those for estimating second-order statistics, such as the cyclic autocorrelation, spectral correlation function, and spectral coherence function. However, as we discussed in the post on Textbook Signals, there are some situations (perhaps only academic; see my question in the Textbook post) for which higher-order cyclostationarity is required. In particular, a probabilistic approach to blind modulation recognition for ideal (textbook) digital QAM, PSK, and CPM requires higher-order cyclostationarity because such signals have similar or identical spectral correlation functions and PSDs. (Other high-SNR non-probabilistic approaches can still work, such as blind constellation extraction.)
Recall that in the post introducing higher-order cyclostationarity, I mentioned that one encounters a bit of a puzzle when attempting to generalize experience with second-order cyclostationarity to higher orders. This is the puzzle of pure sine waves (My Papers ). Let’s look at pure and impure sine waves, and see how they lead to the probabilistic parameters widely known as cyclic cumulants.
Sine Waves Times Sine Waves Equals Sine Waves
Suppose we have a single complex-valued sine wave given by
where is the amplitude, is the frequency, and is the phase. If we consider a quadratic operation on this signal to produce the signal , we will find another sine wave,
That is, this particular quadratic operation leads to a sine wave with frequency phase , and complex amplitude . If we had included a conjugation on the second factor, , the result is a trivial sine wave with zero frequency, amplitude , and phase depending on and .
Next consider a signal consisting of the sum of two sine waves with different frequencies,
The simple (no-conjugate) second-order operation leads to
In other words, the quadratic transformation to leads to sine-wave components with frequencies of , , and . Other quadratic operations, such as those with a single conjugated factor or two conjugated factors, lead to different sine-wave frequencies in , but they all produce at least one sine-wave component.
Impure Sine Waves
Now recall that our working definition of cyclostationarity is that at least one non-trivial sine wave can be generated by subjecting the signal of interest to a quadratic nonlinearity, even when the signal itself contains no such sine-wave component. So how does the previous discussion of `sine waves begetting sine waves’ apply?
A natural extension of second-order theory to higher orders is to simply increase the order of the moment that we considered when defining the cyclic autocorrelation function. The usual second-order lag product
can be generalized to
and then we could consider the third-order or fourth-order versions
It turns out that the odd-order moments are zero for almost all communication signals, so going forward we’ll consider only even-order moments, such as the moment corresponding to the fourth-order lag product. That is, we consider the expected value of the lag product,
We know that this fourth-order moment contains sine-wave components because the second-order lag product contains sine waves, the second-order lag product contains sine waves, and sine waves multiplied by sine waves results in sine waves. In particular, suppose possesses a second-order non-conjugate cycle frequency when the lags are and . It also must possess the non-conjugate second-order cycle frequency of zero (it possesses power), so that the sine wave with frequency in multiplies the sine wave with frequency zero in , producing a sine wave with frequency in the fourth-order lag product (for suitable values of and , such as zero).
This kind of analysis can be performed for other choices of conjugations and lag-vector values . The conclusion is that the sine-wave components in a higher-order lag product are due, at least in part, to the interactions between lower-order sine-wave components, because `sine waves times sine waves equals sine waves.’ So is there a component of any of the fourth-order sine waves that is NOT due to these lower-order interactions?
We’ll call the sine waves in the lag product (the moment sine waves) impure sine waves, because they will often contain products of lower-order components. Any component of the impure sine waves that does not arise from the interactions of the lower-order sine waves will be called a pure sine wave.
Purifying Sine Waves
To purify the impure sine waves in some th-order lag product, we’ll need to identify the impure culprits and subtract them. Let’s start with a simple case where the signal itself does, in fact, contain additive sine-wave components. Then the second-order lag product must contain sine waves, and we can set about purifying the second-order lag product, which simplifies our notation. Then we’ll extend the ideas to higher-order lag products.
Suppose our signal consists of our old friend the rectangular-pulse BPSK signal plus an additive sine-wave component
Recall that the BPSK signal has a symbol rate of and a carrier frequency of (normalized-frequency parameters). If we consider the non-conjugate cycle frequencies of , we see that the BPSK signal contributes the bit-rate harmonics , and that the two sine sine waves contribute . For the conjugate cycle frequencies, the BPSK signal contributes
and the sine waves contribute . Therefore, for both the non-conjugate and conjugate cases, the BPSK and sine wave share cycle frequencies. However, only the BPSK signal produces second-order sine-wave components that are not the result of products of first-order sine waves.
One way to proceed with the purification in this case is to just subtract off the sine waves from and then compute the cyclic autocorrelation or spectral correlation function of the result. But that method doesn’t generalize well to the case where there are no first-order sine waves. So let’s take a more systematic approach.
The lower-order (first-order) sine wave frequencies are , and the sine-wave created frequencies in the second-order lag product are . Considering the non-conjugate second-order cycle frequency of 0, we see that it has a component due to the BPSK signal and a component due to the tones. If the lag product in question is
then the sine-wave induced component has complex amplitude
This is just the sum of the complex amplitudes for each sine wave multiplied by itself. So if we extract the zero-frequency sine wave from the second-order lag product, and then subtract from it the complex amplitude of each unique product of first-order sine waves with resultant frequency zero, we’ll be left with the complex amplitude due solely to the BPSK signal. The lower-order product terms will be removed, and we’ll have purified the second-order sine wave:
The idea is to subtract all unique products of the lower-order sine waves such that the sum of the sine-wave frequencies for each product is equal to the specified second-order cycle frequency. Here we are targeting the second-order cycle frequency of , so we subtract the product of first-order sine waves with frequencies (in ) and (in ), and also the product with frequencies (in ) and (in ).
Suppose the signal actually contained sine-wave components with frequencies . Let the second-order non-conjugate cycle frequency of interest be . Then we seek to subtract off lower-order products for each pair of frequencies and such that
In our example above, the frequency set is and . Then we have the two sets of pairs and . For , we have the single pair , which makes sense because there is only one product of sine waves that has frequency , and that is the sine wave with frequency in the unconjugated (first) factor, and in the conjugated factor.
Partitions of the Index Set
A key phrase in the development above is “unique products of lower-order sine waves.” How do we characterize such unique products? For the general case of an th-order moment, we are thinking about the sine-wave components of the th-order lag product
where denotes an optional conjugation, and there are factors that are conjugated. We want to account for every possible product of lower-order sine waves whose frequencies sum to the cycle frequency of interest . So there could be a sine wave in that multiplies a sine wave in and whose sum is , and maybe a sine wave in that multiplies a sine wave in , or possibly any other partitioning of the involved factors .
So we have to consider all the ways of partitioning the set of factors, and then for each one of these partitions, we have to determine whether or not the lag products for the partitions contain sine-wave components that sum to the desired cycle frequency . Whew. That is indeed complicated.
But the purification idea leads to the idea of partitions of a set of elements, in this case the elements . For simplicity, we can associate these elements with the index set . A partition of is a collection of proper subsets whose union is and whose pairwise intersections are all the empty set. This is just what we need to characterize all the possible ways lower-order sine waves can combine to contribute to the lag-product sine wave.
For , we have and the partitions are
and there are no others (no unique partitions, anyway, because is not unique). For , the partitions are
The number of partitions rapidly increases with the index size . (The number of partitions is given by Bell’s number, which makes use of Stirling numbers of the second kind.) For example, there are partitions for ,
Fortunately, any partitions containing any element that has an odd size, such as a singleton or a three-element set can usually be ignored in the computation of a pure sine wave because the signal under study has no first-order sine-wave components nor any non-zero odd-order moments. So, typically for , the partitions that matter are only
The partitions of can be easily computed recursively.
Getting back to purifying sine waves, we have arrived at the idea that to find the pure sine-wave component of some th-order lag product, we need to subtract from that lag product all products of lower-order sine-wave components that have frequencies that sum to the desired frequency, and this can be done by using the partitions of the index set. To ensure that we don’t end up subtracting off the same sine wave component multiple times, we wish to subtract from the lag product all products of lower-order pure sine-wave components.
We can get all the pure th order sine waves at once by subtracting all the products of pure lower-order sine waves from the sum of all th-order sine waves. That is, let represent the general sine-wave extraction operator,
where the sum is over all frequencies for which the signal contains a sine-wave component with non-zero amplitude,
In other words, extracts all the impure sine waves from its argument.
Our purified sine waves are then given by
where the sum is over all th-order partitions of the index set , and each partition has size with elements having sizes and numbers of optional conjugations .
The Shiryaev-Leonov Moment-Cumulant Formula
In the theory of random variables and processes, the th-order moment function is defined as the expected value of the product of variables. For simplicity, let’s consider random variables for . Then the th-order moment is
The moment also appears in the series expansion of the characteristic function,
which we won’t get into here. However, the cumulants appear in the series expansion of the logarithm of the characteristic function. With the moments, we can directly compute the value by time or ensemble averaging in practice, but how do we compute cumulants (leaving aside for now why we would even want to, but see this post if you can’t help yourself)? Fortunately for us, the relationship between the th-order moment and the th-order cumulant was derived long ago, by Shiryaev and Leonov (My Papers , Downloads).
The th-order cumulant is given by a nonlinear combination of lower-order moments,
where the sum is over all unique partitions of the index set , each having elements . The notation indicates the elements . The reverse formula, that gives the moment in terms of the lower-order cumulants is
which should look familiar. It is the same formula as the pure-sines-waves formula. To see this more clearly, break up the sum over the partitions into one term for and one term for all the other partitions,
The th-order temporal moment function is simply
and the th-order cyclic temporal moment function is a Fourier coefficient of the (poly)periodic moment,
The th-order temporal cumulant function is given by the now-familiar combination of products of lower-order temporal moment functions,
Now the th-order cyclic temporal cumulant function (cyclic cumulant) is just a Fourier coefficient of the temporal cumulant function
The vector of cycle frequencies is the vector of cyclic temporal moment cycle frequencies, and they must sum to the cyclic-cumulant cycle frequency The sum over the cycle-frequency vector must include all for which each of the corresponding cyclic moment functions is non-zero.
We can also express the time-varying moment function as the sum of products of time-varying cumulant functions, using (26):
Rearranging (34), we obtain a sometimes-useful expression relating the moment to the cumulant and the sum of products of all strictly lower-order cumulants,
So that’s it. The cyclic cumulants are really pure sine-wave amplitudes, and they can be computed by combining cyclic-moment amplitudes, which are easy to estimate directly from data. And since cyclic cumulants derive from cumulants, they inherent all the well-known properties of cumulants, such as tolerance to Gaussian contamination. Most importantly for CSP, ‘cumulants accumulate,’ that is, the cumulant for the sum of independent cyclostationary signals is the sum of the cumulants for the individual signals. This property does not hold in general for the moments.
As far as I know, the pure-sine-waves derivation is the only mathematical problem whose solution is a cumulant. All other uses of the cumulant that I know about use it like a found tool. But I could very well be wrong about this, and if so, please leave a comment below.
An Example: Rectangular-Pulse BPSK
The following two videos show estimated fourth-order cyclic temporal cumulant magnitudes for rectangular-pulse BPSK with fixed at zero, and the video frames corresponding to varying . The first video is the cyclic cumulant for and the second is for the BPSK bit rate . In both videos, two conjugations are used .
An Example: SRRC BPSK
The fourth-order cyclic cumulant analysis is repeated for BPSK with square-root raised-cosine pulses in the following two videos.