# CSP Estimators: Cyclic Temporal Moments and Cumulants

In this post we discuss ways of estimating $n$-th order cyclic temporal moment and cumulant functions. Recall that for $n=2$, cyclic moments and cyclic cumulants are usually identical. They differ when the signal contains one or more finite-strength additive sine-wave components. In the common case when such components are absent (as in our recurring numerical example involving rectangular-pulse BPSK), they are equal and they are also equal to the conventional cyclic autocorrelation function provided the delay vector is chosen appropriately.

The more interesting case is when the order $n$ is greater than $2$. Most communication signal models possess odd-order moments and cumulants that are identically zero, so the first non-trivial order $n$ greater than $2$ is $4$. So our estimation task is to estimate $n$-th order temporal moment and cumulant functions for $n \ge 4$ using a sampled-data record of length $T$.

As is often the case in cyclostationary signal processing, there are two quite different situations in which one wishes to perform an estimate: blind and non-blind. In particular, we either have prior knowledge of the cycle frequencies exhibited by the data (non-blind) or we do not (blind). Some people also consider the case where we have partial prior knowledge of cycle frequencies, such as knowing the symbol rate but not the carrier offset frequency, and they term this partially blind processing.

### Signals and Scene Used in this Post

We focus on BPSK with independent and identically distributed bits and a square-root raised-cosine pulse function with roll-off of $1$ (excess bandwidth of $100$%). In other words, textbook BPSK.

There are two scenes. The first one is just a single BPSK signal in noise; we’ll call it BPSK1. The symbol rate is $1/10$ (normalized frequencies are used throughout this post), the carrier offset frequency is $0.01$, and the total signal power is $1$ ($0$ dB). The second scene adds another BPSK signal to the first signal; let’s call it BPSK2. The bit rate, carrier offset frequency, and power level for this interfering BPSK signal are $1/11$, $0.002$, and $1.6 (2$ dB $).$

Estimated PSDs for these two scenes are shown below to set the stage:

### Non-Blind Estimation of Temporal Moments and Cumulants

This kind of estimation is a straightforward application of the formulas that define the temporal moments, temporal cumulants, cyclic temporal moments, and cyclic temporal cumulants. We can simply implement the finite-time versions of the formulas that define these quantities. In particular, since we know all the cycle frequencies that are exhibited by the data under study, we can estimate any $n$th-order cyclic moment function by using a discrete Fourier transform.

Recall that the theoretical $n$th-order cyclic moment function is given by

$\displaystyle R_x^\alpha (\boldsymbol{\tau}; n,m) = \lim_{T\rightarrow\infty} \frac{1}{T} \int_{-T/2}^{T/2} \prod_{j=1}^n x^{(*)_j} (t + \tau_j) e^{-i 2 \pi \alpha t} \, dt \hfill (1)$

So a simple non-parametric estimator for the cyclic moment function is just a non-limiting version of this formula,

$\displaystyle \hat{R}_x^\alpha (\boldsymbol{\tau}; n,m) = \frac{1}{T} \int_{-T/2}^{T/2} \prod_{j=1}^n x^{(*)_j} (t + \tau_j) e^{-i 2 \pi \alpha t} \, dt \hfill (2)$

which is easy to implement in discrete time.

The $n$th-order temporal moment function is simply the Fourier series involving all the non-zero (corresponding to the choices made for $\boldsymbol{\tau}$, $n$, and $m$) $n$th-order cyclic temporal moment functions,

$\displaystyle \hat{R}_x(t, \boldsymbol{\tau}; n,m) = \sum_\alpha \hat{R}_x^\alpha (\boldsymbol{\tau}; n,m) e^{i 2 \pi \alpha t} \hfill (3)$

Given that we know all the impure $n$th-order cycle frequencies $\alpha$, we can estimate the time-varying temporal moment function. Suppose we know the impure cycle frequencies for all orders $1, 2, \ldots, n$ and all choices of the number of conjugated factors $m$. Then we can estimate the temporal moment functions for all orders less than or equal to $n$. This set of temporal moment functions can then be combined using the Shiryaev-Leonov moment-to-cumulant formula to estimate the temporal cumulant function

$\displaystyle \hat{C}_x(t, \boldsymbol{\tau}; n,m) = \sum_{P} (-1)^{p-1} (p-1)! \prod_{j=1}^p \hat{R}_x(t, \boldsymbol{\tau}_j; n_j, m_j) \hfill (4)$

where the set $P$ contains all unique partitions of the index set $\{1, 2, \ldots, n\}$, $\{ \nu_j \}_{j=1}^p$ (see the post on cyclic cumulants for far too many details).

