In this post I introduce the spectral coherence function, or just coherence. It deserves its own post because the coherence is a useful detection statistic for blindly determining significant cycle frequencies of arbitrary data records. See the posts on the strip spectral correlation analyzer and the FFT accumulation method for examples.

Let’s start with reviewing the standard correlation coefficient defined for two random variables and as

where and are the mean values of and , and and are the standard deviations of and . That is,

So the correlation coefficient is the covariance between and divided by the geometric mean of the variances of and .

Now consider the spectral correlation function,

where the expectation operator can be either a time average or the ensemble average used in a stochastic-process formulation. For each finite , we can consider the quantity

which is the correlation between the two random variables and . To form the correlation coefficient, we need the variances for these two random variables,

and

Here we are assuming that the two variables have zero means, which is always true when there is no additive sine-wave component in the data with frequency .

For finite , then, we can form the correlation coefficient

Now, as , the numerator of converges to the spectral correlation function and

and

If the limit of the quotient exists, then the correlation coefficient is given by

We call this function the *spectral coherence*. A similar argument holds for the *conjugate spectral coherence*,

which is a normalization of the conjugate spectral correlation function. Since the coherence is a valid correlation coefficient, its magnitude will be less than or equal to one, and since the involved random variables are complex-valued, in general, so is the spectral correlation function and the coherence. Therefore, the coherence lies in the closed unit disk in the complex plane.

Here is an FSM-based estimate of the coherence for our rectangular-pulse BPSK signal:

The data-block length is samples and the frequency resolution (width of in the FSM) is set to (one percent of the sampling rate, which here is ) for both the numerator spectral correlation function and the denominator PSDs.

In practice, of course, the numerator and denominator of the coherence are estimates corresponding to finite-duration data blocks. The tricky part of estimating the coherence involves good selection of the estimator for the denominator PSDs . If the PSDs are not well resolved, or contain zeros, the coherence quotient can be numerically unstable or erroneous.

I cannot help asking you how to estiamte the denominator PSDs? I supposed that the FAM or SSCA based SCF estimates are obtained. Then I think FFT algorithm is a must to avoid computational complexity when estimating PSDs. I have tried a lot but I cannot present a detail solution to estimate PSDs furthermore.

So could you please give some advice?

Did you look at the discussion of coherence and PSD estimates in the FAM post? Also be sure to look through the comments of that post and the SSCA post; there is relevant material there too.

Hi again Chad,

I’m just wondering what your thoughts are on how to automatically decide whether a PSD is “well resolved” or not. I’m guessing that it depends on the application requirements, but do you have any rules of thumb?

If I’m analyzing recorded signals of various lengths, how do I decide if they’re too short or otherwise might lead to low-confidence analysis results? On the flip side, if I’m recording a signal to be analyzed, how do I decide between a long recording that uses all my storage space (just to be safe) or a shorter recording that would be sufficient?

Thanks for any thoughts!

Clint

I don’t think it will be possible to automatically decide whether a PSD estimate is well-resolved. You’d have to know the true PSD of the signals in the processed data in advance. But there are spectrum features that are inherently unresolvable, such as unit-step functions and additive sine-wave components. They possess features that have spectral widths of zero, and so any spectrum estimate that possesses a finite effective spectral resolution cannot resolve them. For a step-function-like spectrum (such as a square-root raised-cosine QAM signal with zero roll-off), the edges of the spectrum cannot be perfectly estimated. For a sine-wave component, the spectrum for the frequencies near the sine-wave frequency will inevitably reflect the tone’s energy, not just the noise or other continuous-spectrum signal components that lie there.

You might consider building into a power spectrum estimator the ability to detect spectral edges or corners or the like, perhaps using a derivative or similar. Derivatives are tricky when dealing with noisy or random signals though. If you could detect such occurrences in your spectrum, you could then discount high coherence estimates that involve the spectral components near the detected edges/corners.

