The Principal Domain for the Spectral Correlation Function

Let’s talk about that diamond-shaped region in the $(f, \alpha)$ plane we so often see associated with CSP. I’m talking about the principal domain for the discrete-time/discrete-frequency spectral correlation function. Where does it come from? Why do we care? When does it come up?

The Principal Domain in Conventional DSP

When we talk about the Fourier transform of a continuous-time signal $x(t)$ ,

$\displaystyle X(f) = \int_{-\infty}^\infty x(t) e^{-i2 \pi f t} \, dt \hfill (1)$

there is no necessary structure to the transform $X(f)$ . That is, it is not necessarily periodic, constant, bounded (think $x(t) = \cos(t)$ ), of bounded support, etc. It can possess some symmetries for some signals, such as when $x(t)$ is real-valued.

On the other hand, when we consider the discrete Fourier transform (DFT), typically computed by the fast Fourier transform (FFT) algorithm, the transform does have strong structure: it is periodic. The DFT for a sequence $x(k), k = 0, 1, 2, \ldots, N-1$ is given by

$\displaystyle X(j) = \sum_{k=0}^{N-1} x(k) e^{-i 2 \pi (j/N) k}, j = 0, 1, 2, \ldots, N-1 \hfill (2)$

Here I’ve implicitly set the sampling increment equal to unity, as is my custom at the CSP Blog. The Fourier frequencies are $j/N$ and the time instants are $k$ . To relate to a physical analog signal that was sampled at rate $f_s = 1/T_s$ Hz, replace $k$ with $kT_s$ and $j$ with $jf_s$ . In normalized units, the Fourier frequencies range from $0$ to $(N-1)/N$ , and in physical units they range from $0$ to $f_s(N-1)/N$ . But they don’t have to be assigned to that range. Let’s see why.

Removing the restriction on $j$ , we simply have the function

$\displaystyle X(j) = \sum_{0}^{N-1} x(k) e^{-i 2 \pi (j/N) k} \hfill (3)$

defined for all integers $j$ . Let’s look at $X(j+M)$ ,

$\displaystyle X(j+M) = \sum_0^{N-1} x(k) e^{-i 2 \pi (j/N)k} e^{-i 2 \pi (M/N)k} \hfill (4)$

We see immediately that if $M = pN$ , with $p$ any integer, then $\displaystyle e^{-i 2 \pi (pN/N)k} = e^{-i 2 \pi pk} = 1$ . Therefore $X(j)$ is a periodic function of index $j$ with period $N$ .

The periodicity of the DFT means we can view the frequency-domain content of $x(k)$ by looking at any set of $N$ consecutive values of the transform $\displaystyle \{X(m), X(m+1), \ldots, X(m+N-1)\}$ . Mathematicians (and MATLAB through fft.m) choose $m=0$ . Electrical engineers and signal processors typically choose $m= -N/2$ (in MATLAB, follow fft.m with fftshift.m).

We can also view the periodicity of the DFT in terms of the Fourier frequencies themselves–the $j/N$ . Changing the notation slightly, we obtain the alternative expression

$\displaystyle X(f_j) = \sum_0^{N-1} x(k) e^{-i 2 \pi f_j k} \hfill (5)$

where $f_j = j/N$ . This means the DFT is a periodic function of frequency with period equal to $1$ . In signal processing, we like to look at the period that corresponds to the frequency interval $[-0.5, 0.5-1/N]$ , which puts zero frequency in the middle of the interval. In terms of physical frequencies, the interval is $[-f_s/2, f_s/2)$ , where the open right end of the interval means $f_s/2-f_s/N$ . This interval of frequencies is the principal domain in signal processing. You can see it follows from the periodicity of the discrete Fourier transform. We’ll now turn to the principal domain in CSP, which also depends on the periodic property of the DFT.

The Principal Domain in CSP

What is the region in the $(f, \alpha)$ plane for which the spectral correlation function is, or can be, unique? That is, what is the minimal-area region for which each point within the region can be different (not dependent on) all other points within the region, but which is identical to one or more points outside of the region? Where is the spectral correlation function non-redundant?

Let’s approach these questions by looking at the non-conjugate cyclic periodogram,

$\displaystyle I_{x_N}^\alpha(t, f) = \frac{1}{N} X_N(t, f+\alpha/2) X_N^*(t, f-\alpha/2) \hfill (6)$

where $\displaystyle X_N(t, f)$ is the discrete Fourier transform for $N$ samples of $x(k)$ beginning with sample $t$ .

Each of the discrete Fourier transforms $X_N(t, f)$ in (6) are periodic functions of frequency $f$ , with period equal to $1$ . For cycle-frequency $\alpha = 0$ , we have the conventional periodogram,

$\displaystyle I_{x_N}^0(t, f) = \frac{1}{N} \left| X_N(t, f) \right|^2 \hfill (7)$