Now since we are still talking about non-blind estimation, we also know all pure cycle frequencies for all the relevant orders, so we can estimate the cyclic temporal cumulant function (cyclic cumulant) for pure cycle frequency $\alpha$ by extracting the Fourier-series coefficient for $\alpha$ from $\displaystyle \hat{C}_x(t, \boldsymbol{\tau}; n,m).$

Alternatively, we can use the cyclic-moment-to-cyclic-cumulant formula directly, estimating each cyclic moment separately and combining them according to

$\displaystyle \hat{C}_x^\alpha (\boldsymbol{\tau}; n,m) = \sum_{P} \left[ (-1)^{p-1}(p-1)! \sum_{\boldsymbol{\beta}^\dagger \boldsymbol{1}=\alpha} \prod_{j=1}^p \hat{R}_x^{\beta_j} (\boldsymbol{\tau}_j; n_j, m_j) \right]. \hfill (5)$

#### Non-Blind Estimation for a Single BPSK Signal in Noise

Using the method (5) of tediously combining all the proper estimated lower-order cyclic moments on the single-BPSK signal scene yields the following result:

Here we are plotting the warped cyclic-cumulant magnitudes. Each $n$th-order cyclic cumulant is raised to the power $2/n$ (independently of $m$). This warping is discussed in the first post on modulation recognition. In the present context, it renders the plots easier to grasp when multiple cyclic-cumulant orders are shown because the dynamic range of the cyclic cumulants can be quite large–this dynamic range is reduced through the warping. The theoretical cyclic cumulants are relatively easy to compute from our closed-form formula for PSK/QAM:

The error in the estimate is small for this simple case. The error is defined as the absolute value of the difference between the $n$th-order estimate and the $n$th-order theoretical cyclic cumulant divided by the maximum theoretical cyclic-cumulant magnitude for all $(n, m)$ cyclic cumulants for the signal. This measure tends to de-emphasize the errors in the $(n,m,k)$ cyclic cumulants for the larger $k$, which have small values relative to $(n,m,0)$ cyclic cumulants. Here $k$ refers to the harmonic number in the general PSK/QAM cycle-frequency formula $\alpha(n,m,k) = (n-2m)f_c + k/T_0$, where $f_c$ is the carrier frequency and $1/T_0$ is the symbol rate.

#### Non-Blind Estimation for a Two Cochannel BPSK Signals in Noise

Now let’s repeat the non-blind cyclic-cumulant estimation for BPSK1 + BPSK2 + noise:

Notice here that almost all the cyclic-cumulant magnitudes are as expected from the theoretical cyclic cumulants for BPSK, except for those corresponding to $(n, m, k) = (n, n/2, 0)$. That is, those cyclic cumulants corresponding to $(2, 1, 0), (4, 2, 0),$ and $(6, 3, 0)$. The $(2, 1, 0)$ cyclic cumulant is the standard autocorrelation function, and therefore reflects the presence of all signals in the data, including the additive white Gaussian noise. It should have a value equal to the sum of the average power levels of all the signals in the data, which is equal to the average power of the data itself: $P = 1.0 + 1.58 + P_N = 2.6.$ We processed the data by first applying the band-of-interest detector, which automatically finds the occupied band and filters the out-of-band noise away prior to cyclic-cumulant estimation. Therefore, the remaining Gaussian noise in the processed data is approximately $P_N = N_0 B = 0.1 (2/10) = 0.02$. So the value of the $(2, 1, 0)$ cyclic cumulant shown in each of the two plots above is correct.

The $(4, 2, 0)$ and $(6, 3, 0)$ cyclic cumulants reflect the presence of all non-Gaussian signals components in the data, because higher-order cumulants for Gaussian variables and processes are zero. But what are their theoretical values here? The $(4, 2, 0)$ cyclic-cumulant magnitude for our 100-percent EBW unit-power BPSK1 is equal to $2.28$. Therefore, the cyclic-cumulant magnitude for 2-dB BPSK2 is equal to $2.28 (1.58^2) = 5.69$. The sum of the cyclic-cumulant magnitudes is then $2.28 + 5.69 = 7.97$. The warped version of this magnitude is $7.97^{(2/4)} = 2.82$, which is what we see in the graphs. Similarly, the warped $(6, 3, 0)$ cyclic cumulant should be about $4.9$.

Now, what happens when we try to estimate the cyclic cumulants for BPSK1 in the data BPSK1 + BPSK2, but only use the lower-order cycle frequencies for BPSK1? (That is, ignore all the cycle frequencies for BPSK2 in the formula (5) above.) We obtain the warped cyclic cumulants and corresponding relative error plots shown here:

