CSP Estimators: The Time Smoothing Method

The non-blind spectral-correlation estimator called the TSM is favored when one wishes to avoid long FFTs.

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.

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. This amounts to a phase compensation of each cyclic periodogram before it enters the averaging operation.

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 final TSM estimator expression is

\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
Figure 1. TSM-based power spectral density estimate for a rectangular-pulse BPSK signal with bit rate 1/T0=0.1 and carrier frequency offset fc=0.05. The theoretical power spectrum is a squared sinc() function.

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
Figure 2. TSM-based non-conjugate spectral correlation function estimates for a rectangular-pulse BPSK signal with bit rate 1/T0=0.1 and carrier fc=0.05.

and the conjugate-SCF estimates are:

tsm_conj_scfs
Figure 3. TSM-based conjugate spectral correlation estimates for a rectangular-pulse BPSK signal with bit rate 1/T0=0.1 and carrier fc=0.05.

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), or 1/N.

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. The cycle-frequency resolutions for the two estimates are also the same at the reciprocal of the processed data-record length \Delta\alpha \approx 1/(NM) = 1/32768. So the two estimates have comparable time, frequency, and cycle-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.

More on the TSM Phase-Compensation Factor

Some readers question the need for the phase-compensation factor in (5), which is given by

\displaystyle e^{-i 2 \pi \alpha j N}.

In fact it is not needed whenever it is always equal to a constant, such as one. This happens whenever the product \alpha j N is equal to an integer for all values of j. A sufficient condition is that the cycle frequency \alpha is equal to a harmonic of 1/N, that is, \alpha = k/N. However, the phase factor is needed otherwise. Here are two confirming examples (do your own when you implement the TSM).

First, let’s look at the case where the cycle frequency is just the bit rate of our usual BPSK signal, or \alpha = 1/10. Let’s pick the TSM block length N=128 and process a total of NM = 32,768 samples. We’ll plot the SCF estimate magnitude for the correct TSM and for the TSM without the phase-compensating factor:

tsm_and_incorrect_tsm_alpha_0.1
Figure 4. Illustration of the absolute need for the phase-compensating factor in the TSM when the cycle frequency is not a multiple of the reciprocal of the block length N.

In this case, \alpha = 1/10 \neq k/128 for all integers k, and the phase-compensating factor is clearly needed.

On the other hand, consider a BPSK signal with bit rate 1/8 and carrier 0.05. Choose the bit-rate cycle frequency \alpha = 1/8, which is in fact equal to k/128 for k = 16:

tsm_and_incorrect_tsm_alpha_0.125
Figure 5. Illustration of the lack of need for the phase-compensating factor in the TSM in the case where the cycle frequency is a multiple of the reciprocal of the block length, \alpha = k/N.

And so in this case the phase-compensating factor is irrelevant since it is equal to one for all the TSM blocks, and the two TSM estimates are identical.

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.

