## Second-Order Estimator Verification Guide

In this post I provide some tools for the do-it-yourself CSP practitioner. One of the goals of this blog is to help new CSP researchers and students to write their own estimators and algorithms. This post contains some spectral correlation function and cyclic autocorrelation function estimates and numerically evaluated formulas that can be compared to those produced by anybody’s code.

The signal of interest is, of course, our rectangular-pulse BPSK signal with symbol rate $0.1$ (normalized frequency units) and carrier offset $0.05$. You can download a MATLAB script for creating such a signal here.

The formula for the SCF for a textbook BPSK signal is published in several places (The Literature [R47], My Papers [6]) and depends mainly on the Fourier transform of the pulse function used by the textbook signal.

We’ll compare the numerically evaluated spectral correlation formula with estimates produced by my version of the frequency-smoothing method (FSM). The FSM estimates and the theoretical functions are contained in a MATLAB mat file here. (I had to change the extension of the mat file from .mat to .doc to allow posting it to WordPress.) In all the results shown here and that you can download, the processed data-block length is $65536$ samples and the FSM smoothing width is $0.02$ Hz. A rectangular smoothing window is used. For all cycle frequencies except zero (non-conjugate), a zero-padding factor of two is used in the FSM.

For the cyclic autocorrelation, we provide estimates using two methods: inverse Fourier transformation of the spectral correlation estimate and direct averaging of the second-order lag product in the time domain.

## A Gallery of Spectral Correlation

In this post I provide plots of the spectral correlation for a variety of simulated textbook signals and several collected communication signals. The plots show the variety of cycle-frequency patterns that arise from the disparate approaches to digital communication signaling. The distinguishability of these patterns, combined with the inability to distinguish based on the power spectrum, leads to a powerful set of classification (modulation recognition) features (My Papers [16, 25, 26, 28]).

In all cases, the cycle frequencies are blindly estimated by the strip spectral correlation analyzer (The Literature [R3, R4]) and the estimates used by the FSM to compute the spectral correlation function. MATLAB is then used to plot the magnitude of the spectral correlation and conjugate spectral correlation, as specified by the determined non-conjugate and conjugate cycle frequencies.

There are three categories of signal types in this gallery: textbook signals, collected signals, and feature-rich signals. The latter comprises some collected signals (e.g., LTE) and some simulated radar signals. For the first two signal categories, the three-dimensional surface plots I’ve been using will suffice for illustrating the cycle-frequency patterns and the behavior of the spectral correlation function over frequency. But for the last category, the number of cycle frequencies is so large that the three-dimensional surface is difficult to interpret–it is a visual mess. For these signals, I’ll plot the maximum spectral correlation magnitude over spectral frequency $f$ versus the detected cycle frequency $\alpha$ (as in this post).

A complementary gallery of cyclic autocorrelation functions can be found here.

## SCF Estimate Quality: The Resolution Product

The two non-parametric spectral-correlation estimators we’ve looked at so far–the frequency-smoothing and time-smoothing methods–require the choice of key estimator parameters. These are the total duration of the processed data block, $T$, and the spectral resolution $F$.

For the frequency-smoothing method (FSM), an FFT with length equal to the data-block length $T$ is required, and the spectral resolution is equal to the width $F$ of the smoothing function $g(f)$. For the time-smoothing method (TSM), multiple FFTs with lengths $T_{tsm} = T / K$ are required, and the frequency resolution is $1/T_{tsm}$ (in normalized frequency units).

The choice for the block length $T$ is partially guided by practical concerns, such as computational cost and whether the signal is persistent or transient in nature, and partially by the desire to obtain a reliable (low-variance) spectral correlation estimate. The choice for the frequency (spectral) resolution is typically guided by the desire for a reliable estimate.

## The Spectral Coherence Function

In this post I introduce the spectral coherence function, or just coherence. It deserves its own post because the coherence is a useful detection statistic for blindly determining significant cycle frequencies of arbitrary data records.

Let’s start with reviewing the standard correlation coefficient $\rho$ defined for two random variables $X$ and $Y$ as

$\rho = \displaystyle \frac{E[(X - m_X)(Y - m_Y)]}{\sigma_X \sigma_Y}, \hfill (1)$

where $m_X$ and $m_Y$ are the mean values of $X$ and $Y$, and $\sigma_X$ and $\sigma_Y$ are the standard deviations of $X$ and $Y$. That is,

$m_X = E[X] \hfill (2)$

$m_Y = E[Y] \hfill (3)$

$\sigma_X^2 = E[(X-m_X)^2] \hfill (4)$

$\sigma_Y^2 = E[(Y-m_Y)^2] \hfill (5)$

So the correlation coefficient is the covariance between $X$ and $Y$ divided by the geometric mean of the variances of $X$ and $Y$.

## CSP Estimators: The Time Smoothing Method

In a previous post, we introduced the frequency-smoothing method (FSM) of spectral correlation function (SCF) estimation. The FSM convolves a pulse-like smoothing window $g(f)$ with the cyclic periodogram to form an estimate of the SCF. An advantage of the method is that is allows fine control over the spectral resolution of the SCF estimate through the choice of $g(f)$, but the drawbacks are that it requires a Fourier transform as long as the data-record undergoing processing, and the convolution can be expensive. However, the expense of the convolution can be mitigated by using rectangular $g(f)$.

In this post, we introduce the time-smoothing method (TSM) of SCF estimation. Instead of averaging (smoothing) the cyclic periodogram over spectral frequency, multiple cyclic periodograms are averaged over time. When the non-conjugate cycle frequency of zero is used, this method produces an estimate of the power spectral density, and is essentially the Bartlett spectrum estimation method. The TSM can be found in My Papers [6] (Eq. (54)), and other places in the literature.

## CSP Estimators: The Frequency-Smoothing Method

In this post I describe a basic estimator for the spectral correlation function (SCF): the frequency-smoothing method (FSM). The FSM is a way to estimate the SCF for a single value of cycle frequency. Recall from the basic theory of the cyclic autocorrelation and SCF that the SCF is obtained by infinite-time averaging of the cyclic periodogram or by infinitesimal-resolution frequency averaging of the cyclic periodogram. The FSM is merely a finite-time/finite-resolution approximation to the SCF definition.

One place the FSM can be found is in (My Papers [6]), where I introduce time-smoothed and frequency-smoothed higher-order cyclic periodograms as estimators of the cyclic polyspectrum. When the cyclic polyspectrum order is set to $n = 2$, the cyclic polyspectrum becomes the spectral correlation function, so the FSM discussed in this post is just a special case of the more general estimator in [6, Section VI.B].

## Signal Selectivity

In this post I describe and illustrate the most important property of cyclostationary statistics: signal selectivity. The idea is that the cyclostationary parameters for a single signal can be estimated for that signal even when it is corrupted by strong noise and cochannel interferers. Cochannel means that the interferer occupies a frequency band that partially or completely overlaps the frequency band for the signal of interest.

A mixture of signals, whether cochannel or not, is modeled by the simple sum of the signals, as in

$x(t) = s_1(t) + s_2(t) + \ldots + s_K(t) + w(t), \hfill (1)$

where $w(t)$ is additive noise. We can write this more compactly as

$x(t) = \displaystyle \sum_{k=1}^K s_k(t) + w(t). \hfill (2)$

## 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:

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