it would be a good idea to examine this relation in a multiscale fashion, i.e. via wavelet cross-spectrum (magnitude and phase), or wavelet coherence, as different scales may be involved.

Then one may also want to move beyond linear analysis, and employ for example entropy related features.

Dear Ali,

Transfer function analysis in Matlab can be performed by means of "cpsd" function, which estimates the cross power spectral density of two signals. From the results of "cpsd" and the individual PSDs, you may calculate gain, phase and coherence. Coherence can also be estimated directly through "mscohere" function, if necessary.

We have used this analysis in a recent study (DOI: 10.1113/JP274065). Let me know if your needs are similar to what we did.

Best regards,

Luiz

Hi Luiz,

Yes, this should be helpful, I'll try it.

We have same time series as you had in your study except the MSNA, but I guess the method is same. With 2 time series of SBP and IBI, I should be able to follow the method you have used. I already did the sequence method but I would like to compare it with the spectral method and have information on the phase since the gain is not different between conditions based on the sequence method.

Tnx

Ali

Hello Ali,

Yes, no matter if the time series is MSNA or IBI. The procedure is the same.

Even though you have already checked the gain by the sequence method, I would suggest you to calculate gain from TF too. As it can be calculated in different frequency bands, it can add some useful information. Moreover, I would also pay attention to BEI of the sequence method, which carry important information about the effectiveness between the two variables (effectivenes of one changing the another).

Finally, if you are interested in the phase between the two components, I would recommend you the synchrogram or coordigram, which gives you an directly estimate about it. See Chaos 27, 093933 (2017) for a nice reference.

Best,

Luiz

Hi Luiz,

So, using pwelch I now have PSD of SBP and IBI which are in double format.

[PSD_IBI, f] = pwelch(x, ...

[PSD_SBP, f] = pwelch(y, ...

Then, from cpsd I got a complex double output.

[pxy, f] = cpsd(y,x, ...

the angle(pxy) is the phase, correct?

and the abs(pxy) is the magnitude of transfer at each f?

To calculate BRS, the magnitude (abs(pxy)) must be divided by PSD of SBP, correct?

abs(pxy)./PSD_SBP

and then it must be integrated in HF and LF bands (with a 0.5 coherence threshold), correct?

Tnx for your help

Ali

Hello Ali,

Your framework sound pretty good to me. There is only one possible correction to do. The TF gain is estimated as:

TFxy = Pyx/Pxx.

In your case, as you are studying the effects of IBI on SBP, Pxx is PSD_IBI (as you have defined).

The coherence threshold is optional. It is intended to act like a filter for "confident estimatives". It can be calculated by:

Cxy = [abs(Pxy) * abs(Pxy)] / [Pxx * Pyy].

For more details, see the Matlab documentation for "mscohere" and "tfestimate". Although we could use those functions directly to estimate Cxy and TFxy, it would not be efficient because CPSD would be calculated more than one time.

Best,

Luiz

Tnx Luiz,

Actually, I'm looking for BRS defined as how much change in IBI we have for each unit of change in SBP (ms/mmHg), so the input is SBP(y) and output is IBI(x).

Q1:

[Pyx, f] = cpsd(y,x, ...

Is this correct considering y as input and x as output?

Q2:

about the TFyx = Pyx/Pyy

the Pyx is a complex double output. Shouldn't it be abs(Pyx)/Pyy? so the magnitude of cross-spectrum of yx at each f over the PSD of y at that f ?

Q3:

so you recommend using Cyx = [abs(Pyx) * abs(Pyx)] / [Pxx * Pyy] instead of mscohere?

Ok. So let's redefine the framework using the classical terminology of variables to avoid confusion.

For the causal relationship x -> y, we have:

x: SBP

y: IBI

Q1:

Pxx = pwelch(x, ...);

Pyy = pwelch(y, ...);

Pyx = cpsd(y, x, ...);

Txy = Pyx/Pxx;

where Txy is the transfer function estimate for x -> y.

Q2:

Actually, in Matlab you have to code Txy as:

Txy = Pyx./Pxx;

so that each freq. component of Pyx will be divided by the corresponding in Pxx. Txy is also a complex function. After that, you may take the magnitude of Txy.

Q3:

Cxy or Cyx gives exactly the same results as "mscohere". My point is that "mscohere" will calculate the cross-psd even though you have already done it to calculate Txy. Thus, this more efficient to use the individual spectra (Pxx, Pyy and Pyx) to derivate whatever you need, because all PSD estimates will be calculated only once.

Try calculating Txy and Cyx from individual spectra and compare the results obtained from "tfestimate(x, y, ...)" and "mscohere(x, y, ...)". It must be the same.

Best,

Luiz

Hi Luiz, Tnx for this, a great help indeed !

For BRS would you recommend averaging in each frequency band or the gain at peak coherence in each band? Also, would you recommend having coherence of 0.5 as the condition?

Hello Ali,

Well, you may find arguments from both sides about the coherence threshold. I have done some tests in the past and it seemed to me that assuming the coherence threshold would be a better option.

Regarding taking the average or the peak coherence, I have never compared the two options. Intuitivelly, taking only one value to represent the gain seems not to be much robust, even though coherence is the highest. Imagine that within the frequency band of interest there are many values with high coherence. They may have different gains and taking only one will turn the real gain poorly estimated. Assuming a coherence threshold acts as a balance factor because you will average the gain, but only considering those values with enough coherence.

In the study I mentioned before, we decided to average the gain in LF and HF bands, taking only those values where the coherence was higher than 0.5.

Best,

Luiz

Tnx to Luiz here what I could get from the spectral analysis.

Does it look normal? how to deal with the harmonics?

If I perform the PSD based on autoregressive model, there is no such harmonic, but the cpsd is still based on FFT and harmonics are there anyway.

Hi Ali,

The harmonics seems to be present in SDB curves (and maybe in PB). Was SDB data collected under any driven component, let's say paced respiration at 0.1 Hz? If yes, I would guess that the driven component is also creating components at its harmonics, just like in music. But this is only a guess.

Best,

Luiz

Hi Luiz,

SDB = slow deep breathing at 0.1 Hz

PB = paced breathing at 0.23 Hz

is there any physiological mechanism behind harmonics? or its due to the way that we process and analyze the signal, the FFT?

Well, I am not sure about that. Maybe @AlbertoPorta can shed light on this. I know he has a lot of experience with transfer functions.

Best,

Luiz

dear Ali

you can find some suggestions of SAP-HP phase and gain analysis in A. Porta et al, Am. J. Physiol. 300, R378-R386, 2011 and J.C. Milan-Mattos et al, J. Appl. Physiol. 124, 791-804, 2018.

All the best, Alberto

Dears,

@Porta thanks for your response and sharing your amazing studies on this topic.

Does MATLAB have a function for cross power spectral density using the AR method? or it's only with Welch's method?

Thanks

Ali