How to calculate covariance using the wcoherence function

3 views (last 30 days)
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
Guillaume
Guillaume on 6 Mar 2019
@Jakob,
Probably, your question is about a topic that few of us regular contributors know. And possibly, the very few that know about that topic didn't see your question. I doubt the addition of a monetary reward changed anything to that and you were just lucky that Bjorn saw it that time.
This particular question attracted more attention than usual because the monetary reward triggered a discussion off-Answers. I don't think there was any ruffled feathers. Just a bit more clucking than usual.
@Bjorn, Jan,
I think that the consensus of the regular contributors is that we don't want Answers to be a market place. And I believe that Mathworks don't want that (you can pay them for support). However, Mathworks do listen to our input, so Bjorn, I'd be interested to know why you don't like that restriction.
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
Jakob Sievers
Jakob Sievers on 4 Mar 2019
Yes that seems to be a lot better. Though I still come to the conclusion that there seems to be no clear/straightforward connection between spectrograms/wavelets and the covariance of the original dataset. Thats unfortunate but oh well.
In any case you have helped me a lot and so you obviously are rewarded with the 100$ :-)
Please send me a private message with your bitcoin wallet so that I may transfer the money to you.
Thanks again!
Jakob
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!