matlab scripts by pierre for calculating GFP and dissimilarity

[seba]

Hello Cartoolers,
I was curious to see Pierre's Matlab scripts he made (april 2008:
http://groups.google.com/group/cartool/browse_thread/thread/162d69ff4e5b8ae8/4973bf8165ff4e62?lnk=gst&q=dissimilarity#4973bf8165ff4e62)
for calculating GFP and dissimilarity.
However, I am not able to download them (the .zip file found here:
http://groups.google.com/group/cartool/files
seems to be empty).
Am I doing something wrong? Pierre, if you read this...do you still
have these scripts?
Does somebody else have matlab scripts for doing GFP and
dissimilarity?
Thank you very much in advance
Best
Sebastian

[Denis]

Use the right-click, and save the file when prompted, left-clicking
does not lead to a page. Thanks Google Group... :-(

[Pierre]

Hello Sebastian,

I think I still have the individual scripts in my old computer at
home. I will check this and send them to you.

Pierre

On 13 juil, 15:50, seba <ders...@gmx.de> wrote:
> Hello Cartoolers,

> I was curious to see Pierre's Matlab scripts he made (april 2008:http://groups.google.com/group/cartool/browse_thread/thread/162d69ff4...)

[seba]

Yup, tried it, but doesn't work for me. .zip folder is empty or
corrupted.

[Pierre]

Hi Sebastian and everyone,

Since there appears to be a problem with opening the old zip file, and
since I'm unable to transfer any file from my HCUGE PC to this group,
I will simply paste a few of the scripts below. I include the
following:
- openeph.m
- opensef.m
- computeavgref.m
- computedataagainstref.m
- computegfp.m
- computedissandsc.m
- saveeph.m
- savesef.m

All these scripts have a header with a short description of the
variable names and formats used.

All the caveats in my original posting remain! (http://www.google.com/
url?sa=D&q=http://groups.google.com/group/cartool/browse_thread/thread/
162d69ff4e5b8ae8/4973bf8165ff4e62%3Flnk%3Dgst%26q%3Ddissimilarity
%234973bf8165ff4e62&usg=AFQjCNHMeIfzTzPT33feWVpQbhfYne4r4A)

Please signal here any problem with these scripts.

Pierre (PS the e-mail address mentioned in the scripts does no longer
exist: please use pierre....@gmail.com)


function
[numchannels,numtimeframes,samplingrate,thedata]=openeph(openfilename)
% openeph: opens a Cartool evoked potential data file (.ep(h))
%
% inputs: full path and name of the file to open
%
% outputs: number of channels, number of timeframes (and sampling rate
for
% .eph files) as 1-D numeric arrays; data as a 2-D numeric array where
% dimension 1 contains the timeframes, dimension 2 contains the
channels
%
% Cartool: http://brainmapping.unige.ch/Cartool.htm
%
% author of this script: pierre....@medecine.unige.ch


% open filename for reading in text mode
fid=fopen(openfilename,'rt');

% for .eph files
if strcmp(openfilename(end-3:end),'.eph')==1

% read header
eph=textscan(fid,'%s','delimiter','/n');
eph=eph{1};
numchannels=sscanf(eph{1},'%f',1);
numtimeframes=sscanf(eph{1},'%*f %f',1);
samplingrate=sscanf(eph{1},'%*f %*f %f',1);

% prepare for reading data
formatstring='%f';
if numchannels>1
for i=1:numchannels-1
formatstring=[formatstring ' %f'];
end
end

% read data
thedata=zeros(numtimeframes,numchannels);
for i=1:numtimeframes
thedata(i,:)=sscanf(eph{i+1},formatstring);
end

% for .ep files
elseif strcmp(openfilename(end-2:end),'.ep')==1
thedata='';
while ~feof(fid)
thedataline=fgetl(fid);
thedata=strvcat(thedata,thedataline);
end
thedata=str2num(thedata);
numtimeframes=size(thedata,1);
numchannels=size(thedata,2);
samplingrate=0;

else
error('incorrect file type');
end

% close file
fclose(fid);


function
[numchannels,numtimeframes,samplingrate,thedata]=opensef(openfilename)
% opensef: opens a Cartool simple EEG data file (.sef)
%
% inputs: full path and name of the file to open
%
% outputs: number of channels, number of timeframes and sampling rate
as
% 1-D numeric arrays; data as a 2-D numeric array where dimension 1
% contains the timeframes, dimension 2 contains the channels
%
% Cartool: http://brainmapping.unige.ch/Cartool.htm
%
% author of this script: pierre....@medecine.unige.ch


% open filename for reading
fid=fopen(openfilename,'r');

% read fixed part of header
version=strcat(fread(fid,4,'int8=>char')');
numchannels=fread(fid,1,'int32');
numauxchannels=fread(fid,1,'int32');
numtimeframes=fread(fid,1,'int32');
samplingrate=fread(fid,1,'float32');
year=fread(fid,1,'int16');
month=fread(fid,1,'int16');
day=fread(fid,1,'int16');
hour=fread(fid,1,'int16');
minute=fread(fid,1,'int16');
second=fread(fid,1,'int16');
millisecond=fread(fid,1,'int16');

% read variable part of header
channels=strcat(fread(fid,[8,numchannels],'int8=>char')');

% read data
thedata=fread(fid,[numchannels,numtimeframes],'float32')';

% close file
fclose(fid);


function avgref=computeavgref(thedata)
% computeavgref: computes the average reference of a data set
%
% inputs: data as a 2-D numeric array where dimension 1 contains the
% timeframes, dimension 2 contains the channels
%
% outputs: average reference as a 1-D numeric array that contains the
% average reference at each timeframe
%
% Cartool: http://brainmapping.unige.ch/Cartool.htm
%
% author of this script: pierre....@medecine.unige.ch


% define number of channels and time frames
numtimeframes=size(thedata,1);
numchannels=size(thedata,2);

% compute average reference
avgref=zeros(numtimeframes,1);
for i=1:numtimeframes
avgref(i)=sum(thedata(i,:))/numchannels;
end


function thedata=computedataagainstref(thedata,ref)
% computedataagainstref: changes the reference of a data set
%
% inputs: data as a 2-D numeric array where dimension 1 contains the
% timeframes, dimension 2 contains the channels; ref as either a
number
% (1-D numeric array) or the 'avgref' string for average reference
%
% outputs: data as a 2-D numeric array where dimension 1 contains the
% timeframes, dimension 2 contains the channels
%
% Cartool: http://brainmapping.unige.ch/Cartool.htm
%
% author of this script: pierre....@medecine.unige.ch


% define number of channels and time frames
numtimeframes=size(thedata,1);
numchannels=size(thedata,2);

% compute data against selected reference
if ischar(ref)==1&strcmp(ref,'avgref')==1
avgref=computeavgref(thedata);
for i=1:numtimeframes
thedata(i,:)=thedata(i,:)-avgref(i);
end
elseif ischar(ref)==1&strcmp(ref,'avgref')==0
error('Unknown code for reference');
elseif isnumeric(ref)==1&ref>=1&ref<=numchannels
for i=1:numtimeframes
thedata(i,:)=thedata(i,:)-thedata(i,ref);
end
elseif isnumeric(ref)==1&(ref<1|ref>numchannels);
error('Reference specified does not exist');
end


function gfp=computegfp(thedata)
% computegfp: calculates the global field power of a data set
%
% inputs: data as a 2-D numeric array where dimension 1 contains the
% timeframes, dimension 2 contains the channels
%
% outputs: global field power as a 1-D numeric array that contains the
% global field power at each timeframe
%
% Cartool: http://brainmapping.unige.ch/Cartool.htm
%
% author of this script: pierre....@medecine.unige.ch


% define number of channels and time frames
numtimeframes=size(thedata,1);
numchannels=size(thedata,2);

% compute global field power
gfp=zeros(numtimeframes,1);
avgref=computeavgref(thedata);
for i=1:numtimeframes
gfp(i)=sqrt(sum((thedata(i,:)-avgref(i)).^2)/numchannels);
end


function [diss,sc]=computedissandsc(thedata1,thedata2)
% computediss: computes dissimilarity and spatial correlation between
2
% data sets
%
% inputs: data as 2-D numeric arrays where dimension 1 contains the
% timeframes, dimension 2 contains the channels
%
% outputs: diss and sc as 1-D numeric arrays that contain the
% dissimilarity and spatial correlation at each timeframe
%
% Cartool: http://brainmapping.unige.ch/Cartool.htm
%
% author of this script: pierre....@medecine.unige.ch


% verify equal number of channels and time frames
if size(thedata1,1)~=size(thedata2,1)
error('Number of timeframes is different');
elseif size(thedata1,2)~=size(thedata2,2)
error('Number of channels is different');
end

% compute dissimilarity
diss=zeros(size(thedata1,1),1);
for i=1:size(thedata1,1)
diss(i)=sqrt(sum((thedata1(i,:)-thedata2(i,:)).^2)/
size(thedata1,2));
end

% compute spatial correlation
sc=zeros(size(diss,1),1);
for i=1:size(diss,1)
sc(i)=(2-diss(i)^2)/2;
end


function saveeph(savefilename,thedata,varargin)
% saveeph: saves data as a Cartool evoked potential data file (.ep(h/
sd/se))
%
% inputs: full path and name of the file to save, data as a 2-D
numeric
% array where dimension 1 contains the timeframes, dimension 2
contains the
% channels, (optional) samplingrate as 1-D numeric array passed as
% 3rd argument
%
% outputs: .ep(h/sd/se) file
%
% Cartool: http://brainmapping.unige.ch/Cartool.htm
%
% author of this script: pierre....@medecine.unige.ch


if strcmp(savefilename(end-3:end),'.eph')==1
numtimeframes=size(thedata,1);
numchannels=size(thedata,2);
samplingrate=varargin{1};
theheader=[numchannels numtimeframes samplingrate];
dlmwrite(savefilename,theheader,'\t');
dlmwrite(savefilename,thedata,'delimiter','\t','-append');
elseif strcmp(savefilename(end-2:end),'.ep')==1
dlmwrite(savefilename,thedata,'delimiter','\t','-append');
elseif strcmp(savefilename(end-4:end),'.epsd')==1|
strcmp(savefilename(end-4:end),'.epse')==1
if nargin==2
dlmwrite(savefilename,thedata,'delimiter','\t','-append');
elseif nargin==3
numtimeframes=size(thedata,1);
numchannels=size(thedata,2);
samplingrate=varargin{1};
theheader=[numchannels numtimeframes samplingrate];
dlmwrite(savefilename,theheader,'\t');
dlmwrite(savefilename,thedata,'delimiter','\t','-append');
end
end


function savesef(savefilename,thedata,samplingrate)
% savesef: saves data as a Cartool simple EEG data file (.sef)
%
% inputs: full path and name of the file to save, data as a 2-D
numeric
% array where dimension 1 contains the timeframes, dimension 2
contains the
% channels, samplingrate as 1-D numeric array
%
% outputs: .sef file
%
% Cartool: http://brainmapping.unige.ch/Cartool.htm
%
% author of this script: pierre....@medecine.unige.ch


% define fixed part of header
version='SE01';
numchannels=size(thedata,2);
numauxchannels=0;
numtimeframes=size(thedata,1);
year=0;
month=0;
day=0;
hour=0;
minute=0;
second=0;
millisecond=0;

% open savefilename for writing
fid=fopen(savefilename,'w');

%write fixed part of header
fwrite(fid,version,'int8');
fwrite(fid,numchannels,'int32');
fwrite(fid,numauxchannels,'int32');
fwrite(fid,numtimeframes,'int32');
fwrite(fid,samplingrate,'float32');
fwrite(fid,year,'int16');
fwrite(fid,month,'int16');
fwrite(fid,day,'int16');
fwrite(fid,hour,'int16');
fwrite(fid,minute,'int16');
fwrite(fid,second,'int16');
fwrite(fid,millisecond,'int16');

% define and write variable part of header
for i=1:numchannels
fwrite(fid,101,'int8');
currentchannel=uint8(num2str(i));
for j=1:size(currentchannel,2)
fwrite(fid,currentchannel(j),'int8');
end
for k=j+2:8
fwrite(fid,0,'int8');
end
end

% write data
count=fwrite(fid,thedata','float32');

% close file
fclose(fid);

> > does not lead to a page. Thanks Google Group... :-(- Masquer le texte des messages précédents -
>
> - Afficher le texte des messages précédents -

[seba]

Merci Pierre,
looking at your scripts I could make sure that mine were not so wrong
in the end ;-)
I was indeed afraid to have done something wrong, and here I explain
why:
I computed the GFP and Dissimilarity (between successive time points
of the same condition, and not between different conditions) of the
grand-average ERP (average of ALL conditions and ALL participants).
I expected to see dissimilarity peaks at the beginning of each ERP-
component (P1, N170, P2, etc.), and this is also more or less what
happened.
However, I also got huge (actually the biggest overall) dissimilarity
values during the baseline period, before stimulus onset. At this
period of the epoch the GFP is actually close to 0 (see the file
"ERP_GFP_DISS.jpg" that I just uploaded), and the variations in
topography are, I think, just random noise.
So here comes my question: is it normal to have such huge
dissimilarity values during baseline and with GFP and ERP values that
are close to 0?
Thank you for the help.
Sebastian

On Jul 15, 3:11 pm, Pierre <pierre.megev...@gmail.com> wrote:
> Hi Sebastian and everyone,
>
> Since there appears to be a problem with opening the old zip file, and
> since I'm unable to transfer any file from my HCUGE PC to this group,
> I will simply paste a few of the scripts below. I include the
> following:
> - openeph.m
> - opensef.m
> - computeavgref.m
> - computedataagainstref.m
> - computegfp.m
> - computedissandsc.m
> - saveeph.m
> - savesef.m
>
> All these scripts have a header with a short description of the
> variable names and formats used.
>
> All the caveats in my original posting remain! (http://www.google.com/
> url?sa=D&q=http://groups.google.com/group/cartool/browse_thread/thread/
> 162d69ff4e5b8ae8/4973bf8165ff4e62%3Flnk%3Dgst%26q%3Ddissimilarity
> %234973bf8165ff4e62&usg=AFQjCNHMeIfzTzPT33feWVpQbhfYne4r4A)
>
> Please signal here any problem with these scripts.
>
> Pierre (PS the e-mail address mentioned in the scripts does no longer

> exist: please use pierre.megev...@gmail.com)

>
> function
> [numchannels,numtimeframes,samplingrate,thedata]=openeph(openfilename)
> % openeph: opens a Cartool evoked potential data file (.ep(h))
> %
> % inputs: full path and name of the file to open
> %
> % outputs: number of channels, number of timeframes (and sampling rate
> for
> % .eph files) as 1-D numeric arrays; data as a 2-D numeric array where
> % dimension 1 contains the timeframes, dimension 2 contains the
> channels
> %
> % Cartool:http://brainmapping.unige.ch/Cartool.htm
> %

> % author of this script: pierre.megev...@medecine.unige.ch

> % author of this script: pierre.megev...@medecine.unige.ch

> % author of this script: pierre.megev...@medecine.unige.ch

>
> % define number of channels and time frames
> numtimeframes=size(thedata,1);
> numchannels=size(thedata,2);
>
> % compute average reference
> avgref=zeros(numtimeframes,1);
> for i=1:numtimeframes
>     avgref(i)=sum(thedata(i,:))/numchannels;
> end
>
> function thedata=computedataagainstref(thedata,ref)
> % computedataagainstref: changes the reference of a data set
> %
> % inputs: data as a 2-D numeric array where dimension 1 contains the
> % timeframes, dimension 2 contains the channels; ref as either a
> number
> % (1-D numeric array) or the 'avgref' string for average reference
> %
> % outputs: data as a 2-D numeric array where dimension 1 contains the
> % timeframes, dimension 2 contains the channels
> %
> % Cartool:http://brainmapping.unige.ch/Cartool.htm
> %

> % author of this script: pierre.megev...@medecine.unige.ch

> % author of this script: pierre.megev...@medecine.unige.ch

>
> % define number of channels and time frames
> numtimeframes=size(thedata,1);
> numchannels=size(thedata,2);
>
> % compute global field power
> gfp=zeros(numtimeframes,1);
> avgref=computeavgref(thedata);
> for i=1:numtimeframes
>     gfp(i)=sqrt(sum((thedata(i,:)-avgref(i)).^2)/numchannels);
> end
>
> function [diss,sc]=computedissandsc(thedata1,thedata2)
> % computediss: computes dissimilarity and spatial correlation between
> 2
> % data sets
> %
> % inputs: data as 2-D numeric arrays where dimension 1 contains the
> % timeframes, dimension 2 contains the channels
> %
> % outputs: diss and sc as 1-D numeric arrays that contain the
> % dissimilarity and spatial correlation at each timeframe
> %
> % Cartool:http://brainmapping.unige.ch/Cartool.htm
> %

> % author of this script: pierre.megev...@medecine.unige.ch

> % author of this script: pierre.megev...@medecine.unige.ch

> % author of this script: pierre.megev...@medecine.unige.ch

[Pierre]

Hello Sebastian,

From what I remember from looking at mouse evoked potential data
(Cartool can do all sorts :-) it is expected that the timeframe-to-
timeframe DISS value be high when the correlation of the signal at
each electrode between one timeframe and the next is low. This is
indeed the case for the baseline (essentially noise) in grand average
ERPs.

A quick tip: to get the timepoint-to-timepoint dissimilarity of the
grand average, simply open it in Cartool! Using the "add computed
tracks" button (http://brainmapping.unige.ch/cartool/refguide/
eeg.htm#more%20computed%20track) on the command list, you can add (in
addition to the GFP) the DISS curve that will be computed between each
time point and the next. You can then compare it with that you obtain
from matlab.

Pierre

> ...
>
> plus de détails »- Masquer le texte des messages précédents -

[Denis]

Pierre is totally right with the dissimillarity behavior during the
baseline period.
Dissimilarity is the topographical difference between 2 maps (here
successive time frames), once the maps are normalized (by their own
power). As weak ERP maps have lot of noise, but are not mathematically
nulls, the dissimilarity can be computed, but then it will peak a lot.
But this is not really relevant (noise)! That's why all measures have
to be put in perspective with the help of the others (Diss to Gfp
etc...).

And the Diss is normally always displayed by default on the tracks
display.
http://brainmapping.unige.ch/cartool/refguide/terms.htm#tracks

You can also export the Diss track directly (though it's a good
practice to check with your own script!):
http://brainmapping.unige.ch/cartool/refguide/export%20tracks.htm#technical%20points
(last point)


And a last thought:
the peaks of Diss don't match exactly the temporal segments produced
by the segmentation process. The latter is a much more refined method
based on mathematical clustering, you can have a switch of component
wihtout a real significant peak of diss.