Cross spectral power and coherence of 2 signals

[Fonsi]

Hello,

I am trying to get the Cross spectral power and coherence of two signals. I should be able to get this:

"The cross spectral power and coherence of these 2 signals were calculated over a 1024-sample (8.5-minute) window using the fast Fourier transform applied to the 3 overlapping 512-sample subwindows within the 1024 coherence window. In each subwindow, the data were first linearly detrended and windowed using the Hanning (cosine) function before calculation of the Fourier transform. The 1024 coherence window was then advanced by 256 samples (2.1 minutes), and the calculation repeated until the entire N-N interval/EDR series were analyzed.

We used the product of the coherence and the cross-spectral power to weight these 2 effects in order to quantify the degree of the coupling."

I should get this type of graph

The RR signal i detrended before then try to get the cross power with sampling freq = 2. And I put 2 ways to get it. 

coherence too and later the product of both. Im try to get the spectogram and later look at 3d like picture.

I attached the vi.

The coherence maybe is bad fixed.

Thanks for help, guys.

Regards.

[jeanlucpicard]

Hi,

I would do the following:

1. get 1024 samples of your signals

2. detrend your data the way you want

3. calculate the coherence and cross power data

4. calculate the time for this window (eg. the time of the midpoint of time-domain data)

5. calculate the product of coherence and cross-power data

6. move your set by the given number (256 values) and repeat 1-5 until you reach the end of the signal

7. build arrays of the product data and the time points from the iterations (1-5 should be in a for loop, and then you can use the indexing output tunnels)

8. you can build a 3d graph by adding the time array to the x-data, the frequency array t the y-data and the product as the z-data (you might have to transpose the z-array)

For the cross power and coherence I previously used the MathScript Node with crosspsd and coherence_ms functions. (Probably there is some VI-based solution available also, but this worked for me previously.) For these functions you can define a window function of your choice, you need to call both your signals, define a window size (fftsize, that is 512 in your case) and you can set the overlap so you achieve the desired 3 512-sample windows. They can also output the frequency array given you supply proper sampling rate. 

[Coxy, f] = coherence_ms(x, y, win, noverlap, fftsize, fs) is the coherence data with x and y as your input subsets, win is a window function with a length of the window (512 in your case), noverlap you have to calculate, fftsize is 512 and fs is your sampling rate (if you specify it in seconds, you get the frequency array in Hz, by the description it is approx. 2 s)

for the other you don't have to require the frequency range as an output because if you give the same time parameters it should have same frequency spacing

Cxy = crosspsd(x, y, win, noverlap, fftsize, fs) similarly to coherence

The for loop should take subsets of the original full length data with desired length (1024) and an appropriate offset (i*256). You can calculate the number of iterations needed. (TSA Welch has a deeply buried subVI called "tsa_window para.vi" which you can use as a template for that)

[Fonsi]

Hi Jean,

Thanks for your time and help!.

Well, I tried to do the steps you propused. But I dont understand well the point 4.

I detrend the full data previously and then reshape the data in rows of 1024.

The MathScript I believe it is correct. But the point 4, I dont undestand to find the time of window. 

The signal is interpolated previously with 2 samples per second. So 1024 samples would be 512 seconds.

I attached one file, trying to accomplish it. 

Please, if you can explain me I dont know if I do each 216 samples inside of each 1024 sample row.

Regards, Fred.

[jeanlucpicard · Accepted solution]

you made a few mistakes, fr example by reshaping it you cannot do overlapping segments, that's why you should go with a loop and subsetting (and that's why i suggested tsa window para)

you should also make the product output indexing and not concatenating, sou you have 2d output in the end which you can use for 3d graphs

fs in the mathscipt should also be in hz, not in seconds, as i told you last time (it is stated in help)

[Fonsi]

Thank you Jean,

Now I understand it.