# How to calculate covariance using the wcoherence function

13 views (last 30 days)
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).
Cheers
Jakob

Jakob Sievers on 6 Mar 2019
I suppose I might add that I did in fact post this question without any monetary reward first. It only attracted _any_ attention the moment I indicated I might give a reward for a satisfying answer. I did this explicitly stating that there appeared to be a disconnect between the broader relevance of this question and how important it was to me (though that post appears to have been taken down by admin now). I do not intend for the mathworks forum to become a place where monetary rewards is the norm because it will indeed degrade the experience. My point is that I myself was somewhat conflicted at first and as I also mentioned: I wish such a website DID exist so I wouldnt have to cause a stir here.
Secondly I should also add that I had in fact published a similar question a while back and had not received any responses back then either, suggesting that perhaps I was right in my feeling that this question was simply not relevant enough to others.
Either way I apologize for any feathers I may have ruffled.
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 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.

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

Bjorn Gustavsson on 26 Feb 2019
Your choise of frequencies seems a bit odd to me. If you just take the inverse of the first few frequencies you get the corresponding period-times of the low-frequency Fourier-component:
2400 2375.9 2352 2328.4 2305
That is very long period-times compared to the time of your window of 200 samples (10 s). Your entire measurement time is 2400 s, if I got it rights. If you want to separate variations with such periods you'll have to have data from a longer observation-period. Would this be something that gives good enough cross-spectra?
Fin = linspace(0,10,400);
[SA,F,T,PA] = spectrogram(A,400,200,Fin,20);
[SB,F,T,PB] = spectrogram(B,400,200,Fin,20);
XAB = real(SA.*conj(SB));
subplot(3,1,1)
colorbar
subplot(3,1,2)
colorbar
subplot(3,1,3)
colorbar
HTH
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 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.