For a noiseless random signal with cycle frequency , you’ll want to process a data block with length that is many multiples of to ensure you can see that periodicity. For a noisy signal, though, you also need to process enough data to overcome the effects of noise on the spectral correlation estimate. So there is no easy answer. I capture many wideband signals and RF scenes, and am supplied with such by others as well, and I have plenty of storage, since it is cheap and readily available. So I generally try to store long scenes, where the length in seconds greatly exceeds the reciprocal of the smallest non-zero cycle frequency of interest. That turns out to be a small number of seconds, typically.

Hi Chad!

I’m trying to implement the spectral coherence for the task that estimate the symbol rate for low roll-off factor signal.

If the pulse shape is bandlimit(e.g. SRRC with roll-off=0 and bandwith B). There is 0/0 problem when computing the coherence function. Especially for in

What do you do to avoid 0/0? For inband singularity point, I just filter out this NaN points to 0. For outband zeros,I padded with small value(1e-10). But this will amplify candidate cycle frequency ,since the when

It’s the first time I use latex comment. The first equation should be “Especially for in

The coherence is undefined whenever the denominator is zero, which can happen only when you process a noiseless signal. Are you sure you want to process noiseless signals?

/quote Here we are assuming that the two variables have zero means, which is true when there is no additive sine-wave component in the data with frequency

What are your thoughts for when you do have additive sine wave components (e.g. pilot tones)? Of course they can make it easier to detect the signal, but it seems that they would confound the detection of cycle frequencies.

I seem to remember reading somewhere that it may be desirable to remove the additive sinusoid components from the PSD (by filtering the computed PSD?), but I’m not seeing it now.

Thanks!

I don’t think it is a good idea to process data containing finite-strength additive sine-wave components with a CSP estimator such as the SSCA, FAM, FSM, etc. Detecting them and removing them will help avoid false-alarms due to the unresolvability of the spectrum when the signal contains sine-wave components. There is no advantage to such CSP processing anyway–if you can detect the sine-wave, you can estimate its parameters, and nothing that the SCF or cyclic cumulants will tell you about the sine-wave is more information than that. I probably said as much, as you suggest, in other comments or posts, but I can’t find one at the moment.

You’ll need to do more than process the PSD. You should try to detect the sine-wave components and remove them from the data. There are many ways to detect sine waves, but a good non-parametric way is to just use the Fourier transform and some kind of thresholding. But that’s not really CSP, so I don’t talk about it much on the Blog. It would be a reasonable topic for a Signal Processing ToolKit post, though.

Thanks for sharing your CSP knowledge, Chad. It seems I’m overlooking an important detail or two in this topic. First, must expected values be taken at any point?

Second, in eqn 13 & 14 denominators, shifting by $\alpha/2$ and zero filling guarantees div by 0. Currently, I do not take expected values and simply divide the bpsk signal’s SCF by shifted PSD’s and replace div by 0 in result with 0.

My non-conj spectral coherance graph looks like mostly like yours but dosn’t roll off nearly as much as alpha approaches 1. At $\alpha=1$, your non-conj spectral coherance value seems about 0.4. Mine is about 0.8. (I don’t see a way to attach pictures or would do that.) Thanks for any pointers to my misunderstandings!

[Added by Chad]

Thanks for checking out the CSP Blog, Mike, and the comments!

No. Expectations can only be taken when there is an ensemble. This is almost always a fictitious mathematical ensemble. In rare cases one might have an actual ensemble at hand. The CSP Blog is about signal processing, and so any averaging we do in the context of estimation is time (e.g., TSM) or frequency (e.g., FSM) averaging. Here we are concerned with estimating the spectral correlation function, coherence function, cyclic autocorrelation function, cyclic cumulants, etc., using only a single sampled-data signal (“time-series”).

Yes, indeed. In DSP, we only have the Fourier transform over the normalized frequency interval , and so we just have to live with the fact that has increasingly small overlap as increases, and similarly for a PSD estimate. And the dividing-by-zero problem is even worse than that, for it is perfectly possible for a PSD to have zeros. So you must avoid outputting a non-zero coherence estimate for any pair for which the coherence denominator is zero.