The errors are largely constrained to the $(n, n/2, 0)$ cyclic cumulants. This is because the cyclic moment function for $(n, n/2, 0)$ contains contributions from various combinations of BPSK2’s lower-order cycle frequencies that add up to zero, and these are not accounted for in (5) when we use only BPSK1’s lower-order cycle frequencies. If BPSK2 had parameters such that combinations of its lower-order cycle frequencies, or combinations of its and BPSK1’s lower-order cycle frequencies, added to a BPSK1 cycle frequency, we’d see additional large errors in the plot above.

### Blind Estimation of Temporal Moments and Cumulants

In blind estimation we don’t know any of the cycle frequencies exhibited by the data, except of course the pure cycle frequencies $\alpha = 0$ for $(n,m,k) = (n, n/2, 0)$, which all non-Gaussian power signals possess. So we have to jointly or sequentially estimate the cycle frequencies and the cyclic moments/cumulants.

How might we go about blindly estimating the cycle frequencies for some arbitrary data record?

1. For $n=2$, we have already seen that the thresholded spectral coherence function is a good way to blindly estimate both non-conjugate and conjugate cycle frequencies. That’s a good start, but we also already know that key cycle frequencies don’t show up for many signal types until order $n=4$ (such as QPSK). So we have to do something else to get at these higher-order cycle frequencies.
2. We can find impure cycle frequencies by finding the periodic components in the outputs of particular nonlinear operations on the data, such as $(\cdot)^4$. More on this later, in the section below on partially blind cumulant estimation.
3. We can apply synchronized averaging to the outputs of particular nonlinear operations. This is an interesting and under-appreciated approach, both because synchronized averaging is under appreciated and because the output corresponds to the time-varying temporal moment function itself, rather than a cyclic moment or cyclic cumulant. So, let’s look at the application of synchronized averaging to the cyclic-moment and cyclic-cumulant estimation problem. I think along the way this might help solidify some of the concepts in higher-order cyclostationarity.

#### Synchronized Averaging Applied to Time-Varying Moment/Cumulant Estimation

First let us define synchronized averaging (see also The Literature [1] and related work by Gardner). Suppose we have some power signal $x(t)$ which is defined for all time. For period $T_p$, the synchronized averaging operation is

$\displaystyle A_x(t, T_p) = \lim_{K\rightarrow\infty} \frac{1}{2K} \sum_{k=-K}^K x(t + kT_p)r_{T_p} (t + kT_p), \hfill (6)$

for $t \in [0, T_p]$. In (6), $r_{T_p}(t)$ is the rectangle function

