Getting started with Biosemi BDF data

Biosemi makes EEG amplifiers for EEG that have active electrodes, i.e the signal is already pre-amplified at the scalp before it is sent to the amplifier box where it is further amplified and digitized. The active electrodes make the recorded signals less sensitive for environmental noise (e.g. electronics noise in the lab) and for movements of the electrode cables. These Biosemi amplifiers are especially popular in applications with high channel density.

The Biosemi system has a few special characteristics

  • it uses a non-standard referencing scheme (with the DRL and CMS electrode, rather than GND and REF)
  • it stores the data in a 24 bit format
  • it uses a sampling rate that is much higher (from 2 kHz upwards) than most EEG systems (commonly around 512 Hz)
  • the active electrodes are only useable with electrode caps from Biosemi.

The 24 bit file format has the practical consequence that files are slightly smaller than they would have been when stored with 32 bits, but also that reading and converting the 24 bit numerical representation is very slow because 24 bit is not a standard numerical representation on Intel computers. MATLAB allows to read single bits or 8-bit values from a file, which can be used to construct the 24 bit value, but which is very slow. To speed up the reading, FieldTrip uses a mex file that reads the 24 bit values and converts them to a 32-bit representation on the fly.

The high sampling rate (minimally 2kHz, i.e. 2000 Hz) has the consequence that files are much larger than with most acquisition systems. E.g. when compared to a BrainAmp data file sampled at 512 Hz with the same number of electrodes, the Biosemi data file will be approximately 3x as large on disk and 4x as large after having read it in memory. That means that for processing bdf files you typically will want to have a computer with more than the standard amount of RAM. After reading the data in MATLAB memory, a common procedure is to downsample it to reduce the sampling rate to 500 Hz using ft_resampledata. This will make all subsequent analyses run much faster and will facilitate doing the analysis with less RAM.

Most Biosemi electrode caps have a unique channel naming scheme. Also the exact number and positions is different from those in other EEG systems. Consequently, when plotting the channel-level data with a topographic arrangement of the channels, or when plotting the topographies (see the plotting tutorial and the channel layout tutorial), you will have to use a layout that is specific to your Biosemi electrode cap. FieldTrip includes the following template 2D layout files in the fieldtrip/template/layout directory, but you might want to construct your own layout.

  • biosemi16.lay
  • biosemi32.lay
  • biosemi64.lay
  • biosemi128.lay
  • biosemi160.lay
  • biosemi256.lay

Please note that the use of these layout files requires that the channel labeling in the data is consistent with the channel labeling in the layout file. You might want to check the following

cfg = [];
cfg.layout = 'biosemi160.lay'
ft_layoutplot(cfg)

Example analysis pipeline for Biosemi bdf data

The following is an example analysis pipeline that was used for the FieldTrip workshop at CUNY, New York by Stephen and Saskia in 2011.

Modify layout file

