# 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].

Let’s start by reviewing the FSM for power-spectrum measurement. That is, let’s begin by looking at an estimator for the power spectral density (PSD). We’ll use discrete time and discrete frequencies in this post, in contrast to My Papers [6], because almost everybody who might read this post will be interested in sampled data. Moreover, discrete frequencies lead to discrete cycle frequencies in the method, which causes a significant complication in practice. We’ll highlight that complication, and show how it can be mitigated.

### The Frequency-Smoothed Periodogram

The power spectrum estimate is obtained by smoothing the periodogram. The discrete Fourier transform (DFT) of $x(t)$ is defined by

$X(f) = \displaystyle \sum_{t=0}^{N-1} x(t) e^{-i2\pi f t}, \hfill (1)$

where $f = n/N$ for $n = 0,1, \ldots, N-1$. The periodogram is the normalized squared magnitude of the DFT,

$I(f) = \displaystyle\frac{1}{N} \left| X(f) \right|^2. \hfill (2)$

The periodogram is not a particularly good spectrum estimator when the data is random; it is quite useful when the data is nonrandom and contains periodic components. The ML estimator for the frequency of a tone in WGN is the maximum of the Fourier transform for all $f$, which is usefully approximated by the periodogram. For random data (like communication signals), the periodogram is erratic and does not converge to the true PSD even as $N\rightarrow\infty$.

The PSD can be accurately estimated by smoothing the periodogram,

$\hat{S}(f) = g(f) \otimes I(f), \hfill (3)$

where $g(f)$ is some pulse-like function, such as a rectangle. The smoothing implied by this convolution allows the estimate to converge to a biased (distorted) estimate of the true PSD. If $N$ is large and the width of $g(f)$ is small relative to the fluctuations over frequency in the true PSD, then the bias is minimized and the estimator converges to very nearly the true PSD. More explicitly, the FSM for PSD estimation is

$\hat{S}(f) = g(f) \otimes I(f) \hfill (4)$

$= \displaystyle \sum_\nu g(f-\nu) I(\nu) \hfill (5)$

$= \displaystyle\frac{1}{N} \displaystyle \sum_\nu g(f-\nu) \left| X(\nu) \right|^2 \hfill (6)$

$= \displaystyle\frac{1}{N} \displaystyle \sum_{n=0}^{N-1} g(f - n/N) \left| X(n/N) \right|^2, \hfill (7)$

for $f = m/N$. So to implement this FSM, we need only use a FFT algorithm and a general convolution algorithm capable of handling complex inputs. A FSM-based PSD estimate for our rectangular-pulse BPSK signal is shown here:

This estimate is formed from a data record with length $32,768$ samples and a rectangular smoothing window $g(f)$ with width $164$ frequency bins, which is about $0.005$ (normalized) Hz.

The frequency-smoothed periodogram method of spectrum estimation is referred to as the Daniell method after P. J. Daniell (The Literature [R42]).

### The Frequency-Smoothed Cyclic Periodogram

For the spectral correlation function, we modify the FSM for the PSD
by replacing the periodogram with the cyclic periodogram, which is defined
by

$I^\alpha(f) = \displaystyle\frac{1}{N} X(f + \alpha/2) X^*(f - \alpha/2), \hfill (8)$

where $\alpha$ is a cycle frequency of interest. Since $X(f)$ is a discrete-frequency function, the arguments $f \pm \alpha/2$ are valid only if $\pm \alpha/2 = k/N$. That is, the discrete-frequency function $X(f)$ can be shifted to the left and right only by multiples of the discrete frequency $1/N$. It may happen that the true cycle frequency does not satisfy this constraint. In practice, the closest discrete shifts are chosen for $\pm \alpha/2$. When the two frequencies do not correspond to multiples of $1/N$, then zero-padding the data prior to FSM estimation is advisable; this will allow a closer match to discrete frequencies. In general, the $+\alpha/2$ and $-\alpha/2$ shifts don’t have to be negatives of each other. The idea is to choose the shifts $a_1$ and $a_2$ such that $\left|\alpha - |a_1/N| - |a_2/N|\right|$ is minimized (My Papers [16]).

To estimate the SCF, then, we convolve the cyclic periodogram with a smoothing window as before,

$\hat{S}^\alpha (f) = g(f) \otimes I^\alpha(f) \hfill (9)$

$= \displaystyle\frac{1}{N} \displaystyle \sum_{n=0}^{N-1} g(f - n/N) X(n/N + a_1/N) X^*(n/N - a_2/N) \hfill (10)$

for $f = m/N$.  All that is required to implement this algorithm is an FFT, complex multiplication, and a general convolution algorithm. The conjugate SCF is estimated by frequency smoothing the conjugate cyclic periodogram given by

$I_*^\alpha(f) = \displaystyle\frac{1}{N} X(f + \alpha/2) X(\alpha/2-f). \hfill (11)$

Applying this algorithm to our rectangular-pulse BPSK signal, we obtain the following estimates for the non-conjugate SCF:

The conjugate SCF estimates are:

To obtain these estimates, we used a data-record length of $32,768$ samples and a rectangular $g(f)$ with width $164$ frequency points, or about $0.005$ Hz. Recall the sampling frequency is set to unity here, the signal’s bit rate is $1/T_0 = 0.1$, and the carrier frequency is $f_c = 0.05$. The $32,768$-sample data record is zero-padded by a factor of two prior to forming the cyclic periodogram.

### Choice for the Frequency-Smoothing Window $g(f)$

We used the rectangular pulse $g(f)$ in our estimator examples for a good reason: computational cost. For arbitrary $g(f)$, a general-purpose convolution routine can be used to find the SCF estimate. Computational cost can be reduced by selecting only a subset of all possible output frequencies, and computing the convolution just for those particular spectral frequencies.

Arbitrary output frequencies can be accommodated, however, by using a rectangular smoothing window $g(f)$ because the convolution can be efficiently performed by a running sum over the cyclic periodogram. The estimates for successive spectral frequencies can be obtained by subtracting the cyclic-periodogram value that is no longer covered by the support of $g(f)$ (left edge of the window) and adding the one that is newly included in the support.

In the final analysis, the FSM offers good control over spectral resolution, but requires an FFT operation on the entire data record. When the signal of interest has low signal-to-interference-and-noise (SINR) ratio, the required data-record length can be very large. For non-real-time CSP applications running on modern general-purpose computers, this isn’t usually a problem. But when attempting to transition CSP to hardware, long complex-valued data vectors are problematic. For such cases, we can turn to the time-smoothing method (TSM) of SCF estimation, which we discuss in another post.

## 25 thoughts on “CSP Estimators: The Frequency-Smoothing Method”

1. AN says:

Hey! This is a great blog! Learning a lot here. Can you please provide MATLAB codes to generate these plots?
Thanks!

• Thanks AN. The CSP Blog is a self-help blog. I don’t give out much code. Let me know if you get stuck, though, trying to write your own. That’s the only way you’ll really understand CSP.