Coherence MATLAB Resources
March 27, 2019
MATLAB® offers operators supporting several statistical estimators using Welch’s method. These include the following:
Estimate ordinary coherence.
Estimate power spectral density.
Estimate cross power spectral density.
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
fprintf(“Cxy integrity check OK!\n”);
fprintf(“Cxy integrity check FAILED!!!\n”);
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:
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
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:
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.
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.