Coherence of two EEG signals
in [DSP]
|
Prev: FFT / IFFT
Next: TI 320C6713
From: gilmar on 6 Jan 2006 17:05 I am interested in comparing segments of minute-long EEG records from two different brain regions. I want to determine the coherence of the two signals before an 'event' and after the 'event', with the hypothesis that the two signals would be more coherent in a given frequency range after the event than before it (perhaps reflecting synchronization of neuronal activity). I want a 0.5Hz resolution so I'm using 2-s windows to calculate FFTs. I'm trying to program this in VB, while also using Matlab to generate the results in order to confirm that my VB program is working correctly. However, here is my problem. I am using 2 test sine waves (4Hz and 4+8Hz) to test the program, but I do not get the coherence graph I expect in either VB or Matlab, and I'm running around in circles trying to determine which definition/equation/function of coherence is the one I need. For these two signals, I'm expecting a coherence graph where the value is 0 at all frequencies except 4Hz, where it would be 1. Instead, using mscohere in Matlab7, I get all values at 1, except for a few spurious lower values at non-meaningful frequencies. When I add noise to my test sine waves, 4Hz remains at coherence 1 and the rest are below 1, but not 0. If I calculate the cross- and auto-spectra and plug them into the following equation (in Matlab): Coh = sqrt((abs(SCab).^2)./(SCaa.*SCbb)); (where SCab = cross-spectral density and SCaa and SCbb are the auto-spectra), I get a similarly strange graph. UNLESS I change the equation to Coh = sqrt((abs(SCab).^2)/(SCaa.*SCbb)); (no dot before the division--i.e. changing it from array division to matrix division). Then I get the graph I expect (all frequencies but 4hz have coherence=0), but the coherence value at 4hz is around 0.4 instead of 1. I do not know why changing the division causes this. Also, if I calculate the cross- and auto-CORRELATIONS then use the latter equation, where SCab now represents the FFT of the cross-correlation and SCaa and SCbb are the FFTs of the auto-correlations, the graph is also what I expect, and the value at 4hz is 1.0. This is great, except I need to know why this approach worked and not the others so that I can 1) program the equations and loops, etc. correctly in VB and 2)understand and interpret my results correctly when I finally apply the analysis to EEG data. If someone could please explain why these different methods of calculating coherence produce different results and which one(s) is correct, I would greatly appreciate the help. Thanks, Marieke
From: Rune Allnor on 6 Jan 2006 18:27 gilmar wrote: > I am interested in comparing segments of minute-long EEG records from two > different brain regions. I want to determine the coherence of the two > signals before an 'event' and after the 'event', with the hypothesis that > the two signals would be more coherent in a given frequency range after > the event than before it (perhaps reflecting synchronization of neuronal > activity). I want a 0.5Hz resolution so I'm using 2-s windows to > calculate FFTs. I'm trying to program this in VB, while also using Matlab > to generate the results in order to confirm that my VB program is working > correctly. The VB + matlab combo is a very good way to do things. I like it. > However, here is my problem. I am using 2 test sine waves (4Hz and 4+8Hz) > to test the program, but I do not get the coherence graph I expect in > either VB or Matlab, and I'm running around in circles trying to determine > which definition/equation/function of coherence is the one I need. For > these two signals, I'm expecting a coherence graph where the value is 0 at > all frequencies except 4Hz, where it would be 1. Instead, using mscohere > in Matlab7, I get all values at 1, except for a few spurious lower values > at non-meaningful frequencies. When I add noise to my test sine waves, > 4Hz remains at coherence 1 and the rest are below 1, but not 0. First, estimating coherence with noise-free synthetic signals is not a good thing to do. I have seen the coherence defined as c(f) = |Sxy(f)|^2/Sxx(f)Syy(f) (1) where Sxy(f) is the cross spectrum between x(t) and y(t), and Sxx(f) and Syy(x) are the respective autospectra. According to eq. (1) one might ask how it is possible at all to get a result NOT equal 1 for all frequencies, regardless of whether there is noise present or not. I think the key is in how the spectra are computed, as your results quoted below seem to indicate. > If I calculate the cross- and auto-spectra and plug them into the > following equation (in Matlab): > Coh = sqrt((abs(SCab).^2)./(SCaa.*SCbb)); > (where SCab = cross-spectral density and SCaa and SCbb are the > auto-spectra), I get a similarly strange graph. Seems to be the square root of my eq. (1) above. > UNLESS I change the > equation to Coh = sqrt((abs(SCab).^2)/(SCaa.*SCbb)); (no dot before the > division--i.e. changing it from array division to matrix division). Then I > get the graph I expect (all frequencies but 4hz have coherence=0), but the > coherence value at 4hz is around 0.4 instead of 1. > I do not know why changing the division causes this. This is interesting. Could you please post the dimensions of the cross- and autospectra at the time you perfrm the above divisions? > Also, if I calculate the cross- and auto-CORRELATIONS then use the latter > equation, where SCab now represents the FFT of the cross-correlation and > SCaa and SCbb are the FFTs of the auto-correlations, the graph is also > what I expect, and the value at 4hz is 1.0. This is great, except I need > to know why this approach worked and not the others so that I can 1) > program the equations and loops, etc. correctly in VB and 2)understand and > interpret my results correctly when I finally apply the analysis to EEG > data. First: I really like the way you insist on understanding what's going on; it is a very welcome break from a few too many posts asking "I have this homework for tomorrow, can somebody post some code that solves the problem." As for the question about why the time-domain correlation + FFT works, it is because that's the definition of auto- and cross spectra. Most of the time, we (at least I) make the short-cut that Sxx'(f) = |X(f)|^2. (2) where X(f) is the N point Discrete Fourier Transform of some signal x(t) of length N. This is not true. Formally, Sxx(f) is defined as Sxx(f) = DFT{rxx(t)} (3) which in time domain is of length 2N-1. A different expression than (2). I have just started looking into this, but it seems the key is in that the spectra somehow are become different when they are computed by the two methods. I can't quite work it out over the keyboard right now, but I find the problem intriguing and I hope to spend a couple of hours with it in the next few days. > If someone could please explain why these different methods of calculating > coherence produce different results and which one(s) is correct, I would > greatly appreciate the help. You see a dramatic difference between the two approaches. The correlate in time domain + FFT gives the correct result straight out, and also complies to the formal definitions, so I think you should use that approach. I'll have to work a bit on explaining why the FFT + correlate approach doesn't work. It is an intriguing question, so I will come back to it in a few days. Rune
From: Rune Allnor on 6 Jan 2006 21:43 Rune Allnor wrote: > First, estimating coherence with noise-free synthetic signals is not a > good > thing to do. I have seen the coherence defined as > > c(f) = |Sxy(f)|^2/Sxx(f)Syy(f) > (1) > > where Sxy(f) is the cross spectrum between x(t) and y(t), and Sxx(f) > and Syy(x) are the respective autospectra. According to eq. (1) one > might ask how it is possible at all to get a result NOT equal 1 for all > > frequencies, regardless of whether there is noise present or not. > > I think the key is in how the spectra are computed, as your results > quoted below seem to indicate. --- > As for the question about why the time-domain correlation + FFT > works, it is because that's the definition of auto- and cross spectra. > Most of the time, we (at least I) make the short-cut that > > Sxx'(f) = |X(f)|^2. (2) > > where X(f) is the N point Discrete Fourier Transform of some signal > x(t) of length N. This is not true. Formally, Sxx(f) is defined as > > Sxx(f) = DFT{rxx(t)} (3) > > which in time domain is of length 2N-1. A different expression than > (2). > I have just started looking into this, but it seems the key is in that > the spectra somehow are become different when they are computed > by the two methods. I can't quite work it out over the keyboard right > now, > but I find the problem intriguing and I hope to spend a couple of hours > > with it in the next few days. > > > If someone could please explain why these different methods of calculating > > coherence produce different results and which one(s) is correct, I would > > greatly appreciate the help. > > You see a dramatic difference between the two approaches. The > correlate in time domain + FFT gives the correct result straight out, > and also complies to the formal definitions, so I think you should > use that approach. I'll have to work a bit on explaining why the > FFT + correlate approach doesn't work. > > It is an intriguing question, so I will come back to it in a few days. This question is more than intriguing, it is almost obsessing. Let's have a look at the basic equations. The Wiener-Kintchine theorem goes as follows: Define the cross correlation between x(n) and y(n) as Rxy(m) = sum[n,-inf,inf] x(n+m)y'(n) where sum[n,-inf,inf] means "the sum over n from -infinite to infinite". Then, the cross spectrum Sxy(w) is defined as Sxy(w) = sum[m,-inf,inf] sum[n,-inf,inf] x(n+m)y'(n) exp(-jmw) Substitute m = q-n to find Sxy(w) = sum[q,-inf,inf] sum[n,-inf,inf] x(q)y'(n) exp(-j(q-n)w) = sum[q,-inf,inf]x(q)exp(-jqw)sum[n,-inf,inf] y'(n) exp(jnw) = X(w)Y'(w). where X(w) = sum[n,-inf,inf]x(n)exp(-jnw) Y(w) = sum[n,-inf,inf]y(n)exp(-jnw) ans x' is the conjugate of x. The same line of arguments yield Sxx(w) = |X(w)|^2 Syy(w) = |Y(w)|^2 For a given w, represent the complex-valued spectrum coefficients on polar form, X(w)= R_x exp(jphi_x) Y(w)= R_y exp(jphi_y) and insert in (1) above: c(w) = |R_xR_y|^2/R_x^2R_y^2 = 1. So the question is not why the spectrum representation does not work. The question is why the time domain correlation + FFT works at all. I suspect this must have to do with the properties of the DFT. In the derivation above, I have used -infinite and +infinite as summation limits in both the autocorrelation sequence and the Fourier transform. My gut feeling right now (which may have as much to do with WAY too much coffee the last few days...) is that something funny happens due to windowing caused by the truncated summations. This will be another weekend of insomnia... Rune
From: Rune Allnor on 6 Jan 2006 22:01 Rune Allnor wrote: > gilmar wrote: --- > First, estimating coherence with noise-free synthetic signals is not a > good > thing to do. I have seen the coherence defined as > > c(f) = |Sxy(f)|^2/Sxx(f)Syy(f) > (1) > > where Sxy(f) is the cross spectrum between x(t) and y(t), and Sxx(f) > and Syy(x) are the respective autospectra. According to eq. (1) one > might ask how it is possible at all to get a result NOT equal 1 for all > > frequencies, regardless of whether there is noise present or not. > > I think the key is in how the spectra are computed, as your results > quoted below seem to indicate. --- > As for the question about why the time-domain correlation + FFT > works, it is because that's the definition of auto- and cross spectra. > Most of the time, we (at least I) make the short-cut that > > Sxx'(f) = |X(f)|^2. (2) > > where X(f) is the N point Discrete Fourier Transform of some signal > x(t) of length N. This is not true. Formally, Sxx(f) is defined as > > Sxx(f) = DFT{rxx(t)} (3) > > which in time domain is of length 2N-1. A different expression than > (2). > I have just started looking into this, but it seems the key is in that > the spectra somehow are become different when they are computed > by the two methods. I can't quite work it out over the keyboard right > now, > but I find the problem intriguing and I hope to spend a couple of hours > > with it in the next few days. > > > If someone could please explain why these different methods of calculating > > coherence produce different results and which one(s) is correct, I would > > greatly appreciate the help. > > You see a dramatic difference between the two approaches. The > correlate in time domain + FFT gives the correct result straight out, > and also complies to the formal definitions, so I think you should > use that approach. I'll have to work a bit on explaining why the > FFT + correlate approach doesn't work. > > It is an intriguing question, so I will come back to it in a few days. This question is more than intriguing, it is almost obsessing. Let's have a look at the basic equations. The Wiener-Kintchine theorem goes as follows: Define the cross correlation between x(n) and y(n) as Rxy(m) = sum[n,-inf,inf] x(n+m)y'(n) where sum[n,-inf,inf] means "the sum over n from -infinite to infinite". Then, the cross spectrum Sxy(w) is defined as Sxy(w) = sum[m,-inf,inf] sum[n,-inf,inf] x(n+m)y'(n) exp(-jmw) Substitute m = q-n to find Sxy(w) = sum[q,-inf,inf] sum[n,-inf,inf] x(q)y'(n) exp(-j(q-n)w) = sum[q,-inf,inf]x(q)exp(-jqw)sum[n,-inf,inf] y'(n) exp(jnw) = X(w)Y'(w). where X(w) = sum[n,-inf,inf]x(n)exp(-jnw) Y(w) = sum[n,-inf,inf]y(n)exp(-jnw) ans x' is the conjugate of x. The same line of arguments yield Sxx(w) = |X(w)|^2 Syy(w) = |Y(w)|^2 For a given w, represent the complex-valued spectrum coefficients on polar form, X(w)= R_x exp(jphi_x) Y(w)= R_y exp(jphi_y) and insert in (1) above: c(w) = |R_xR_y|^2/R_x^2R_y^2 = 1. So the question is not why the spectrum representation does not work. The question is why the time domain correlation + FFT works at all. I suspect this must have to do with the properties of the DFT. In the derivation above, I have used -infinite and +infinite as summation limits in both the autocorrelation sequence and the Fourier transform. My gut feeling right now (which may have as much to do with WAY too much coffee the last few days...) is that something funny happens due to windowing caused by the truncated summations. This will be another weekend of insomnia... Rune
From: Rune Allnor on 6 Jan 2006 23:30
Rune Allnor wrote [duplicate posts]: Sorry for the duplicate(s). Google had a hickup and didn't seem to react to my first couple of attempts to post. Rune |