CSP Estimators: Cyclic Temporal Moments and Cumulants

How do we efficiently estimate higher-order cyclic cumulants? The basic answer is first estimate cyclic moments, then combine using the moments-to-cumulants formula.

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. That is, the two-dimensional delay vector $\boldsymbol{\tau} = [\tau_1\ \ \tau_2]$ is set equal to $[\tau/2\ \ -\tau/2]$ .

The more interesting case is when the order $n$ is greater than two. Most communication signal models possess odd-order moments and cumulants that are identically zero, so the first non-trivial order $n$ greater than two is four. 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 call this partially blind processing.

Signals and Scenes 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.0$ (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 frequency offset 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 frequency offset, 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:

psd_bpsk_noisy — Figure 1. Estimated PSD for BPSK1 in noise. This BPSK signal is a square-root raised-cosine signal with excess bandwidth (filter roll-off) of 100% (1).

psd_bpsk_interf — Figure 2. Estimated PSDs for BPSK1 plus BPSK2 in noise, as well as PSD estimates for noise-free BPSK1 and BPSK2.

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. Let me elaborate on this idea.

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 integral can be replaced by a sum, and the cycle-frequency parameter $\alpha$ restricted to values on a regular grid. This is easily implemented by the discrete Fourier transform, which is efficiently computed by the fast Fourier transform.

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 nth-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:

sig1_hoccs — Figure 3. **Estimated** warped cyclic-cumulant magnitudes using the cyclic-moment-to-cyclic-cumulant formula (5). Warping a cyclic-cumulant magnitude means raising it to the power $2/n$ , where $n$ is the cyclic-cumulant order. The corresponding theoretical warped cyclic-cumulant magnitudes are shown in Figure 4.

Figure 3. **Estimated** warped cyclic-cumulant magnitudes using the cyclic-moment-to-cyclic-cumulant formula (5). Warping a cyclic-cumulant magnitude means raising it to the power $2/n$ , where $n$ is the cyclic-cumulant order. The corresponding theoretical warped cyclic-cumulant magnitudes are shown in Figure 4.

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:

ideal_bpsk_hoccs — Figure 4. **Theoretical** warped cyclic-cumulant magnitudes for the BPSK1 signal. Compare to the estimates in Figure 3.

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 offset and $1/T_0$ is the symbol rate.

error_sig1_hoccs — Figure 5. Normalized (relative) error in the cyclic-cumulant estimates corresponding to Figure 3. Note that the vertical scale is about one hundred times smaller than those in Figures 3 and 4, indicating small errors in general.

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:

joint_bpsk_hoccs_soi — Figure 6. Estimated warped cyclic-cumulant magnitudes for the case of BPSK1 + BPSK2 + noise. The cyclic-moment-to-cyclic-cumulant formula (5) is used. The cyclic-cumulant cycle frequencies are those of BPSK1, but the set of lower-order cyclic moments used in (5) considers all impure (moment) cycle frequencies exhibited by BPSK1, BPSK2, and the noise.

joint_bpsk_hoccs_interf — Figure 7. Estimated warped cyclic-cumulant magnitudes for the case of BPSK1 + BPSK2 + noise. The cyclic-cumulant cycle frequencies are those of BPSK2, but the set of lower-order cyclic moments used in (5) considers all impure (moment) cycle frequencies exhibited by BPSK1, BPSK2, and the 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 processed 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 = P_{BPSK1} + P_{BPSK2} + P_N = 1.0 + 1.58 + P_N = 2.58 + P_N.$ 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. (Note that the warping is inconsequential for $n=2$ because $2/n = 2/2 = 1$ .)

The $(4, 2, 0)$ and $(6, 3, 0)$ cyclic cumulants reflect the presence of all non-Gaussian signal 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 in Figures 8 and 9.

solo_bpsk_hoccs_soi — Figure 8. Estimated warped cyclic-cumulant magnitudes for BPSK1 using input data equal to BPSK1 + BPSK2 + noise. The cyclic-cumulant cycle frequencies are those of BPSK1, and the set of lower-order cyclic moments used in (5) considers only impure (moment) cycle frequencies exhibited by BPSK1.

error_soi_solo_hoccs — Figure 9. Normalized (relative) error in the cyclic-cumulant estimates corresponding to Figure 8.

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 Figure 9.

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 (1) the cycle frequencies, and (2) the cyclic moments/cumulants.

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

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.
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.
We can apply synchronized averaging to the outputs of particular nonlinear operations. This is an interesting 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/2 \leq t < T_p/2 \\ 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 (which is 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 us to select the parameter $T_p$ , which is presumably related to the cycle frequencies exhibited by $x(t)$ , which we are saying we don’t know (we’re in the blind processing context). 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)$ . A 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. If $A_x(t, T_p, K)$ represents an actual periodic component in $x(t)$ , then subtracting its periodized version $Z_x$ from $x(t)$ will produce a signal with smaller power than $x(t)$ . Before applying to our BPSK data records, let’s look at a simple case to fix the ideas.

The data record for our simple illustrative case 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, as nearly always at the CSP Blog). The period is therefore $600$ samples, but the signal is also periodic with period $1200$ samples, $1800$ samples, or generally $k600$ samples for all integers $k >0$ .

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$ :

sa_power_tone_input — Figure 10. Measured power in the error (difference) between the original signal $x(t)$ and the periodized extracted period given by $Z_x(t, T_p, K)$ . For most candidate periods $T_p$ , there is no periodic component in $x(t)$ , and so the synchronized averaging operation produces a low-power ouput. For $T_p = 600$ , however, there is a periodic component, and the low power in the error reflects the fact that the periodic component is successfully subtracted from the original data.

Figure 10. Measured power in the error (difference) between the original signal $x(t)$ and the periodized extracted period given by $Z_x(t, T_p, K)$ . For most candidate periods $T_p$ , there is no periodic component in $x(t)$ , and so the synchronized averaging operation produces a low-power ouput. For $T_p = 600$ , however, there is a periodic component, and the low power in the error reflects the fact that the periodic component is successfully subtracted from the original data.

sa_period_tone_input — Figure 11. Plot of the period estimated by the synchronized averager when $T_p = 600$ . This is simply one period of the sine wave in $x(t)$ , with residual noise in the period estimate.

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.

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 first find the various-order temporal moment functions using synchronized averaging and then combine them using the moment-to-cumulant formula (4). For the synchronized-averaging results here, I always take the delays $\tau_j$ to be zero. This is fine for most signal types. An exception is those PSK/QAM signals that employ rectangular pulses. (See the gallery of cyclic correlations for where the cyclic-autocorrelation magnitudes peak as functions of the delay $\tau$ )

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

bpsk_2_0_tmf_residual — Figure 12. Residual power in the synchronized-averager output as a function of the candidate period $T_p$ .

Figure 12. Residual power in the synchronized-averager output as a function of the candidate period $T_p$ .

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\}$ in units of samples. 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:

bpsk_2_0_tmf — Figure 13. Estimated temporal moment function using synchronized averaging for the candidate period $T_p = 500$ , which is actually $10$ periods of the temporal moment function.

Figure 13. Estimated temporal moment function using synchronized averaging for the candidate period $T_p = 500$ , which is actually $10$ periods of the temporal moment function.

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 of the signal in Figure 13 leads to the cyclic moments for $(n,m) = (2, 0)$ :

