The Cyclic Autocorrelation Function

The cyclic autocorrelation function is the amplitude of a Fourier-series component of the time-varying autocorrelation for a cyclostationary signal.

In this post, I introduce the cyclic autocorrelation function (CAF). The easiest way to do this is to first review the conventional autocorrelation function. Suppose we have a complex-valued signal x(t) defined on a suitable probability space. Then the mean value of x(t) is given by

M_x(t, \tau) = E[x(t + \tau)]. \hfill (1)

For stationary signals, and many cyclostationary signals, this mean value is independent of the lag parameter \tau, so that

\displaystyle M_x(t, \tau_1) = M_x(t, \tau_2) = M_x(t, 0) = M_x. \hfill (2)

The autocorrelation function is the correlation between the random variables corresponding to two time instants of the random signal, or

\displaystyle R_x(t_1, t_2) = E[x(t_1)x^*(t_2)]. \hfill (3)

To see how the autocorrelation varies with some particular central time t, we can use a more convenient parameterization of the two time instants t_1 and t_2, such as

\displaystyle R_x(t, \tau) = E[x(t+\tau/2)x^*(t-\tau/2)], \hfill (4)

where

\displaystyle t_1 = t+\tau/2

\displaystyle t_2 = t-\tau/2

and

\displaystyle t = (t_1 + t_2)/2

\displaystyle \tau = t_1 - t_2

So time t represents the center point of the two time instants involved in the correlation and \tau is their separation. If the autocorrelation depends only on the separation between the two time instants \tau, and not their center point t, the signal is stationary of order two, or just stationary, and we have

\displaystyle R_x(t_1, \tau) = R_x(t_2, \tau) = R_x(0, \tau) = R_x(\tau). \hfill (5)

This is a convenient way of saying that the autocorrelation E[x(t_1)x^*(t_2)] depends only on the difference between t1 and t2, and not on their midpoint.

For nonstationary signals, on the other hand, the autocorrelation does depend on central time t. For the special case of nonstationary signals called cyclostationary signals, the autocorrelation is either a periodic function or an almost periodic function. In either case, it can be represented by a Fourier series

\displaystyle R_x(t, \tau) = \sum_\alpha R_x^\alpha (\tau) e^{i 2 \pi \alpha t}, \hfill (6)

where R_x^\alpha(\tau) is a Fourier-series coefficient called the cyclic autocorrelation function. The Fourier frequencies \alpha are called cycle frequencies (CFs). The cyclic autocorrelation functions are obtained in the usual way for Fourier coefficients,

\displaystyle  R_x^\alpha(\tau) = \lim_{T\rightarrow\infty} \frac{1}{T} \int_{-T/2}^{T/2} R_x(t,\tau) e^{-i 2 \pi \alpha t}\,dt. \hfill (7)

If the signal is a cycloergodic signal (or we are using fraction-of-time probability), then the cyclic autocorrelation function can be obtained directly from a sample path of the random process (the observed signal itself),

\displaystyle R_x^\alpha(\tau) = \lim_{T\rightarrow\infty} \frac{1}{T} \int_{-T/2}^{T/2} x(t+\tau/2) x^* (t-\tau/2) e^{-i 2 \pi \alpha t} \, dt. \hfill (8)

For many cyclostationary signals, such as BPSK, the conjugate autocorrelation function is also non-zero (and is also useful in practice). This function is defined by

\displaystyle R_{x^*}(t, \tau) = E[x(t+\tau/2)x(t-\tau/2)], \hfill (9)

and is represented by its own Fourier series

\displaystyle R_{x^*}(t, \tau) = \sum_\alpha R_{x^*}^\alpha (\tau) e^{i 2 \pi \alpha t}, \hfill (10)

where R_{x^*}^\alpha(\tau) is the conjugate cyclic autocorrelation function and here \alpha is a conjugate cycle frequency. I explain in detail why we need two autocorrelation functions in the post on conjugation configurations. The problem is worse when we look at higher-order moments and cumulants, where we need n/2 + 1 functions to properly characterize a signal at order n.

Symmetric versus Asymmetric Cyclic Autocorrelation Functions

The conventional autocorrelation and cyclic autocorrelation functions are symmetric in the delay variable \tau in that it appears as +\tau/2 in one factor of the delay product and as -\tau/2 in the other:

\displaystyle x(t+\tau/2)x^*(t-\tau/2).

When attempting to estimate the cyclic autocorrelation function, using discrete-time data, the lag variable \tau and the time index variable t take on integer values, so that \tau/2 does not correspond to a known value of x(t) when \tau is odd. To work around this practical problem, we can operate in the frequency domain, directly estimating the spectral correlation function and inverse Fourier transforming it to obtain the cyclic autocorrelation. But we can also employ the asymmetric cyclic autocorrelation, which is friendly to discrete-time data, and then scale it appropriately to obtain the symmetric version. Let’s go through the supporting analysis here.

The asymmetric cyclic autocorrelation is

\displaystyle \tilde{R}_x^\alpha(\tau) = \lim_{T\rightarrow\infty} \frac{1}{T} \int_{-T/2}^{T/2} x(t) x^*(t-\tau) e^{-i2\pi \alpha t} \, dt, \hfill (11)

but it could be defined in other ways, including

\displaystyle \tilde{R}_x^\alpha(\tau) = \lim_{T\rightarrow\infty} \frac{1}{T} \int_{-T/2}^{T/2} x(t+\tau) x^*(t) e^{-i2\pi \alpha t} \, dt. \hfill (12)