% load biosemi160 into q (which has sensor locations in polar coordinates
% add fiducials then save

% cfg=[];
% cfg.layout='biosemi160.lay';
% q=ft_prepare_layout(cfg);

zeroline = 92/40*50;

q{161,1} = 'nasion';
q{161,2} = zeroline; %ten percent under FPz
q{161,3} = q{104,3};

q{162,1} = 'inion';
q{162,2} = zeroline; %ten percent under Oz
q{162,3} = q{23,3};

q{163,1} = 'left';
q{163,2} = -zeroline; %ten percent under T7
q{163,3} = q{137,3};

q{164,1} = 'right';
q{164,2} = zeroline; %ten percent under T8
q{164,3} = q{64,3};

save('LayoutNY_fiducials',q);



% construct 3D electrode positions

load LayoutNY_fiducials.mat

ph = cell2mat(q(:,2));
th = cell2mat(q(:,3));

x = sin(ph*pi/180) .* cos(th*pi/180);
y = sin(ph*pi/180) .* sin(th*pi/180);
z = cos(ph*pi/180);

plot3(x, y, z, '.');
elec.label = q(:,1);
elec.pnt = [x y z];

% scale to get into mm
elec.pnt = 100*elec.pnt;

vol = ft_read_vol('headmodel/standard_vol.mat');

%realign electrodes to headmodel
cfg = [];
cfg.method = 'interactive';
cfg.elec = elec;
cfg.headshape = vol.bnd(1); %1 = skin
elec = ft_electroderealign(cfg);

%save new electrodes
save elec160.mat elec

Preprocessing

clear all

cond = {'Pitch20' 'Pitch40', 'Timbre490', 'Timbre510'};

for c = 1:4 % loop over the 4 conditions
  for block = 1:4 % loop over the 4 blocks within each condition
    
    % create the trial definition
    cfg=[];
    cfg.filename              = ['data/RON_', cond{c}, '_Block', num2str(block), '_Sub03.bdf'];
    cfg.trialfun              = 'trialfunNY';
    cfg.trialdef.eventtype    = 'STATUS';
    cfg.trialdef.eventvalue   = [121 122 103 104 111 112]; % stimulus triggers (standard before (2), deviant (2), standard after(2))
    cfg.trialdef.eventcorrect = [1   2   1   2   1   2  ]; % correct response triggers
    cfg.trialdef.prestim      = 0.2; % latency in seconds
    cfg.trialdef.poststim     = 1;   % latency in seconds
    
    cfg = ft_definetrial(cfg);
    
    % remember condition and block
    trl = cfg.trl;
    trl(:,6) = c;
    trl(:,8) = block;
    
    % define the relevant trials
    for i=1:length(trl)
      if ismember(trl(i,4), [121 122])
        trl(i,7) = 1; % standard before
      elseif ismember(trl(i,4), [103 104])
        trl(i,7) = 2; % deviant
      elseif ismember(trl(i,4), [111 112])
        trl(i,7) = 3; % standard after
      end
    end
    
    % read and preprocess the data using the trial definition
    cfg=[];
    cfg.dataset       = ['data/RON_', cond{c}, '_Block', num2str(block), '_Sub03.bdf'];
    cfg.trl           = trl;
    cfg.reref         = 'yes';
    cfg.refchannel    = 'EXG5';
    
    data = ft_preprocessing(cfg);
    
    % save the data to disk
    save(['analysis/prepro/data_', cond{c}, '_block', num2str(block)], 'data')
  end
end

Artifact rejection

clear all
cond = {'Pitch20' 'Pitch40', 'Timbre490', 'Timbre510'};

% concatenate all trials/blocks
ipart=1;
for icond = 1:4
  for iblock = 1:4
    disp(['Loading analysis/prepro/data_', cond{icond}, '_block', num2str(iblock)]);
    load(['analysis/prepro/data_', cond{icond}, '_block', num2str(iblock)], 'data');
    datapart(ipart) = data;
    ipart=ipart+1;
  end
end

cfg=[];
data = ft_appenddata(cfg, datapart(1), datapart(2),datapart(3),datapart(4),datapart(5),datapart(6),datapart(7),datapart(8),datapart(9),datapart(10),datapart(11),datapart(12),datapart(13),datapart(14),datapart(15),datapart(16));
clear datapart*

% detect eog artifacts using ICA
cfg=[];
cfg.method  = 'runica';
cfg.channel = 1:160; % EEG channels only
datacomp = ft_componentanalysis(cfg, data);
save('analysis/ica/datacomp', 'datacomp')

% plot the components to detect the artifacts
figure
k=1; f=1;
for icomp=1:length(datacomp.topo)
  if k>20
    k=1;
    figure
  end
  cfg=[];
  cfg.layout = 'biosemi160.lay';
  cfg.xlim   = [icomp icomp];
  
  subplot(4,5,k);
  ft_topoplotER(cfg, datacomp);
  title(icomp);
  
  k = k+1;
end


% remove components that reflect eog artifacts
cfg=[];
cfg.component = [12 49]; % note the exact numbers will vary per run
data = ft_rejectcomponent(cfg, datacomp);
save('analysis/data_clean', 'data')



%% manual artifact rejection

% remove the mean
cfg=[];
cfg.demean = 'yes';
data = ft_preprocessing(cfg, data);

% shuffle the trials for non-biased artifact rejection
shuffle = randperm(length(data.trial));
datashuff = data;
for i=1:length(data.trial)
  datashuff.trial{i} = data.trial{shuffle(i)};
  datashuff.time{i}  = data.time{shuffle(i)};
  datashuff.trialinfo(i,:) = data.trialinfo(shuffle(i),:);
end

% browse for artifact rejection
cfg = [];
cfg.channel = 'EEG';
cfg.continuous = 'no';
cfg = ft_databrowser(cfg, data);
data = ft_rejectartifact(cfg,data);

% visual artifact rejection in summary mode
cfg = [];
cfg.method = 'summary';
data = ft_rejectvisual(cfg, data);

% the following trials were removed: 43, 87, 174, 261, 307, 333, 343, 398, 557, 583, 593, 606, 666, 670, 674, 741, 742, 743, 744, 745, 899
% the following channels were removed: A7, A13, A14, A25, A26, A27, A28, B8, B9, B12, B19, B32, D1, D21, D23, D25, D26, D31, D32, E1, E17, E21, E32


save('analysis/data_clean', 'data')

Timelock Analysis

clear all

cond = {'Pitch20' 'Pitch40', 'Timbre490', 'Timbre510'};

load('analysis/data_clean', 'data')

for icond = 1:4 % loop over the 4 conditions    
  for istim = 1:3 % loop over standard-before, deviant, standard-after
      
    % compute timelocked averages for each condition
    cfg=[];
    cfg.lpfilter   = 'yes';
    cfg.lpfreq     = 40;
    cfg.trials   = find(data.trialinfo(:,3) == icond & data.trialinfo(:,4) == istim);
    timelock{icond,istim} = ft_timelockanalysis(cfg, data);
      
    % baseline correct
    cfg=[];
    cfg.baseline = [-0.2 0];
    timelock{icond,istim} = ft_timelockbaseline(cfg, timelock{icond,istim});
  end  
end

save('analysis/timelock/timelock_all','timelock');

% plot the ERPs over all sensors
figure
for i = 1:4
  subplot(2,2,i);
  ft_singleplotER([],timelock{i,1},timelock{i,2},timelock{i,3});
  title(cond{i});
  if i==1
    legend('standard-before', 'deviant', 'standard-after', 'Location', 'NorthWest')
  end
end
print(gcf, '-dpng', 'figures/fig1_ERP')

% plot ERP in interactive mode, only for standard-before
cfg = [];
cfg.layout      = 'biosemi160lay.mat';
cfg.interactive = 'yes';
figure; ft_multiplotER(cfg,timelock{1,1});

% compute the contrasts
for icond = 1:4
  % standard before vs deviant
  stb_vs_dev{icond} = timelock{1,1};
  stb_vs_dev{icond}.avg = timelock{icond,2}.avg - timelock{icond,1}.avg;
  
  % standard before vs standard after
  stb_vs_sta{icond} = timelock{1,1};
  stb_vs_sta{icond}.avg = timelock{icond,3}.avg - timelock{icond,1 }.avg;
  
end

% plot the contrasts
figure
for icond = 1:4
  cfg=[];
  cfg.xlim        = [0.5 7];
  cfg.zlim        = 'maxabs';
  cfg.interactive = 'yes';
  cfg.layout      = 'biosemi160lay.mat';
  cfg.colorbar    = 'yes';
  
  subplot(2,2,icond);
  ft_topoplotER(cfg, stb_vs_dev{icond});
end
print(gcf, '-dpng', 'figures/fig2_ERP')

figure
for icond = 1:4
  cfg=[];
  cfg.xlim        = [0.5 7];
  cfg.zlim        = 'maxabs';
  cfg.interactive = 'yes';
  cfg.layout      = 'biosemi160lay.mat';
  cfg.colorbar    = 'yes';
  
  subplot(2,2,icond);
  ft_topoplotER(cfg, stb_vs_sta{icond});
end
print(gcf, '-dpng', 'figures/fig3_ERP')



%% now collapse

% compute averages for pitch and timbre
for istim = 1:3
  cfg=[];
  cfg.lpfilter   = 'yes';
  cfg.lpfreq     = 40;
  cfg.keeptrials = 'yes';
  cfg.trials     = find((data.trialinfo(:,3) == 1 | data.trialinfo(:,3) == 2 ) & data.trialinfo(:,4) == istim);
  timelock_pitch{istim}  = ft_timelockanalysis(cfg, data);
  
  cfg.trials     = find((data.trialinfo(:,3) == 3 | data.trialinfo(:,3) == 4 ) & data.trialinfo(:,4) == istim);
  timelock_timbre{istim} = ft_timelockanalysis(cfg, data);
  
  cfg=[];
  cfg.baseline   = [-0.2 0];
  timelock_pitch{istim}  = ft_timelockbaseline(cfg, timelock_pitch{istim});
  timelock_timbre{istim} = ft_timelockbaseline(cfg, timelock_timbre{istim});
end
save('analysis/timelock/timelock_avg','timelock_*');


% plot contrasts
figure
subplot(1,2,1);
ft_singleplotER([],timelock_pitch{1},timelock_pitch{2},timelock_pitch{3});
title('Pitch');

subplot(1,2,2);
ft_singleplotER([],timelock_timbre{1},timelock_timbre{2},timelock_timbre{3});
title('Timbre');

print(gcf, '-dpng', 'figures/fig4_ERP')

Statistics

cfg=[];
cfg.layout = 'biosemi160lay.mat'; %in meters
%cfg.layout        = 'elec160.mat'; %in mm
cfg.neighbourdist = .1; 
cfg.neighbours    = ft_neighbourselection(cfg,timelock_pitch{1});

cfg.latency       = [0 1];
cfg.parameter     = 'trial';
cfg.method        = 'montecarlo';
cfg.design        = [1:size(timelock_pitch{2}.trial,1) 1:size(timelock_pitch{1}.trial,1);
  ones(1,size(timelock_pitch{2}.trial,1)), ones(1,size(timelock_pitch{1}.trial,1))*2];

cfg.numrandomization = 1000;
cfg.correctm      = 'cluster';
cfg.correcttail   = 'prob';
cfg.ivar          = 2;
cfg.uvar          = 1;

cfg.statistic     = 'indepsamplesT';

stat = ft_timelockstatistics(cfg, timelock_pitch{2}, timelock_pitch{1});

save('analysis/timelock/stat','stat');



%% plot

% cfg=[];
% cfg.layout = 'biosemi160lay.mat';
% ft_clusterplot(cfg, stat)


% find relevant clusters
ipos = find([stat.posclusters.prob]<=0.05);
ineg = find([stat.negclusters.prob]<=0.05);

% loop over all sig positive clusters
for i=ipos
  
  cfg=[];
  cfg.highlight = 'on';
  cfg.zparam    = 'stat';
  cfg.layout    = 'biosemi160lay.mat';
  cfg.style     = 'straight';
  cfg.gridscale = 500;
  
  % find the significant time range for this cluster
  tmp=[];
  for t = 1:length(stat.time)
    if ~isempty(find(any(stat.posclusterslabelmat(:,t)==ipos)))
      tmp = [tmp t];
    end
  end
  cfg.xlim      = [stat.time(tmp(1)) stat.time(tmp(end))];
  
  % find the channels belonging to this cluster
  cfg.highlightchannel = [];
  for c = 1:length(stat.label)
    if ~isempty(find(any(stat.posclusterslabelmat(c,:)==ipos)))
      cfg.highlightchannel = [cfg.highlightchannel c];
    end
  end
  
  figure
  ft_topoplotER(cfg, stat);
  title('positive cluster')
  print(gcf, '-dpng', ['figures/fig5_STAT_pos', num2str(i)])
end

% loop over all sig negative clusters
for i=ineg
  
  cfg=[];
  cfg.highlight = 'on';
  cfg.zparam    = 'stat';
  cfg.layout    = 'biosemi160lay.mat';
  cfg.style     = 'straight';
  cfg.gridscale = 500;
  
  % find the significant time range for this cluster
  tmp=[];
  for t = 1:length(stat.time)
    if ~isempty(find(any(stat.negclusterslabelmat(:,t)==ineg)))
      tmp = [tmp t];
    end
  end
  cfg.xlim      = [stat.time(tmp(1)) stat.time(tmp(end))];
  
  % find the channels belonging to this cluster
  cfg.highlightchannel = [];
  for c = 1:length(stat.label)
    if ~isempty(find(any(stat.negclusterslabelmat(c,:)==ineg)))
      cfg.highlightchannel = [cfg.highlightchannel c];
    end
  end
  
  figure
  ft_topoplotER(cfg, stat);
  title('negative cluster')
  print(gcf, '-dpng', ['figures/fig6_STAT_neg', num2str(i)])
end



load analysis/timelock/timelock_all.mat
load headmodel/standard_mri.mat

elec = ft_read_sens('elec160.mat');
vol  = ft_read_vol('headmodel/standard_vol.mat');

  vol.r = [86 88 92 100];
  vol.o = [0 0 40];
  figure, ft_plot_vol(vol)
ft_plot_vol([],vol);

% prepare data
cfg = [];
cfg.demean = 'yes';
cfg.baselinewindow = [-inf 0];
cfg.reref = 'yes';
cfg.refchannel = 'all'; % thereby averaging systematic leadfield error over sensors
for icond = 1 : 4
    for istim = 1 : 3
        timelock{icond,istim} = ft_preprocessing(cfg, timelock{icond,istim});
    end
end

% do some averages
before = timelock{1,1};
before.trial{1} = (timelock{1,1}.trial{1} + timelock{2,1}.trial{1} + timelock{3,1}.trial{1} + timelock{4,1}.trial{1} ) ./ 4;
deviant = timelock{1,1};
deviant.trial{1} = (timelock{1,2}.trial{1} + timelock{2,2}.trial{1} + timelock{3,2}.trial{1} + timelock{4,2}.trial{1} ) ./ 4;
after = timelock{1,1};
after.trial{1} = (timelock{1,3}.trial{1} + timelock{2,3}.trial{1} + timelock{3,3}.trial{1} + timelock{4,3}.trial{1} ) ./ 4;

% plot ERP's
cfg = [];
cfg.layout = 'biosemi160lay.mat';
cfg.interactive = 'yes';
figure; ft_multiplotER(cfg, before, deviant, after);
legend;

% dipole fitting
cfg = [];
cfg.vol = vol;
cfg.elec = elec;
cfg.model = 'regional';
cfg.numdipoles = 2;
cfg.latency = [0.080 0.1]; %n100
cfg.gridsearch = 'yes';
cfg.resolution = 20; % mm
cfg.nonlinear = 'yes';
cfg.symmetry = 'x';
source = ft_dipolefitting(cfg, before);

% topoplot ERP
cfg = [];
cfg.layout = 'biosemi160lay.mat';
cfg.zparam = 'Vdata';
figure; ft_topoplotER(cfg, source);

% topoplot dipole model
cfg.zparam = 'Vmodel';
figure; ft_topoplotER(cfg, source);

% topoplot difference data with model
source.Vdifference = source.Vdata ./ source.Vmodel;
cfg.zparam = 'Vdifference';
cfg.zlim = [0 2];
figure; ft_topoplotER(cfg, source);

% plot first dipole position on mni
cfg = [];
cfg.method = 'ortho';
cfg.location = source.dip.pos(1,:);
figure; ft_sourceplot(cfg, mri);

% plot second dipole position on mni
cfg.location = source.dip.pos(2,:);
figure; ft_sourceplot(cfg, mri);

% use dipole positions to extract timecourse of whole trial
cfg = source.cfg;
cfg.dip.pos = source.dip.pos;
cfg.gridsearch = 'no';
cfg.nonlinear = 'no';
cfg.latency = [-inf inf];
cfg.vol = vol;
cfg.elec = elec;
source2 = ft_dipolefitting(cfg, before);

% plot timecourses of both dipoles
figure; hold; 
plot(source2.time, source2.dip.mom(1:3,:), '-')
plot(source2.time, source2.dip.mom(4:6,:), ':')

% use Singular Value Decomposition to extract component from all
% three orientations of first dipole. u = identity matrix, s = scaling, v =
% signal strength
[u1, s1, v1] = svd(source2.dip.mom(1:3,:));
[u2, s2, v2] = svd(source2.dip.mom(4:6,:));
figure; hold;
plot(source2.time, v1(:,1)); 
plot(source2.time, v2(:,1),':'); 

% use pythagoras to make one dimensional dipole from three orientations
figure; hold;
plot(source2.time, sqrt(sum(source2.dip.mom(1:3,:).^2)));
plot(source2.time, sqrt(sum(source2.dip.mom(4:6,:).^2)),':');

Appendix: the BDF fileformat

BDF is a 24 bit version of the popular 16 bit EDF format, which was used on previous BioSemi models with 16 bit converters. BDF is almost the same as EDF Although initially the EDF format was mainly used in sleep research, BDF/EDF is now quickly gaining popularity in other EEG applications, ECG body surface potential mapping as well as EMG.

The EDF format was designed and published in 1992 by Bob Kemp, Alpo Värri, Agostinho C. Rosa, Kim D. Nielsen and John Gade as “A simple format for exchange of digitized polygraphic recordings” in Electroencephalography and Clinical Neurophysiology, 82 (1992) 391-393. The original EDF specifications can be found here

Each BDF/EDF file starts with a header followed by the number of Data records indicated in the header.

Length in bytes BDF Header EDF Header Description
8 bytes Byte 1: “255” (non ascii) Byte 1: “0” (ASCII) Identification code
Bytes 2-8 : “BIOSEMI” (ASCII) Bytes 2-8 : ” ”(ASCII)
80 bytes User text input (ASCII) Local subject identification
80 bytes User text input (ASCII) Local recording identification
8 bytes dd.mm.yy (ASCII) Startdate of recording
8 bytes hh.mm.ss (ASCII) Starttime of recording
8 bytes (ASCII) Number of bytes in header record
44 bytes “24BIT” (ASCII) “BIOSEMI” (ASCII) Version of data format
8 bytes (ASCII) Number of data records ”-1” if unknown
8 bytes e.g.: “1” (ASCII) Duration of a data record, in seconds
4 bytes e.g.: “257” or “128”(ASCII) Number of channels (N) in data record
N x 16 bytes e.g.: “Fp1”, “Fpz”, “Fp2”, etc (ASCII) Labels of the channels
N x 80 bytes e.g.: “active electrode”, “respiration belt” (ASCII) Transducer type
N x 8 bytes e.g.: “uV”, “Ohm” (ASCII) Physical dimension of channels
N x 8 bytes e.g.: ”-262144” (ASCII) e.g.: ”-32768” (ASCII) Physical minimum in units of physical dimension
N x 8 bytes e.g.: “262143” (ASCII) e.g.: “32767” (ASCII) Physical maximum in units of physical dimension
N x 8 bytes e.g.: ”-8388608” (ASCII) e.g.: ”-32768” (ASCII) Digital minimum
N x 8 bytes e.g.: “8388607” (ASCII) e.g.: “32767” (ASCII) Digital maximum
N x 80 bytes e.g.: “HP:DC; LP:410” e.g.: “HP:0,16; LP:500” Prefiltering
N x 8 bytes For example: “2048” (ASCII) Number of samples in each data record (Sample-rate if Duration of data record = “1”)
N x 32 bytes (ASCII) Reserved

Bold Text Number of samples in each data record (Sample-rate if Duration of data record = “1”) N x 32 bytes (ASCII)

Total header length (for BDF and EDF) is: {(N+1)*256} bytes, where N is number of channels (including the status channel). The “gain” of a specific channel can be calculated by: (Physical max - Physical min) / (Digital max - Digital min). The result is the LSB value in the specified Physical dimension of channels. (31,25nV / 1uV in the BDF/EDF example Header from above). The last 10 fields are defined for each fields separately. Each channel can be different.

getting_started/bdf.txt · Last modified: 2013/11/12 20:36 by 37.251.53.205

You are here: startgetting_startedbdf
This DokuWiki features an Anymorphic Webdesign theme, customised by Eelke Spaak and Stephen Whitmarsh.
www.chimeric.de Valid CSS Driven by DokuWiki do yourself a favour and use a real browser - get firefox!! Recent changes RSS feed Valid XHTML 1.0