Coherence corticomuscular with mscohere

Hi

I'm trying to calculate cortico-muscular coherence (EEG, EMG) with 'mscohere' in various frequency bands, eg 8-13 Hz. I have blocks of 30 sec of a task interspaced with blocks of 30 sec of rest. I want to know how coherence changes as a function of time during the task, so across the 30 sec of the task, not rest. I plan to separate the 30 sec EMG and EEG signals into 30 1 sec bin and then use for each of 30 bins:

[cxy,fc]  =  mscohere(dataEMG, dataEEG, [], [], [8:0.1:13], 1000)
mean(cxy) 

Is that correct to calculate coherence in the 8-13 Hz across the 30 sec?

Thanks

help mscohere

 MSCOHERE   Magnitude Squared Coherence Estimate.
    Cxy = MSCOHERE(X,Y) estimates the magnitude squared coherence estimate
    of the system with input X and output Y using Welch's averaged,
    modified periodogram method. Coherence is a function of frequency with
    values between 0 and 1 that indicate how well the input X corresponds
    to the output Y at each frequency.  The magnitude squared coherence,
    Cxy, is given by Cxy = (abs(Pxy).^2)./(Pxx.*Pyy) where Pxx and Pyy are
    the power spectral density (PSD) estimate of X and Y, respectively; and
    Pxy is the Cross-PSD (CPSD) estimate of X and Y. See "help pwelch" and
    "help cpsd" for complete details.
 
    By default, X and Y are divided into eight sections with 50% overlap,
    each section is windowed with a Hamming window, and eight modified
    periodograms are computed and averaged.  See "help pwelch" and "help
    cpsd" for complete details.
 
    X and Y may be either vectors or two-dimensional matrices.  If both are
    matrices and have the same size, MSCOHERE operates column-wise:
    Cxy(:,n) = MSCOHERE(X(:,n),Y(:,n)).  If one is a matrix and the other
    is a vector, the vector is converted to a column vector and expanded
    internally so both inputs have the same number of columns.
 
    Cxy = MSCOHERE(X,Y,WINDOW), when WINDOW is a vector, divides each
    column of X and Y into overlapping sections of length equal to the
    length of WINDOW, and then windows each section with the vector
    specified in WINDOW. If WINDOW is an integer, each column of X and Y
    are divided into sections of length WINDOW, and each section is
    windowed with a Hamming window of that length.  If WINDOW is omitted or
    specified as empty, a Hamming window is used to obtain eight sections
    of X and Y.
     
    Cxy = MSCOHERE(X,Y,WINDOW,NOVERLAP) uses NOVERLAP samples of overlap
    from section to section. NOVERLAP must be an integer smaller than
    WINDOW, if WINDOW is an integer; or smaller than the length of WINDOW,
    if WINDOW is a vector. If NOVERLAP is omitted or specified as empty,
    it is set to obtain a 50% overlap.
 
    When WINDOW and NOVERLAP are not specified, MSCOHERE divides X into
    eight sections with 50% overlap and windows each section with a Hamming
    window. MSCOHERE computes and averages the periodogram of each section
    to produce the estimate.
 
    [Cxy,W] = MSCOHERE(X,Y,WINDOW,NOVERLAP,NFFT) specifies the number of
    FFT points, NFFT, used to calculate the PSD and CPSD estimates. For
    real X and Y, Cxy has length (NFFT/2+1) if NFFT is even, and (NFFT+1)/2
    if NFFT is odd. For complex X or Y, Cxy always has length NFFT. If NFFT
    is omitted or specified as empty, it is set to  either 256 or the next
    power of two greater than the length of each section of X (or Y),
    whichever is larger.
 
    If NFFT is greater than the length of each section, the data is
    zero-padded. If NFFT is less than the section length, the segment is
    "wrapped" (using DATAWRAP) to make the length equal to NFFT. This
    produces the correct FFT when NFFT is smaller than the section length.
 
    W is the vector of normalized frequencies at which Cxy is estimated.
    W has units of radians/sample.  For real signals, W spans the interval
    [0,pi] when NFFT is even and [0,pi) when NFFT is odd.  For complex
    signals, W always spans the interval [0,2*pi).
 
    [Cxy,W] = MSCOHERE(X,Y,WINDOW,NOVERLAP,W) computes the two-sided
    coherence estimate at the normalized angular frequencies contained in
    the vector W.  W must have at least two elements.
 
    [Cxy,F] = MSCOHERE(X,Y,WINDOW,NOVERLAP,NFFT,Fs) returns the coherence
    estimate as a function of physical frequency.  Fs is the sampling
    frequency specified in hertz.  If Fs is empty, it defaults to 1 Hz.
 
    F is the vector of frequencies (in hertz) at which Cxy is estimated.
    For real signals, F spans the interval [0,Fs/2] when NFFT is even and
    [0,Fs/2) when NFFT is odd.  For complex signals, F always spans the
    interval [0,Fs).
 
    [Cxy,F] = MSCOHERE(X,Y,WINDOW,NOVERLAP,F,Fs) computes the coherence
    estimate at the physical frequencies contained in the vector F.  F must
    be expressed in hertz and have at least two elements.
 
    [...] = MSCOHERE(...,'mimo') returns the multiple coherence function
    for a multiple-input/multiple-output system. The multiple coherence
    function of the ith signal in Y is contained in the ith column of Cxy.
    If X is a vector, multiple coherence is equivalent to ordinary
    magnitude squared coherence. If X and Y have a different number of
    columns and both have more than one column, this option is used by
    default.
 
    [...] = MSCOHERE(...,FREQRANGE) returns the coherence estimate computed
    over the specified range of frequencies based upon the value of
    FREQRANGE:
 
       'onesided' - returns the one-sided coherence estimate of real input
          signals X and Y. If NFFT is even, Cxy has length NFFT/2+1 and is
          computed over the interval [0,pi]. If NFFT is odd, Cxy has length
          (NFFT+1)/2 and is computed over the interval [0,pi). When Fs is
          optionally specified, the intervals become [0,Fs/2] and [0,Fs/2)
          for even and odd NFFT, respectively.
 
       'twosided' - returns the two-sided coherence estimate for either
          real or complex input X and Y. Cxy has length NFFT and is
          computed over the interval [0,2*pi). When Fs is specified, the
          interval becomes [0,Fs).
 
       'centered' - returns the centered two-sided two-sided coherence
          estimate for either real or complex X and Y.  Cxy has length NFFT
          and is computed over the interval (-pi, pi] for even length NFFT
          and (-pi, pi) for odd length NFFT. When Fs is specified, the
          intervals become (-Fs/2, Fs/2] and (-Fs/2, Fs/2) for even and odd
          NFFT, respectively.
 
       FREQRANGE may be placed in any position in the input argument list
       after NOVERLAP.  The default value of FREQRANGE is 'onesided' when X
       and Y are both real and 'twosided' when either X or Y is complex.
 
    MSCOHERE(...) with no output arguments plots the coherence estimate (in
    decibels per unit frequency) in the current figure window.
 
    % Example 1:
    %   Compute and plot the coherence estimate between two colored noise 
    %   sequences x and y.
    h = fir1(30,0.2,rectwin(31));   % Window-based FIR filter design
    h1 = ones(1,10)/sqrt(10);       
    r = randn(16384,1);             
    x = filter(h1,1,r);  % Filter the data sequence
    y = filter(h,1,x);   % Filter the data sequence  
    noverlap = 512; nfft = 1024;
    mscohere(x,y,hann(nfft),noverlap,nfft); % Plot estimate
 
    % Example 2:
    %   Compute and plot the multiple coherence estimate between x and y.  
    h1 = fir1(30,0.3,rectwin(31));
    h2 = fir1(30,0.5,rectwin(31));
    x = randn(16384,2);
    y = filter(h1,1,x(:,1)) + filter(h2,1,x(:,2));
    noverlap = 512; nfft = 1024;
    mscohere(x,y,hann(nfft),noverlap,nfft,'mimo'); % Plot estimate
 
    See also TFESTIMATE, CPSD, PWELCH, PERIODOGRAM.

    Documentation for mscohere
       doc mscohere

    Other uses of mscohere

       gpuArray/mscohere