bpsk_2_0_ctmf — Figure 14. Comparison of the Fourier components of the synchronized-averager output in Figure 13 with direct estimates of the cyclic moment function magnitudes. They should be equal if the period selected by the synchronized averaging is a multiple of the true cyclic-moment period.

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 the temporal cumulant 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.

Figures 15 through 22 show temporal moment functions and their Fourier components for other choices of $(n,m)$ . Comparison with directly estimated cyclic moment functions confirms the efficacy of synchronized averaging.

bpsk_2_1_tmf — Figure 15. The temporal moment function obtained by synchronized averaging for the noise BPSK1 signal and $(n,m) = (2,1)$ .

Figure 15. The temporal moment function obtained by synchronized averaging for the noise BPSK1 signal and $(n,m) = (2,1)$ .

bpsk_2_1_ctmf — Figure 16. Fourier components of the temporal moment function in Figure 15 along with direct cyclic-moment estimates for comparison.

bpsk_4_0_tmf — Figure 17. The temporal moment function obtained by synchronized averaging for the noise BPSK1 signal and $(n,m) = (4, 0)$ .

Figure 17. The temporal moment function obtained by synchronized averaging for the noise BPSK1 signal and $(n,m) = (4, 0)$ .

bpsk_4_0_ctmf — Figure 18. Fourier components of the temporal moment function in Figure 17 along with direct cyclic-moment estimates for comparison.

bpsk_4_1_tmf — Figure 19. The temporal moment function obtained by synchronized averaging for the noise BPSK1 signal and $(n,m) = (4, 1)$ .

Figure 19. The temporal moment function obtained by synchronized averaging for the noise BPSK1 signal and $(n,m) = (4, 1)$ .

bpsk_4_1_ctmf — Figure 20. Fourier components of the temporal moment function in Figure 19 along with direct cyclic-moment estimates for comparison.

bpsk_4_2_tmf — Figure 21. The temporal moment function obtained by synchronized averaging for the noise BPSK1 signal and $(n,m) = (4, 2)$ .

Figure 21. The temporal moment function obtained by synchronized averaging for the noise BPSK1 signal and $(n,m) = (4, 2)$ .

bpsk_4_2_ctmf — Figure 22. Fourier components of the temporal moment function in Figure 21 along with direct cyclic-moment estimates for comparison.

The cyclic cumulants obtained through synchronized averaging, the moment-to-cumulant formula, and Fourier analysis are shown in Figures 23-35, along with directly estimated cyclic cumulants (using the cyclic-moment-to-cyclic-cumulant formula (5)).

bpsk_4_0_ctcf — Figure 23. Cyclic-cumulant magnitudes for synchronized-averaging-based estimation, cyclic-moment-to-cyclic-cumulant estimation, and from theory for BPSK1 and $(n,m) = (4,0).$

bpsk_4_1_ctcf — Figure 24. Cyclic-cumulant magnitudes for synchronized-averaging-based estimation, cyclic-moment-to-cyclic-cumulant estimation, and from theory for BPSK1 and $(n,m) = (4,1).$

bpsk_4_2_ctcf — Figure 25. Cyclic-cumulant magnitudes for synchronized-averaging-based estimation, cyclic-moment-to-cyclic-cumulant estimation, and from theory for BPSK1 and $(n,m) = (4,2).$

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 machinery to the data that consists of BPSK1+BPSK2+noise. 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$ .

bpsk_2_0_ctmf_interf — Figure 26. Fourier components of the temporal moment function obtained by applying synchronized averaging to BPSK1+BPSK2+noise, along with direct cyclic-moment estimates for comparison, for $(n,m) = (2,0).$

Figure 26. Fourier components of the temporal moment function obtained by applying synchronized averaging to BPSK1+BPSK2+noise, along with direct cyclic-moment estimates for comparison, for $(n,m) = (2,0).$

bpsk_2_1_ctmf_interf — Figure 27. Fourier components of the temporal moment function obtained by applying synchronized averaging to BPSK1+BPSK2+noise, along with direct cyclic-moment estimates for comparison, for $(n,m) = (2,1).$

bpsk_4_0_ctmf_interf — Figure 28. Fourier components of the temporal moment function obtained by applying synchronized averaging to BPSK1+BPSK2+noise, along with direct cyclic-moment estimates for comparison, for $(n,m) = (4,0).$

bpsk_4_1_ctmf_interf — Figure 29. Fourier components of the temporal moment function obtained by applying synchronized averaging to BPSK1+BPSK2+noise, along with direct cyclic-moment estimates for comparison, for $(n,m) = (4,1).$

bpsk_4_2_ctmf_interf — Figure 30. Fourier components of the temporal moment function obtained by applying synchronized averaging to BPSK1+BPSK2+noise, along with direct cyclic-moment estimates for comparison, for $(n,m) = (4,2).$

And the cyclic cumulants:

bpsk_4_0_ctcf_interf — Figure 31. Cyclic-cumulant magnitudes for synchronized-averaging-based estimation, cyclic-moment-to-cyclic-cumulant estimation, and from theory for BPSK1+BPSK2+noise and $(n,m) = (4,0).$

bpsk_4_1_ctcf_interf — Figure 32. Cyclic-cumulant magnitudes for synchronized-averaging-based estimation, cyclic-moment-to-cyclic-cumulant estimation, and from theory for BPSK1+BPSK2+noise and $(n,m) = (4,1).$

bpsk_4_2_ctcf_interf — Figure 33. Cyclic-cumulant magnitudes for synchronized-averaging-based estimation, cyclic-moment-to-cyclic-cumulant estimation, and from theory for BPSK1+BPSK2+noise and $(n,m) = (4,2).$

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 a cycle 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. See also the post on synchronization for more discussion about the problem of blind fourth-order cycle-frequency estimation.

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

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. View all posts by Chad Spooner

49 thoughts on “CSP Estimators: Cyclic Temporal Moments and Cumulants”

Ramesh Rao says:

March 9, 2018 at 1:16 am

I would like some help regarding estimation of fourth order cyclic cumulants. Can someone post the code for that. Much appreciated!

Reply
1. Chad Spooner says:
  
  March 9, 2018 at 11:04 am
  
  Please see the following post for understanding how to get help from the CSP Blog:
  http://cyclostationary.blog/2017/04/13/blog-notes-and-how-to-obtain-help-with-your-csp-work
  
  Reply
quyangdl says:

April 4, 2018 at 7:55 pm

Hi, Dr Spooner

In the figure, “warped measured cc mags for 0-dB bpsk1 in BPSK1+BPSK2 DATA using all valid lower-order CFs”, the C(2,0,0) is 1, which is equal to the power of bpsk1.

But if we want the cc mags of data(bpsk1+bpsk2) together, the C(2,0,0), C(2,1,0) and C(2,1,1) will be same and be equal to the average power of the data itself. The second order cyclic moment CM(2,0,0) will be equal to E^{\alpha_{1}}[s^2]+E^{\alpha_{2}}[s^2] , where s=bpsk1+bpsk2, \alpha_{1} is cyclic frequency of bpsk1, and \alpha_{2} is cyclic frequency of bpsk2. Meanwhile, The second order cyclic moment is equal to the second order cyclic cumulant. In addition, the fourth order CCs will be same, such as C(4,0,0) = C(4,2,0) =C(4,1,0).