$\displaystyle r_{T_p}(t) = \left\{ \begin{array}{ll} 1, & -T_p \leq t < T_p \\ 0, & \mbox{\rm otherwise} \end{array} \right. \hfill (7)$

So, for example, $A_x(0, T_p)$ is the average of $\{x(kT_p)\}$ over all integers $k$. If there is a periodic component of $x(t)$ that has period $T_p$ or a period that is a submultiple of $T_p$, it will survive the averaging operation in (6). If not, the average will be zero.

If the periodic component in $x(t)$ is a sine wave, then you’re probably thinking why not just Fourier transform $x(t)$ and look for the largest peak? And that’s fine, and we do that all the time in CSP and signal processing in general, including in a later part of this post. But if the periodic component is “sophisticated” or “rich”, meaning that its Fourier series has many components with similar Fourier-coefficient magnitudes, Fourier analysis can be tricky in terms of deciding what is a significant peak and what isn’t (a thresholding problem).

In practice, we can implement (6) in discrete time and for finite $K$ determined by the length of the available data record.

So how does this apply to blind estimation of temporal moments and cumulants for a cyclostationary signal?  Recall that an $n$th-order temporal moment function is just the (almost) periodic component of a homogeneous nonlinear transform of the signal,

$\displaystyle R_x(t, \boldsymbol{\tau}; n,m) = E^{\{\alpha\}} \left[ \prod_{j=1}^n x^{(*)_j} (t + \tau_j) \right]. \hfill (8)$

If we apply the synchronized averaging operation to the delay product in (8), ideally we should obtain $R_x(t, \boldsymbol{\tau}; n,m)$. That is, we obtain the whole time-varying moment function itself in one shot.

However … and you knew there would be complications … the synchronized averaging operation requires selection of the parameter $T_p$, which is presumably related to the cycle frequencies exhibited by $x(t)$, which we are saying we don’t know. So how do we choose $T_p$?

What follows is an approach that builds on the fundamental tool of synchronized averaging, but is still fully blind. Let’s call the finite-$K$ synchronized average estimate $A_x(t, T_p, K)$. The periodized version of this signal is built by concatenating multiple instances of $A_x(t, T_p, K)$ until the length of the result is just greater than the length $T$ of the data record $x(t)$. The periodized waveform is denoted by $Z_x(t, T_p, K)$.

The power of the error between $x(t)$ and $Z_x(t, T_p, K)$ can be used to assess whether the signal $A_x(t, T_p, K)$ captures an actual periodic component in $x(t)$ or not. Before applying to our BPSK data records, let’s look at a simple case to fix the ideas. The data record is a single sine wave in white Gaussian noise; the sine wave and the noise each have unit power and the frequency of the sine wave is $1/600$ (normalized times and frequencies are used). The period is therefore $600$ samples, but the signal is also periodic with period $1200$ samples, $1800$ samples, etc., for all multiples of $600$.

We apply the synchronized averager to $65536$ samples of the sine-wave data record, and we store the power in the error signal as well as the obtained period for the value of $T_p$ that achieves the minimum error power. Here is the plot of the error power versus candidate $T_p$:

Notice that the error power (or “residual” power) reaches a minimum of $1.0$ at $T_p = 600$. This indicates that the full power of the signal component is removed by the periodized waveform using the $T_p = 600$ period because we know the noise power is unity. In this way, the periodicity of the input is blindly estimated.

Of course things get more complex when the input contains a sophisticated periodic component, made up of many different (Fourier-series) sine-wave components, but the basic idea is the same: look for a minimum in the error-power function to estimate the period.

So now let’s apply the synchronized averaging operation to the noisy BPSK1 signal. To find the various $n$th-order temporal cumulant functions, we can find the various $n$th-order temporal moment functions using synchronized averaging and then combine them using the moment-to-cumulant formula (4).

First up is $(n,m) = (2, 0)$:

For this input, what is the period? It’s not quite as simple as in the sine-wave illustration above. For $(n,m) = (2,0)$, the BPSK1 cycle frequencies are $2f_c + k/T_0 = 0.02 + 0.1k$ for $k \in \{-1, 0, 1\}$. That set of numbers turns out to be $\{-0.08, 0.02, 0.12\}$, which correspond to periodicities with periods $\{12.5, 50, 8.3333\}$. Each of these periods divides the candidate period $T_p = 50$. The synchronized averager shows a nice local minimum at $T_p =50$, and of course there are other local minima at multiples of $50$.  The algorithm picks $T_p = 500$, and produces the following period:

Now, the plot of the period so obtained is the actual time-varying moment function. We almost never directly plot this function, preferring to focus on its Fourier-series components, the cyclic moments. Taking the discrete Fourier transform leads to the cyclic moments for $(n,m) = (2, 0)$:

In this plot of the CTMFs we have also included the CTMF values using the theoretical formulas for BPSK (no estimation), and we see good agreement.

We can continue this kind of analysis, finding the time-varying temporal moment functions for all needed $(n,m)$, and then simply combining them according to the moment-to-cumulant formula. The result is a blind estimate of the temporal cumulant function, and the cyclic cumulants are obtained by Fourier transforming and then picking peaks with high local SNR. It is important to note that the plots of the CTMFs are for our edification–they show agreement between theory and estimate–but we don’t need to do that Fourier analysis until we obtain a temporal cumulant if we are looking for pure cycle frequencies.

And here are the temporal cumulants:

We conclude that the synchronized averaging approach to blind estimation of cyclic cumulants can work, provided the final temporal cumulant function can be analyzed into its significant Fourier components through some thresholding operation.

Now let’s apply this to the data that consists of BPSK1+BPSK2. Recall that BPSK1 has unit power and BPSK2 has power level $1.6$, or $2$ dB. Also, the bit rate for BPSK1 is $1/10$ and that for BPSK2 is $1/11$. The carrier offset for BPSK1 is $0.01$ and that for BPSK2 is $0.002$.

And the cyclic cumulants:

### Partially Blind Estimation of Temporal Cumulants

This post is already overly long; thanks to those that have made it to the end. I’m running out of steam, but we have the lingering question of “What to do about partially blind estimation of cyclic cumulants?” A typical example is the case of a digital QAM signal in noise. The symbol rate is known, but the carrier offset and modulation type (exact constellation) are not. For such signals, the carrier frequency is needed for specification of cycle frequencies starting with order $n=4$.

We can definitely apply synchronized averaging to the $(n,m) = (4,0)$ delay product and look for the peak in the Fourier transform of the result. We can also realize that we can simply look at the Fourier transform of the lag product itself. The peak will occur for frequency equal to the quadrupled carrier. Other signal types may have to be handled in ways that are specific to their cycle-frequency structure;  the synchronized averaging method is universally applicable.

As always, please note in the comments any questions, errors, or suggestions for improvement related to this post.