Answer 1

[Star Strider]

The first output of mscohere is a vector, and is intended to show the magnitude-squared coherence of the two signals. It seems to me that calculating the mean of it destroys the information it provides. I would keep it (and others like it, since I assume this is part of a set of experiments) as vectors (along with any other necessary information), store them, and then process them. How you process them depends on the information you want from them.

[Pat]

Hi

thanks for answering

I'm not really inerested in differences within the 8-13 range, that's why I was thinking of using mean.

Also, I will look at others frequency ranges, lke Theta (4-8 Hz), Beta (13-30) and Gamma (30-50).

But, I wanted to know if the way I was usign mscohere was correct.

best

[Star Strider]

My pleasure!

The fifth argument to mscohere is the size of the fft. A vector is not going to work for that, so in that sense, no. That has to be a scalar, preferably a power of 2 (512, 1024, 4096, etc.).

I suggest looking at several different EEG bands, to see what works best. (My neurophysiology studies were not the sort you’re doing, so I’m not certain what the best approach would be. You need to do a PubMed search on similar research, if you’ve not already done so.) Obviously, leads over the motor cortex are going to be most important if you’re studying them with respect to EMG.

I get the impression that you’ve not completely understood what mscohere does and the output it produces. The result is a vector of amplitudes of the magnitude-squared coherence estimate, with (if you request it, and I recommend that) a matching vector of frequencies. Taking the mean of the first output destroys all that information, making the mscohere call essentially pointless. That’s the reason I suggest storing the vectors in a cell array or some other convenient array to compare later.

