CSP Estimators: The Strip Spectral Correlation Analyzer

In this post I present a very useful blind cycle-frequency estimator known in the literature as the strip spectral correlation analyzer (SSCA) (The Literature [R3-R5]). We’ve covered the basics of the frequency-smoothing method (FSM) and the time-smoothing method (TSM) of estimating the spectral correlation function (SCF) in previous posts. The TSM and FSM are efficient estimators of the SCF when it is desired to estimate it for one or a few cycle frequencies (CFs). The SSCA, on the other hand, is efficient when we want to estimate the SCF for all CFs.

The SSCA is a kind of time-smoothing method and a detailed examination of the method and its properties can be found in The Literature [R3-R5]. Here I’ll briefly describe the estimator mathematically, then try to explain why the rather odd-looking estimator form produces estimates of the spectral correlation function. Finally, I’ll provide a few processing examples. In future posts I’ll contrast the SSCA performance and computational cost with those for the FSM and TSM.

The block diagram for the SSCA looks like this:

ssca_block_diagram

Essentially, the input data is first channelized, and then each of the narrowband outputs of the channelizer is multiplied by the original input signal and Fourier transformed. The channelizer is usually implemented by a sliding FFT. So the method clearly harnesses the computational power and efficiency of the FFT, leading to its overall low cost compared to other exhaustive-CF spectral correlation estimators (such as the FSM and TSM).

Mathematically, a particular cross (between arbitrary signals x(t) and y(t)) spectral correlation function point estimate is given by

\displaystyle S = S_{xy_T}^{f_k +q\Delta\alpha} (n, \frac{f_k}{2} - \frac{q\Delta\alpha}{2})_{\Delta t} \displaystyle = \sum_{r=-N/2}^{N/2 -1} X_T (r, f_k) y^* (r) g(n-r) e^{-i2\pi q r/N}, \hfill (1)

where the demodulates are defined as

\displaystyle X_T(n,f) = \sum_{r=-N^\prime/2}^{N^\prime/2 - 1} a(r) x(n-r) e^{-i 2 \pi f(n-r)}. \hfill (2)

Let’s try to show that this expression can be manipulated into a form such that we can more easily see that the spectral correlation function is estimated. The first step is to substitute the expression for the demodulate into the SSCA point estimate and simplify,

\displaystyle S = \sum_r \left[ \sum_{n=N^\prime/2}^{N^\prime/2-1} a(n) x(r-n) e^{-i2\pi f_k (r-n)} \right] y^*(r) g(n-r) e^{-i2 \pi qr/N} \hfill (3)

\displaystyle = \sum_{r,n} (a(n)g(n-r)) (x(r-n)y^*(r)) e^{-i2\pi f_k (r-n)} e^{-i2\pi qr/N} \hfill (4)

\displaystyle = \sum_n a(n) e^{i2\pi f_k n} \left( \sum_r g(n-r) [x(r-n) y^*(r)] e^{-i2\pi (f_k + q/N)r} \right). \hfill (5)

Define \alpha_0 = f_k + q/N. Then

\displaystyle S = \sum_n a(n) e^{i2\pi f_k n} \left( \sum_r g(n-r) [x(r-n)y^*(r)] e^{-i2\pi \alpha_0 r} \right). \hfill (6)

Now, substitute n = -n^\prime, rename n^\prime to n, substitute r = r^\prime -n/2, and finally rename r^\prime to r:

\displaystyle S = \sum_n a(-n) e^{-i2\pi f_k n} \left( \sum_r g(-n-r) [x(r+n) y^*(r)] e^{-i2\pi \alpha_0 r} \right) \hfill (7)

\displaystyle = \sum_n a(-n) e^{-i2\pi f_k n} \left( \sum_r g(-n-r+n/2) [x(r+n/2)y^*(r-n/2)] e^{-i2\pi \alpha_0 (r-n/2)} \right). \hfill (8)

