How to calculate covariance using the wcoherence function

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

[Jan]

I dislike offering money in the forum and I'm happy, that this does not belong to the usual behavior here.

[Jakob Sievers]

Agreed. yet I dislike even more being stuck with a fairly case-specific problem and not receiving any help or knowing where else to turn to. In reality I would have loved to be able to post this on a website which caters specifically to the subset of people who are willing to pay for help with small tasks like this. I didnt know of such a matlab-centric site and so chose to take a chance here.

[John D'Errico]

This is not a web site where services/solutions are offered for pay.

[Jakob Sievers]

Before adding a monetary reward I took care to read the Comunity Guidelines which do not state that this is not allowed.

[John D'Errico]

Then I'll see about getting that added to the guidelines in the next meeting.

[Jakob Sievers]

Good. I didn't intend to offend anyone.

[Guillaume]

Mathworks probably needs to tidy up their guidelines page. I wasn't aware of the above community guideline page, which doesn't really match the Matlab Central Terms of Use which you can find in the footer of every page.

The terms of use are clear, under 2.iv:

"Content that you submit must be offered free of charge. You may not use the Site to sell or market your products or services to others."

[Bjorn Gustavsson]

Just to add my 5 p to this discussion of principles: Neither the OP nor I have broken that term. I have offered the advice (for what they're worth) free of charge, he asked a question where he indicated he might pay for replies he would/might find useful. If Mathworks deem it preferable to prohibit posters to throw their money away using their forum that's their absolute right - but something about such a restriction rubs me the wrong way.

The value I've extracted from the matlab-news-group and file-exchange the last 20 odd years by far outweigh what any sensible person would consider paying for my advice, and I would never be cheap or vain enough to peddle them for money, nor do I care too much about the community brownie-points either (as an illustration: John d'Errico was a matlab-guru already 20 years ago and no amount of additional brownie-points will change that).

[Jan]

Beyond any guidelines, we have Kant's idea:

Act only according to that maxim whereby you can, at the same time, will that it should become a universal law.

Then we have two aspects:

1. Someone offers money for a solution. If this is usual in the forum, say 90% of the questoins will cause a transer of money, the main idea of the forum would degrade soon: The voluntary mutual assistance and the interaction as a community. We would get advertising, envy and competition. The forum would not be a list of the best solutions anymore, but of the most aggressively advertised questions and answers.
2. Someone offers money for a solution as a very rare exception. This will not influence the nature and efficiency of the forum. Maybe it has a positive effect on this single thread.

I've voted for the question, because I really like the discussion about offering money in the forum. I apperciate, that this will be a rare exception. In my opinion the OP Jakob did not offend anybody and did not break any rules or guidelines.

Some or many of the frequent contributors earn money with professional programming. It is very useful and likable, that this forum is not a colletions of cries to "buy me!", but that commercial interests are not discussed in public here.

[Jakob Sievers]

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]

@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]

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.

Answer 1

[Bjorn Gustavsson]

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

[Jakob Sievers]

So, more broadly, it is not possible to calculate a cross-spectrum or scalogram, using wavelets, that is directly related to covariance? I mean, does a function like that not exist in Matlab? Maybe my mistake was just to look at the wcoherence function?

cheersJakob

[Bjorn Gustavsson]

Well why wouldnt you calculate a cross-spectrum like:

fta = fft(A);
ftb = fft(B);
Xspec = fta*conj(ftb)

Or with help of the spectrogram function:

[SA,F,T] = spectrogram(A,WINDOW,NOVERLAP,F,Fs);
[SB,F,T] = spectrogram(B,WINDOW,NOVERLAP,F,Fs);
XAB = SA.*conj(SB);

That should give you something to calculate time-variations for covariances...

[Jakob Sievers]

Thanks for your reply. I still can't quite seem to get this to work. When I try the fft solution I receive a memory warning (my datasets, A and B, are 48000 long, measured at 20hz, and presumably 48000^2 values is too big for my system) and when I try the spectrogram it shows the scales linearly from 0 to 10 (naturally) but since the majority of the covariance exists on scales lower than 0.5hz the resulting plot is mostly zeros. There doesnt appear to be an option to assign a range of frequencies of interest in the spectrogram function, but I might be wrong?

I suspect it might be beneficial if I take a step back and explain what I am trying to do in a bit more detail:

