The CSP Blog recently received a comment from a signal processor that needed a small amount of debugging help with their python spectral correlation estimator code.
The code uses a form of the time-smoothing method and aims to compute and plot the spectral correlation estimate as well as the corresponding coherence estimate. What is cool about this code is that it is clear, well-organized, on github, and is written using Jupyter Notebook. Moreover, there is a Google Colab function so that anyone can run the code from a chrome browser and see the results, even a python newbie like me. Tres moderne.
The researcher is Fabio Casagrande Hirono, and he is happy to share his code with the world. We found a little bug relating to the calculation of the denominator of the coherence function and after that things worked well.
The plotting of the obtained spectral correlation and coherence estimates is a bit different than my style. I use MATLAB’s waterfall.m, but Fabio uses matplotlib, which is likely to be of high interest to a lot of CSP Blog readers. Here is a sample:
Cool colors! And unlike my 3D plots, the hidden-line effect is not completely opaque–you can see the PSD curve lurking behind the bit-rate feature.
So go check it out. But remember that if you really want to learn CSP, or any kind of signal processing, writing your own code is indispensable. But now you have another tool in the toolkit to help you along the way.
Thanks Fabio!
Chad,
I investigated Fabio’s Python code for analyzing cyclostationary data you posted, cycloStationaryAnalFabio.p.
The code works nicely for when the symbol or bit rate matches any of the cyclic frequency vector elements used in the search. That is, he uses T_bit = 10 with alpha_vec = [0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ] for example. However, when the symbol rate is not matching any of the elements in the cyclic frequency vector (say T_bit = 7), I cannot see any cyclic frequency. When the cyclic search reaches 0.5% resolution, then the cyclic frequency barely begins to be displayed. This resolution is impractical.
Was wondering if you have a suggestion to rectify this problem since you have already reviewed the code. I also found a minor bug in his 3D plotting routine which I sent to his github.
Thanks.
Fabio’s code is a variation on what I call the time-smoothing method of spectral correlation estimation. As such, it is defined as an estimator of the spectral correlation function for a single value of cycle frequency. If that value is exactly equal to a cycle frequency exhibited by the input signal, you’ll get a good estimate. If that cycle-frequency value is within the basic cycle-frequency resolution of a cycle frequency exhibited by the input signal, you’ll get a reasonable estimate. What is that basic cycle-frequency resolution? About the reciprocal of the length of the processed data record.
You can, if you choose, modify Fabio’s selection of a handful of cycle frequencies for which to estimate the spectral correlation function. If you modify it to the set
, and your input signal has 
samples per bit, you’ll get good estimates. If you modify it to the set 
and you retain the original signal with 
samples per bit, your visited cycle frequencies will be many cycle-frequency resolution widths away from the true ones (because the data-record length here is 
), and you should not expect to get good estimates of the spectral correlation for 
. You can continue to increase the set of visited cycle frequencies until they are spaced by the cycle-frequency resolution, or smaller, and then you’ll see the spectral correlation function for the signal’s cycle frequencies 
in a plot.
However, using the time-smoothing method here, or the one in my TSM post, as a method of estimating the spectral correlation function over all cycle frequencies (modulo the cycle-frequency resolution) and spectral frequencies is not advisable. This is because the computational cost of that kind of “TSM in a loop over cycle frequency” is orders of magnitude greater than that for better alternatives like the SSCA and FAM.
Does that help?
Yes and thanks.