Assume that a(-n) = a(n). The inner sum is just an estimate of the cyclic autocorrelation, with data tapering (windowing) and with limits on the sum that depend on the lag n. This is the cyclic correlogram (The Literature [R5]).  Thus, we have the following relation,

\displaystyle S = \sum_n a(n) e^{-i2\pi (f_k - \alpha_0/2)n} \tilde{R}_{xy_{\Delta t}}^{\alpha_0} (n) \hfill (9)

\displaystyle = S_{xy_{\Delta t}}^{\alpha_0} (f_k - \alpha_0/2)_{1/\Delta \! f}. \hfill (10)

Here we can identify the (f, \alpha) pair from our usual parameterization of the spectral correlation function.  First, we have \alpha = \alpha_0 = f_k + q/N = k/N^\prime + q/N. The spectral frequency f is the argument of S(\cdot), and is given by f = f_k - \alpha_0/2 = k/(2N^\prime) - q\Delta\alpha/2 = k/(2N^\prime) - q/(2N).

Therefore, the point estimate is confirmed. We see that the time-smoothed SSCA estimate is approximately equal to a Fourier-transformed cyclic autocorrelation estimate. The requirement for reliability on this estimate is that the length of time used to estimate the cyclic correlogram, \Delta t = N, is sufficiently large and much larger than the approximate span of the function, which is contained in the window a(n), which itself has width 1/\Delta\!f = N^\prime.  That is, \Delta t \Delta\!f \gg 1 or N \gg N^\prime once again.

The demodulates are formed by short-duration (N^\prime) FFTs, which are often called the channelizer FFTs, and the final point-estimate outputs are formed by long (N) FFTs, which are often called the strip FFTs, because each one of those FFTs produces a set of point estimates that lie on strip parallel to the line \alpha = -2f in the (f, \alpha) plane.

Numerical Example: Rectangular-Pulse BPSK

Let’s turn to our usual example: the textbook rectangular-pulse BPSK signal. This signal has normalized bit rate of 1/10 and normalized carrier frequency of 0.05. The non-conjugate cycle frequencies are harmonics of the bit rate, and so are equal to \alpha = k/T_{bit} = k/10, and the conjugate cycle frequencies are the non-conjugate cycle frequencies plus the doubled carrier \alpha = 2f_c + k/T_{bit} = 0.1 + k/10.

We process 65536 samples of the rectangular-pulse signal using a straightforward C-language implementation of the SSCA. The number of strips, N^\prime, is set to 64. This SSCA processing results in 65536 * 64 = 4194304 SCF point estimates for the non-conjugate SCF and another 4194304 for the conjugate SCF. However, symmetries allow us to look at the non-conjugate function for only the non-negative cycle frequencies. Nevertheless, there are a lot of point estimates in the example. Plotting all of them is a difficult exercise in visualization. Here I plot only the (f, \alpha) pairs that correspond to the 500 largest values of spectral correlation magnitude:

scf_nc_top_500

And the 500 largest values of conjugate spectral correlation magnitude:

scf_c_top_500

These two plots show that the algorithm is producing relatively large SCF values for the signal’s true cycle frequencies. There are also some relatively strong SCF magnitudes for cycle frequencies near zero (non-conjugate) and the doubled carrier 0.1 (conjugate) that are not actual cycle frequencies of the textbook rectangular-pulse BPSK signal. These high values are due to the general and important phenomenon of cycle leakage, which we’ll tackle in a future post. For now, think of cycle leakage as akin to spectral leakage in a normal spectrum analysis setting: the PSD estimate for frequencies near to a powerful sine wave signal will be inflated due to proximity to the sine wave and the imperfect frequency response of the estimator.

So, if the SSCA is applied to a new signal, whose cycle frequencies are not yet known, it is clear from the two graphs above that it could be used to blindly find the significant cycle frequencies, simply by thresholding the spectral correlation magnitude. However, there is a better way to sort through the large number of SSCA-produced spectral correlation point estimates, and that is by using the spectral coherence.

Computing and Using the Coherence in the SSCA