Am I right?

Thanks

Andy

Reply
1. Chad Spooner says:
  
  April 5, 2018 at 10:16 am
  
  Hey Andy. Thanks for looking at that estimation post.
  
  In the figure, “warped measured cc mags for 0-dB bpsk1 in BPSK1+BPSK2 DATA using all valid lower-order CFs”, the C(2,0,0) is 1, which is equal to the power of bpsk1.
  
  Yes.
  
  But if we want the cc mags of data(bpsk1+bpsk2) together, the C(2,0,0), C(2,1,0) and C(2,1,1) will be same and be equal to the average power of the data itself.
  
  No. Perhaps my notation is confusing you here. C(n, m, k) means the cyclic cumulant of order n, using m conjugations, and harmonic number k. This applies to QAM/PSK/PAM/CPM/OQPSK/MSK etc. for which the cycle frequencies conform to the formula
  
  alpha(n,m,k) = (n-2m)fc + k/T0
  
  where fc is the carrier offset and 1/T0 is the symbol rate.
  
  Now, the two signals, BPSK1 and BPSK2, have different carrier offsets and different symbol rates. So if we say “C(2, 0, 0)” we have to say which signal this is relative to–what are fc and T0 for this cyclic cumulant? For the BPSK1 + BPSK2 data, C(2, 0, 0) for BPSK1 corresponds to the cycle frequency alpha = 2fc_1 + 0/T0_1 = 2fc_1 = 2(0.01) = 0.02.
  
  Cyclic cumulant C(2, 1, 0) is a special case because it doesn’t depend on the carrier offsets or the symbol rates, it corresponds to alpha = 0 = (2 – 2(1))fc + 0/T0 = 0. So, C(2, 1, 0) is in fact the average power of the data (= P1 + P2 + PNoise). C(2, 0, 0) for BPSK1 will be the doubled-carrier cyclic cumulant for BPSK1 as if BPSK2 was not present. C(2, 0, 0) for BPSK2 will be the doubled-carrier cyclic cumulant for BPSK2 as if BPSK1 was not present.
  
  The second order cyclic moment CM(2,0,0) will be equal to E^{\alpha_{1}}[s^2]+E^{\alpha_{2}}[s^2] , where s=bpsk1+bpsk2, \alpha_{1} is cyclic frequency of bpsk1, and \alpha_{2} is cyclic frequency of bpsk2.
  
  I think you’re mixing up some notation and concepts here. The notation E^alpha[x(t)] means the sine-wave extraction operator applied to x(t). If y = E^alpha[x(t)], then if x(t) is cyclostationary, y is a function of time, y(t). It extracts all the periodic components of its
  input x(t). On the other hand CM(2, 0, 0) is a cyclic moment, with order 2, no conjugations, and harmonic of 0. It is a number, not a function of time (although it is a function of delay tau, which is suppressed here because I usually consider tau = 0).
  
  Hope that helps.
  
  Reply
Charan says:

April 17, 2018 at 7:41 am

Hi, Dr. Spooner

I have been trying to do the estimation of second order temporal moment the same way as you explained in the blind estimation part of this post with the BPSK you have given in the downloads section of this blog “MATLAB script for rectangular-pulse BPSK example” (I changed fc to 0.01 for getting the same results as in this post). But, while calculating Synchronous averaging and plotting the power difference between original Vs the extracted sine wave, I could not get the little dips at odd multiples of 25, which were obtained by you as shown in this figure linked below

https://cyclostationary.files.wordpress.com/2017/09/bpsk_2_0_tmf_residual.jpg

I have written the Synchronous averaging code myself and also used the built-in “tsa” function in MATLAB. Nonetheless, I was not getting the small dips when I considered Tp = 500 or 50. And hence I was getting only one peak in the cyclic moment of (2,0) at 0.02. I have taken \tau as zero and just tried to calculate R_x(2,0) = TSA[(y_of_t).*(y_of_t)]. Can you help me rectify the wrong I have done in the calculation.

Thanks,
Charan

Reply
1. Chad Spooner says:
  
  April 17, 2018 at 9:23 am
  
  Thanks for the comment Charan, and for studying that post on estimating cyclic moments and cumulants.
  
  I didn’t know about MATLAB’s tsa, so thanks for that tip!
  
  The good news is that I don’t think you’ve done anything wrong. You won’t be able to replicate my results unless you use the same kind of signal. Recall that I used BPSK with square-root raised-cosine pulses ( $r=1$ ) in that post. You used my code that creates a rectangular-pulse BPSK signal.
  
  So should the result be the same? No, because when you look at squaring the rectangular-pulse BPSK signal (and “squaring” is implied by $(n,m) = (2,0)$ ), you get a single sine-wave at the squarer output. It has frequency equal to the doubled carrier frequency which is $0.02$ . But when you square the square-root raised-cosine BPSK signal, you get the doubled carrier plus two other sine-wave components. What about other quadratic nonlinearities? What if we don’t restrict the delays in the second-order lag product to zero?
  
  I redid the synchronized averaging using a rectangular-pulse BPSK signal with the same parameters as you mention in your comment above. Here is the result (which should match what you’ve been seeing) for delays of zero in the lag product:
  
  We see only residual-power dips for periods corresponding to harmonics of the doubled carrier ( $k*(1/0.02) = k*50$ ).
  
  Then I set the delay of one of the factors in the lag product to $5$ and kept the other at $0$ . The synchronized averaging result is:
  
  See if you can replicate.
  
  And try to understand why the smaller dips happen at multiples of $25$ samples. Hint: What are the other $(2,0)$ cycle frequencies for this signal?
  
  Reply
Andy says:

June 22, 2018 at 1:59 pm

Hi, Dr.Spooner

I followed your instruction. I could got the synchronized averaging result of sine wave. However, I could not get the the synchronized averaging results of BPSK sine, which are shown after “First up is (n,m) = (2, 0):”

My procedure is as following,
1) Generate a BPSK signal x(t).
2) Set a Tp vector.
3) According to different Tp, select different length (Tp length) signals s(t) from x(t).
4) perform average operation for different s(t).^2. r(t) = 1/N(\sum s(t).^2 )
5) combine multiple r(t) into an new signal f(t), which has same length as x(t).
6) compare the power difference between f(t) and x(t). Where f(t) is already in power format.

Did I miss anything?

Thanks.

Best,

Andy

