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

where for . The periodogram is the normalized squared magnitude of the DFT,

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

The PSD can be accurately estimated by smoothing the periodogram,

where is some unit-area 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 is large and the width of 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

for . 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 samples and a rectangular smoothing window with width frequency bins, which is about (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

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

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

for . All that is required to implement this algorithm is an FFT, complex multiplication, and a general convolution algorithm (such as MATLAB’s conv.m or filter.m). The conjugate SCF is estimated by frequency smoothing the conjugate cyclic periodogram given by

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 samples and a rectangular with width frequency points, or about Hz. Recall the sampling frequency is set to unity here, the signal’s bit rate is , and the carrier frequency is . The -sample data record is zero-padded by a factor of two prior to forming the cyclic periodogram.

**Choice for the Frequency-Smoothing Window **

We used the rectangular pulse in our estimator examples for a good reason: computational cost. For arbitrary , a general-purpose convolution routine can be used to find the SCF estimate, but the convolutions can be expensive. 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 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 (left edge of the window) and adding the one that is newly included in the support (right edge of the window).

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 ratio (SINR), 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.

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.

Hi Chad

I think you chose complex-valued waveform of BPSK signal to estimate the PSD and cyclic spectrum, right?

And how can I judge the simulation results of the FSM priodogram and cyclic priodogram of BPSK signal correct or not?

The value range of n in Equation (10) can not reach 0 and N. The reason is that for certain cycle frequency, “n+a1=0” , right?

Yes, the BPSK signal I use throughout the CSP Blog is complex-valued.

Do you mean how can you judge whether my plots are correct? Or your plots? My plots are consistent with theory in the cases where I know the theory, such as PSK and QAM signals. So you can use my plots in the BPSK posts (here and here) or the gallery posts (here and here) to compare with yours.

Yes, when in the non-conjugate SCF, you cannot use all the indices in the sum over unless you are willing to view the FFT output as a periodic function (with period ). But generally you don’t want to do that, and the edges of the cyclic periodogram will contain zeros. The larger becomes, the more zeros you’ll see at the ends (this also corresponds to the diamond shape of the support for the SCF).

But why the estimation performance of FSM method decreased as I increased the length of data record? The simulation parameters I chose were the same as yours.

It isn’t possible to substantively answer this question. Not enough information!

Hi Chad

According to the equation (25) in Gardner’s paper “Measurement of Spectral Correlation”, I think the result of equation (3) should be divided by the frequency bins number of smoothing window. Do as what I advised, the FSM-based PSD matches well with the theoretical PSD for BPSK signal.

What about your opinion?

My (3) is a bit more general than Gardner’s (25), which smooths using a unit-area rectangle. My can be any pulse-like unit-area window function, but I do use a rectangle a lot too. The factor you say is missing is actually embedded in the definition of : It has unit area. Therefore if it is a rectangle, it has height equal to the reciprocal of the width. I’ve added “unit area” to the post to clarify.

The Figure 1 in this post is similar to what your code make_rect_bpsk.m. But you did not use the convolution after the FFT, but the plot looks same.

So my question is related to this: To estimate the SCF, we should take the FFT of non conjugate CAF and then we need to convolve this with a smoothing window g(f). Does it give any extra benefit of using the g(f)? What is a good size of g(f)?

Is the FFT of non conjugate CAF called the cyclic Periodogram?

I don’t quite understand this. The plot in Figure 1 looks the same as what?

No, to estimate the SCF, take the FFT of the

data, shift it up and down by , multiply, divide by the length of the data, then convolve the result with a smoothing window. See the posts on the frequency-smoothing and time-smoothing methods of spectral correlation estimation for details and examples.Yes. The main point of the Periodogram post is that if you don’t average (“smooth”) the periodogram in either time or frequency, you will have a highly erratic function of frequency. Note also that the reliability (variance) of a spectral correlation estimate is inversely proportional to the product of the data-block length and the frequency resolution. For the frequency-smoothing method, the frequency resolution is the width of . So, to generate a reliable (low variance) estimate, you want to use large data-block lengths and large (coarse) frequency resolution (wide ).

A good width for is typically 1-5% of the sampling rate.

No. The cyclic periodogram is defined in (3) and (5) in the Periodogram post. The Fourier transform of the non-conjugate CAF is the non-conjugate spectral correlation function.

Thanks for the reply. I meant, the Figure 1 in this post is similar to what your code make_rect_bpsk.m. generates.

Thanks for the reply. I need some more clarification. You said “To estimate the SCF, take the FFT of the data, shift it up and down by \alpha/2, multiply, divide by the length of the data”.

1. Shift it up and down by alpha/2: Is it not giving back the same matrix if I shift it up by alpha/2 and then down by the same amount. Also is this shifting the same what fftshift does? If not same, then what this shift up and down by alpha/2 is doing here?

2. Then you said “multiply”: what do I multiply with what? I understood the divide by the length of the data.

3. Is the cyclic periodogram calculated as: FT[x(t)*exp(j*2*pi*alpha/2*t)] * conj{FT[x(t)*exp(-j*2*pi*alpha/2*t)]}

Well, you shift the FFT up by to get and shift it down by to get . Those are the two quantities you need to multiply together to get the cyclic periodogram.

No. MATLAB’s fftshift.m swaps the second half of a vector for the first half: fftshift([1 2 3 4]) = [3 4 1 2]. This is identical to circshift([1 2 3 4], 2). The shifting we need to do in CSP depends on the value of the cycle frequency .

The mathematics of CSP says that you can estimate the spectral correlation function by averaging the cyclic periodogram, which is the product of two frequency-shifted Fourier transforms. More physically, the spectral correlation function is the temporal correlation between two narrowband spectral components of a signal, after those components have been isolated and frequency-shifted to zero frequency. And those two frequency components have frequencies specified by and in the usual parameterization.

See answer to Question 1 above.

I’m not sure about the plus and minus signs in the exp() functions, and you’ll need to divide by length(x(t)), but, yes, this is one way to obtain the cyclic periodogram. I typically don’t use this method, because if you’re careful, you can do the frequency shifting in the frequency domain using only things like circshift.m, which is much cheaper than creating the complex sine waves and multiplying them by the signal x(t) and taking two (not one) Fourier transforms.

I find an interesting phenomenon. To meet the condition that the approximated value, obtained by integrating PSD values over frequency and multiplying that sum by the frequency increment, must be equal to the estimate of power. Therefore, different PSD estimating methods have different frequency increment and different PSD value. For example, periodogram based PSD values are much greater than that based on TSM and FSM method.

How do you explain that phenomenon?

Yes, the different methods can have different frequency-bin increments as well as different spectral resolutions. This means that for any particular given frequency , the PSD estimates from the various methods will likely have different values. But for properly set up measurements (the measurement spectral resolution is well-matched to the data’s spectral characteristics, the data-block length is sufficient to allow averaging away erratic components), the values will be quite similar, as you can see from my comparisons between the FSM and TSM outputs.

And the sum of the PSD estimate values, multiplied by the frequency increment, should be very close for each method.

Yes, the periodogram is highly erratic, and is a terrible estimator of the PSD, so there will be periodogram values that are very much greater than the true PSD value and periodogram values that are very much less than the true PSD value–that’s why we need the averaging.

Convolving the periodogram of x(t) with g(f) would obtain PSD estimate with length of N+164-1. So what is the frequency increment between the adjacent frequency bins of that PSD estimates? Is the frequency increment fs/(N+164-1)?

This is a convolution question, not a CSP question! The spacing between the samples does not change as a result of convolution.

The posts are very helpful for the beginners like me. It would be better if you supply more detailed unit labels on the figures. I think the unit of the nomalized frequency is [cycle/sample], and of the PSD is [dB/(rad/sample)]. I am not sure whether the units I added here is correct? Or maybe it is a traditional label manner in the CPS society.

Thanks for visiting the CSP Blog and leaving a comment Smoon! I appreciate it.

I agree with you about normalized frequency–strictly speaking I should label axes with [Cycles/Sample] or the like. But when the actual physical sampling rate really is , then the unit is also Cycles/Second, or Hz. Mostly I’m following convention, which is to use ‘Hz’ for both physical frequency and normalized frequency. To translate any frequency axis from normalized frequency to physical frequency, just multiply by the actual sampling rate.

And yes, the physical unit of a power spectrum is Watts/Hz. Normally expressed in decibels, but often in decibels relative to some basic power level, such as one Watt (dBW) or 1 milli Watt (dBm). In almost all of my posts, I’m unconcerned with the actual power level of the data–even when I capture data using my lab equipment I don’t keep track of the relationship between the integers I obtain from a sampler and the power level of the electromagnetic wave incident on the attached antenna. So I rather lazily just use decibels and the relationship to actual physical power is suppressed.

What is the ‘CPS society’?