18 thoughts on “CSP Estimators: The Time Smoothing Method”

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

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

  2. I am having trouble understanding and implementing the delay needed to estimate the cyclic periodogram from the FFT. To start with, I am not completely sure whether I understand the need for it. I take that the idea behind this is to actually estimate the band-pass zero centered signals from the FFT bins in the spectrogram, but I am not sure of why this is needed (the covariance should cancel out this delay as it is the same for the two signals). Furthermore, I’ve been testing my implementation and with/without the delay compensation results seem fairly similar (if any, the version without the delay seems slightly better). It may be that the signal length I am using is big enough to “average out” the delay effects, but I have not been able to find a case where the delay compensed option is significantly better (I take I may have an error on my own implementation).

    Any insight that you can provide?

    Thanks and really incredible job with this blog.

    1. Thanks for your interest Pablo!

      What is the TSM block length N in samples, and what are the cycle frequencies \{\alpha\} that you are passing to the TSM (normalized frequencies)?

      Because this phase-compensating factor has historically been met with confusion and resistance, I’m working on an addendum to the TSM post… I’m hoping to be much more clear about it than I’ve been in the past. Sorry about that!

      1. I’ve been playing around with different values and signals. Using your BPSK signal and the values given in the post (N = 256, etc), I’ve been able to reconstruct something similar to your plots, both with no phase-compensation and with my maybe incorrect one. With respect to the cycle frequencies, I am generating the complete square matrix (every combination of spectral correlation). For the delay compensation I am using normalized frequencies (alpha = [0, 1/256, 2/256, … 255/256]) and multiples of N for D, i.e: [0, 256, 512, …] up to the last segment. I understand that this delay does not change between spectral components (FFT bins). Not sure if this gives you an idea of where my mistake is. If you have any reference I can use to better understand the phase compensation, I am willing to read it.

        Even with my likely ill-formed SCF I have been able to see some interesting features of different modulations, this is a really interesting topic.

        Thanks a lot!

        1. Let’s say you are processing a total of T samples, and your TSM block size is N samples. Recalling the basics of cycle-frequency resolution, cycle frequencies separated by about 1/T will be distinguishable. So a full spectral correlation analysis must consider the cycle frequencies k/T, k = 0, 1, 2, \ldots T-1. If you consider instead the set of cycle frequencies k/N, k = 0, 1, 2, \ldots N-1, you will be missing a lot of cycle frequencies (more and more as T/N increases with N fixed). In either case, the phase-compensating factor will be one for all cycle frequencies and all block indices if both T and N are dyadic (powers of two).

          More to the point, what happens with your two versions (one with my phase correction, one without) when you use my rectangular-pulse BPSK signal as an input and you specify the single cycle frequency 1/10 (which is the bit rate of the rectangular-pulse BPSK signal)? Here the cycle frequency is not the reciprocal of a dyadic number. For my TSM implementations, only the one with the phase-compensator works.

          1. I now see where my mistake was. I thought I was using N=256, but actually, I was using N=100, so 1/10 was one of the values generated. This way, the delay compensation plays no role, and so the estimation held up without delay compensation. Changing to N=256, the delay compensation is indeed needed and works as expected.

            Thanks a lot for your support. I will keep on digging in these concepts!

  3. Quick question: 1) “The spectral resolution is 0.005 Hz (164 points in g(f))”.. would it actually be 1/164 which is 0.0061 or do you take the RMS bandwidth which in the case of rectangular window would match the support?

    Thanks

    1. Masoud: Thanks for visiting the CSP Blog and leaving a comment!

      In the FSM post example, I mentioned that the processing block length is 32768 samples. That means each FFT bin would represent 1/32768 Hz (normalized frequencies here). And therefore the samples in the cyclic periodogram are separated by that same amount, 1/32768 Hz. If the smoothing window width is 164 frequency samples, the spectral resolution of the smoothed cyclic periodogram is therefore 164/32768 = 0.005 Hz. I also mention using zero-padding in the FSM post. To maintain an effective spectral resolution of 0.005 Hz, the number of points in the smoothing window in that case is 328 because the FFT bin width is then 1/65536.

  4. Another question about the phase correction term. It seems that I am having issues to align the cyclic periodogram whenever the expected cyclic frequency does not exactly fall on a frequency bin of the FFT. I am using a QPSK signal at 10Hz carrier and 1.25Hz symbol rate. My FFT has 600points which gives me a resolution of 0.1Hz. with my implementation, the cyclic frequencies at k*2.5Hz are all properly estimated while the frequencies at (2*k-1)*1.25Hz are all zero.
    Looking at the real and imaginary parts of each cyclic periodograms clearly shows an energy peak on the two closest FFT bins (e.g. at 1.2Hz and1.3Hz, 3.7Hz and 3.8Hz) however the phase of these peaks are not aligned, and the average in [equ(5)] makes those peaks vanish. On the other hand, averaging the absolute value of the cyclic periodograms (instead of the real and imaginary parts individually as per [equ (5)]) clearly show peaks at the expected location.
    For probably the same reason, using an FFT of 602 samples (instead of 600) makes all spectral correlation peaks disappear (except at \alpha=0 of course). Using the average of the absolute value (instead of the average of the real/imaginary parts) shows strong peaks at the expected locations.
    Do you think this a normal behavior, or do you suspect that I have an issue with my implementation? Any advice would be highly appreciated.
    Thank you for your all the dedication on the work you have been doing with this blog and papers.

    1. Thanks for stopping by, Arnaud! All signs point to an error in the phase-correction term that you apply. You didn’t include a description of that term, so I can’t say for sure that it is wrong. Your description implies that the sampling rate is 60 Hz, and your QPSK signal has 48 samples per symbol. So, the even harmonics of the symbol rate are cycle frequencies that contain an integer number of FFT bins (for TSM FFT length of 600), and the odd ones don’t. Since you get good SCF estimation (I assume “the cyclic frequencies … are all properly estimated …” means that the SCF is properly estimated, not the CFs) for even harmonics of the symbol rate, I suspect that your applied phase-compensation factor is always one.

      No, this is not normal behavior. Exactly what is your phase-compensation factor?

      1. Chad, thank you very much for your reply.

        You are right the sampling frequency is indeed 60Hz.

        Given the cyclic periodogram I^\alpha (j, f), where \alpha is the non-normalized cyclic frequency, and f the non-normalized frequency [f=fs*(0:(Nfft-1))/Nfft] (fs is the sampling frequency and Nfft is the size of the FFT vector).
        I obtain the periodogram by multiplying the column vector FFT by its conjugate transpose (non-conjugate periodogram). The resulting “matrix” [X(f1)*conj(transpose(X(f2)))]_m,n has a frequency and cyclic frequency given by f_m,n=f1_m/2+f2_n/2 and \alpha_m,n=f1_m-f2_n.

        The phase correction term I am using is exp(-2*i*pi*\alpha/fs*offset(j)), where offset(j) is the starting point (i.e. index value) of the jth cyclic periodogram block. Since \alpha is not normalized in my implementation, I need to divide it by the sampling frequency to make it consistent with equ5.

        1. I obtain the periodogram by multiplying the column vector FFT by its conjugate transpose (non-conjugate periodogram). The resulting “matrix” [X(f1)*conj(transpose(X(f2)))]_m,n has a frequency and cyclic frequency given by f_m,n=f1_m/2+f2_n/2 and \alpha_m,n=f1_m-f2_n.

          So you are using the outer product of the FFT with itself to obtain a matrix that contains several cyclic periodograms. The main diagonal is the conventional periodogram. The off diagonals are the cyclic periodograms. Agree?

          The phase correction term I am using is exp(-2*i*pi*\alpha/fs*offset(j)), where offset(j) is the starting point (i.e. index value) of the jth cyclic periodogram block.

          Your description of the phase compensation term is ambiguous. The matrix you obtain contains multiple cyclic periodogram values, and so each value of the matrix will require its own compensation term, which is related to the particular value of cycle frequency for the value. That is, there is no single phase compensation term for you. I can’t tell from your description whether you are applying the correct term to each individual matrix element or are applying a single compensation term to the entire matrix.

          I recommend that you first code and debug a single-cycle-frequency TSM algorithm. Make sure that it works for arbitrary TSM block lengths and arbitrary cycle frequencies. That is, the cycle frequencies that are implied by your outer product are a small subset of the possible cycle frequencies.

          1. “So you are using the outer product of the FFT with itself to obtain a matrix that contains several cyclic periodograms. The main diagonal is the conventional periodogram. The off diagonals are the cyclic periodograms. Agree?”

            Yes, this is exactly what I have been doing.

            “I can’t tell from your description whether you are applying the correct term to each individual matrix element or are applying a single compensation term to the entire matrix.”

            I am using a different compensation term for each individual matrix element. If I_m,n(j) is the jth “periodograms matrix” obtained from the outer product of the jth FFT vector, the compensation term I am using for the element _m,n (i.e. row m, column n) is exp(-2*i*pi*\alpha_m,n/fs*offset(j))

            “I recommend that you first code and debug a single-cycle-frequency TSM algorithm. Make sure that it works for arbitrary TSM block lengths and arbitrary cycle frequencies.”

            Thank you for the recommendation. I will try to implement the TSM for a single frequency and follow up if I have some more issue.
            Thank you again for the helpful feedback.

  5. Your blog and information is fantastic. I’m confused by something here though, and it’s pretty basic. Apologies.

    In (1) you define the DFT in a way I’ve never seen that seems odd to me. The DFT definition is usually agnostic to sampling rate, and even if you normalize such that your sampling interval is 1 the DFT equation should be the same. Except your definition does not normalize the numerator of the exponent term by N. Even for normalized sampling rates shouldn’t the numerator be (-i 2 pi f t / N)?

    1. Okay, I think I see what’s going on. It seems to me there are *two* normalizations here:

      1. Ignore actual sample rate and pretend our data is sampled on 1 second intervals.
      2. The DFT is defined such that the output frequency isn’t indexed like a normal DFT but defined on the interval [0, (N-1)/N], or simplified a bit, [0, 1).

      Of course, that might work well for showing the math and equations. If I’m using an actual computer to compute the DFT I’m going to have sampled indices on number 2 above. But I think I can get there. I’m just running into issues with the phase normalization term. When I try to apply it “as is” (with alpha being the baud in Hz divided by the actual sample rate) it’s not working. So I’m guessing there’s something in this DFT definition here that will explain my troubles.

      1. Thanks for stopping by the CSP Blog cubingthesphere! And for the comments.

        In (1) you define the DFT in a way I’ve never seen that seems odd to me. The DFT definition is usually agnostic to sampling rate, and even if you normalize such that your sampling interval is 1 the DFT equation should be the same. Except your definition does not normalize the numerator of the exponent term by N. Even for normalized sampling rates shouldn’t the numerator be (-i 2 pi f t / N)?

        Well, I am indeed hiding the exact values of f in (1). The usual way of defining f in the DFT is as MATLAB does: f = k/N for k = 0, 1, \ldots, N-1:

        But engineers like to see f = k/N for k=-N/2, -N/2+1, \ldots, N/2-1, so MATLAB provides fftshift.m.

        In other words, in my (1), time (t and u) are integers, which is consistent with a unit sample rate, but f is a value like k/N, the exact nature of k is up to you. Mathematicians typically like the MATLAB definition, electrical engineers like the alternative definition (supplied by fftshift.m acting on the output of fft.m in MATLAB).

        1. Ignore actual sample rate and pretend our data is sampled on 1 second intervals.

        Yes, and this is pervasive at the CSP Blog. See my lengthy recent post on the sampling rate for a discussion (scroll down to near the end if you don’t want the tutorial).

        2. The DFT is defined such that the output frequency isn’t indexed like a normal DFT but defined on the interval [0, (N-1)/N], or simplified a bit, [0, 1).

        Well, I’m hiding the indexing. Perhaps I shouldn’t have done that. The TSM will work for whatever you use for f, as long as it is consistent from TSM block to TSM block.

        If I’m using an actual computer to compute the DFT I’m going to have sampled indices on number 2 above.

        Well… see above. The choice is yours. ‘To fftshift or not to fftshift, that is the question!’

        I’m just running into issues with the phase normalization term. When I try to apply it “as is” (with alpha being the baud in Hz divided by the actual sample rate) it’s not working. So I’m guessing there’s something in this DFT definition here that will explain my troubles.

        I doubt that the definition of f is your problem. Here are some TSM estimates for the rectangular-pulse BPSK signal using the TSM except I just removed all fftshifts:

        What do you think?

        1. Thanks for the reply!

          You’re right, that wasn’t my issue. I noticed an error with which leg I was conjugating. Well, or I could change the sign on the phase term. Either way that seemed to be my problem. Thanks for the explanation.

Leave a Reply to cubingthesphereCancel reply

Discover more from Cyclostationary Signal Processing

Subscribe now to keep reading and get access to the full archive.

Continue reading