Recall from the post on the coherence that it is simply a normalized version of the spectral correlation, in exactly the same way that the correlation coefficient from statistics is a normalized version of the correlation. The non-conjugate coherence is given by

\rho = C_x^\alpha(f) = \displaystyle\frac{S_x^\alpha(f)}{[S_x^0(f+\alpha/2)S_x^0(f-\alpha/2)]^{1/2}}, \hfill (11)

and the conjugate coherence is

C_{x^*}^\alpha(f) = \displaystyle\frac{S_{x^*}^\alpha(f)}{[S_x^0(f+\alpha/2)S_x^0(\alpha/2-f)]^{1/2}}. \hfill (12)

The SSCA produces the values needed for the numerators of the coherences; all we need is the power spectrum values in the denominators. In principle, these can be obtained from the N^\prime values of the power spectrum automatically produced during SSCA computation. However, it is often better to use a more finely sampled PSD estimate, such as can be obtained through use of the frequency-smoothing method or time-smoothing method of spectral correlation estimation, with the cycle frequency set to zero.

Once the coherence is computed in the SSCA implementation, it can be thresholded to yield the significant cycle frequencies. Moreover, a cycle-frequency binning operation can be applied to those cycle frequencies that pass threshold, limiting the output of the algorithm to only unique values of cycle frequency. For our rectangular-pulse BPSK signal, the top 500 non-conjugate coherence values correspond to the following (f, \alpha) pairs:

coh_nc_top_500

and for the conjugate coherence, we obtain

coh_c_top_500

The false cycle frequencies do not pass through the thresholding process when the coherence is used!

After coherence-based thresholding and cycle-frequency binning, we obtain the following output for the non-conjugate cycle frequencies:

Frequency         Cycle Frequency       Coherence        SCF

ssca_nc_normal

and for the conjugate cycle frequencies:

Frequency           Cycle Frequency         Coherence          SCF

ssca_c_normal

Final Example: Captured CDMA-EVDO

There is an EVDO downlink signal at 885 MHz where I live, and I have recorded several seconds of it for analysis purposes. A typical PSD for the captured signal is shown here:

evdo_psd

After applying the SSCA to this EVDO data, we obtain the non-conjugate and conjugate cycle frequency sets shown in the following graphs:

evdo_nc_cfs

evdo_c_cfs

We conclude that the signal has five strong non-conjugate cyclic features for cycle frequencies that are harmonics of 1.2300 MHz (not 1.2288 MHz), and no significant conjugate cyclic features. The CDMA chip rate of 1.2288 MHz is also present, but has smaller spectral correlation and spectral coherence than the feature at 1.2300 MHz. This kind of analysis illustrates the power of the SSCA: Efficient (fast) blind estimation of all second-order cycle frequencies.

