Coherence of two EEG signals

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

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

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

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

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

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

"Rune Allnor" <allnor@tele.ntnu.no> wrote in message
news:1136590023.538910.188560@g47g2000cwa.googlegroups.com...
>
> 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)
>
>
This is magnitude squared coherence (MSC).

Flora

>
Rune, thanks for your responses.  I'm glad to have provided an intriguing
problem, even if some of the mystery turns out to be due to my limited
knowledge of coherence and Matlab.  Here are a couple more things:

>> 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.
>
>This is interesting. Could you please post the dimensions of the cross-
>and autospectra at the time you perfrm the above divisions?

The value of 0.4 would have been 1.  I was testing my programs with two
sets of test signals.  A 4Hz and 8 Hz signal which I produced with a
wave-generator and then recorded through our system and then separately, a
4 and 8 hz sine wave that I made within Matlab.  So the value of 0.4
happened because I was using the wave-generator signals which appeared to
have different amplitudes.  Using the Matlab sine waves, both the
correlate+FFT and cross/auto-spectra versions of this equation, Coh =
sqrt((abs(SCab).^2)/(SCaa.*SCbb)), gave the same value of 1 at 4Hz.

As for the dimensions of the spectra, all were 1025 x 1 before array
division and the output was also 1025 x 1.  But when I changed it to
matrix division, the output dimensions of Coh were 1025 x 1025.  Although
the graph looked nice with a single peak of 1 at 4Hz, there were actually
1025 graphs drawn, all flatlines at 0, except for one graph, which had the
solitary peak at 1.  It turned out to be the 9th column in the 1025 x 1025
matrix.  This makes some intuitive sense to me, given that the 9th
position of the FFT outputs at 0.5Hz resolution is the 4hz signal.  So, is
this Coh output really just another spectrum, displaying only the
frequencies common to both signals and having a value between 1 and 0?  If
this isn't coherence, what is it?  Certainly, the array division version of
the equation gives me a similar output as the mscohere function.  I've done
a bit more reading about coherence from some of the links in this newsgroup
and I've played around with windowing and more varied signals.  It seems
like the typical graph says just as much about where two signals do not
cohere as where they do, so the coherence graph of a 4hz wave and a 4+8hz
wave with noise added will have a value of 1 at 4hz, but 0 at 8hz, with
all other frequencies being somewhere in the middle.

Am I on the right track to think of these outputs like this?:

The matrix division version (output of a single peak at 4Hz) is telling me
what frequencies are common to both signals and that's it.
The array division version (coherence) is telling me what proportion of
one signal can be explained by the 2nd signal (via a transfer function). 
So a pure sine wave of 4hz would cohere with a 4+8hz wave at all
frequencies because a 4+8hz wave could be produced from a 4hz wave and a
specific transfer function.  But with noise (or other signals) added in,
the 4hz component can be explained, but not the 8hz given all the other
signals, so coherence at 4hz is 1 and at 8hz is 0.  I'm not sure I fully
believe this last part, but it seems to fit the graphs.

Marieke
  

gilmar wrote:
> >
> Rune, thanks for your responses.  I'm glad to have provided an intriguing
> problem, even if some of the mystery turns out to be due to my limited
> knowledge of coherence and Matlab.  Here are a couple more things:
>
> >> 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.
> >
> >This is interesting. Could you please post the dimensions of the cross-
> >and autospectra at the time you perfrm the above divisions?
>
> The value of 0.4 would have been 1.  I was testing my programs with two
> sets of test signals.  A 4Hz and 8 Hz signal which I produced with a
> wave-generator and then recorded through our system and then separately, a
> 4 and 8 hz sine wave that I made within Matlab.  So the value of 0.4
> happened because I was using the wave-generator signals which appeared to
> have different amplitudes.  Using the Matlab sine waves, both the
> correlate+FFT and cross/auto-spectra versions of this equation, Coh =
> sqrt((abs(SCab).^2)/(SCaa.*SCbb)), gave the same value of 1 at 4Hz.
>
> As for the dimensions of the spectra, all were 1025 x 1 before array
> division and the output was also 1025 x 1.  But when I changed it to
> matrix division, the output dimensions of Coh were 1025 x 1025.

I can't see why that happens without seeing your code, but a working
coherence computation should take two vectors as input and produce
one vector as output. The length of the output vector may differ
from the inputs, but it should be a vector, not a matrix.

> Although
> the graph looked nice with a single peak of 1 at 4Hz, there were actually
> 1025 graphs drawn, all flatlines at 0, except for one graph, which had the
> solitary peak at 1.  It turned out to be the 9th column in the 1025 x 1025
> matrix.  This makes some intuitive sense to me, given that the 9th
> position of the FFT outputs at 0.5Hz resolution is the 4hz signal.  So, is
> this Coh output really just another spectrum, displaying only the
> frequencies common to both signals and having a value between 1 and 0?  If
> this isn't coherence, what is it?

I don't know. There are too many important details missing in your post
to attempt an answer.

For now, if you get a vector, not a matrix, that shows the expected
behaviour by first computing time domain correlations and computing
the spectra of those, that is the method you should use.

Rune

Actually, FFT is rather crude and inappropriate technique for analyzing
something as complex as EEG signal...

If I were you I would look closely at HHT (Hilbert-Huang Transform),
which was developed (and patented) by NASA.

Go to http://www.fuentek.com/technologies/hht.htm