# The Spectral Coherence Function

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.

Let’s start with reviewing the standard correlation coefficient $\rho$ defined for two random variables $X$ and $Y$ as

$\rho = \displaystyle \frac{E[(X - m_X)(Y - m_Y)]}{\sigma_X \sigma_Y}, \hfill (1)$

where $m_X$ and $m_Y$ are the mean values of $X$ and $Y$, and $\sigma_X$ and $\sigma_Y$ are the standard deviations of $X$ and $Y$. That is,

$m_X = E[X] \hfill (2)$

$m_Y = E[Y] \hfill (3)$

$\sigma_X^2 = E[(X-m_X)^2] \hfill (4)$

$\sigma_Y^2 = E[(Y-m_Y)^2] \hfill (5)$

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

Now consider the spectral correlation function,

$S_x^\alpha(f) = \displaystyle \lim_{T\rightarrow\infty} \frac{1}{T} E[ X_T(t, f+\alpha/2) X_T^* (t, f-\alpha/2)], \hfill (6)$

where the expectation operator $E[\cdot]$ can be either a time average or the ensemble average used in a stochastic-process formulation. For each finite $T$, we can consider the quantity

$\displaystyle \frac{1}{T} E[X_T(t, f+\alpha/2) X_T^* (t, f-\alpha/2)], \hfill (7)$

which is the correlation between the two random variables $\displaystyle \frac{1}{T^{1/2}} X_T(t, f+\alpha/2)$ and $\displaystyle \frac{1}{T^{1/2}} X_T(t, f-\alpha/2)$. To form the correlation coefficient, we need the variances for these two random variables,

$\sigma_1^2(T) = E[\displaystyle\frac{1}{T} \left| X_T(t, f+\alpha/2) \right|^2] \hfill (8)$

and

$\sigma_2^2(T) = E[\displaystyle\frac{1}{T} \left| X_T(t, f-\alpha/2) \right|^2]. \hfill (9)$

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 $f\pm \alpha/2$.

For finite $T$, then, we can form the correlation coefficient

$\rho(T) = \displaystyle \frac{\frac{1}{T} E[X_T(t, f+\alpha/2) X_T^* (t, f-\alpha/2)]}{[\sigma_1^2(T) \sigma_2^2(T)]^{1/2}}. \hfill (10)$

Now, as $T\rightarrow \infty$, the numerator of $\rho(T)$ converges to the spectral correlation function $S_x^\alpha(f)$ and

$\displaystyle \lim_{T\rightarrow\infty} \sigma_1^2(T) = S_x^0(f+\alpha/2) \hfill (11)$

and

$\displaystyle \lim_{T\rightarrow\infty} \sigma_2^2(T) = S_x^0(f-\alpha/2). \hfill (12)$

If the limit of the quotient exists, then the correlation coefficient 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 (13)$

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

$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 (14)$

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 $65536$ samples and the frequency resolution (width of $g(f)$ in the FSM) is set to $0.01$ (one percent of the sampling rate, which here is $1.0$) 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 $S_x^0(f\pm\alpha/2)$. If the PSDs are not well resolved, or contain zeros, the coherence quotient can be numerically unstable or erroneous.

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.

## 9 thoughts on “The Spectral Coherence Function”

1. Chen says:

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.

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

2. Clint says:

If the PSDs are not well resolved, or contain zeros, the coherence quotient can be numerically unstable or erroneous.

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

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

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.

how do I decide if they’re too short or otherwise might lead to low-confidence analysis results?

For a noiseless random signal with cycle frequency $\alpha$, you’ll want to process a data block with length that is many multiples of $1/\alpha$ 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.

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 $f_0 >B-\frac{1}{2alpha}$ in $\alpha S^\alpha(f_0)$
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 $\hat{\alpha} > B$ ,since the $S^0(f) = 0$ when $f>B$

1. It’s the first time I use latex comment. The first equation should be “Especially for $|f_0| >B-\frac{\alpha}{2}$ in $S^\alpha(f_0)$

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

4. Clint says:

/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 $f\pm \alpha/2$

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!

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