[Eeglablist] why is the result of coherence not between 0~1? (3)

[wu xiang]

hi all:

I have found the bug, it's for carelessness when do averaging. Thanks
for the suggestions from philip grieve,

Sarnthein Johannes and Andre Achim. 

Now I want to post my file calculating coherence here, which do welch average using matlab function cpsd.m and pwelch.m (the matlab version is 7.01). Hope it be useful to somebody need it. 

---------------------------------------------------------------------------
function [coh,F] = coherence(x,y,nfft,fs)
%
% coherence		Return coherence of signals from two electrodes calculateing using welch average. 
% Usage	[coh,F] = coherence(x,y,fs)
% Input 
%	x,y -- Signals from two electrodes. 
%			if they are 2-D matrixs, rows are trials (i.e. trial 1:n), columns are sample points.
%				the returned coherence is multi-trial coherence, which is usually used in eeg analysis.
%			if they are 1-D vectors(row), they are sample points. 
%				the returned coherence is single trial coherence, as cohere.m or mscohere.m of matlab.
%	nfft -- number of FFT points. It should be power of 2.
%	fs -- sampling frequency.
% Output
%	coh -- returned coherence
%	F -- phisical frequency (Hz). coh is a function of F.


% get trial number to m, sample points to n
[m,n] = size(x); % or [m,n] = size(y);


%parameters of matlab power spectra density (psd) function, you can modify here if need
nwindow = []; % default 8 window
noverlap = []; % default 50% ovelap
% nfft, as input parameter

% the length of returned Psd or F is nfft/2+1 when nfft is even, or (nfft+1)/2 when nfft is odd. 
% the range of returned f is [0,fs/2].
% the length of returned psd or F should > fs/2+1, so that every frequency at least has a psd value.
len_Psd = nfft/2+1; % length of Psd returned by matlab function cpsd,pwelch. 

Sxy = zeros( len_Psd,1 ); % initialize cross power spectra
Sxx = zeros( len_Psd,1 ); % initralize auto power spectra
Syy = zeros( len_Psd,1 );


for i=1:m
	%get cross power spectra density using matlab function cpsd
	[Pxy,F] = cpsd( x(i,:),y(i,:),nwindow,noverlap,nfft,fs );  
	Pxy = Pxy.*fs; %convert power spectra density to power spectra 
	Sxy = Sxy + Pxy;

	%get auto power spectra density using matlab function pwelch
	[Pxx,F] = pwelch( x(i,:),nwindow,noverlap,nfft,fs ); 
	Pxx = Pxx.*fs; 
	Sxx = Sxx + Pxx;

	[Pyy,F] = pwelch( y(i,:),nwindow,noverlap,nfft,fs ); 
	Pyy = Pyy.*fs; 
	Syy = Syy + Pyy;
end

% average 
Sxy = Sxy/m;
Sxx = Sxx/m;
Syy = Syy/m; % the mistake I have made is that, I have writen Syy = Sxx/m :)

% get coherence
coh = abs(Sxy).^2 ./ (Sxx .* Syy);

----------------------------------end-------------------------------------



-- 
wu xiang
Cognitive neuroscience lab of life school of university of science and technology of china
Hefei, anhui, P.R.China. 230027
web-page:http://202.38.64.10/~rwfwu
You can mail to: rwfwu at ustc.edu or rwfwuwx at gmail.com