34 thoughts on “CSP Estimators: The Strip Spectral Correlation Analyzer

  1. aapocketz says:

    From what I understand the FAM and SCCA algorithm generates a tilted bifrequency plane, that being the value of alpha=0 lie across the primary diagonal of the matrix output, and there is only a single value for alpha = +/- Fs. This makes some sense, due to the fact that at least in the FAM I am conjugate multiplying each frequency bin column against all the other columns. Does that sound correct? Is there a proper technique for correcting this shift, so alpha=0 lies upon the x-axis? Would that cause points to lie upon Fs/(2*N) frequency spacing?

    In a related question, when computing the spectral coherence (SOF), I can take the value of the SCF at alpha=0 which is essentially a PSD. I can take that and do circular shifts conjugate multiplies to get NxN matrix. In theory I can use that matrix to be the denominator in the SOF calculation, but since the FAM outputs in that diagonal fashion, how would I deal with that?

    Also, is it a true statement that the SOF magnitude will always be 1 when alpha=0? Is the SCF or SOF magnitude always symmetric about the alpha=0 axis?

    • From what I understand the FAM and SCCA algorithm generates a tilted bifrequency plane, that being the value of alpha=0 lie across the primary diagonal of the matrix output, and there is only a single value for alpha = +/- Fs.

      The difficulty here is that we must agree what “the matrix output” means to proceed. We can agree right now, however, that the SSCA proceeds strip-by-strip. Each strip FFT is the FFT of the product of the original complex-valued signal with one of the channelizer outputs. Each of these strip FFTs will contain a contribution to the PSD (alpha = 0). And each of the strip FFTs will contain N points, each of which corresponds to N different values of cycle frequency. So, I don’t agree with “there is only a single value for alpha = +/- Fs.”

      Is there a proper technique for correcting this shift, so alpha=0 lies upon the x-axis?

      What I do is create a table of SCF values as the SSCA proceeds with its strip FFTs. I end up with a data structure that has four elements: spectral frequency f, cycle frequency alpha, spectral correlation S, and coherence C. Then I can extract whatever “slices” I want from this structure.

      when computing the spectral coherence (SOF), I can take the value of the SCF at alpha=0 which is essentially a PSD. I can take that and do circular shifts conjugate multiplies to get NxN matrix. In theory I can use that matrix to be the denominator in the SOF calculation,

      Yes, you can use the SSCA-generated PSD values (those elements of the output that correspond to alpha = 0 in the non-conjugate SSCA output). However, there are only M of these values, one for each strip, and M << N, where N is the total number of samples processed (also the length of each strip FFT). So the spectral resolution of that effective PSD estimate is low (coarse) relative to the cycle resolution. This causes problems when the coherence quotient is formed, because the numerator and denominator have different resolutions. One way out of this is to form a stand-alone PSD estimate and use that for the denominator of the coherence. You can control the resolution of the stand-alone measurement.

      but since the FAM outputs in that diagonal fashion, how would I deal with that?

      I’m not sure about the FAM. I’ll need to look into it and do a FAM post!

      Also, is it a true statement that the SOF magnitude will always be 1 when alpha=0?

      Yes, the non-conjugate spectral coherence is always 1 for alpha=0, theoretically, provided that there is no point where the PSD is equal to zero. A good coherence estimator should always output very nearly 1 for this case. If it isn’t, you know you’ve got problems with your estimator, so definitely pay attention to this case.

      Is the SCF or SOF magnitude always symmetric about the alpha=0 axis?

      You mean, is |S^(-a)(f)| = |S^(a)(f)|? Yes, for the non-conjugate SCF and coherence functions. But it is not true for the conjugate functions.

      • aapocketz says:

        Thanks for the reply, I am new to this so still processing the concepts.

        > And each of the strip FFTs will contain N points, each of which corresponds to N different values of cycle frequency. So, I don’t agree with “there is only a single value for alpha = +/- Fs.”

        I am hung up on this concept. Once you do the strip FFTs, you get a resultant matrix which must then be mapped back into Sx in a particular way. Via this mapping, there appears to be only a single non-zero value at alpha=Fs, at F=0 and at alpha=-Fs, at f=0.
        This mapping is mentioned in this reference http://www.dtic.mil/dtic/tr/fulltext/u2/a289815.pdf (24), and described in detail in section 3.3.3. It appears to essentially be a rotation mapping, where the region of support appears “diamond shaped” in the bifrequency plane. In the SSCA, the strips are diagonal slices. The highest alpha value in that region is alpha = Fs, which appears at f=0. Likewise at f=Fs/2 there is only a single value at alpha=0.

        Does that make sense, or did I misinterpret things? I think each strip FFT does correspond to N different values of cycle frequency as you said, but each strip works at an offset of each other, so only one strip actually outputs alpha=-Fs, and the next strip outputs alpha=-Fs+1/Np/2, etc.

        • I think you are understanding just fine; I misinterpreted your earlier comment. Yes, that report by April has the diamond-shaped support region correct (this is independent of the SSCA of course), and you can see the diamond outline in the figures in my original post. I had forgotten about that report; thanks for reminding me aapocketz.

  2. Vincent Voogt says:

    Hi there,

    Thank you for the great informative posts related to this subject. I am a bit confused however regarding the frequency shifting of Sx by an amount of gamma=q*delta_alpha (see eq. 13 in aapocketz’ reference: http://www.dtic.mil/dtic/tr/fulltext/u2/a289815.pdf), such that the total frequency shifted expression can be evaluated using the FFT algorithm. Basically I just cannot mathematically figure out why multiplication of the TS cylic periodogram with exp[j*2*pi*gamma*m*Ts] relates to a shift in the value of alpha (one would expect an frequency shift such that Sx(n,f+gamma) ).

    Thankfully your approach is the other way around, proofing that the SSCA provides us with an evaluation of the Fourier transform of the cyclic correlogram. What confuses me in your derivation though, it that in your third equation (where you substitute the complex demodulate X_T in spectral correlation function point estimate) you end up with an inner sum expression with iterative variable n. This variable is also used as the argument of the function S(n,f) and the window function g(n-r) (which is outside of the inner summation). Isn’t there something going wrong in the subsequent derivation due to this conflict in variable definition?

    I hope you can help me,

    Kind regards,

    Vincent

    • I’m on vacation this week, but I will look into your question next week and provide a substantive reply. I may very well have made a math or notation error! Thanks for the careful read…

    • In the first paragraph of your comment, I believe you are talking about the complex-valued factor that we use in the TSM. I’ve previously tried to explain the origin of this factor in the TSM post. Can you take a look at that post and let me know if it doesn’t help? We can take it from there…

      In the second paragraph, you point out an error in the SSCA development. Yes, there is a problem with the derivation, and some bad choices for notation! I will have to fix this. I don’t believe the basic conclusion will change, though. Thanks very much for pointing me to a flaw in the blog.

  3. William S says:

    Hi Chad,

    I’ve been following your blog recently and I really appreciate the detailed contents on cyclostationary. Currently I’m trying to recreate the simulation for SSCA, but I could not get converged coherently (not normalizing to 1).

    I did try to use FSM method to estimate the PSD, but they seemingly on different scale altogether, although, with scaling that can be avoided. The result still looks that it didn’t normalizes properly.

    With full FSM estimation, I can get it to normalizes properly.

    Do you have any pointers on where I should look into?

    • Can you post a few more details about what is going wrong? It may be helpful to others to see what you are struggling with. You and I may need to take this to email to solve it, but for now, can you let us all know a little more?

      We’ll get to the bottom of it eventually.

        • It’s hard to tell from the FSM plots, but it looks like the coherence (non-conj) for alpha = 0 is not unity across the band. This is the easiest thing to debug, perhaps, since we know that the coherence should be exactly one for every frequency. (You are just dividing the PSD by itself.) So if that is not coming out ‘1’, then you probably have an indexing error. Either that, or you are using two different PSD estimates, one for the numerator of the coherence, and one for the denominator, and they don’t match.

        • Something looks wrong with the FSM PSD in that plot. Can you send me the data files for both the red and blue dots?
          Also, what were the parameters used in the FSM PSD estimate? (block length, width of smoothing window)

          • William S says:

            I forgot to close this topic thread and offer thanks to Chad for helping me out.
            Looking forward to hearing more cs blog in the future.

          • Calvin L says:

            Hi William S and Chad,

            May I know how you managed to resolve your problem of the Coherence function not normalizing properly? I am having the same problem with magnitude scaling.

            Currently, I am using matlab’s pwelch function to estimate the PSD, and then scaling the pwelch output to match the PSD estimate from SSCA at alpha=0.

            I am sure this is incorrect, but this method of finding coherence using pwelch allows me to better detect spectral peaks much better than SSCA alone.

            I have also tried estimating PSD using a frequency-smoothed periodogram, using a rectangular smoothing window of different lengths. However, I still obtain PSD estimates with a different scale from the SSCA.

            This is the cyclic domain profile (maximum coherence across all alphas) of my test signal: http://imgur.com/gTUe8mY.

            May I have some insights as to how I can solve this problem? Thank you!

          • Thanks for the comment Calvin. I also let William know you left it so he might reply as well.

            I’m not sure if William had enough time to complete his coherence debugging to his satisfaction. Of course, my several implementations of the SSCA produce good coherences.

            I think a possible strategy for you is to make sure that your various PSD estimates (SSCA, FSM, pwelch) have the correct shape and scaling independent of each other. That is, process a signal with an exactly known (analytically) PSD, such as white Gaussian noise with unit power (variance). Make sure each PSD is close to one across frequency. Then move on to other signals with known PSDs, such as PSK/QAM with IID symbols.

            Also, the degree to which a spectral coherence estimate is accurate depends on the length of the processed data record and the spectral resolutions used in the SSCA and in the off-to-the-side PSD estimate, so there is some art to this task. Coherence can also go awry if the PSD contains discontinuities (e.g., step functions, impulses) or very rapidly fluctuating spectral regions.

          • What is the signal that you used to generate the estimates in your imgur link? I see that the maximum coherence for alpha=0 is about 2. Since the coherence for alpha=0 (non-conjugate) should be \hat{S}_x^0(f) / \hat{S}_x^0(f) = 1, you clearly have a factor-of-two difference between the psd-on-the-side estimate (denominator) and the SSCA output (numerator).

          • Calvin L says:

            Thanks for the input Chad! You gave me a couple of ideas as to how to start debugging the coherence. One possible area I’ll look into is a power correction for pwelch and my PSD estimates. Once I get something working again I’ll update this thread.

            As for my input signal, I’m using a srrc 16-QAM signal. Thanks for the help.

          • If you used a square-root raised-cosine 16-QAM signal with IID symbols to produce the cyclic domain profile shown in your imgur link above, we’ve got bigger problems than a PSD scaling factor. That signal has three non-conjugate cycle frequencies: -f_{sym}, 0, and +f_{sym}, and no conjugate cycle frequencies. Your plot shows at least 10 significant cycle frequencies. Seems closer to a rectangular-pulse QAM pattern.

          • Calvin L says:

            Hi Chad,

            First, regarding my previous CDP plot:

            Seems closer to a rectangular-pulse QAM pattern.

            You are right, I made a careless mistake in my signal generation. These are the corrected plots for a SRRC 16QAM signal. CDP: http://imgur.com/WetD7O9, Coherence: http://imgur.com/msOwflK.

            Next, regarding the problem with normalising, it seems that the best coherence plots come from FSM PSD estimates that strongly mirrors the SSCA PSD (Sx at alpha=0). This means that a SSCA with coarser spectral resolution (i.e. smaller df) requires a corresponding FSM with a larger smoothing window. A FSM that is not sufficiently smoothed will produce spurious peaks in the coherence plot along the alpha=0 axis. These peaks are the cause of the larger than 1 magnitudes at alpha=0 on the CDP. A PSD that is not sufficiently smoothed will produce peaks.
            Coherence: http://imgur.com/lwcp4vv
            CDP: http://imgur.com/YrL65pM

            Conversely, a PSD that is oversmoothed will cause dips below magnitude 1 along alpha=0 of the coherence plot. However, this does not affect the CDP plot.
            Coherence: http://imgur.com/3KGP0u7
            CDP: http://imgur.com/d1AUtuQ

            My conclusion from all these is that the FSM PSD should mirror the SSCA PSD as closely as possible. This makes me question the need for a separate PSD estimate. Why should I be using a fine PSD estimate if ultimately I am smoothing it to match the SSCA PSD, which has a much lower spectral resolution. Would it not be easier just to interpolate the SSCA PSD to match spectral resolution (block length of alpha) and use that in my coherence calculations?

            I would appreciate any help as I am quite puzzled by this seemingly contradictory situation. Thank you!

          • I view the interpolation of the SSCA-produced PSD values as a “separate PSD estimate.” That can definitely work. I would caution you to not make any solid conclusions, though, based on this one rather easy test case with a smooth underlying true PSD.

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