Staying with definition (11), and using the change of variables t_0 = t - \tau/2, we have

\displaystyle \tilde{R}_x^\alpha(\tau) = \lim_{T\rightarrow\infty} \frac{1}{T} \int_{-T/2-\tau/2}^{T/2 - \tau/2} x(t_0+\tau/2) x^*(t_0-\tau/2) e^{-i2\pi \alpha (t_0 + \tau/2)} \, dt_0, \hfill (13)

or

\displaystyle \tilde{R}_x^\alpha(\tau) = e^{-i\pi\alpha\tau} R_x^\alpha(\tau) \hfill (14)

so that we can compute the asymmetric version and then scale it by a simple complex exponential to obtain the symmetric version.

And that is it! The definitions in this post are the basic definitions for the (second-order) probabilistic parameters of cyclostationary signals in the time-domain. In other posts, I have much to say about their utility, their estimation,  their connection to the frequency-domain parameters, and their generalization to higher-order parameters.

Support the CSP Blog here.

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.

75 thoughts on “The Cyclic Autocorrelation Function”

  2. Hi Dr. Spooner,

    Thank you for making this blog. Are there certain qualities of a signal that are more noticeable using the Cyclic Autocorrelation Function (CAF) versus the Spectral Correlation Function (SCF)? In other words, are there situations where the usage of the CAF is favorable over the SCF?

    Thank you,
    Abdul

    Reply

      1. Haha, my apologies. I initially expected my comment to be viewable immediately. I don’t normally use wordpress and didn’t realize that my comment was under moderation until the second time I tried to comment.

        Reply

  3. Hi Dr. Spooner,

    Thank you for this blog. Are there situations where it would be advantageous to utilize the Cyclic Autocorrelation Function over the Spectral Correlation Function?

    Thanks

    Reply

    1. I’ve found the CAF to be useful when doing CSP for OFDM signals (see the work of Octavia Dobre).

      I think the general answer to your question is that the advantages of the two functions depend on how much prior information you have about your signals. When you don’t know anything, and you want to do RFSA, the SCF is quite useful because of FFT-based estimators like the SSCA and FAM. But if you know quite a lot, such as a cycle frequency and some CAF lags of interest, then just estimating the CAF over a limited range of lags can be quite inexpensive.

      Reply

      1. Thank you for your response. I’ll take a look at that article and see how they used the CAF for OFDM signal detection. I’ve been meaning to learn about OFDM signals, so this paper will be a good motivator.

        So in general, if we know certain characteristics of our signal (e.g. cycle frequencies, lags of interest, etc.), then the CAF would be a computationally inexpensive means of calculating Second Order Cyclostationary Features. Makes sense, thanks for the clarification!

        Reply

  4. Hi Dr. Spooner. Great blog. I know for ergodic signals, time average equals to mean value. So there should be two integral symbols in equation 8. Why do you delete one?
    Thank you.

    Reply

    1. Thanks Liang. We refer to “cycloergodicity” here, which means that ensemble averages in the stochastic framework equal the output of the sine-wave extraction operator (the expected value in the fraction-of-time probability framework):

      E[F(X(t))] = E^{\{\alpha\}}[F(x(t))]

      where X(t) is the random process and x(t) is a sample path thereof. When cycloergodicity holds, we can obtain the cyclic autocorrelation function directly from the second-order lag product for almost all sample paths (almost all == with probability one).

      Agree?

      Reply

      1. What does F mean? And alpha? Can you give me some reference papers for the proved equation. I know little about cycloergodicity. Thanks.

        Reply

        1. F is just some functional, like F(x(t)) = x(t+\tau/2)x^*(t-\tau/2). \alpha denotes a cycle frequency. Or, more generally, the frequency of a finite-strength additive sine-wave component in F(x(t)). For information on cycloergodicity and the fraction-of-time probability framework, see The Literature [R8, R67, R68].

          Reply

      2. So for equation 8, can I calculate CAF with fft[x*conj(x)]? Because I find it has same form as DFT operation.

        Reply

        1. Well, you can estimate the cyclic autocorrelation by picking out one element of the FFT of the lag product. But notice that the cyclic autocorrelation function is a limiting version of a time average. Because there is no guarantee that the cycle frequency \alpha is exactly equal to a Fourier frequency in the DFT, you’ll get better results by just computing the single DFT directly. That works fine if you know the cycle frequency \alpha in advance. If you don’t, you have to search for the cycle frequencies, and that is best done using the spectral correlation function and the strip spectral correlation analyzer.

          Reply

  5. Hello Sir,
    Thank you so much for your valuable information. I just wanted to ask how to cite this tutorial?
    Thanks in advance

    Reply

    1. Hey Sara. Do you mean to reference the CSP Blog, or a post within it, in the reference list of a published paper?

      I’ve referenced it in a journal paper and a couple of conference papers. I just use the URL. So if you want to reference the CSP Blog in general, you could include an entry in your reference list that looks like this:

      [1] C. M. Spooner, The Cyclostationary Signal Processing Blog, https://cyclostationary.blog.

      or even

      [1] https://cyclostationary.blog.

      If you wanted to reference a particular post within the CSP Blog, you could get the URL from your browser. For the post you’ve commented on, it would be

      [2] https://cyclostationary.blog/2015/09/28/the-cyclic-autocorrelation.

      or

      [2] C. M. Spooner, The Cyclostationary Signal Processing Blog, https://cyclostationary.blog/2015/09/28/the-cyclic-autocorrelation.

      Does that answer your question?

      I would appreciate very much this kind of citation. It will help spread the word on the Blog, which I consider a valid alternative to published papers in conventional journals and conference proceedings.

      Reply

  6. Dear Dr. Spooner,

    Many thanks for such an extensive blog on cyclostationary.

    I am new to this cyclostationary analysis, though I have been reading through many your posts. In particular the one on cyclic autocorrelation.

    You wrote in one of your replies that for OFDM signal CAF be used especially when the properties of the signal are known.

    I also read several publications from Octavia Dobre, as you suggested. One of her recent publications (see the link below) is about identification of GSM and LTE signals.

    https://ieeexplore.ieee.org/document/7151426

    I have downloaded the sample LTE signal from your “Data Sets” and used it for testing the algorithm from the paper above. The algorithm seemed to be straightforward, however I didn’t obtain the expected results. Basically I computed the CCF at CF alpha. The CF for LTE is known, i.e. 2 kHz (corresponds to 0.5 ms time slot).

    I would like to know if the LTE sample file in your data set is based on OTA measurement or a simulated signal?

    Do you think I should better use SCF instead of CAF for this purpose?

    Thanks in advance.

    Best regards,
    Johann

    Reply

  7. Hi Dr. Spooner,

    Thank you for your prompt reply. I see, I will stick with the CAF then.

    I computed the CAF with a zero delay only, as described in the paper. Maybe I need to include a range of delay values, as you suggested?

    So far I haven’t validated my CAF estimator with a simpler simulated signal. I guess I should do that first.

    I will keep you updated. Thanks again for your help.

    Best regards,
    Johann

    Reply

    1. Thank You Dr Spooner for such a informative and guiding blog. My questio is to Mr. Johann if he can his codes with me I am also trying to work On LTE signals.

      Reply

      1. Rinkoo:

        I’ve approved your comment, and good luck with Johann. However, my advice is to write your own code. You will be able to understand it much better, and apply it properly.

        Reply

  8. Dear Dr. Spooner,

    I managed to get similar results to the paper using your LTE and GSM data sets. There was a bug in my original code.

    The corresponding CAF plots can be found on the links below. It’s my first time using imgur, so I don’t know for how long the links will be available.

    https://imgur.com/a/gs8a0wi

    https://imgur.com/a/9Anibdx

    For your information, I am still using a zero delay. The CFs for LTE and GSM are equal to the multiple of 2 kHz and 1733 Hz, respectively.

    Could you elaborate why you think using a range of delay values might be better?
    Thanks in advance.

    Best regards,
    Johann

    Reply

  9. Hi Dr. Spooner,

    Many thanks for your reply and the nice plots.

    I will try to replicate your plots.

    Best regards,
    Johann

    Reply

  10. Hi sir, nice blog for beginners to study more about cyclostationary processing. I am new to this cyclostationary processing. I tried to write a matlab code to find cyclic auto correlation function, Here are my observations please help me that i understood this concept or not.

    I generated one sine wave of 1000 samples with freq 10 and sampling freq 1000 (f/fs = 0.01)
    alpha = cyclic frequency , Tau = time shift , N = 1000

    x = sin(2*pi*f/fs*(0:N-1))

    for alpha = 0:10
    carrier = exp(-1i*2*pi*alpha*(1/fs)*(0:N-1)) ;

    for tau = 1:980
    M = N-tau;
    Rx(alpha+1,tau) = 1/M*(sum( x(1:N-tau).*x(tau+1:N).*carrier(1:M) )) ;
    end

    plot(abs(fft(Rx(alpha+1,:)))) ;
    figure

    end

    I have written the above matlab code , For each alpha i am getting the peaks , i am expecting only one peak at alpha = 0 but at different alpha’s also i am getting peaks with less amplitude values compare to alpha = 0 .

    So please suggest me am i doing correct or not?
    any way thanks a lot for ur nice blog

    Reply

    1. Thanks for stopping by the CSP Blog and leaving a comment Hari. I really appreciate it.

      You should expect the cyclic autocorrelation to be non-zero even when the cycle frequency you use does not correspond to a cycle frequency of the signal.

      A couple things stand out in your code. The first is the use of a relatively short data-record length N. Try increasing it and you’ll see that most of the peaks get very small. Second, you are not quite implementing the cyclic autocorrelation as I define it in the CSP Blog (which is the most common definition) in that you are shifting up and down by \tau, not by \tau/2. When you do shift up and down by \tau, you are then dividing by N-\tau and not N-2\tau. Finally, you are not considering the actual cycle frequencies for your sinusoidal signal! Real-valued sine waves have cycle frequencies equal to zero, 2f_c and -2f_c, where f_c here is your f. So your doubled-carrier cycle frequency would be 20, but you go up to 10 only.

      But why bother with a sine-wave input? I suggest you use my rectangular-pulse BPSK signal and my second-order verification guide. That way, you know exactly what to expect.

      Reply

      1. Thanks a lot for ur quick reply sir. From ur above suggestion , i am not getting this statement

        Finally, you are not considering the actual cycle frequencies for your sinusoidal signal! Real-valued sine waves have cycle frequencies equal to zero, 2f_c and -2f_c, where f_c here is your f. So your doubled-carrier cycle frequency would be 20, but you go up to 10 only.

        could u please elaborate it ?

        Reply

        1. If you have a sine wave x(t) = e^{i 2 \pi f_c t} and you form x(t+\tau/2)x^*(t-\tau/2) you get a factor of e^{i 2 \pi (2 f_c) t}. This indicates that the signal has a cycle frequency of \alpha = 2f_c. Now do the same thing for a real-valued sine wave \sin(2 \pi f_c t) using Euler’s formula and you’ll see that you get components with frequencies 2f_c and -2f_c. Sine waves are trivially cyclostationary in that they are non-random, but a quadratic functional yields sine-wave components with frequencies not found in the sine wave itself. For a more detailed explanation of all this, see the post on cyclic cumulants.

          Reply

          1. Many thanks for ur reply sir, After analyzing ur above suggestions ,i have some doubts ,please calrify it

            A couple things stand out in your code. The first is the use of a relatively short data-record length N. Try increasing it and you’ll see that most of the peaks get very small. Second, you are not quite implementing the cyclic autocorrelation as I define it in the CSP Blog (which is the most common definition) in that you are shifting up and down by \tau, not by \tau/2. When you do shift up and down by \tau, you are then dividing by N-\tau and not N-2\tau.

            1) After increasing the data record length N , I got very small peaks ,thanks a lot for that suggestion

            2)As u mentioned that i am not implemented cyclic autocorrelation as defined in CSP blog in that \tau shifting up and down.
            Is any wrong with this? if so where will it get affected if i am not using \tau/2.
            As per my knowledge with respect to central time (t) i am using the tau delayed sample to perform correlation.
            To perform \tau/2 based cyclic auto correlation in matlab ,it throws error since \tau/2 is fraction(if tau=1).

            3)Regarding the division factor ( N – \tau). since i am using N- tau elements for that operation ,so i am dividing with N-\tau. but u suggested to divide with (N – 2tau) .Is any reason behind that?

            So please give ur time to clarify this doubts

          2. i am not implemented cyclic autocorrelation as defined in CSP blog in that \tau shifting up and down.
            Is any wrong with this? if so where will it get affected if i am not using \tau/2.

            Well, by using that particular definition of a cyclic autocorrelation, your results will no longer be consistent with the literature and with the CSP Blog. In particular, if you Fourier transform your cyclic autocorrelation, you will not get the spectral correlation function as commonly defined. By shifting up and down by \tau (instead of the usual definition of shifting up and down by \tau/2), you effectively stretch the function in the lag variable–you are scaling time.

            To perform \tau/2 based cyclic auto correlation in matlab ,it throws error since \tau/2 is fraction(if tau=1).

            Yes, that is true. One way around this is to use the asymmetric definition that uses the following delay product:

            x(t+\tau)x^*(t)

            or

            x(t) x^*(t-\tau)

            The relationship between the conventional (symmetric in \tau) definition and the asymmetric definition is that they are the same function except for a complex exponential scaling factor related to \alpha and $\latex \tau/2$. See the comments on the cyclic autocorrelation posts here and here and here for further hints or, better, derive it yourself.

            Regarding the division factor ( N – \tau). since i am using N- tau elements for that operation ,so i am dividing with N-\tau. but u suggested to divide with (N – 2tau) .Is any reason behind that?

            OK, maybe I am wrong, but if you shift one x(t) up by \tau and the other by -\tau, the overlapping region has width N-2\tau.

            Actually, I just divide the sum by N no matter what \tau is because the range of \tau for which the function is significant is usually much much smaller than N.

  11. Hi Chad,

    I’ve tried to implement the CAF in python following pseudo-code in Professor Napolitano’s newest book page 209. I am using a BPSK like you recommend. Below is an excerpt of the CAF computation loop. The non-conjugate CAF this produces is comparable to the one presented in the BPSK CAF post. But the conjugate CAF is not. I did notice in the pseudo-code in the book uses exp(i*2*pi*alpha*N) but the formula in the CAF post uses exp(-1*i*2*pi*alpha*N). Meaning it is up shifting the frequency rather than down shifting like the formula suggests. Could this be an error? Or maybe my implementation is wrong. Any help or comment on this will be much appreciated. Thanks.

    # x_of_t is the BPSK time series

    N = np.arange(0, x_of_t.shape[0])

    ## Cyclic Autocorrelation function

    fs = 1 # Sampling Freq
    Ts = 1/fs # Define sampling period (s)

    alphaMax = 1.0 #Maximum cyclic frequency
    dAlpha = 0.1 #step in cyclic frequency

    alpha = np.arange(0, alphaMax, dAlpha) # Define alpha

    # Preallocation
    RxTauAlpha = np.zeros((alpha.shape[0], x_of_t.shape[0]), dtype=np.complex128)

    RxTauAlpha_conj = np.zeros((alpha.shape[0], x_of_t.shape[0]), dtype=np.complex128)

    cnt = 0

    for alpa in alpha:
    RxTauAlpha[cnt, :] = np.correlate(x_of_t, x_of_t * np.exp(1j * 2 * np.pi * alpa * N), mode=’same’)
    RxTauAlpha_conj[cnt, :] = np.correlate(x_of_t, np.conj(x_of_t) * np.exp(1j * 2 * np.pi * alpa * N), mode=’same’)
    cnt+=1

    Reply

    1. Thanks for the comment Joel. Maybe the following responses will help; let me know.

      I am using a BPSK like you recommend.

      Are you using my BPSK signal? Or your own? I ask because I want to know the carrier frequency.

      But the conjugate CAF is not.

      In what way does your conjugate CAF not match mine?

      I did notice in the pseudo-code in the book uses exp(i*2*pi*alpha*N) but the formula in the CAF post uses exp(-1*i*2*pi*alpha*N). Meaning it is up shifting the frequency rather than down shifting like the formula suggests. Could this be an error?

      I don’t think so. He is multiplying by e^{i 2 \pi \alpha t} instead of e^{-i2\pi\alpha t} because he is using MATLAB’s xcorr.m, which, like your np.correlate, conjugates its second argument. So at that point the sign of the exponent for e is reversed. My formulas are more straightforward in that they do not hide any details behind a function call, so you see the correlation itself in detail, which requires the exponent to be negative.

      alphaMax = 1.0 #Maximum cyclic frequency
      dAlpha = 0.1 #step in cyclic frequency

      alpha = np.arange(0, alphaMax, dAlpha) # Define alpha

      You are estimating the non-conjugate CAF and the conjugate CAF for the same set of cycle frequencies, which is equal to \{0, f_{bit}, 2f_{bit}, \ldots\}. This is sufficient for the non-conjugate CAF, but leaves out cycle frequencies for the conjugate CAF. Their symmetries are different, and in general one wants to compute the conjugate CAF (or conjugate SCF) for cycle frequencies stretching from -1 to 1.
      That’s why I asked about your carrier frequency at the top of this reply.

      Reply

      1. Hi Chad,

        Yes I am using your BPSK signal. Thanks for clarifying that it was no an error. Below is the google colab link where you can view, run, comment and edit the code I’ve used. Sticking to Prof. Napolitano’s method for now, hope that does no annoy you too much.
        I am saying the conjugate CAF is wrong because the non-conjugate CAF at alpha=0 (figure 1 in symmetry post) and conjugate CAF at alpha=0.1 are the same. I am trying to get my head around what is wrong. Thanks.

        https://colab.research.google.com/drive/1jdVJsyzSMlwEzKnm7KGftxDsStiIp3SE

        Reply

        1. I don’t think you are doing anything wrong, per se, and I hope I didn’t come across as annoyed in my previous reply.

          Mostly we just have to interpret your results correctly, especially if/when comparing to my results in the Symmetries post.

          In particular, you are using numpy’s correlate function np.correlate, which is an asymmetrical correlation (fine!), but is a slightly different asymmetric correlation than the one I consider in this post on the cyclic autocorrelation. If you want to get to the symmetric cyclic autocorrelation, you’ll have to do something similar to what I explain above in the section called “Symmetric versus Asymmetric Cyclic Autocorrelation Functions.” I tried to make the modifications to your collaboration worksheet, but I’m not great at python. Check out my additions to your code and the section above and adjust as needed. You know you are on the right track when the magnitudes of the functions match mine. The particular correlation you use and the asymmetric compensation will determine whether the real and imag parts of your estimates match mine or not.

          Reply

        2. Hi Joel,

          I’d love to look at your code as a reference for a project I’m working on at the moment. The colab link is down, and I was wondering if you’d moved it and if so could update me with the new repo?

          Thanks in advance,
          Billy

          Reply

      1. Happens to me all the time! At work the other day, your website was mentioned by another engineer as the definitive resource for learning cyclostationary analysis. You probably don’t hear much about the word of mouth traffic you get to your blog but people are talking about it.

        Reply

        1. Very nice to hear! I get a lot of “thanks for the blog,” which I very much appreciate, but of course I don’t hear what people are saying to each other.

          I haven’t forgotten your suggestion about mutual guest-authoring, by the way… Just busy.

          Reply

          1. I always like to say, “better busy than bored”. You do a lot of good work so I’m not surprised you don’t have much time. Anytime you want to do something like that, I’m all ears.

  13. Hello,

    I implemented in Matlab cyclic autocorrelation using the A. Napolitano book. The code is presented below. I have one question, because I am not an expert in signal processing. N stands for the number of DFT points. It is better to use N equal to length of the signal and plot the map with only some of the alphas? Or chose some smaller N ?

    function [R]=cyclic_cross_correlogram(y,x,taumax,N);
    R=zeros(N,2*taumax+1);
    for tau=-taumax:taumax
    R(1:N,tau+taumax+1)=fftshift(fft(wshift(1,y,tau).*conj(x),N))/N;
    end
    end

    Reply

    1. Hey Piotr! Thanks for stopping by and asking a question.

      I don’t know the answer to your question as you have written it. What are you trying to do with the cyclic cross correlation? Search for cycle frequencies? Or do you have some cycle frequencies that you already know and you want to compute the cyclic cross correlation for those cycle frequencies?

      If you use N samples in your code above, but then only plot or examine some subset of size M < N, you could miss seeing a significant cycle frequency. The cycle-frequency resolution of your measurement is approximately 1/N, and you are computing estimates for cycle frequency alpha on a grid with grid spacing of 1/N, so you won’t miss any cycle frequencies in the calculation itself, for any N. However, for smaller values of N, the grid spacing is larger (1/N becomes larger), so you will have a less precise estimate of the cycle frequencies that do give rise to large cyclic cross correlations than if you used a larger value of N.

      Does that answer your question?

      Reply

      1. I appreciate that you respond me so quickly. Thank you for this great blog, I ve just started so probably I will have more questions related to different posts.

        I am mostly interested in the case with unknown cycle frequencies and want to search them. What is your advice? When I plotted it for all N frequency points it was difficult to find any useful information from the plot.

        Reply

  14. Hi Chad, I try to achieve the CAF of UFMC signal, I have proved that my CAF function is right
    when I use your BPSK signal. But I can not get correct the X-axis(cycle frequencies axis) and Y-axis(time delay) for UFMCsignal, Could you recommend the parameter settings of cycle frequencies and time delay?

    Reply

  15. Dear Dr. Spooner,

    Thank you for making this great blog,
    I’m new for Cyclostationary Signal and now i have to find the autocorrelation of Wide-Sense Cyclostationary Signals which are realistic for communication signals used in wireline.
    Do you have any suggestion to find the autocorrelation function of the PAM signals such as PAM2, PAM4? or related papers
    Thank you very much!

    Reply

    1. Hey Yugi! Thanks for checking out the CSP Blog. PAM signals are a special case of general complex-valued digital QAM signals, and so the formulas for the spectral correlation function, cyclic autocorrelation, and their degenerate special cases of the power spectrum and autocorrelation are well known and widely published. In The Literature, see Gardner’s work [R1], [R47], and [R145]. The latter is tutorial. In the CSP Blog, I’ve summarized the higher-order parameters, which degenerate to the second-order parameters you’re looking for, in the posts on digital QAM/PSK and the cyclic polyspectrum.

      Reply

  16. Hi Dr. Spooner. Great blog. very helpful, i appreciate.

    Please could you explain how to calculate the cyclic frequencies for the WBFM modulation, for both type of input signal ( sinusoidal signal and aleatory signal), for witch freq should give us the picks.

    Thanks a lot

    Reply

    1. A frequency modulated signal with a sinusoidal message is a signal with a purely discrete spectrum–it is nothing more than the sum of sinusoids. A complex-valued sinusoid has no non-conjugate cycle frequencies and a single conjugate cycle frequency equal to twice the sine-wave frequency. A frequency modulated signal with a random (aleatory) message is typically not a cyclostationary signal. The degree to which it possesses cycle frequencies depends on the modulation index, but for all practical purposes it does not possess exploitable cyclostationarity. See The Literature [R1] and [R46].

      Reply

  17. Hello Dr. Spooner,

    I am having trouble grasping a good elementary understanding what is going on.

    So my understanding is we take the autocorrelation of a signal, and then we take something similar to a Fourier Transform of the output of the autocorrelation.

    So if I have a BPSK signal for example, and Ii take the auto correlation, I will get an periodic waveform based on the BPSK chip length. Then I take the Fourier transform of different (cyclic?) frequencies; does this pull out and highlight the chip width at the output? So the shorter the chip width, the higher rate the periodic output of the autocorrelation function will be. So in that case, we should see more energy in the higher cyclic frequency?

    Thanks!

    Reply

    1. Thanks for stopping by Vince!

      So my understanding is we take the autocorrelation of a signal,

      Yes, the autocorrelation function is time-varying, and so it has a Fourier series representation as all periodic signals do. You can look at its Fourier frequencies (which we rename cycle frequencies in the CSP context) and its Fourier coefficients (which we rename as the cyclic autocorrelation in the CSP context). But there is a lot packed into “we take the autocorrelation.” On paper, we can “take” the autocorrelation of a random process using a probability model and the expectation operator. When we have a signal at hand, though, we don’t have a process, so we can’t apply the expectation. That’s why I bring up the cycloergodicity idea in the post. If time-averages and ensemble averages are equal (with probability one as they say), then you can find the Fourier coefficients of the time-varying autocorrelation by directly operating on the time-series at hand. If you want, you can compute them all (use all the different real numbers for the alphas), and then construct the time-varying autocorrelation through the Fourier-series expression. But usually we don’t want to do that; we’re just interested in the cycle frequencies and the cyclic autocorrelation functions, or, much more often, the spectral correlation functions.

      I do discuss how to estimate the time-varying autocorrelation and time-varying higher-order moments and cumulants in another post.

      The elementary understanding I’m trying to get across is: You can find the Fourier-series coefficients of the time-varying autocorrelation by directly operating on data using a finite-T version of (8).

      So if I have a BPSK signal for example, and Ii take the auto correlation, I will get an periodic waveform based on the BPSK chip length.

      Yes, as long as we remember what “take the autocorrelation” means.

      Then I take the Fourier transform of different (cyclic?) frequencies; does this pull out and highlight the chip width at the output? So the shorter the chip width, the higher rate the periodic output of the autocorrelation function will be. So in that case, we should see more energy in the higher cyclic frequency?

      Yes, if you Fourier transform the time-varying autocorrelation, you’ll see peaks at Fourier frequencies equal to the cycle frequencies, which for BPSK will be harmonics of the bit rate (I only use “chip” in the context of DSSS signals). (There is a complication here because you also need to consider the conjugate autocorrelation, and BPSK gives rise to a different set of cycle frequencies for the conjugate compared to the non-conjugate.) And yes, the shorter the bit interval is, the larger the bit rate is, and since one of the cycle frequencies is the bit rate, it will become larger as the bit interval becomes smaller.

      Reply

  18. Hi Spooner,

    In my code I use the FAM method to estimate the SCD matrix. But after the second FFT. I have no idea about how to arrange the result of second FFT to SCD matrix. And I try some paper but only symmetry in Alpha. Could you please give me some reference about how to symmetry in both alpha and frequency.

    Best Regards.

    Reply

    1. Hey Caro! Thanks for stopping by the CSP Blog.

      What have you tried? Have you tried to use (7) and (8) in Step 5 of the FFT Accumulation Post?

      Also, I don’t quite understand your questions and comments because of your English grammar problems. Can you more carefully rephrase?

      Reply

  19. Hi, Dr Spooner, I want to ask: what do cyclic frequencies(CFs) represent, cause it is a little bit abstract.

    Reply

    1. I’m not sure what you mean by “represent” here. But … the cycle frequencies are the Fourier-series frequencies for the representation of the time-varying autocorrelation as a Fourier series. So each cycle frequency, from this point of view, is a frequency of a sine-wave component in the periodic (or almost periodic) time-varying autocorrelation function for the cyclostationary signal.

      The cyclic autocorrelation function is the inverse Fourier transform of the spectral correlation function. The spectral correlation function S_x^\alpha(f) is interpreted as the density of correlation between two narrowband spectral components of a signal or of some data. Only certain pairs of such narrowband spectral components are correlated, and the separations between the center frequencies of these pairs are exactly equal to the cycle frequencies. So each cycle frequency, from this point of view, is the separation between narrowband signal components that contain redundant (correlated) information.

      None of these mathematical interpretations mention the connection of the cycle frequencies to parameters of actual communication signals. That would be a question less like “what do the CFs represent” and more like “how are CFs and modulation parameters connected?” For that latter question, see some of my other posts on the cyclostationarity of various broad signal classes, such as QAM, PSK, and DSSS.

      Reply

  21. I feel like I’m overlooking something with respect to Eq. 7. Why let the limits of integration get large? It seems like the Fourier series coefficients (CAFs) should be found based on integrating over one period of the autocorrelation.

    Reply

    1. Why let the limits of integration get large?

      You’re right that for a periodic signal, you can always get at the Fourier coefficient for Fourier frequency \alpha by integrating over a single period. However, I say

      the autocorrelation is either a periodic function or an almost periodic function.

      and an almost-periodic signal here is the sum of two or more periodic signals with periods that have ratios that are irrational. These are called incommensurate numbers, and so the two or more periodic components have incommensurate periods. For such signals we have to integrate over all time to get at the Fourier coefficient.

      So we can deal with both cases if we just integrate over all time as in (7) and (8).

      Reply

      1. Ahh. That makes sense. For now I’m just deriving the analytical expression for the CAF and SCF of a baseband BPSK signal. Trying to overcome years of mental atrophy.

        FWIW, I much prefer the probabilistic paradigm over fraction of time. With all due respect to Dr. Gardner as a pioneer in the field (and his very helpful website too), I grew up on probability theory. Thank you for you straightforward presentations.

        Reply

  22. What is the proof for equation (8)? How were you able to drop the expectation E[…] from the definition of R_x(t, tau) in equation (4)?

    Reply

    1. Good morning, Aaron. You were up early today. 🙂

      The expectation operator (i.e. averaging over the ensemble) is absent in equation 8 because of the postulated condition that the equation holds for a cycloergodic signal. If ergodicity is unfamiliar, you will want to read up on it a bit to understand a response to your question. If ergodicity is familiar, maybe you could take a stab at a proof and let the readership know if you find an error or have a more specific question. Writing your own proofs, like writing your own code, is a great way to internalize a new topic.

      Just my $0.02.

      Reply

      1. You might try: R.A. Boyles, W.A. Gardner, Cycloergodic properties of discrete parameter non-stationary stochastic processes,
        IEEE Trans. Inform. Theory IT-29 (1983) 105–114

        I haven’t read that one yet, but it was cited in “Cyclostationarity, Half a century of research.” That paper rivals Chad’s Literature page for crazy number of citations.

        Reply

      2. We gotta make an all-out stand
        Hey yo, I’m gonna need a right-hand man
        …
        I cannot be everywhere at once, people!
        I’m in dire need of assistance

        -Right Hand Man
        Lin-Manuel Miranda
        (Hamilton)

        Thanks John! I agree.

        It follows from the definition of a cycloergodic process, as John says. There, E[g(X(t))] = E^{\{\alpha\}}[g(x(t))] almost surely, where the latter is the fraction-of-time expectation, also known as the multiple-sine-wave extraction operator.

        But … we could also sidestep the invocation of cycloergodicity by staying completely within the non-stochastic time-series framework and just stating that (8) is the definition of the cyclic autocorrelation for some signal (time-series) x(t).

        Reply

        1. Doh! I must have been up too early myself. I thought through my bleary eyes that Aaron was responding to me rather than starting a new question.

          My apologies to you both for getting in the middle.

          Reply

  23. I’m trying follow the article, and at equation 8, it looks to me like if the complex exponential were removed (for example if alpha = 0), that the resulting single integral on the RHS of the equation is just the autocorrelation function(ACF) of a WSS process. And that the integration is part of the convolution/correlation.

    And if that is right, and it might not be, cause I’m super rusty on my signal processing fundamentals, I’m having trouble with what I see above it in equations 6 and 7, which seem to indicate, that the integration is part of what I’d call the fourier integral, the integrand of which is the ACF.

    It feels like there must be an integration falling out somewhere that I’m missing. And I’m thinking about eq.7 being that a fourier coefficient is just a correlation of a time series with the basis complex exponential at whatever given frequency, is there some linearity or some property that allows us to ignore one of these integral terms.

    Reply

    1. Hey Joe! Thanks for checking out the CSP Blog and leaving a thoughtful comment.

      equation 8, it looks to me like if the complex exponential were removed (for example if alpha = 0), that the resulting single integral on the RHS of the equation is just the autocorrelation function(ACF) of a WSS process.

      Not quite. The autocorrelation of a process is defined in terms of an expected value. Equation (8) contains no expected value. I think what you are struggling with is the concepts of ergodicity and cycloergodicity. I’ve updated the text near (8) to provide a link to a rather long and difficult post I recently did on random processes, stationarity, cyclostationarity, and ergodicity.

      The idea is that we can start with our usual stochastic process (random process) formalism, which leads to the periodically time-variant autocorrelation function for a cyclostationary process, then invoke cycloergodicity, which allows us to operate with time averages directly on the sample paths of the process, and we’ll get results consistent with the process parameters.

      So, yeah, if you apply the expectation to the delay product X(t+\tau/2)X^*(t-\tau/2), you’ll get a periodic function, which we can then express as a Fourier series, and there we see a Fourier-type integral. But we can invoke cycloergodicity to obtain that same Fourier-series coefficient (the cyclic autocorrelation) by operating on a sample path (actually, by operating on almost any path).

      Where it gets really weird is when you decide to abandon the conventional stochastic process theory and embrace a kind of probability theory that uses only a single sample path–an infinite-duration model of the data you actually want to process. In that case, you can use the theory of fraction-of-time probability to define all probabilistic quantities and you’ll never need ergodicity or cycloergodicity because there is no stochastic process to deal with. See The Literature [R1] and Napolitano.

      So we have a couple ways of justifying (8).

      What do you think?

      Reply

  24. Hi Chad,
    I am struggling a bit, I think as a few others have commented as well, on utilising this on a discrete set of data i.e. your BPSK signal from the first post.
    I used one of the other codes on here and I am able to get the same plots for the ACF as you are for the original BPSK signal from the first post you put up. The problem occurs when I try to implement Equation (14) from https://cyclostationary.blog/2015/09/28/the-spectral-correlation-function/. Where your SCF output is symmetric, almost an X ( https://cyclostationary.blog/2015/09/28/the-spectral-correlation-function-for-rectangular-pulse-bpsk/). My version is shifted , so that its a skewed X.
    I calculate (14) in this post with the following piece of code, using the BPSK from the first page with noise removed:
    tau_lim = 30;
    taus = -tau_lim:tau_lim;
    alphas = -0.5:0.1:0.5;
    Rx = zeros(length(alphas), length(taus));
    fx=-0.5:0.01:0.5;
    for tau = taus
    unshift = y_of_t(1 : end – abs(tau));
    shifted = y_of_t(1 + abs(tau) : end);
    len = length(shifted);
    if tau > 0
    comp_exp = exp(-1i * 2*pi * alphas’* (tau : len + tau – 1));
    else
    comp_exp = exp(-1i * 2*pi * alphas’* (0 : len – 1));
    end
    result = shifted .*unshift .* comp_exp;
    sum(result)/len;
    Rx(:, 1+(tau+tau_lim)) = sum(result, 2) / len .* exp(-i*pi*alphas’*tau) ;
    end
    I then convert it to the SCF using (14) (https://cyclostationary.blog/2015/09/28/the-spectral-correlation-function/) :
    %Create the FFT function exp(-i*2*pi*fx*taus)
    Fx=exp(-1i * 2*pi * fx’* (taus));
    %Create a container for the SCF
    Sx=zeros(length(fx),length(alphas));
    %Convert from ACF to SCF
    for i=1:length(alphas)
    % Rx = sqrt(Rx .* conj(Rx));
    Sx(:,i)=sum(Rx(i,:).*Fx,2);
    end
    Sx=Sx’;
    And then print:
    Sx = sqrt(Sx .* conj(Sx));
    figure(1)
    hold on
    for i = 1:length(alphas)
    hp = plot3(alphas(i)*ones(size(fx)), fx, Sx(i, :));
    set(hp, ‘linewidth’, 2);
    end
    view(3);
    grid on;
    xlabel(‘alpha’);
    ylabel(‘f’);
    zlabel(‘Lin. Magnitude’);
    If you run this, you will see its a skewed X.
    I think my lack of understanding comes from the fact that I don’t understand what an Asymmetrical Autocorrelation Function is, and why you need to scale it by a complex exponential ? If I remove the exp(-i*pi*alphas’*tau) portion of my code, I see no difference in the ACF, however the X skews more in the SCF. So I think that is my issue.
    Thanks

    Reply

    1. Hey Sandford! Welcome to the CSP Blog and thanks for the comment.

      We’ll get you going soon.

      I can’t run your code because I get an error early on:

      But I can see problems with that first section, where you are estimating the cyclic autocorrelation function:

      Recall that the normal (symmetric) cyclic autocorrelation functions shifts the data x(t) up and down (or left and right if visualizing time along the x-axis) by \tau/2. But in discrete time, for odd values of the delay \tau, a delay of \tau/2 is not available. We can only shift the data by integer numbers of samples (unless we resample). So the asymmetric version of the cyclic autocorrelation function keeps one factor unshifted (x(t)) and shifts the other by \tau:

      When we analyze this asymmetric function, we can relate it to the desired symmetric one by (14),

      In your code, you are shifting one of the factors by |\tau|. So that is already a problem if you then consider the resulting function (estimate) a function of \tau: it must be a different function than desired. You need to deal with negative as well as positive \tau.

      I suggest that you use a circular shift, which is sensitive to the sign of \tau. If you want, you can take care to zero out the appropriate end of the to-be-shifted vector so that the shift is no longer circular, but you can still use the convenient circshift.m. If the range of \tau is small compared to the length of the data to be processed (and it is here), then such edge effects are tiny anyway.

      Also, regarding the asymmetric compensation line of your code (last line of the tau loop), notice that (14) says that to obtain the asymmetric version of the cyclic autocorrelation from the symmetric one, multiply by e^{-i\pi \alpha \tau}. But to obtain the symmetric one from the asymmetric one, you’ll need to multiply by e^{i\pi \alpha \tau}.

      Go ahead and reread that section on Symmetric Versus Asymmetric Cyclic Autocorrelation Functions and let me know where it might be confusing–I can clarify.

      Reply

Leave a Comment, Ask a Question, or Point out an Error