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

The basic idea is to segment the provided data record into M contiguous blocks of N samples each, compute the cyclic periodogram for each block, and average the results. Since we will likely use the FFT to compute the Fourier transform, we will be viewing each N-sample block as if its time samples correspond to t = 0, 1, \ldots, N-1, and so the cyclic polyspectrum formula of My Papers [6] will have to be slightly modified to take into account the actual temporal start time for each block.

So let’s consider the Fourier transform (DFT) of a block of data that is shifted from the origin by some amount of time u,

X(u, f) = \displaystyle\sum_{t=0}^{N-1} x(t+u) e^{-i 2 \pi f t}. \hfill (1)

The periodogram and cyclic periodogram are then functions of time offset u as well,

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

and

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

and similarly for the conjugate cyclic periodogram. The TSM estimate of the SCF is simply the average value of the cyclic periodogram over all available values of u,

\hat{S}_x^\alpha (f) = h(u) \otimes I^\alpha(u, f), \hfill (4)

where h(u) is some pulse-like temporal window. In practice, the FFT is used to create each cyclic periodogram, so their relative phases are no longer taken into account. According to our Fourier transform result for a delayed signal, however, we can easily take this into account by multiplying each cyclic periodogram by e^{-i 2 \pi \alpha D}, where D represents the left edge (starting point) of the subblock. For blocks having length N samples, then, the value of D for the jth block is simply jN. Our TSM estimator is then

\hat{S}_x^\alpha(f) = \displaystyle\frac{1}{M} \displaystyle\sum_{j=0}^{M-1} \left[\tilde{I}^\alpha(jN, f) e^{-i 2 \pi \alpha j N} \right], \hfill (5)

where \tilde{I}^\alpha(jN, f) is just the cyclic periodogram created from the jth block of N samples using the FFT. Notice that when the cycle frequency is set to zero, the SCF estimate is an estimate of the PSD, and the TSM just averages M periodograms, as in the Bartlett spectrum estimation method. Here is the TSM (Bartlett) PSD estimate for our rectangular-pulse BPSK signal:

tsm_psd

For this PSD estimate the data-record length is 32,768 samples and the TSM block length is 256 samples, leading to M = 32768/256 = 128 blocks. Recall that the bit rate for the BPSK signal is 1/T_0 = 1/10 and the carrier frequency is f_c = 0.05 (in normalized frequency units).

The TSM PSD estimate matches the FSM PSD estimate in the FSM post.

The TSM-based spectral correlation function estimates for the BPSK signal’s non-conjugate cycle frequencies are shown below:

tsm_nonconj_scfs

and the conjugate-SCF estimates are:

tsm_conj_scfs

Again, these TSM estimates match quite well with the FSM estimates.

The reason the TSM and FSM estimates match so well is that the temporal and spectral resolution parameters of the estimates are similar. For both methods, the temporal resolution is equal to the data-record length (32,768 samples). For the FSM, the spectral resolution of the estimates is equal to the width of the frequency-smoothing window g(f), and for the TSM, the spectral resolution is equal to the intrinsic spectral resolution of each cyclic periodogram, which is equal to the reciprocal of the TSM block length (in normalized units).

For the FSM results in the FSM post, the spectral resolution is 0.005 Hz (164 points in g(f)), and for the TSM results in this post, the spectral resolution is 1/256 = 0.0039 Hz. So the two estimates have comparable time and frequency resolution parameters, and so produce similar results. The relationship between estimator quality and the temporal, spectral, and cycle-frequency resolutions is discussed in this post.

14 thoughts on “CSP Estimators: The Time Smoothing Method

  1. aapocketz says:

    When using the TSM, you must calculate the conjugate cyclic periodogram by conjugate multiplication of the Fourier transform at offset +alpha/2 and -alpha/2.

    Functionally, that seems to mean that you are multiplying the Fourier transform that has been circularly shifted left by 1 bin, by a conjugate circularly shifted right by 1 bin, to get the cyclic value at alpha = 2*Fs/N. You would circular shift by 2 bins to get alpha at 4*Fs/N, and so on all the way up to N bins, the maximum alpha shift at Fs.

    Does that sound right? That means that including negative alpha values, the most alpha points will be N, and the max alpha value will be Fs, and the min alpha will therefore be 2*Fs/N?

    Is there a way to get the alpha resolution (2*Fs/N) as fine as the frequency resolution (Fs/N) via some other mechanism? Do the FAM and SCCA algorithms have the same limitations on alpha?

    Thank you for your blog and papers and time!

    • When using the TSM, you must calculate the conjugate cyclic periodogram by conjugate multiplication of the Fourier transform at offset +alpha/2 and -alpha/2.

      The conjugate cyclic periodogram is used to estimate the conjugate spectral correlation function and, unfortunately, it involves no conjugations. That is, the conjugate cyclic periodogram is X(f+a/2)X(a/2-f), whereas the (non-conjugate) cyclic periodogram goes like X(f+a/2)X*(f-a/2). So that is just terminology. I’m assuming that the rest of the comment is concerned with the non-conjugate cyclic periodogram and spectral correlation function.

      Functionally, that seems to mean that you are multiplying the Fourier transform that has been circularly shifted left by 1 bin, by a conjugate circularly shifted right by 1 bin, to get the cyclic value at alpha = 2*Fs/N.

      Right.

      That means that including negative alpha values, the most alpha points will be N, and the max alpha value will be Fs, and the min alpha will therefore be 2*Fs/N?

      You must use shifts that are equal to integer numbers of the FFT bins, yes, but there is nothing preventing you from using different shifts for the two involved shifted Fourier transforms. For example, you could form X(f + 1/N)X*(f). So that is a shift of one FFT bin for X(f + a/2) and no shift for X(f – a/2). This product ends up providing an estimate of S^b(g), where b = 1/N and g = f-1/N. So the minimum cycle frequency that you can look at (besides zero) is 1/N (or in your terminology Fs/N). Agree?

      Is there a way to get the alpha resolution (2*Fs/N) as fine as the frequency resolution (Fs/N) via some other mechanism?

      I’m working on a post on the three different kinds of resolutions involved in CSP (temporal, spectral, and cycle), but for now I think my comment above shows that the cycle resolution is Fs/N. Yes, this is equal to the native spectral resolution of the DFT using N points, but be careful because that is not the spectral resolution of the spectral correlation measurement in general. For the TSM, where we have N total samples and the FFTs have size M << N, the spectral resolution is 1/M (Fs/M in physical Hz). For the FSM, it is the width of the spectral smoothing window, which usually contains many of the N FFT points. For the SSCA, it is the reciprocal of the number of channels in the front-end channelizer.

      It turns out the cycle resolution is always about equal to the reciprocal of the total number of samples processed, whereas the spectral resolution varies by estimator and is strongly affected by estimator parameters.

      Great comments, aapocketz, thanks much.

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 )

Twitter picture

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

Facebook photo

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

Google+ photo

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

Connecting to %s