The Spectral Correlation Function

Spectral correlation is perhaps the most widely used characterization of the cyclostationarity property. The main reason is that the computational efficiency of the FFT can be harnessed to characterize the cyclostationarity of a given signal or data set in an efficient manner. And not just efficient, but with a reasonable total computational cost, so that one doesn’t have to wait too long for the result.

Just as the normal power spectrum is actually the power spectral density, or more accurately, the spectral density of time-averaged power (variance), the spectral correlation function is the spectral density of time-averaged correlation (covariance). What does this mean? Consider the following schematic showing two narrowband spectral components of an arbitrary signal:

scf_schematic

The sequence of shaded rectangles on the left are meant to imply a time-series corresponding to the output of a bandpass filter centered at f-A/2 with bandwidth B. Similarly, the sequence of shaded rectangles on the right imply a time-series corresponding to the output of a bandpass filter centered at f+A/2 with bandwidth B.

Let’s call the first time-series Y_1(t, f-A/2) and the second one Y_2(t, f+A/2). Since these time-series, or signals, are bandpass in general, if we attempt to measure their correlation we will get a small value if A \neq 0. However, if we downconvert each of them to baseband (zero frequency), we will obtain lowpass signals, and there is the possibility that these new signals are correlated to some degree.

As a limiting case, suppose each of the Y_j(t, \cdot) signals were noiseless sine waves. By the construction of the figure, their frequencies must be different if A \neq 0, and so their correlation will be zero. However, if each sine wave is perfectly downconverted to baseband, the resulting signals are simply complex-valued constants, and the correlation between the two constant time-series is high.

In general, the narrowband time-series Y_1(t, \cdot) and Y_2(t, \cdot) are not simple sine waves, but complicated random processes. But the correlation between separated spectral components–spectral correlation–is still a highly useful characterization of the signal for a large class of interesting signals.

The spectral components (the individual downconverted narrowband spectral components of the signal) are most often obtained through the use of the Fourier transform. As the transform of length T slides along the signal, if produces a number of downconverted spectral components with approximate bandwidth 1/T.  The two involved time-series of interest are then renamed as X_T(t, f-A/2) and X_T(t, f+A/2). A measure of the spectral correlation is given by the limiting average of the cyclic periodogram, which is defined by

\displaystyle I_T^A (t, f) = \frac{1}{T} X_T(t, f-A/2) X_T^*(t, f+A/2), \hfill (1)

as the amount of processed data increases without bound, and then the spectral resolution (B = 1/T) is allowed to decrease to zero,

\displaystyle S_x^A (f) = \lim_{T\rightarrow\infty} \lim_{U\rightarrow\infty} \displaystyle\frac{1}{U} \int_{-U/2}^{U/2} I_T^A(t, f) \, dt. \hfill (2)

The limit spectral correlation function we just wrote down is a time-smoothed (time-averaged) cyclic periodogram. But the limit function can also be obtained by frequency smoothing the cyclic periodogram

\displaystyle S_x^A(f) = \lim_{\Delta\rightarrow 0} \lim_{T\rightarrow\infty} g_\Delta(f) \otimes I_T^A(t, f), \hfill (3)

where g_\Delta(f) is a unit-area pulse-like smoothing kernel (such as a rectangle). In (3), the symbol \otimes denotes convolution.

The Significance of the Frequency A

The spectral correlation function (SCF) S_x^A(f) is typically zero for almost all real numbers A. Those A for which the SCF is not identically zero are called cycle frequencies (CFs). The set of SCF CFs is exactly the same as the set of cycle frequencies for the cyclic autocorrelation function (CAF)! That is, the separation between correlated narrowband signal components of x(t) is the same as a frequency of a sine wave that can be generated by a quadratic nonlinearity (for example, a squarer or a delay-and-multiply device) operating on the original time-series data x(t).

The Cyclic Wiener Relationship

It can be shown that the Fourier transform of the CAF is equal to the SCF (The Literature [R1], My Papers [5,6]):