I work with turbulence and stumbled across a paper by Hudgins et al., 1993 titled "Wavelet transforms and atmospheric turbulence". In it the authors show this neat plot:

which shows the wavelet cross scalogram and the actual covariances at a given time above (presumably derived using e.g. the trapz function along the individual columns). I am trying to re-create this with my own dataset. I have attached an example dataset ("AB.mat") with 48000x2 entries reflecting 40 minutes of 20hz observations in case this helps to illustrate my situation.

The paper clearly describes how they arrive at this result, yet I have been unable to translate it into a working matlab example. Undoubtedly because I am new to wavelets. In case you do not have access to the paper through official channels I have added a photo of the central part of their brief section on wavelet treatment here:

I really hope you can help me clarify what I am doing wrong.

Sincerely

Jakob Sievers

[Bjorn Gustavsson]

OK, I think I made a typo for the cross-spectrum over the entire data-set. The moultiplication should be element-wise:

fta = fft(A);
ftb = fft(B);
Xspec = fta.*conj(ftb);

That said, you certainly can tell the spectrogram function what frequencies to use:

[S,F,T] = spectrogram(X,WINDOW,NOVERLAP,F,Fs);

There you can give it an array for F - granted you have to use a long enough window to capture variations at your frequencies of interest. If you want something like 0:(1/4):20, you should have a window over at least ~4 s. If you try spectrogram first for your data then you should get a feel for how a standard short-time FT works - if that doesn't give you enough resolution at higher frequencies then you proceed to wavelets or multi-taper versions. I'm "very pedestrian" so I have to aproach this on a step-by-step path.

HTH

[Jakob Sievers]

Hi again

The first solution just yields a single vector with the cospectrum. I am familiar with this method as that is what I already use. My goal is to expand into a 2D "space" where the second axis is time. Perhaps you meant to transpose one of the vectors? I.e.:

fta = fft(A);
ftb = fft(B);
XAB = real(fta.*conj(ftb'));

Though that doesnt seem to work either. At least not in the sense that if I perform a trapz along the columns I dont get results which sum up to the regular covariance of A and B. I.e.:

cc_cov=cov(A,B);cc_cov=cc_cov(2);
% and
cc_trapz=sum(trapz(XAB));
% do not produce the same result

As for the spectrogram solution I am still having trouble getting anything meaningful out of it. Perhaps this is because I am not setting my window the right way. This is an example of where I am at:

You'll notice that the final covariance is again not even close. Though I am beginning to wonder if it is even possible that there should be a connection between a spectrogram/wavelet illustration and the "raw" covariance of two signals. I am beginning to think that spectrograms and wavelet transforms are a way to bring out certain features or aspects of a dataset but that these results arent necessarily linked to the raw covariance. Perhaps that is what I am getting wrong here?

In any case, the plot was achieved with the script below. Looking forward to hearing your opinion on all of this.

cc_cov=cov(A,B);cc_cov=cc_cov(2);
    freq=10.^linspace(log10(hz/(length(A))),log10(hz/2),1000);
    [SA,F,T,PA] = spectrogram(A,200,[],freq,20); 
    [SB,F,T,PB] = spectrogram(B,200,[],freq,20); 
    XAB = real(SA.*conj(SB)); 
    figure
    subplot(3,1,2:3)
    imagesc(T,log10(F),XAB)
    ylabel('frequency')
    xlabel('time')
    ax=axis;
    subplot(3,1,1)
    grid on
    hold on
    Z = trapz(F,XAB,1);
    cc_spec=sum(Z);
    plot(T,Z)
    set(gca,'xlim',ax(1:2))
    plot(ax(1:2),[0 0],'-k')
    ylabel('trapz(F,XAB)')
    gc=gca;
    gc.XLabel.String='Time (s)';
    gc.Title.String = {'Wavelet cross spectrum';['cov(x,y)=',num2str(cc_cov),'  |  cov_{spectrogram}=',num2str(cc_spec)]};

[Bjorn Gustavsson]

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)
pcolor(T,F,asinh(XAB)),shading flat
colorbar
subplot(3,1,2)
pcolor(T,F(5:end),XAB(5:end,:)),shading flat
colorbar
subplot(3,1,3)
pcolor(T,F(1:15),XAB(1:15,:)),shading flat  
colorbar

HTH

[Jakob Sievers]

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]

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.