[Pat]

Hi Star

Maybe I should have asked how to compute coherence on a 30 sec blocks of a task, in different frequency bands (5-8, 8-13, 13-30, 30-50 Hz), with good enought temporal resolution to see what's changing during that 30 sec, maybe close to 1 sec but it could close to 2 and that should be fine too. That's why I divided my 30 sec into 1 sec bin and I though of computing coherence on these bins, like I did for other metrics. But maybe here it's different...

I see in the doc:

[cxy,f] = mscohere(x,y,window,noverlap,f,fs) returns the magnitude-squared coherence estimate at the frequencies specified in f.

f seems to be the frequencies of interest, that's why I wrote it like I did above.

[Star Strider]

For EEG and EMG analyses, the sampling frequency is usually at least 1 kHz, so 1 second is a fair amount of data. It would be relatively straightforward to create elliptic filters for the bands-of-interest using the bandpass function (use the 'ImpulseResponse','iir' name-value pair for best results) and getting both outputs from the first call to it. (After that, use the second output, that being the digital-filter object, with filtfilt to filter the other signals.)

Probably the most efficient way to segment the signals would be to use the buffer function to divide the signal into separate one-second intervals (the second argument to buffer would be the sampling frequency for 1 second segments) and use that as the input to filtfilt along with the digital-filter object. The filtered outputs would have the same dimensions as the inputs, so you would then present chosen columns of the filtered output matrices to mscohere and then store those results (both ‘cxy’ and ‘f’) for each result.

That is essentially how I would do it.

[Pat]

Hi Star Strider

My EEG/EMG are already processsed, I just need to calculate coherence to get a sense of how it changes in time within various frequency bands.

[Star Strider]

O.K.

I would use buffer to create the appropriate matrices (ideally the same sizes) and then give one of them to mscohere as the EMG matrix (perhaps ‘x’) and the other as the EEG matrix (‘y’). If they are the same sizes, the result will be a matrix of the magnitude-squared coherence vectors at the times corresponding to the columns in the argument matrices, and as a function of frequency, so keep at least one frequency vector as well (since they should all share it).

A waterfall plot might be an appropriate way of presenting the results graphcially, with frequency on one axis and time on the other.