\displaystyle S_x^\alpha(f) = \int_{-\infty}^\infty R_x^\alpha(\tau) e^{-i 2 \pi f \tau}\, d\tau, \hfill (4)

which is called the cyclic Wiener relationship. The Wiener relationship (sometimes called the Wiener-Khintchine theorem)  is a name given to the familiar Fourier transform relation between the conventional power spectral density and the autocorrelation

\displaystyle S_x^0 (f) = \int_{-\infty}^\infty R_x^0(\tau) e^{-i 2 \pi f \tau} \, d\tau, \hfill (5)

where S_x^0(f) is the conventional power spectrum and R_x^0(\tau) is the conventional autocorrelation function.

Conjugate Spectral Correlation

This post has defined the non-conjugate spectral correlation function, which is the correlation between X_T(t, f-\alpha/2) and X_T(t, f+\alpha/2). (The correlation between random variables X and Y is defined as E[XY^*]; that is, the standard correlation includes a conjugation.)

The conjugate SCF is defined as the Fourier transform of the conjugate cyclic autocorrelation function,

\displaystyle S_{x^*}^\alpha (f) = \int_{-\infty}^\infty R_{x^*}^\alpha (\tau) e^{-i 2 \pi f \tau}\, d\tau . \hfill (6)

From this definition, it can be shown that the conjugate SCF is the density of time-averaged correlation between X_T(t, f+\alpha/2) and X_T^*(t, \alpha/2-f),

\displaystyle S_{x^*}^\alpha (f) = \lim_{T\rightarrow\infty} \lim_{U\rightarrow\infty} \displaystyle\frac{1}{U} \int_{-U/2}^{U/2} J_T^\alpha(t, f) \, dt, \hfill (7)

where J_T^\alpha(t, f) is the conjugate cyclic periodogram

J_T^\alpha(t,f) = \displaystyle \frac{1}{T} X_T(t, f+\alpha/2) X_T(t, \alpha/2-f). \hfill (8)

The detailed explanation for why we need two kinds of spectral correlation functions (and, correspondingly, two kinds of cyclic autocorrelation functions) can be found in this post.

Illustrations

The SCF below is estimated (more on that estimation in another post) from a simulated BPSK signal having bit rate of 333.3 kHz and carrier frequency of 100 kHz. A small amount of noise is added to the signal prior to SCF estimation. The (non-conjugate) SCF shows the power spectrum for \alpha = 0 and the bit-rate SCF for \alpha = 333.3 kHz. The conjugate SCF plot shows the prominent feature for the doubled-carrier cycle frequency \alpha = 200.0 kHz, and features offset from the doubled-carrier feature by \pm 333.3 kHz. More on the spectral correlation of the BPSK signal can be found here and here. The spectral correlation surfaces for a variety of communication signals can be found in this gallery post.

A closely related function called the spectral coherence function is useful for blindly detecting cycle frequencies exhibited by arbitrary data sets.

ww_1

Now consider a similar signal: QPSK with rectangular pulses. Let’s switch to normalized frequencies here for convenience. The signal has a symbol rate of 1/10, a carrier frequency of 0.05, unit power, and a small amount of additive white Gaussian noise. A power spectrum estimate is shown in the following figure:

spec_corr_example_psds

Consider also four distinct narrowband (NB) components of this QPSK signal as shown in the figure. The center frequencies are 0.0, 0.1, 0.08, and 0.18. We know that this signal has non-conjugate cycle frequencies that are equal to harmonics of the symbol rate, or k/10 for k = 0, \pm 1, \pm 2, \ldots. This means that the NB components with separations k/10 are correlated. So if we extract such NB components “by hand” and calculate their correlation coefficients as a function of relative delay, we should see large results for the pairs (0.0, 0.1) and (0.08, 0.18), and small results for all other pairs drawn from the four frequencies.

So let’s do that. We apply a simple Fourier-based ideal filter (ideal meaning rectangular pass band) with center frequency f_0 \in \{0.0, 0.1, 0.08, 0.18\}, frequency shift to complex baseband (zero center frequency) and decimate. The results are our narrowband signal components. Are they correlated when they should be and uncorrelated when they should be?