which must be periodic with period $1$ . That simply means that the periodogram is unique (need only be considered) on an interval of length $1$ , and we signal processors typically choose the interval $[-0.5, 0.5-1/N]$ as discussed earlier.

For spectral frequency $f = 0$ , we have the cyclic periodogram slice given by

$\displaystyle I_{x_N}^\alpha(t, 0) = \frac{1}{N} X_N(t, \alpha/2) X_N^*(t, -\alpha/2) \hfill (8)$

which is periodic with period $2$ ,

$\displaystyle I_{x_N}^{\alpha+2}(t, 0) = \frac{1}{N} X_N(t, \alpha/2 + 1) X_N^*(t, -\alpha/2 - 1) = \frac{1}{N} X_N(t, \alpha/2) X_N^*(t, -\alpha/2). \hfill (9)$

So we know the non-redundant set of frequencies can be specified by the interval $[-0.5, 0.5-1/N]$ , and the non-redundant set of cycle frequencies can be specified by the interval $[-1.0, 1.0-1/N]$ . That is, the maximum and minimum frequencies to consider are $0.5$ and $-0.5$ , and the maximum and minimum cycle frequencies to consider are $1.0$ and $-1.0$ , when $N$ is large. And this does suggest the diamond shape we’ve seen in previous CSP Blog posts and in many papers–the one with vertices $(-0.5, 0), (0, 1.0), (0.5, 0),$ and $(0, -1.0)$ .

The diamond shape of the spectral correlation principal domain can be inferred by considering the act of shifting the $N$ -point FFT to the left and to the right by the amount $\lfloor N\alpha/2 \rfloor$ , where $\alpha$ is the cycle frequency in normalized Hz. The amount of overlap is $N$ for $\alpha =0$ (the periodogram), and linearly decreases to zero for $\alpha = 1$ . For any valid $\alpha$ , the overlap is about $1 - \alpha$ in Hz.

The periodicity means that we can write the following relationship between the principal domain centered at the origin and all other domains:

$\displaystyle S_x^\alpha(f) = S_x^{\alpha+j-k}(f+(j+k)/2) \hfill (10)$

where $j$ and $k$ are arbitrary integers.

So each point in the principal domain diamond translates to a corresponding point in each of an infinite number of non-principal-domain diamonds. The redundancy of a point within the principal domain with another point within the principal domain is possible too, and depends on the nature of the signal. See the CSP Blog post on symmetry relations for second-order cyclostationarity for more details. For complex-valued data, we can safely restrict our attention to $\alpha \ge 0$ and $f \in [-0.5, 0.5)$ for the non-conjugate spectral correlation function, but must consider the entire principal-domain diamond for conjugate spectral correlation.

Figure 1. The principal domain for the spectral correlation function and all the other non-principal domains. Each blue diamond is redundant with the principal-domain pink diamond.

Time- and Frequency-Domain Reconciliation

But there is a catch or stumbling block. The kinds of results we get depends on whether we proceed with estimation in the time domain or in the frequency domain. When we estimate the cyclic autocorrelation function directly in the time domain we can obtain results that appear to contradict the principal-domain results.

Staying with our usual rectangular-pulse BPSK signal, with bit rate $1/10$ and carrier offset $0.05$ , we can straightforwardly compute the cyclic autocorrelation function by a finite-time average,

$\displaystyle \hat{R}^\alpha(\tau) = e^{i \pi \alpha \tau} \frac{1}{N} \sum_{n=0}^{N-1} x(n) x^*(n-\tau) e^{-i 2 \pi \alpha n} \hfill (11)$

where we have adopted the asymmetric version of the cyclic autocorrelation because we are in discrete time. Now substitute for the cycle frequency $\alpha$ a shifted-by-one cycle frequency $\alpha_* = \alpha + 1$ ,

$\displaystyle \hat{R}^{\alpha_*} (\tau) = e^{i \pi \alpha \tau + i \pi\tau} \frac{1}{N} \sum_{n=0}^{N-1} x(n)x^*(n-\tau) e^{-i 2 \pi \alpha n} \hfill (12)$

$\displaystyle \hat{R}^{\alpha_*} (\tau) = e^{i \pi\tau} \hat{R}^\alpha(\tau) \hfill (13)$

This means, for example, that the magnitude of the non-conjugate cyclic autocorrelation function for $\alpha = -0.1$ is the same as for $\alpha = 0.9$ . We already know that the magnitude of the non-conjugate cyclic autocorrelation for $\alpha = 0.1$ is the same as that for $\alpha=-0.1$ , so that means that the cyclic autocorrelation magnitude for $\alpha = 0.1$ is the same as that for $\alpha = 0.9$ !

