How to calculate covariance using the wcoherence function

14 views (last 30 days)
Jakob Sievers
Jakob Sievers on 18 Feb 2019
QUESTION:
Hi everyone
Is it possible to calculate the covariance of two signals, not by the standard cov function, but by wcoherence instead? Below is a code example:
% make data
rng default;
t = 0:0.001:2;
x = cos(2*pi*10*t).*(t>=0.5 & t<1.1)+cos(2*pi*50*t).*(t>= 0.2 & t< 1.4)+0.25*randn(size(t));
y = sin(2*pi*10*t).*(t>=0.6 & t<1.2)+sin(2*pi*50*t).*(t>= 0.4 & t<1.6)+ 0.35*randn(size(t));
% calculate cross-spectrum
[~,wcs,f,coi] = wcoherence(x,y,1/0.001);
wcs=real(wcs);
% covariances
cc_cov=cov(x,y);cc_cov=cc_cov(2);
cc_wcoh=trapz(t,trapz(f,wcs)); %????!
% plot results
figure
h = pcolor(t,log10(f),wcs);
h.EdgeColor = 'none';
ax = gca;
hold on;
plot(ax,t,log10(coi),'k--','linewidth',2);
colorbar;
ax.XLabel.String='Time (s)';
ax.YLabel.String='Logarithmic frequency';
ax.Title.String = {'Wavelet cross spectrum';['cov(x,y)=',num2str(cc_cov),' | cov_{wavelet}=',num2str(cc_wcoh)]};
Now, obviously my understanding of how to calculate the covariance from the wcoherence output (cov_wav) is flawed.
Can anyone help me get this right?
Essentially what I am trying to do is to calculate the cross spectrum of two signals and determine which regions in the resulting scalogram to include in an estimate of a covariance, as opposed to including everything (as is the case with the cov function).
Thanks in advance!
Cheers
Jakob
  12 Comments
Bjorn Gustavsson
Bjorn Gustavsson on 7 Mar 2019
The main objection is that it would restrict "how loud someone can shout for help" and there is something repugnant with that. If you aren't allowed to shout "a kingdom for a horse" we will never know that you're 5 minutes from dying from dehydration...
That's why I don't like the idea of a restriction, it has not till now been difficult to judge whether a help-me-Ill-pay-plea comes from an undergraduate that's been too lazy to learn or from someone with a real (more interesting) problem. It seems as the community culture of persistent tutting at those that do have been good enough to keep this at a very small fraction of the questions. As long as there is a restriction againg "selling" there will be no market, and I think that is enough to keep this at an unproblematic level. I see myself as a regular contributor, even though I no longer participate as actively as during the news-group-days, and though I dont want this forum to convert to a buying-selling consultation afair I certainly am no part of a consensus for a ban.
Finally, you're completely correct in that I only saw this question and not the previous ones - and I realised that I actually knew something about the topic.

Sign in to comment.

Answers (1)

Bjorn Gustavsson
Bjorn Gustavsson on 19 Feb 2019
Well, to me it looks as if the wcoherence function returns a normalized cross-spectra, i.e. something like
WC = FT(A).*FT(B)/sqrd(|FT(A)|*|FT(B)|)
Since the covariance are not normalized like that (the correlation-matrix is), you cant go from WC - COV. Perhaps you can select periods to include in covariance clculations from time-periods with interesting coherence.
HTH
  8 Comments
Bjorn Gustavsson
Bjorn Gustavsson on 5 Mar 2019
Well dont jump to that conclusion, the power-spectra is the Fourier-transform of the autocorrelation function, likewise for the cross-covariance and cross-spectra:
x = 0:100;
K = exp(-(x-20).^2/5^2) + exp(-(x-70).^2/12^2);
a = randn(5000,1);
b = randn(5000,1);
A = conv2(a,K','full');
B = conv2(b,K','full');
clf
subplot(2,1,1)
plot(fftshift(ifft((fft(A).*conj(fft(A)))))/numel(A),'r')
hold on
plot(xcorr(A,numel(A)/2,'unbiased'))
subplot(2,1,2)
xcAB = xcorr([A,B],numel(A)/2,'unbiased');
plot(fftshift(ifft((fft(A).*conj(fft(B)))))/numel(A),'r')
hold on
plot(xcAB(:,2))
The "only" difference between this and what you get with the spectrogram/wavelet spectra is that those methods cuts the data into shorter (possibly overlapping segments) and applies the same functions to those data-segements. Thus making it possible to see temporal variations in spectral content/covariance.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!