Reply
1. Chad Spooner says:
  
  June 22, 2018 at 3:13 pm
  
  Your description is a bit vague, so there could be errors lurking in places I can’t see yet… But I would not do Steps 3 and 4 the way you’re doing them. I would apply the exact same code I had created for the case of a sine-wave input to every other signal. In other words, take your BPSK signal s(t), square it, then use that as the input to your synchronized averager.
  
  Step 6 is not clear to me here. You want to find the variance of the difference between f(t) and x(t). But I don’t understand what “f(t) is already in power format” means.
  
  Overall, if you want to leave a more detailed comment, I believe we will be able to find the problem(s).
  
  Reply
  1. Andy says:
    
    June 25, 2018 at 2:09 pm
    
    Hi, Dr.Spooner
    
    Thanks for you help. I got the result. However, I am not quite understand that if we want to show up the symbol duration, why we need time delay tau in the product function x(t)x(t+tau)? I read some of your papers, everywhere has time delay tau, such as spectral correction function, cyclic cumulant function and so on.
    
    How could we realize the delay tau? When we could set the tau is equal to 0?
    
    Thanks
    
    Reply
    1. Chad Spooner says:
      
      June 25, 2018 at 2:49 pm
      
      Glad you got a good result for BPSK!
      
      Well, for the square-root raised-cosine pulses used by BPSK1 and BPSK2 here in this post, we don’t need any non-zero delays $\tau_j$ . So in all those results I posted, the delay vector was all zeros. I’ve updated the post to make that explicit.
      
      If you did want to use non-zero delays, you can approximate the delay in MATLAB using circshift.m. By including leading and trailing zeros, and then applying circshift.m, you can also avoid the edge effects implied by circularly shifting. If the delays are small relative to the length of the data block to be processed, the circular shifting isn’t a big error factor.
      
      Yes, I do typically use a delay vector of zero in cyclic moments and cumulants, but not in the spectral correlation function. When we compute an estimate of the spectral correlation function for a particular value of cycle frequency $\alpha$ , we can inverse transform that estimate (over frequency $f$ ) to obtain the cyclic autocorrelation for a large number of delays. So, in effect, the spectral correlation function includes many values of delay. Similar statements apply for the cyclic polyspectrum.
      
      Reply
      1. Andy says:
        
        July 2, 2018 at 5:46 pm
        
        Hi, Dr.Spooner
        
        For the synchronized averaging method, If we set the normalized carrier frequency f, whose reciprocal is not an integer. In this case, we could not set the period T_p as a float number in MATLAB. The T_p will not be exactly equal to the 1/f.
        
        If use integer T_p vector, I got different size dips in T_p, 2T_p, 3T_p… If we choose a carrier frequency f value, let 1/f is an integer, I will get multiple same dips at 1/f, 2/f, 3/f ….
        
        So, how could we deal with this case? Or do I miss anything ?
        
        Thanks,
        
        Best.
        
        Andy
      2. Chad Spooner says:
        
        July 3, 2018 at 11:48 am
        
        I don’t think you are missing anything. The different depths of the dips arises because the different periods $T_p, 2T_p, \ldots$ have different distances to the true periods $1/f, 2/f, \ldots$ .
        
        Oversampling will help. Also, you could fashion an iterative algorithm that finds the minima, estimates the carrier offset, shifts the data by that estimate, then repeats.
Andy says:

October 2, 2018 at 6:58 am

Hi, Dr.Spooner

I have a question for blind calculate CTCM of a BPSK/QPSK via synchronized averaging DFT.

Could I believe that synchronized averaging method is kind of perform general sine-wave extraction operator, like Equation (20) on the link https://cyclostationary.blog/2016/02/14/cyclic-temporal-cumulants/ ?