The result of applying (11) to our BPSK signal is shown in Figure 2. By comparing to the cyclic autocorrelation function gallery (or Figure 4 below), we can see that all is well for cycle frequencies $\{0, 0.1, 0.2, 0.3, 0.4, 0.5\}$ , but after $0.5$ , the feature magnitudes repeat in reverse order: the feature for $0.6$ is the same as that for $0.4$ , the one for $0.9$ is the same as that for $0.1$ , etc. So this seems to imply that the range of cycle frequency is not $[-1.0, 1.0)$ , but is more like $[-0.5, 0.5)$ .

What happened? This is a question I’ve gotten many many times.

My answer is that when we start with (11), we are explicitly constructing and using a sine wave of the form $\displaystyle e^{-i 2 \pi \alpha n}$ , and the unique sine waves of this form have frequencies on the interval $[-0.5, 0.5)$ . That is, we cannot unambiguously represent all sine waves in this discrete-time formulation of the cyclic-autocorrelation estimation problem. We can only get results from (11) for sine waves with frequencies that are in the signal-processing principal domain. If we try to use $\alpha = 0.6$ , for example, we see that that sine wave is the same as the sine wave with frequency $\alpha = -0.4$ . It’s a result of the sampling theorem (which we’ll get to in the Signal Processing ToolKit series soon). Essentially, if we stick with this time-domain approach (equation (11)), we are limited in the values of the cycle frequency we can operate with relative to those that make up the full spectral-correlation principal domain.

One way out of this problem is to increase the sampling rate of the data. If we doubled the sampling rate of our BPSK signal, we’d have a bit rate of $1/20$ , and we could still ‘see’ the cyclic autocorrelation estimates for cycle frequencies up to $10/20 = 0.5$ , but this means ten harmonics of the symbol rate, not just five. The signal has an infinite number of cycle frequencies, but after eight or nine, the magnitudes of the cyclic autocorrelation are too small to be useful anyway, so upsampling helps as much as we need it to.

Figure 2. Cyclic autocorrelation estimates for a rectangular-pulse BPSK signal with bit rate 0.1 and carrier offset of 0.05 using a direct time-domain estimator like (11). Compare with theory in the PSK verification post, the cyclic autocorrelation gallery post, or Figure 4 below to see that these surfaces are not completely correct. They are only correct for a limited range of the cycle frequencies that are unique to the spectral correlation function’s principal domain.

But another way out is to abandon the direct time-domain estimation method (11) and estimate the cyclic autocorrelation by first estimating the spectral correlation function and then inverse Fourier transforming. Because the cyclic periodogram requires only the multiplication of two frequency-shifted Fourier transforms, each with the usual periodicity of one, we sidestep the need to explicitly generate a sine wave with frequency related to a cycle frequency, unlike (11). If we estimate the spectral correlation function we get the results in Figure 3, and if we then inverse transform each slice (fixed cycle frequency, variable spectral frequency), we obtain the desired cyclic autocorrelation estimates shown in Figure 4.

Figure 3. The spectral correlation function estimated for true cycle frequencies for the rectangular-pulse BPSK signal described in this post. Here the spectral correlation function is estimated using the frequency-smoothing method. Note the visible diamond-shaped region of support–the frequency-smoothing method produces spectral correlation estimates on the principal-domain diamond.

Figure 4. Cyclic autocorrelation estimates obtained by inverse Fourier transforming each slice of the surface shown in Figure 3 (using the complex-valued version, not just the magnitudes). These estimates conform to theory.

So the principal domain for the spectral correlation function is the pink diamond in Figure 1, although we cannot achieve cyclic autocorrelation estimates for each correct slice in the spectral correlation function if we insist on direct time-domain estimation of the cyclic autocorrelation. If we focus on the frequency domain, we know that we will not miss any signal’s cycle frequencies if we cover the principal domain with spectral correlation function point estimates with appropriate resolution.

Are you buying what I’m selling? Let me, and all the CSP Blog readers, know in the comments.

Author: Chad Spooner

I'm a signal processing researcher specializing in cyclostationary signal processing (CSP) for communication signals. I hope to use this blog to help others with their cyclo-projects and to learn more about how CSP is being used and extended worldwide. View all posts by Chad Spooner

2 thoughts on “The Principal Domain for the Spectral Correlation Function”

Kaushal kumar says:

October 19, 2022 at 7:05 am

Dr. Spooner

In the lines “So we know the non-redundant set of frequencies can be specified by the interval [-0.5, 0.5-1/N], and the non-redundant set of frequencies can be specified by the interval [-1.0, 1.0-1/N].”,just after equation 9 it should be “So we know the non-redundant set of frequencies can be specified by the interval [-0.5, 0.5-1/N], and the non-redundant set of cyclic frequencies can be specified by the interval [-1.0, 1.0-1/N].” Am I correct

Loading...

1. Chad Spooner says:
  
  October 19, 2022 at 7:11 am
  
  Welcome to the CSP Blog Kaushal!
  
  Yes, you are correct and I have fixed that sentence to include “cycle” prior to the second occurrence of “frequencies.”
  
  Thanks much!
  
  Loading...