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. That is, the two-dimensional delay vector
is set equal to
.
The more interesting case is when the order is greater than two. Most communication signal models possess odd-order moments and cumulants that are identically zero, so the first non-trivial order
greater than two is four. 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 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 (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 frequency offset 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 frequency offset, 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. Let me elaborate on this idea.
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 integral can be replaced by a sum, and the cycle-frequency parameter 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 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 nth-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 offset 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 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:
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. (Note that the warping is inconsequential for
because
.)
The and
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
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 in Figures 8 and 9.


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 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 for
, 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
, 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 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 (which is 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 us to select the parameter , which is presumably related to the cycle frequencies exhibited by
, which we are saying we don’t know (we’re in the blind processing context). 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
. A 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. If
represents an actual periodic component in
, then subtracting its periodized version
from
will produce a signal with smaller power than
. 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 (normalized times and frequencies are used, as nearly always at the CSP Blog). The period is therefore
samples, but the signal is also periodic with period
samples,
samples, or generally
samples for all integers
.
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.
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 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
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
)
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
in units of samples. 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 of the signal in Figure 13 leads to the cyclic moments for :

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 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 . Comparison with directly estimated cyclic moment functions confirms the efficacy of synchronized averaging.








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



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 , 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 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.
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:
http://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
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
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.
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
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
and we are using a sampling rate of one. So we might estimate the period as
,
, or
when things are going well. The corresponding frequency estimate is
.
The error in the frequency estimate is
or
Supposing the estimate is good, we have
The normalized error is
So the resolution depends on the period length
. This is why oversampling helps in synchronized averaging.
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
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,
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.
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
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.
No. I do not assume they have the same carrier phase.
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.
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.
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
and
, 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.
Dr.Spooner Thanks for your replies. Very helpful.
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
Welcome to the CSP Blog Jiang’an!
I don’t think your equation here is consistent with (5). The equation you wrote here is an oversimplification. You must obey (5).
The multiplications in (5) are simple–note that there is no time
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
,
and
. So you have to find those in the
FFT values.
It looks like you are mistaking the Fourier transform of a nonlinear function of
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
. Remember the formula for the CFs for PSK/QAM:
. Once you know/estimate
and
you know all the CFs.
That’s up to you.
Yes, all CSP algorithms have performance that increases with increasing averaging time
. See the resolution product post. Yes, all signal processing algorithms have performances that decrease with decreasing SNR.
Does that help?
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
You are treating
and
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.
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
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.
Hello,Mr.Spooner
Thank you for your kind attention and reply.
it benefite me a lot.
Best regards,
Jiang’an
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
What do you think? Let me ask you: What is the meaning of
in synchronized averaging?
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
This is essentially an arithmetic problem. Also, “CF” is not the same thing as “Tp.”
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
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.
Because
is a cycle frequency for the noise for
and
.
Look at papers in The Literature or my dissertation or the post on the cyclostationarity of PSK/QAM.
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
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
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?
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
You haven’t offended me, but you also did not answer my questions.
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.
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