In this post we discuss ways of estimating -th order cyclic temporal moment and cumulant functions. Recall that for , 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 is greater than . Most communication signal models possess odd-order moments and cumulants that are identically zero, so the first non-trivial order greater than is . So our estimation task is to estimate -th order temporal moment and cumulant functions for using a sampled-data record of length .

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 (excess bandwidth of %). 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 (normalized frequencies are used throughout this post), the carrier offset frequency is , and the total signal power is ( 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 , , and 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 th-order cyclic moment function by using a discrete Fourier transform.

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

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

which is easy to implement in discrete time.

The th-order temporal moment function is simply the Fourier series involving all the non-zero (corresponding to the choices made for , , and ) th-order cyclic temporal moment functions,

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

where the set contains all unique partitions of the index set , (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 by extracting the Fourier-series coefficient for from

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

#### 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 th-order cyclic cumulant is raised to the power (independently of ). 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 th-order estimate and the th-order theoretical cyclic cumulant divided by the maximum theoretical cyclic-cumulant magnitude for all cyclic cumulants for the signal. This measure tends to de-emphasize the errors in the cyclic cumulants for the larger , which have small values relative to cyclic cumulants. Here refers to the harmonic number in the general PSK/QAM cycle-frequency formula , where is the carrier frequency and 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 . That is, those cyclic cumulants corresponding to and . The 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: 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 . So the value of the cyclic cumulant shown in each of the two plots above is correct.

The and 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 cyclic-cumulant magnitude for our 100-percent EBW unit-power BPSK1 is equal to . Therefore, the cyclic-cumulant magnitude for 2-dB BPSK2 is equal to . The sum of the cyclic-cumulant magnitudes is then . The warped version of this magnitude is , which is what we see in the graphs. Similarly, the warped cyclic cumulant should be about .

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 cyclic cumulants. This is because the cyclic moment function for 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 for , 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?

- For , 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 (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 . 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 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 which is defined for all time. For period , the synchronized averaging operation is

for . In (6), is the rectangle function

So, for example, is the average of over all integers . If there is a periodic component of that has period or a period that is a submultiple of , it will survive the averaging operation in (6). If not, the average will be zero.

If the periodic component in is a sine wave, then you’re probably thinking why not just Fourier transform 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 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 th-order temporal moment function is just the (almost) periodic component of a homogeneous nonlinear transform of the signal,

If we apply the synchronized averaging operation to the delay product in (8), ideally we should obtain . 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 , which is presumably related to the cycle frequencies exhibited by , which we are saying we don’t know. So how do we choose ?

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

The power of the error between and can be used to assess whether the signal captures an actual periodic component in 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 (normalized times and frequencies are used). The period is therefore samples, but the signal is also periodic with period samples, samples, etc., for all multiples of .

We apply the synchronized averager to 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 that achieves the minimum error power. Here is the plot of the error power versus candidate :

Notice that the error power (or “residual” power) reaches a minimum of at . This indicates that the full power of the signal component is removed by the periodized waveform using the 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 th-order temporal cumulant functions, we can find the various th-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 to be zero. This is fine for most signal types. An exception is those PSK/QAM signals that employ rectangular pulses.

First up is :

For this input, what **is** the period? It’s not quite as simple as in the sine-wave illustration above. For , the BPSK1 cycle frequencies are for . That set of numbers turns out to be , which correspond to periodicities with periods . Each of these periods divides the candidate period . The synchronized averager shows a nice local minimum at , and of course there are other local minima at multiples of . The algorithm picks , 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 :

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 , 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 , or dB. Also, the bit rate for BPSK1 is and that for BPSK2 is . The carrier offset for BPSK1 is and that for BPSK2 is .

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 .

We can definitely apply synchronized averaging to the 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.

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

Please see the following post for understanding how to get help from the CSP Blog:

https://cyclostationary.blog/2017/04/13/blog-notes-and-how-to-obtain-help-with-your-csp-work

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

Hey Andy. Thanks for looking at that estimation post.

Yes.

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.

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.

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

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

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 () 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 ), you get a single sine-wave at the squarer output. It has frequency equal to the doubled carrier frequency which is . 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 ().

Then I set the delay of one of the factors in the lag product to and kept the other at . The synchronized averaging result is:

See if you can replicate.

And try to understand why the smaller dips happen at multiples of samples. Hint: What are the other cycle frequencies for this signal?

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

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).

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

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 . 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 , we can inverse transform that estimate (over frequency ) 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.

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

I don’t think you are missing anything. The different depths of the dips arises because the different periods have different distances to the true periods .

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.

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.~~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.

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

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 that are integer multiples of the sampling increment.

Well, for QPSK, if we look at the lag product, we see that there are three cycle frequencies: and (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.

If the synchronized averaging works out and we’ve identified the fundamental period for the periodic temporal variation of the 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.