Assume that we have a QPSK signal S(t), If we get the correct Tp(actually carrier frequency), we can get its synchronized averaging result ~~via Tp, which should be Fs/(4Fc）, and then we can get its CTCM result through taking a FFT operation for the~~ .

~~I am right?~~

~~Thanks.~~

Reply
1. Chad Spooner says:
  
  October 2, 2018 at 8:32 am
  
  Hey Andy.
  
  Your comment is garbled and difficult to read. You’ve crossed out the last few lines (perhaps accidentally).
  Can you rewrite it and repost? Also, what do you mean by “CTCM”? I discuss the CTMF and the CTCF in various posts. These are the cyclic temporal moment function and cyclic temporal cumulant function, or cyclic moment and cyclic cumulant, respectively.
  
  Reply
Andy says:

October 2, 2018 at 8:55 am

Hi, Dr.Spooner

I have a question for blind calculate cyclic moment of a BPSK/QPSK via synchronized averaging DFT.

Could I believe that synchronized averaging method is kind of perform general sine-wave extraction operator, like Equation (20) on the link https://cyclostationary.blog/2016/02/14/cyclic-temporal-cumulants/ ?

Assume that, we have a QPSK signal S(t), If we get the correct Tp(actually carrier frequency) via synchronized averaging method, we can get its synchronized averaging result according to the Tp, which should be Fs/(4Fc), and then we can get its cyclic moment after taking a FFT operation for the synchronized averaging result.

Am I right?

Thanks.
Best,

Andy

Reply
1. Chad Spooner says:
  
  October 3, 2018 at 8:15 am
  
  Could I believe that synchronized averaging method is kind of perform general sine-wave extraction operator, like Equation (20)
  
  Yes, SA can do that for you, but you have to be careful because with DSP we have a sampled time-series, and so without resampling, you are constrained to check for candidate periods $T_p$ that are integer multiples of the sampling increment.
  
  Assume that, we have a QPSK signal S(t), If we get the correct Tp(actually carrier frequency) via synchronized averaging method, we can get its synchronized averaging result according to the Tp, which should be Fs/(4Fc)
  
  Well, for QPSK, if we look at the $(n, m) = (4, 0)$ lag product, we see that there are three cycle frequencies: $\alpha = 4f_c, \alpha = 4f_c+f_{sym},$ and $\alpha = 4f_c - f_{sym}$ (for SRRC QPSK). If those three cycle frequencies are commensurate, they will all correspond to some harmonic of a common period that might be found using synchronized averaging. Otherwise, depending on your sampling rate, you might obtain a candidate period from synchronized averaging that corresponds to just one of them.
  
  and then we can get its cyclic moment after taking a FFT operation for the synchronized averaging result.
  
  If the synchronized averaging works out and we’ve identified the fundamental period for the periodic temporal variation of the $(4,0)$ lag product, then yes, if you create a signal by periodizing the obtained period, you can analyze its sine-wave components using the Fourier transform.
  
  Reply
Andy says:

December 11, 2018 at 12:17 pm

Hi, Dr.Spooner

Does synchronized averaging method has better frequency resolution than FFT?

It seems that synchronized averaging method require more data to get good frequency resolution, which enable us to find the multiple peaks(if more signals). However, if we use FFT with longer FFT size, we still can get good frequency resolution.

So, what is the big difference/ benefit between synchronized averaging method and FFT?

Thanks,

Best,
Andy

Reply
1. Chad Spooner says:
  
  December 11, 2018 at 5:01 pm
  
  In my context of CSP, synchronized averaging is a time-domain process. Its benefit is that it can estimate the period of a periodic signal component even when the period waveform is highly complex. In that case, the periodic signal has a rich harmonic structure–sometimes no single element is particularly strong. Fourier methods can then have a hard time detecting the various harmonics. So Fourier analysis is particularly good for finding strong sine-wave components in a data set, and synchronized averaging is good for finding periods of periodic signals that do not possess a single dominant sine-wave component.
  
  If we want to talk about frequency resolution, we can talk about estimating the fundamental frequency of a periodic waveform with synchronized averaging. Suppose the period is $T_p$ and we are using a sampling rate of one. So we might estimate the period as $\hat{T}_p = T_p-1$ , $T_p$ , or $T_p + 1$ when things are going well. The corresponding frequency estimate is $1/\hat{T}_p$ .
  
  The error in the frequency estimate is
  
  $\displaystyle e = \left[ \frac{1}{T_p} - \frac{1}{\hat{T}_p} \right]$
  
  or
  
  $\displaystyle e = \frac{|\hat{T}_p - T_p |}{T_p \hat{T}_p}$
  
  Supposing the estimate is good, we have
  
  $\displaystyle e = \frac{1}{T_p\hat{T}_p} \approx \frac{1}{T_p^2}$
  
  The normalized error is
  
  $\displaystyle e_n = \frac{e}{1/T_p} = 1/T_p$ .
  
  So the resolution depends on the period length $T_p$ . This is why oversampling helps in synchronized averaging.
  
  Reply
Andy says:

February 21, 2019 at 7:54 am

Hi, Dr. Spooner

The high order cyclic cumulant value is sensitivity with signal power. Is any way to remove the power information from the cyclic cumulant?

Thanks

Reply
1. Chad Spooner says:
  
  February 21, 2019 at 1:51 pm
  
  Strictly speaking, no. Unless you multiply the data by zero, then all the cyclic cumulants are zero, and there is no sensitivity to the power.
  
  But maybe you mean normalize the cyclic cumulant or data so that the cyclic cumulant corresponds to a signal power of, say, one, no matter what the original power was. In that case, you can normalize the data by the square root of the signal power, if you know it. If you don’t know it, you can estimate it. If you don’t want to do that, and you know the modulation type, you can compute the cyclic cumulants for the original data, then determine the signal power by comparing, say, $C_x^0(0; 2, 1)$ with the corresponding cumulant from the known unit-power modulated signal. That will give you an estimate of the signal power and you can normalize after that.
  
  If you don’t know the modulation type, or if the signal is cochannel with other signals, you’ll have to extend these ideas.
  
  Reply
Andy says:

October 4, 2019 at 10:39 am

Hi, Dr. Spooner

In the figure, the cyclic cumulant result of BPSK1 +BPSK2, do you assume that two signals are coming from same channel? and they have same carrier phase shift? But if in real-time system, BPSK1 and BPSK2 may come from different/same channel with different phase shift. In this case, the cyclic cumulant result of BPSK1+BPSK2 will not be unique. The result will be changed by two random phase shift. Please correct me if I am wrong.

So, If I want to classify BPSK1+BPSK2 with other communication signal, I have to use different order cumulant instead of cyclic cumulant, right? Becasue cyclic cumulant is sensitive to phase shift, in real-time transmission, we cannot assume two signals have same phase shift. However, cumulant is statistical value, which may not sensitive to phase shift.

Thanks,
Best,
Andy

Reply
1. Chad Spooner says:
  
  October 5, 2019 at 8:48 am
  
  In the figure, the cyclic cumulant result of BPSK1 +BPSK2, do you assume that two signals are coming from same channel?
  
  No. I do not assume that they propagate through the same physical channel, nor do I assume that they are emitted from the same antenna.
  
  and they have same carrier phase shift?
  
  No. I do not assume they have the same carrier phase.
  
  But if in real-time system, BPSK1 and BPSK2 may come from different/same channel with different phase shift. In this case, the cyclic cumulant result of BPSK1+BPSK2 will not be unique. The result will be changed by two random phase shift. Please correct me if I am wrong.
  
  The plots are associated with the case in which the two BPSK signals have slightly different symbol rates and carrier frequency offsets. Their relative delay (symbol-clock phase) and carrier phase do not matter to the displayed results, where I am only showing the magnitudes of the estimated cyclic cumulants. If the symbol-clock phases and/or carrier phases were changed, the phases of the complex-valued cyclic cumulants would also change. But not the magnitudes. The estimation processes do not require knowledge of the phases; they are carried along in the various estimates of the complex-valued cyclic moments and are carried along in the synchronized-averaging process.
  
  If I want to classify BPSK1+BPSK2 with other communication signal, I have to use different order cumulant instead of cyclic cumulant, right?
  
  I don’t understand the grammar of “classify BPSK1+BPSK2 with other communication signal”. But you can successfully process the BPSK1 plus BPSK2 signals I use in this post using cyclic cumulants because the signals are statistically independent, and they have distinct modulation parameters (symbol rate and carrier offset). The phases do matter more, however, in the special case in which BPSK1 and BPSK2 have identical symbol rates and carrier frequencies, but that is a topic I’ve not yet introduced on the CSP Blog.
  
  Becasue cyclic cumulant is sensitive to phase shift, in real-time transmission, we cannot assume two signals have same phase shift. However, cumulant is statistical value, which may not sensitive to phase shift.
  
  No need to assume anything about the phases provided the two signals have distinct modulation parameters, in which case almost all the cyclic cumulants for each signal can be extracted from the data that consists of their noisy sum. The cyclic cumulants that cannot be extracted for each signal individually are those that correspond to $(n, m) = (n, n/2)$ and $\alpha = 0$ , as noted in the post. But those are just a few cyclic cumulants relative to all those that make up the full set of cyclic cumulants.
  
  I think your ideas about this might be more focused if you had a particular modulation classification approach or algorithm in mind. Then you can check for yourself whether the features for one signal can be extracted in the presence of the other signal.
  
  Reply
  1. Andy says:
    
    October 6, 2019 at 6:32 pm
    
    Dr.Spooner Thanks for your replies. Very helpful.
    
    Reply
Jiang’an Dai says:

May 29, 2021 at 11:42 pm

Hello，Mr.Spooner
Thank you for your long post about how to estimate Temporal Moments(M)/Temporal Cumulants(C)/Cyclic Temporal Moments(CM)/Cyclic Temporal Cumulants(CC).
In conclusion, two ways have been presented.
The first way is（CM->CC）：
1、estimate CM (equation 2) 2、CM-CC formula (equation 5).
The second way is（M->C->CC）：
1、synchronized averaging to get M(equation 6) 2、M-C formula (equation 4) 3、extract fourier coefficient CC.

My goal is to blind estimate high-order CC(n=6,8) such as CC60,CC80 with T-point data.
And i set all tao=0 for simplicity. So the only variable i care about in CC is cyclic frequence α.

Let’s first talk about CC40,then turn to higher-order CC(60,80)
For odd-order CMs equal zero，CC40=CM40-3*CM20^2.

The following is my navie method of estimation.
1、Using T-point data to get CM20,CM40 (equation 2).
This can be done by T-point fft easily.
2、Using CM-CC formula to compute CC (equation 5).
And the difficulty appears in this step because we can’t operate simple multiplication （CM20*CM20）.
CM20(β1)、CM20(β2)、CC40（α）must match the condition that α=β1+β2.

In practice, i implementation this by two layer for-loop in Matlab.
For 1000-point data(T=1000,100 symbols)，it takes 16s to compute CC40.That’s barely acceptable.

However,when i want to estimate CC60(CC60 = CM60-15CM20*CM40+30CM20^3), the maximum power of multiplication (CM20^3) is 3.
That means if we do as above ，we must implementation this by three layer for-loop.
That will increase the computation cost by 1000 times.
And CC80 by 1,000,000 times.That’s really unacceptable for practical application.
I’m wondering if there are any mistake about my thought of the first way。

For the second way，the key step is synchronized averaging to get time-varying moment.
It‘s a little complicated.So i divide it to four step.
1、synchronized averaging to get Ax(t,Tp）(equation 6)
2、periodic extension to get zx(t,Tp,K) with the same size of x(t)
3、compute the power of the error between x(t) and zx(t,Tp,K), determine Tp
4、using Tp get time-varying moment estimated value zx(t,Tp,K)
And i have some question about the step 3.
If we want to determine Tp，it associates with peak-picking operation.So how can we choose appropriate threshold？

It seems the performance of synchronized averaging is related to the data length T.
If T is not large enough(for example T=1000), will the performance deteriorate?

If environment is much worse (SNR=-20),will the performance deteriorate?

If for some time-varying moment with no cyclostationarity but higher time-varying moment with cyclostationarity （M20 of QPSK with no cyclostationarity and M40 with cyclostationarity ）, then we can’t find peek in the error plot between x(t) and zx(t,Tp,K).In this situation ,we can’t find an appropriate Tp，then the zx(t,Tp,K)=0?

And final i want to know if synchronized averaging method can have acceptable computational complexity for practical application compared with the first way.

Thanks,
Jiang’an

Reply
1. Chad Spooner says:
  
  May 31, 2021 at 8:11 am
  
  Welcome to the CSP Blog Jiang’an!
  
  For odd-order CMs equal zero，CC40=CM40-3*CM20^2.
  
  I don’t think your equation here is consistent with (5). The equation you wrote here is an oversimplification. You must obey (5).
  
  1、Using T-point data to get CM20,CM40 (equation 2).
  This can be done by T-point fft easily.
  2、Using CM-CC formula to compute CC (equation 5).
  And the difficulty appears in this step because we can’t operate simple multiplication （CM20*CM20）.
  CM20(β1)、CM20(β2)、CC40（α）must match the condition that α=β1+β2.
  
  The multiplications in (5) are simple–note that there is no time $t$ in this equation. It is simply a combination of cyclic-moment Fourier coefficients. I think the problem with your description is that you don’t acknowledge this fact. You say that you can use a ‘T-point FFT’ to get the ‘CM20’ and ‘CM40’ values. Consider the ‘CM20’ values. For many signal types (e.g., QPSK), CM20 is zero. So you have to find a way to throw away all the T FFT points–none of them will be a valid cyclic moment. For SRRC BPSK, CM20 will have three components corresponding to the three conjugate CFs $2f_c$ , $2f_c - f_{bit}$ and $2f_c + f_{bit}$ . So you have to find those in the $T$ FFT values.
  
  It looks like you are mistaking the Fourier transform of a nonlinear function of $x(t)$ with cyclic moments. They are not the same thing.
  
  Also, you don’t need to do completely blind searches for cycle frequencies for the larger values of $n$ . Remember the formula for the CFs for PSK/QAM: $\alpha(n,m,k) = (n-2m)f_c + kf_{symbol}$ . Once you know/estimate $f_c$ and $f_{symbol}$ you know all the CFs.
  
  And i have some question about the step 3.
  If we want to determine Tp，it associates with peak-picking operation.So how can we choose appropriate threshold？
  
  That’s up to you.
  
  It seems the performance of synchronized averaging is related to the data length T.
  If T is not large enough(for example T=1000), will the performance deteriorate?
  If environment is much worse (SNR=-20),will the performance deteriorate?
  
  Yes, all CSP algorithms have performance that increases with increasing averaging time $T$ . See the resolution product post. Yes, all signal processing algorithms have performances that decrease with decreasing SNR.
  
  Does that help?
  
  Reply
Jiang’an Dai says:

May 31, 2021 at 10:28 am

Hello，Mr.Spooner
Thank you very much for your kind attention and reply.
My objective is using T-point data to estimate cyclic cumulant(CC40/60/80). And we don’t the modulation type prior.
I set all time delays(τ1…τn) to 0 for simplicity, then the only variables in CC is cyclic frequency.
So the equation become CC40(α)=CM40(α)-sum[3*CM20(β1)*CM20(β2)], α=β1+β2.

CM20(β1) and CM20(β2) are all T-point.
When compute the term “sum[3*CM20(β1)*CM20(β2)]” , I traverse all T-point optional β1\β2 and find appropriate(α=β1+β2)
to conduct CM20(β1)*CM20(β2).This will take T*T complex-number multiplication. If T is large(65536), the computation time will be very long. And for CC60/CC80, it’s nearly unacceptable for practice application. That‘s the biggest problem bother me.

The cyclic moments and the Fourier transform of lag product are not the same thing.But they have proportional relationship with T. And i estimate cyclic moments by FFT of lag product in practice.（equation 2）

When i don’t know the modulation type prior, how can i throw away all the FFT points corresponding to the theoretical zero terms of CM.

Please correct me if I am wrong.
Thanks,
Jiang’an

Reply
1. Chad Spooner says:
  
  May 31, 2021 at 12:00 pm
  
  You are treating $\beta_1$ and $\beta_2$ as if any possible pair of these values are cycle frequencies for the data. I agree that if you do that any method of using cyclic cumulants will be very costly. In reality only a tiny fraction of possible pairs correspond to actual cycle frequencies exhibited by the data you are processing.
  
  So the answer is to not do that. I suggest that you first figure out how to effectively and efficiently estimate the cyclic cumulants of any order if you know the cycle frequencies of the data. Then work on a strategy for blindly estimating those cycle frequencies, or the underlying key modulation parameters, using an efficient algorithm. Such an algorithm is unlikely to involve sixth- or eighth-order cumulants. Please see My Papers [26]. Also My Papers [25, 28, 43] are relevant.
  
  Reply
  1. Jiang’an Dai says:
    
    June 1, 2021 at 4:09 am
    
    Hello，Mr.Spooner
    Thank you for your kind attention and reply.
    
    In your opinion，if we want to estimate cyclic cumulants effectively and efficiently，it’s better to estimate the cycle frequency first.
    And in your paper[26]，you suggest estimate second-order CFs by thresholding spectral coherence.
    However, from your post ‘The Spectral Coherence Function’，it seems not easy to blindly estimate the cycle frequencies and modulation parameter (carry frequency and symbol rate).
    And this method to estimate cyclic cumulants seems a little complex and depend heavily on the performance of cycle frequency estimate.
    Or that is to say that if we don’t have the prior knowledge of cycle frequency ,the it’s a difficult task to estimate high-order cyclic cumulants.
    
    Please correct me if I am wrong.
    Thanks,
    Jiang’an
    
    Reply
    1. Chad Spooner says:
      
      June 1, 2021 at 5:59 am
      
      if we don’t have the prior knowledge of cycle frequency ,the it’s a difficult task to estimate high-order cyclic cumulants.
      
      It isn’t all that difficult, but it isn’t as easy as you were implying with your approach to computation of the cyclic cumulants.
      To apply CSP to signal analysis, parameter estimation, detection, and modulation recognition requires both understanding of the details of CSP and understanding of the probability structure of mathematical models of RF modulation types.
      
      Reply
      1. Jiang’an Dai says:
        
        June 1, 2021 at 8:48 am
        
        Hello，Mr.Spooner
        Thank you for your kind attention and reply.
        it benefite me a lot.
        
        Best regards,
        Jiang’an
Jiang’an Dai says:

June 2, 2021 at 4:51 am

Hello，Mr.Spooner
Do we need to estimate different Tp for every lower-order temporal moment when using synchronized averaging to estimate higher-order temporal cumulant？
For example，C60=M60-15M20*M40+30M20^3. Do we need to estimate different Tp for M20/M40/M60 respectively？

Thanks
Best regards,
Jiang’an

Reply
1. Chad Spooner says:
  
  June 2, 2021 at 6:24 am
  
  What do you think? Let me ask you: What is the meaning of $T_p$ in synchronized averaging?
  
  Reply
  1. Jiang’an Dai says:
    
    June 2, 2021 at 7:31 am
    
    Hello，Mr.Spooner
    In my opinion,Tp is the estimated period of periodic component in temporal moment.
    And the CF of higher-order temporal moment always have a proportional relation with lower-order temporal moment.
    So i want to know should i estimate Tp of higher-order temporal moment if we know the Tp of lower-order temporal moment.
    
    Thanks
    Best regards,
    Jiang’an
    
    Reply
    1. Chad Spooner says:
      
      June 3, 2021 at 6:25 am
      
      This is essentially an arithmetic problem. Also, “CF” is not the same thing as “Tp.”
      
      Reply
Jiang’an Dai says:

June 3, 2021 at 1:05 am

Hello，Mr.Spooner
“For the synchronized-averaging results here, I always take the delays \tau_j to be zero.
This is fine for most signal types. An exception is those PSK/QAM signals that employ rectangular pulses.
(See the gallery of cyclic correlations for where the cyclic-autocorrelation magnitudes peak as functions of the delay \tau)”

Is that because the cyclic moments equal zero for PSK/QAM signals that employ rectangular pulses when τ=0？
And why they equals zeros？

In figure 14, https://i2.wp.com/cyclostationary.blog/wp-content/uploads/2017/09/bpsk_2_0_ctmf.jpg?resize=733%2C458&ssl=1
the peak value of CTMF at CF=2*fc equals 1.
But in figure 14, https://i0.wp.com/cyclostationary.blog/wp-content/uploads/2017/09/bpsk_2_1_ctmf.jpg?resize=699%2C437&ssl=1
why is the peak value of CTMF at CF=2*fc larger than 1？
And there are same question of the theoretical values in other pictures like figure18\20\22\23\24\25.
In a word,how do we compute the theoretical values of CTMF/CTCF?

Thanks
Best regards,
Jiang’an

Reply
1. Chad Spooner says:
  
  June 3, 2021 at 6:24 am
  
  Is that because the cyclic moments equal zero for PSK/QAM signals that employ rectangular pulses when τ=0？
  And why they equals zeros？
  
  Because that is the way the math works out. I suggest you try to write down the equations for yourself, using a rectangular pulse in PSK and QAM.
  
  why is the peak value of CTMF at CF=2*fc larger than 1？
  
  Because $\alpha = 0$ is a cycle frequency for the noise for $(n, m, k) = (2, 1, 0)$ and $N_0 \neq 0$ .
  
  how do we compute the theoretical values of CTMF/CTCF?
  
  Look at papers in The Literature or my dissertation or the post on the cyclostationarity of PSK/QAM.
  
  Reply
  1. Jiang’an Dai says:
    
    June 3, 2021 at 6:17 pm
    
    Hello，Mr.Spooner
    When i try to compute the cyclic cumulant of PSK/QAM，i find there is an item correlated with the shape pulse.And i don’t know how to compute the integration of it and get the theoretical values of cyclic cumulant/moment.
    
    Thanks
    Best regards,
    Jiang’an
    
    Reply
Jiang’an Dai says:

June 5, 2021 at 1:34 am

Hello，Mr.Spooner
I have programed to estimate time-varying moment/cumulant using synchronized averaging.And i get some result as same as yours like figure 10,12.
But when SNR is low(snr=0) and the order of time-varying moment/cumulant n is high(n=8), the plot of residual power in 12 is similar to random noise sequence.That make it hard to get right Tp.

I really want to using higher-oder CC as discriminating feature to implement blind modulation recognition.
However, i find it’s not easy to estimate the CC especially for high-order.
Both the non-blind and blind estimation method have many limitation.
For non-blind estimation ,many prior knowledge is not suitable for blind modulation recognition.
And for blind estimation, synchronized averaging method seems not work well for higher-oder CC estimation in low SNR.
Now I’m a little pessimistic about using higher-oder CC as discriminating feature to implement blind modulation recognition.
It seems high-order CC is suitable for signal property analysing rather than practical application.

And I want to appologize to you for asking many naive questions and thank you for your kind attention and reply.

Best regards,
Jiang’an

Reply
1. Chad Spooner says:
  
  June 5, 2021 at 6:17 am
  
  How long have you been studying cyclostationary signal processing and cyclostationary random process theory? How long have you tried to construct a blind modulation recognition system based on exploitation of cyclostationarity?
  
  Reply
  1. Jiang’an Dai says:
    
    June 5, 2021 at 7:47 am
    
    Hello，Mr.Spooner
    First and most important,if i have said anything offencive,please accepct my sincerely aplogize.
    I don’t mean to do that abosulutly.
    
    I have read some works about cyclostationarity written by you\Professor Gardner\Giannakis\Napolitano\Dobre.
    And i try using high-order cyclic cumulant to implement blind modulation recognition without prior knowledge just like FAM/SSCA algorithm for second-order cyclostationarity signal processing.
    But i can’t find a feasible way with acceptable computation cost.
    When we only have some modulation data without prior knowledge, how can we get a quick and robust estimation of
    higher-order CC?And how can we implement blind modulation recognition using higher-order CC robustly?
    
    Finally,I want to repeat that if i have said anything offencive,please accepct my sincerely aplogize.
    
    Best regards,
    Jiang’an
    
    Reply
    1. Chad Spooner says:
      
      June 5, 2021 at 9:11 am
      
      You haven’t offended me, but you also did not answer my questions.
      
      blind modulation recognition without prior knowledge just like FAM/SSCA algorithm for second-order cyclostationarity signal processing. But i can’t find a feasible way with acceptable computation cost.
      
      I had previously suggested you look to My Papers [26] for a possible path to doing this, but you’ve rejected that as too difficult (“it seems not easy”).
      
      The creation, reduction to practice, and evaluation of complex signal-processing algorithms is not a “turn the crank” kind of process like machine learning. It takes a lot of mental effort. You have to develop both skill and intuition. To develop intuition means serious engagement with the prior art. To develop skill means serious engagement with the mathematics of both probability theory and communication-signal models.
      
      I can process, completely blindly, a BPSK signal with SNR of 0 dB using maximum cumulant orders of 2, 4, 6, and 8 in 0.6, 0.8, 1.0, and 14.0 seconds, respectively, getting the correct answer, correct parameter estimates, and no false-alarm decisions in all four cases. On a six-year-old laptop using no parallelization. Whether 1.0 seconds is “practical” or not is up to you I suppose. It’s no secret that CSP is computationally costly.
      
      Reply
      1. Jiang’an Dai says:
        
        June 5, 2021 at 6:00 pm
        
        Hello，Mr.Spooner
        Thank you for your kind suggestion and reply.
        I can‘t agree more that complex signal-processing algorithms takes a lot of mental effort. Acturally i have studied in modulation recognition and CSP for more than three years.
        
        For your papers [26], i have some questions.
        1、I‘m not sure the purpose of first two step of the classification algorithm.
        Is that to convert the signal to baseband so we can elimenate the effort of carrier offset frequency?
        2、Using spectral coherence，we can get CF of the signal.But if we want to compute the theory value of CC, there is other parameter we need to know such as the amplitude a and roll-off factor of pulse shaping filter. So how can we get these parameter without prior knowledge.
        
        I’m glad to know that blind estimating of CC can be implement in order of 1 second because it’s
        absolutely acceptable for me. But i want to know which algorithm can get this effect. Is it the synchronized averaging method in this post or the algorithm in you paper[26]?
        
        Best regards,
        Jiang’an
Qin says:

December 24, 2023 at 11:27 pm

Hi, Dr. Spooner
I am new to this field, and I want to thank you for your post. I have learned a lot from it. I have been trying to replicate the method of using synchronous averaging mentioned in your post, but I have encountered some problems. I applied the synchronized averaging operation to the noisy BPSK1 signal, and I was able to achieve similar results to your post when I set fsym=0.1, fc=0.01, and (n, m)=(2,0) as you did. However, when I changed fsym to 1/11, I could not obtain the expected result. I expected to see three spectral lines appearing at 0.02+k/11, where k={-1, 0, 1}, but I only observed one spectral line at the 0.02 position. In fact, I could only achieve the expected results when I set fsym =1/10, 1/20, 1/30…… I am not sure if I have made a mistake somewhere.
Thank you for your help.
Sincerely,
Qin

Reply
jingya ren says:

December 25, 2023 at 11:52 pm

Dear Spooner,
I am new to this field, and I want to thank you for your post. I have learned a lot from it. I have been trying to replicate the method of using synchronous averaging mentioned in your post, but I have encountered some problems. I applied the synchronized averaging operation to the noisy BPSK1 signal, and I was able to achieve similar results to your post when I set fsym=0.1, fc=0.01, and (n, m)=(2,0) as you did. However, when I changed fsym to 1/11, I could not obtain the expected result. I expected to see three spectral lines appearing at 0.02+k/11, where k={-1, 0, 1}, but I only observed one spectral line at the 0.02 position. In fact, I could only achieve the expected results when I set fsym =1/10, 1/20, 1/30…… I am not sure if I have made a mistake somewhere.
Thank you for your help.
Sincerely,
Qin

Reply
1. Chad Spooner says:
  
  December 26, 2023 at 2:27 pm
  
  Qin:
  
  Thanks for reading the CSP Blog and the comment! I really appreciate it.
  
  One of the weaknesses of the synchronized-averaging method is that the true period of the periodic component (or the individual periods in the case of an almost cyclostationary signal) needs to be included in the search window for the period $T_p$ . What was the period for my example in the post where the BPSK signal had symbol rate 1/10 and carrier offset 0.01? For $(n,m) = (2,1)$ is is simply 10, and for $(n,m) = (2,0)$ it is 50 (the period for the doubled carrier is 1/0.02 = 50 and the periods for the doubled carrier plus-and-minus the symbol-rate cycle frequencies $0.02 \pm 1/10$ are 8.333333 and 12.5, which divide 50, so 50 is the common period).
  
  When you switch the symbol rate from 1/10 to 1/11, the period of the periodic component in the $(n,m)$ temporal moment and cumulant functions changes as well. Now the period associated with the doubled carrier is still 50, but the periods associated with the other two cycle frequencies are $|1/(0.02 + 1/11)| = 9.016393442622951$ and $|1/(0.02 - 1/11)| = 14.102564102564102$ . Using MATLAB’s rat.m, I find
  
  [n, d] = rat(14.102564102564102/9.016393442622951) ==> n = 61, d = 39,
  
  [n, d] = rat(50/9.016393442622951) ==> n = 61, d = 11,
  
  [n, d] = rat(50/14.102564102564102) ==> n = 39, d = 11.
  
  So we have
  
  11*50 = 61*9.016393442622951 = 39*14.102564102564102 = 550 samples
  
  If a value for $T_p$ equalling 550 is not included in the search window for $T_p$ (or one of its multiples), then the synchronized averaging algorithm cannot find the period–it will find a period within the search window that corresponds to one or two of the individual periods (arising from the individual cycle frequencies), such as $T_p =$ 50. In that latter case, all you’d see in your output would be the doubled-carrier cycle frequency 0.02 = 1/50.
  
  So the question is: what is your search window?
  
  Using the same search window as I used in the post ([5, 2000]), my code finds a period of 1650 = 3*550 (as in the post, the SA algorithm finds a slightly lower residual error for one of the multiples of the period rather than the period itself, which is just fine). Here are the plots so you can have some confidence in my response:
  
  Reply
  1. jingya ren says:
    
    December 30, 2023 at 12:04 am
    
    Thank you, Professor Spooner, for your guidance. I have successfully obtained the desired outcome. I apologize for the multiple comments I made, as I mistakenly believed that my initial comment did not go through.
    
    Reply
    1. Chad Spooner says:
      
      December 30, 2023 at 8:08 am
      
      Was the search interval for $T_p$ the problem?
      
      Reply
      1. jingya ren says:
        
        January 2, 2024 at 5:05 am
        
        Hi, Dr. Spooner
        Yes, the search interval for Tp is indeed the issue. Once again, I would like to express my gratitude. However, I have discovered another problem. When I increase the search range for Tp significantly, the plot of the error power versus candidate Tp seems to become skewed. If I select the minimum value as Tp, the chosen Tp tends to be biased towards larger values. Additionally, if the length of the signal is relatively short, this skewness appears to be more pronounced. In such cases, how can I accurately search for reasonable Tp values?
        
        Sincerely,
        Qin
      2. Chad Spooner says:
        
        January 2, 2024 at 1:24 pm
        
        When I increase the search range for Tp significantly, the plot of the error power versus candidate Tp seems to become skewed.
        
        Yes, and you can see the slight decrease in the upper values of the error-power graph in my plots as well. This is understandable. Think of the limiting case where we pick $T_p = N$ , where $N$ is the length of the data to be analyzed.
        
        If I select the minimum value as Tp, the chosen Tp tends to be biased towards larger values.
        
        Yes, this happens for the same reason as the gradual decrease in the value of the error as $T_p$ increases.
        
        if the length of the signal is relatively short, this skewness appears to be more pronounced.
        
        Yes, if you have a shorter signal, but keep the same search range for $T_p$ , the decrease in the error value for non-period-related $T_p$ will be worse, because you’ll approach that limiting case $T_p = N$ more quickly.
        
        What can you do? Well, this is where my help ends and your algorithm development begins. One caution is that with synchronized averaging or Fourier methods, you cannot reliably and accurately find a period or tone frequency unless your data record is long in the sense that it contains many periods of the periodicity you seek. (See the Resolution post for related discussions.) So trying to use short data records to do synchronized averaging in the context of CSP will be as difficult as doing normal CSP with short data records (again, “short” here is with respect to the period or periods of the involved cyclic moments and cyclic cumulants, which are the periodic things we seek).
        
        One thing to consider to avoid estimating the period as a multiple of the true minimum period (say, selecting $KT_p$ instead of just $T_p$ ) is to do segmented maximization. Find the min over multiple successive segments, then see if those minima are harmonics of a single number–that would be your $T_p$ estimate.