Here are the correlation-coefficient results:

spec_corr_example_cc

Here the signals y_j(t) arise from the frequencies \{0.0, 0.1, 0.08, 0.18\}. So the spectral correlation concept is verified here. One can also simply plot the decimated shifted narrowband components and assess correlation visually:

spec_corr_example_time

Estimators

I’ve written several posts on estimators for the spectral correlation function; they are listed below. I think of them as falling into two categories: exhaustive and focused. For exhaustive spectral correlation estimators, the goal is to estimate the function over its entire (non-redundant) domain of definition as efficiently as possible. For focused estimators, the goal is to estimate the spectral correlation function for one or a small number of cycle frequencies with high accuracy and selectable frequency resolution.

Exhaustive Estimators

Strip Spectral Correlation Analyzer (SSCA)

FFT Accumulation Method (FAM)

Focused Estimators

The Frequency-Smoothing Method (FSM)

The Time-Smoothing Method (TSM)

 

73 thoughts on “The Spectral Correlation Function

    • Do you mean you want to plot the spectral correlation function? Or the cyclic autocorrelation?
      I’m a bit confused because you mention the CAF but the comment is in response to the SCF post.

      For either function, the cycle frequencies are, in the case of the post, the known cycle
      frequencies for the BPSK signal with bit rate of 333.3 kHz and carrier offset frequency of 100 kHz.
      This particular signal has a bandlimited pulse shape (not rectangular in the time domain,
      as we usually use on this site), so the non-conjugate cycle frequencies are +333.3, 0, -333.3 kHz.
      The conjugate cycle frequencies are 2×100 = 200 kHz, 200 – 333.3 kHz, and 200 + 333.3 kHz.

      In your estimator, you need to use normalized frequencies at some point.
      In my case, the sampling rate was 1 MHz, so that the normalized bit rate is 1/3, and
      the normalized carrier offset is 0.1. That leads to normalized non-conjugate cycle
      frequencies of -1/3, 0, 1/3 and normalized conjugate cycle frequencies of 0.2, 0.2-1/3,
      and 0.2+1/3.

      For the CAF, you can first create estimates of the SCF for each cycle frequency,
      then inverse Fourier transform them. The SCF estimates will correspond to all spectral
      frequencies between negative half the sampling rate and positive half the sampling rate,
      and transforming will then provide you with a CAF estimate over many distinct values of
      the lag variable. Using that method, you don’t have to choose any particular lags,
      you get them all.

  1. I would like to know the formula of the threshold or the average to make the decision about spectrum sensing ? (likely as energy detector algorithm for spectrum sensing )

    • I don’t have a formula for the threshold for detection based on spectral correlation. I typically choose a threshold using empirical methods, which can take into account all the real-world deviations from a typical simple AWGN model. However, if you use the spectral coherence function, there is a way to compute a threshold that is based on some older work on cross spectral density measurements. I have a post in the works on how to compute and apply that thresholding formula. For now, look at The Literature [G. Carter, R64].

  2. Ian Frasch says:

    Hey Chad, I just want to say that this post really helped me understand what the spectral correlation function is conceptually (the correlation between different FFT frequency bins as they change over time). That first diagram was especially helpful. Thanks.

    • Eq (8) is correct as is. Notice that this is the conjugate cyclic periodogram, not the normal one (which is Eq (1)). You might want to read the post on conjugation configurations, which also explains why we need to consider a variety of numbers of conjugated terms when we study higher-order statistics for complex-valued signals. Yes, X_T(t, f) is the complex envelope.

      • Gabriel. says:

        Hi Dr Spooner,

        Second bracket, right of eq (8) has (α/2-f) which is different from (f- α/2) as in eq (1). Please is this reversal correct? If it’s correct, is it because is the conjugate cyclic periodogram?

        Kind regards,
        Gabriel.

        • Yes, the reversal is correct. And you are right, it is because that quantity corresponds to the conjugate cyclic periodogram. Probably the easiest way to see this is to look at my post on the cyclic polyspectrum, and take the general nth-order case to the specific case of n=2 and no conjugations.

          I haven’t forgotten about your previous comment.

  3. Samith says:

    Dear Chad, Thank you very much for your explanation. I followed several resources to understand the concept of spectral correlation, but those weren’t so helpful. Your post is the ideal.

  4. Mohamed says:

    Hey Chad,

    Thanks for this valuable blog.
    I have some difficulty to understand the physical meaning of the cyclic spectrum.
    Let’s pick a specific point (f1,alpha1,a1) in the 3D spectral correlation figure, where f1 is the spectral frequency value, alpha1 is the cyclic frequency value, and a1 is the magnitude value. And Let’s say f1=10KHz and alpha1=100KHz. What these values are representing related to the original signal?

    • a_1 represents the complex-valued correlation between two narrowband frequency components. The first is at f_1 + \alpha_1/2 and the second is at f_1 - \alpha_1/2 (assuming we’re talking about the non-conjugate spectral correlation function for now). But it is really the density of correlation, meaning that the ideal spectral correlation function considers the two narrowband frequency components to have infinitesimal width.

      So a_1 represents the “idealized” correlation between the time-varying fluctuations in a very narrow band centered at 60 kHz and the time-varying fluctuations in a second very narrowband centered at -40 kHz.

      Does that help?

      • Mohamed says:

        For a specific alpha, the correlation has to be between 2 narrowband frequency components or it can be between 2 shifted versions of the whole spectrum?

        • For a specific \alpha and f, the correlation is between two narrowband components. But the spectral correlation function is a (typically continuous) function of f and discrete function of \alpha, so if you compute the entire function, you get something that you might call a set of correlations between all possible shifts of the entire spectrum.

          If you fix \alpha, and let f vary over [-f_s/2, f_s/2], then you’ll get a set of all the possible correlations between narrowband spectral components whose separation is \alpha.

          If you fix f, and let \alpha vary over [-f_s, f_s], then you’ll get a set of all the possible correlations between narrowband spectral components whose center frequencies have average value f.

  5. aaron says:

    Hi, Chad,
    Thanks for your post. I have a question regarding to “first time-series Y_1” and “second time-series Y_2”. Are they both the sequence of shaded rectangles on the right? Or Y_1 should be the left ones with center frequency f-A/2? Thanks!

    • Aaron:

      Thanks for finding a typo in the Spectral Correlation post! Yes, Y_1 is supposed be the sequence of rectangles centered at f - A/2, but I had written Y_1(t, f+A/2) and Y_2(t, f+A/2). So I’ve changed the first + to a -.

  6. Deepu says:

    Hi Chad,

    Thanks for this wonderful blog. I am new to this cyclostationary signal analysis. I just wanted to know in case of QPSK signal whether taking a conjugate SCF will help us in estimating an unknown carrier frequency/Offset. Or have you covered parameter estimation using SOC in anywhere in detail?

    • Thanks for writing Deepu.

      No, the conjugate spectral correlation function cannot be used to obtain a high-accuracy estimate of the carrier offset for QPSK (it can for BPSK). That is because QPSK signals, and other symmetric-constellation digital QAM/PSK signals (but not BPSK), do not possess conjugate cyclostationarity. You have to use higher-order cyclostationarity to create a high-accuracy estimator of the carrier offset.

      For parameter estimation, I’ve only covered cycle-frequency estimation and time-difference-of-arrival estimation so far.

      • Deepu says:

        Thanks for the reply Dr Chad Spooner. But from your post regarding the “Conjugation Configuration” i have read that for second order quadrature signals possess conjugate cyclostationarity. I am quoting the content of that post.

        “”This means that in the end, for second-order, we always need to consider the “no conjugations” case z(t+\tau_1)z(t+\tau_2) as well as the “one conjugation” case z(t+\tau_1)z^*(t+\tau_2), which provide us with the conjugate and non-conjugate cyclic autocorrelation and spectral correlation functions, respectively.”

        Am I confused with something else here please help?

        • We always need to consider both the non-conjugate and conjugate cyclic autocorrelation and spectral correlation functions for any signal. However, the conjugate cyclic autocorrelation (and therefore the conjugate spectral correlation function) can be zero for a particular signal type. So in the cases of BPSK, MSK, OOK, and GMSK, the conjugate spectral correlation function is not zero for all possible cycle frequencies, but for QPSK, 16QAM, 8PSK, etc., it is zero. This situation is one of the motivating factors for studying and using higher-order cyclostationarity.

          Does that help?

  7. Ramesh says:

    When I obtain SCF for BPSK and QPSK(or any QAM), I get 4 and 2 peaks in the plot respectively? Can you explain why it is so?

    • Hey Ramesh. Can you elaborate on what you have done? For example, are you using real-valued or complex-valued BPSK, QPSK, QAM? Which “plot” are you referring to?

      I believe I will be able to explain once I get a clear question.

      In the meantime, have you read my posts on QAM/PSK? They might help answer your questions. Here they are:

      QAM and PSK

      SRRC Pulses

      Gallery

  8. Ramesh says:

    And in the estimation of cyclic cumulants in other post. We get peaks at different values of \alpha. Can you explain why?

  9. Hi, Dr.Spooner

    For a complex valued QPSK signal, there is no peaks at its alpha domain (alpha = +/-2fc) because of real component and image component cancel each out. So, if we just take real part of received signal, is possible to use SCF to blindly estimate its carrier frequency ?

    Thanks,

    Best,

    Andy

    • I don’t think so, as the ability to generate a carrier-related cycle frequency is dependent on the signal’s constellation, and that doesn’t change if you take the real part of the complex-valued signal representation. Let’s try a mathematical argument.

      The baseband QAM/PSK signals are given by

      y(t) = \sum_k a_k p(t - kT_0)

      which neglects the symbol-clock phase t_0 and the carrier phase $\phi_0$ that I often include in the model. Not important here.

      The complex-valued transmitted signal representation is

      x(t) = y(t) e^{i 2 \pi f_c t}

      again neglecting carrier phase, and noting that f_c >> 1/T_0. The actual real-valued transmitted-signal model is just the real part of x(t), which is what you’re thinking of:

      x_r(t) = \Re [x(t)]

      or

      x_r(t) = \frac{1}{2} \left[ x(t) + x^*(t)\right]

      Now you are suggesting that if we square x_r(t), we might get at a cycle frequency that is related to the carrier (like the doubled-carrier cycle frequency of BPSK). So let’s take a look at x_r^2(t).

      x_r^2(t) = \frac{1}{4} \left[ x(t) + x^*(t)\right] \left[x(t) + x^*(t)\right]

      or

      x_r^2(t) = \frac{1}{4} \left[ x^2(t) (x^*(t))^2 + 2x(t) x^*(t) \right] = z(t)

      The first-order cycle frequencies for x_r^2(t) can be found by applying the cyclic autocorrelation formula directly to the squared signal, resulting in terms like

      \int x^2(t) e^{-i 2 \pi \alpha t} dt.

      So we see, then, that the cycle frequencies for x_r^2(t) are those for x^2(t), (x^*(t))^2, and
      x(t) x^*(t). This is just the analysis in the conjugation configuration post.

      For BPSK, which has a constellation $a_k \in \{\pm 1\}$, the first-order cycle frequencies are:

      Term 1: \alpha_1 \in \{2f_c, 2f_c\pm 1/T_0\} (SRRC pulses)

      Term 2: \alpha_2 \in \{-2f_c, -2f_c \pm 1/T_0\}

      Term 3: \alpha_3 \in \{0, \pm 1/T_0\}

      And for QPSK,

      Term 1: \alpha_1 = \emptyset

      Term 2: \alpha_2 = \emptyset

      Term 3: \alpha_3 \in \{0, \pm 1/T_0\}

      This just means that if you take the FFT of x_r^2(t), you’ll see the usual second-order cycle frequencies for the underlying signal type.

      If you apply a second-order analysis to x_r^2(t), well, then you’re back to a fourth-order transformation of the original signal, which you are seeking to avoid.

      Did I err?

  10. Gabriel. says:

    Hi Dr Spooner,

    Thank you for your posts. They are very useful.

    I have a question on this topic.
    Since the conjugate spectral correlation function is zero for QPSK, is it the same with offset QPSK (OQPSK) or staggered QPSK (SQPSK)?

    Are there peaks at (+/- 2fc) and (2fc +/- 1/T) where 1/T is the symbol rate for OQPSK or SQPSK?

    Kind regards,

    Gabriel.

    • Good question. First, staggered QPSK (SQPSK) and offset QPSK (OQPSK) refer to the same modulation (which you probably already know, but other readers might not). It consists of two BPSK signals summed in phase quadrature (one baseband binary PAM multiplies a \cos(\cdot), the other multiplies a \sin(\cdot)). This is also true of QPSK. However, for OQPSK/SQPSK, the second binary PAM is delayed relative to the first by half a symbol width.

      The second thing to know here is that MSK is exactly a form of SQPSK/OQPSK. It is an OQPSK signal with pulses in the PAM signals that are half-sines. That is, they are shaped like half a period of a sine wave. So, MSK is exactly an OQPSK signal. This is useful because if we analyze generic OQPSK signals, we get the result for MSK (and GMSK) too. The tricky part of all this is that an MSK signal is a FSK signal, and it has an underlying bit rate (the rate of the PAM signal inside the \cos(\cdot)), whereas we usually talk about the symbol rate of an OQPSK signal. The bit rate of the MSK signal, when viewed as an FSK signal, is twice the symbol rate of the MSK signal, when viewed as an OQPSK signal. Or, the symbol rate is half the bit rate.

      All of this is fairly well known. The staggering of the second PAM signal causes the second-order cycle-frequency pattern to deviate from that for BPSK or QPSK. Basically, every other feature is cancelled. So in the non-conjugate domain, we don’t see the symbol-rate feature, but do see the second harmonic of the symbol rate (if the pulse supports it). In the conjugate domain, we do not see the doubled carrier, but we do see the two features on either side of the doubled carrier.

      So consider two signals. The first is an MSK signal with bit rate 1/5 and carrier 0.1. The second is an OQPSK signal (SRRC pulses) with symbol rate 1/10 and carrier 0.1. So the symbol rate of the OQPSK signal is half the bit rate of the MSK signal, and therefore we should see similar cycle frequencies for the two signals:

      Cycle frequencies for MSK and OQPSK

  11. Gabriel. says:

    Hi Prof. Spooner,

    Thank you so much for taking your time to give such detailed explanations and illustrations. This blog is highly commendable. It is difficult to get articles with such clarity. I came across a table in an article that presented cyclic features of some modulation types including OQPSK/SQPSK. It wasn’t easy to get clear explanations from many articles.

    Yes , just as you have analyzed, the cyclic features of OQPSK/SQPSK were presented as similar to that of MSK in that table. Does that put OQPSK/SQPSK in a position of being analyzed by the second order cyclostationary structure? In other words, can the conjugate SCF be sufficiently used to analyze it with respect to the carrier frequencies? Something not possible with the normal QPSK. I hope that is right?

    Once again, many thanks for your efforts in this blog.

    King regards,

    Gabriel.

    • Yes, absolutely. If you look at the graph from the previous comment I made, although the doubled-carrier feature is missing for MSK/GMSK/OQPSK, there are enough strong conjugate features that are related to the doubled carrier such that it is easy to get an estimate of the carrier provided you properly identify each feature. So, MSK/GMSK/SQPSK is easy to recognize and characterize using only second-order features.

      If you really want a mind-bender, check out \pi/4 DQPSK.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.