That explains the appearance of my coherence plot in the post. As the cycle frequency increases in magnitude, the overlap between the two shifted transforms decreases. In continuous-time mathematical analysis, this doesn’t happen. The coherence is perfectly well defined for all and , except for the case where the denominator PSDs have a zero.

As for why your coherence is not decaying as rapidly as mine, what are the signal and noise power levels?

Thanks, Chad, for your helpful reply, and seems I’m in better shape than I thought! I converted your BPSK sample code to python as a starting point. Barring bugs (!) I create the BPSK signal of random 1/-1 (duplicated to samples/sym), mix to carrier, then add noise to achieve 15 dB s/n.

And I now see if I drop to 10 dB s/n it falls much more inline with yours! I’ll get a grasp of the basics yet. 🙂

As a short follow up to my questions, for my results described in the previous note I used the SCF FSM subroutine that I used to generate correct looking results for previous blog lessons.

Dear Chad,

I am trying to use the SCohF to blindly detect the symbol rate of a set of signals.

I am using here the code from J. Antoni — I do have my code for the symbol rate estimation and I have read your critics on the algorithm, but I wanted to be sure that what I observe is not due to a bug in my code (yeah, pretty low self confidence ;)).

The window size for the STFT is 16 -> frequency resolution = 18.75 kHz in the example below.

For clean signals it seems to work pretty well, take for instance the SCohF for a 4-FSK signal (20ksym/s, sampling rate = 200 ksps, generated with a signal generator, SNR = 60 dB):

https://ibb.co/1R11CXr

In the picture I show, for each alpha value on the X-axis, the maximum of the mean of the absolute values returned by Antoni’s code (from the ‘S’ matrix in his code).

The peak is crystal clear, and gives the correct result.

Given the properties of CSP techniques, I would expect the SCohF to be pretty insensitive to AWGN, therefore I’ve taken the signal and lowered the SNR to 10 dB (using Matlab’s awgn() function, so there should be no issues there). In some of your posts you’ve worked with far, far worse SNRs, e.g., I recall a ROC curve where you do signal detection at -11 dB.

However, running again the code with the same parameters on the corrupted signal, this is what I obtain:

https://ibb.co/xLD9n8w

All the peaks are gone… 😦

Am I missing anything?

Thanks in advance for your help!

Best,

Rob

The dominant peak is at 20 kHz, which is good, but the maximum coherence is tiny. This implies that either (1) the SNR is not 60 dB, (2) something is going wrong in your estimator, or (3) the particular kind of FSK signal you are using has very weak cyclostationarity. There are many varieties of FSK, and their statistical properties depend on how the symbol-clock phase and carrier phase are related, and vary, over time. I suggest you use a simpler, more well-understood signal, such as my usual rectangular-pulse BPSK signal (and set its power to unity [0 dB]). That way you have lots of solid reference points to compare to what comes out of your code.

Since the peaks have completely disappeared, I would check this step. What was the exact command you used? Did you look at a PSD?

Dear Chad,

Thanks for the quick reply!

(1) I haven’t recorded the data myself, I have taken this value from the metadata that came with the signal dataset.

(2) the estimator is the one from Antoni’s paper. I simply take his code, load my signal, and run it as:

[S, alpha, ~, ~] = Fast_SC(x, Nw, alpha_max, Fs, opt);

where:

– x = my signal

– Nw = 16

– alpha_max = 100e6

– Fs = 200e6

– opt.coh = 1 (to compute the SCohF)

Then I plot mean(abs(S(:, 2 : end)))

My code does the same, but in a non-optimized way and on selected alpha values (to reduce the computational load), and the result is virtually the same (in the AWGN case, I get garbage out).

(3) this might be the answer. In my dataset I do have multiple xFSK signals, and everything’s nice for all of them except this one… (I do few other mistakes, but these are due to the fact that I pick the wrong peak — either an harmonic or a sub-harmonic –, thus it is a mistake in my detection algorithm and not an issue with the CSP output)

(4) the signal actually disappeared from the PSD…

https://ibb.co/bbXjB5B

which is weird. The SNR reduction is simply made by invoking the Matlab function:

x = awgn(x, 10, ‘measured’);

Thanks again for your help and have a nice day!

Rob