Coherence MATLAB Resources

March 27, 2019

MATLAB® offers operators supporting several statistical estimators using Welch’s method. These include the following:

   

mscohere(…)

Estimate ordinary coherence.

https://www.mathworks.com/help/signal/ref/mscohere.html

   

pwelch(…)

Estimate power spectral density.

https://www.mathworks.com/help/signal/ref/pwelch.html

   

cpsd(…)

Estimate cross power spectral density.

https://www.mathworks.com/help/signal/ref/cpsd.html

Example Code from MATLAB

Complex coherence can be calculated in MATLAB using CSD and PSD as demonstrated below.

clear                                      % clear all variables
format long g;                     % set format of text output
N = 16;                                  % length of input data sequence x
Fs = 1;                                  % sample rate (used for scaling)
percent = 50;                      % percent overlap
Nfft = N * (1-percent/100);
x= [0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0]’;
y= [0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0]’;
%   0        3          7              11  13  15
%—————————————————————————-
% smoothing window
%—————————————————————————-
window = hanning(Nfft);     % Hanning       window
%—————————————————————————-
% calculate estimates
%—————————————————————————-
[Sxy,freq]=cpsd(               x,y,window,Noverlap,Nfft,Fs);   % estimate Sxy
[Sxx,freq]=pwelch(          x,    window,Noverlap,Nfft,Fs);   % estimate Sxx
[Syy,freq]=pwelch(             y,window,Noverlap,Nfft,Fs);    % estimate Syy
[Cms,freq]=mscohere(   x,y,window,Noverlap,Nfft,Fs);    % estimate ordinary coherence
% |      |         |                     |  |  |                |              |       |_____sample rate
% |      |         |                     |  |  |                |              |________estimate length
% |      |         |                     |  |  |                |______________num. overlapping positions
% |      |         |                     |  |  |____________________ window type
% |      |         |                     | |______________________data sequence y
% |      |         |                     |_______________________data sequence x
% |      |         |_______________________________ estimation operation
% |     |____________________________________x-axis sequence
% |______________________________________ estimate sequence
Cxy = Sxy ./ sqrt( Sxx .* Syy );                                                 % estimate complex coherence
%—————————————————————————-
% perform integrity check
%—————————————————————————-
if norm(Cms-abs(Cxy).^2)<1e-12
fprintf(“Cxy integrity check OK!\n”);
else
fprintf(“Cxy integrity check FAILED!!!\n”);
end

A Note of Caution

As the name implies, estimators yield estimates. In general, these estimates contain some error with respect to the true value.

It would be difficult to find a more striking example of this situation than the K=1 Welch estimate of coherence. To understand this, compare the rectangular window K=1 Welch estimate of coherence with the true coherence of uncorrelated white noise signals:

If:
a.   x[n] and y[n] are zero-mean, white noise sequences
b.  x[n] and y[n] are uncorrelated
c.  the estimate Ĉxy (w) is the K=1 Welch estimate of x and y

Then:

(1)   \begin{equation*} \hat{C}_{xy}(\omega)=1 \end{equation*}

Proof:

(2)   \begin{equation*} |C(\omega)|=\left|\frac{S_{xy}}{\sqrt{S_{xx}S_{yy}}}\right| \end{equation*}

(3)   \begin{equation*} =\left|\frac{(\tilde{F}x)(\tilde{F}y)^*}{\sqrt{|\tilde{F}x|^2|\tilde{F}y|^2}}\right| \end{equation*}

(4)   \begin{equation*} =1 \end{equation*}

Thus, the K=1 Welch estimate of coherence is always 1!  This is definitely an estimate; but the question is, “Is it a good estimate?”

To answer this question, let’s compute the true (non-estimated) coherence:

(5)   \begin{equation*} |C_{xy}(\omega)|\triangleq|\tilde{F}[R_{xy}(m)]|\text{ by definition of coherence }C_{xy}(\omega) \end{equation*}

(6)   \begin{equation*} \triangleq|\tilde{F}\mathbf{E}[x(n)y(n+m)]|\text{ by definition of cross-correlation }R_{xy} \end{equation*}

(7)   \begin{equation*} =|\tilde{F}(\mathbf{E}[x(n)]\mathbf{E}[y(n+m)])|\text{ because x[n] and y[n] are uncorrelated} \end{equation*}

(8)   \begin{equation*} =|\tilde{F}(0\times0)|\text{because x[n] and y[n] have mean=0} \end{equation*}

(9)   \begin{equation*} =0\text{ because the Fourier transform of 0 is 0} \end{equation*}

So, we see that the true value of the coherence is Cxy(ω)= 0 (the smallest possible value), but the Welch estimate is Ĉxy(ω)= 1 (the largest possible value). It is likely that, in most cases, Welch will yield a very “good” estimate. However, as the above example demonstrates, this is not always the case and some caution is in order.

References

Bendat, J.S. & Piersol, A.G. (2011). Random data: Analysis and measurement procedures. (4th ed.). John Wiley and Sons.

Carter, G.C. & Arnold, C.R. (1971). Some practical considerations of coherence estimation. (Report No. ADA954033). Naval Underwater Systems Center. http://www.dtic.mil/docs/citations/ADA954033.

Carter, G.C. (1972). Estimation of the magnitude-squared coherence function (spectrum). (Report No. AD0743945). Naval Underwater Systems Center. https://apps.dtic.mil/dtic/tr/fulltext/u2/743945.pdf.

Carter, G.C., Knapp, C.H., & Nuttall, A.H. (1973). Estimation of the magnitude-squared coherence function via overlapped fast Fourier transform processing. IEEE Transactions on Audio and Electroacoustics, 21(4), 337 – 344. https://ieeexplore.ieee.org/document/1162496.

Goldman, S. (1999). Vibration spectrum analysis: A practical approach (2nd ed.). (pp. 184–186). Industrial Press.

Greenhoe, D.J. (2019). A book concerning statistical signal processing.

Liang, Z. & Lee, G.C. (2015). Random vibration: Mechanical, structural, and earthquake engineering applications (1st ed.). (pp. 200-201, 363-365). CRC Press.

Otnes, R.K. & Enochson, L. (1978). Applied time series analysis: Basic techniques. (Vol. 1). John Wiley and Sons.

Shin, K. & Hammond, J. (2008). Fundamentals of signal processing for sound and vibration engineers (pp. 284–287, 349–350, 385–386). John Wiley and Sons.