fiff_read_evoked

 [data] = fiff_read_evoked(fname,setno)

 Read one evoked data set

0001 function [data] = fiff_read_evoked(fname,setno)
0002 %
0003 % [data] = fiff_read_evoked(fname,setno)
0004 %
0005 % Read one evoked data set
0006 %
0007 
0008 
0009 %
0010 %   Author : Matti Hamalainen, MGH Martinos Center
0011 %   License : BSD 3-clause
0012 %
0013 %
0014 %   Revision 1.11  2009/03/30 11:37:37  msh
0015 %   Added copying of measurement info blocks from the original like in mne_browse_raw
0016 %
0017 %   Revision 1.10  2008/03/28 13:42:01  msh
0018 %   Accommodate evoked files from xplotter, which may have channel information
0019 %   local to the evoked block
0020 %
0021 %   Revision 1.9  2007/11/13 10:55:32  msh
0022 %   Specify the second argument to all calls to the exist function
0023 %
0024 %   Revision 1.8  2006/10/14 11:47:32  msh
0025 %   Omitted listing of nepoch from the diagnostic output
0026 %
0027 %   Revision 1.7  2006/09/21 19:27:16  msh
0028 %   fiff_read_evoked did not tolerate missing comment
0029 %
0030 %   Revision 1.6  2006/05/03 19:03:19  msh
0031 %   Eliminated the use of cast function for Matlab 6.5 compatibility
0032 %
0033 %   Revision 1.5  2006/04/23 15:29:40  msh
0034 %   Added MGH to the copyright
0035 %
0036 %   Revision 1.4  2006/04/18 20:44:46  msh
0037 %   Added reading of forward solution.
0038 %   Use length instead of size when appropriate
0039 %
0040 %   Revision 1.3  2006/04/12 13:13:51  msh
0041 %   Added fiff_find_evoked.m
0042 %   Use aspect_kind field name instead of aspect_type
0043 %
0044 %   Revision 1.2  2006/04/12 10:29:02  msh
0045 %   Made evoked data writing compatible with the structures returned in reading.
0046 %
0047 %   Revision 1.1  2006/04/10 23:26:54  msh
0048 %   Added fiff reading routines
0049 %
0050 %
0051 
0052 global FIFF;
0053 if isempty(FIFF)
0054     FIFF = fiff_define_constants();
0055 end
0056 
0057 me='MNE:fiff_read_evoked';
0058 
0059 if nargin == 1
0060     setno = 1;
0061 elseif nargin ~= 2
0062     error(me,'Incorrect number of arguments');
0063 end
0064 if setno <= 0
0065     error(me,'Data set selector must be positive');
0066 end
0067 %
0068 %   Open the file
0069 %
0070 fprintf(1,'Reading %s ...\n',fname);
0071 [ fid, tree ] = fiff_open(fname);
0072 %
0073 %   Read the measurement info
0074 %
0075 [ info, meas ] = fiff_read_meas_info(fid,tree);
0076 info.filename = fname;
0077 %
0078 %   Locate the data of interest
0079 %
0080 processed = fiff_dir_tree_find(meas,FIFF.FIFFB_PROCESSED_DATA);
0081 if length(processed) == 0
0082     fclose(fid);
0083     error(me,'Could not find processed data');
0084 end
0085 %
0086 evoked = fiff_dir_tree_find(meas,FIFF.FIFFB_EVOKED);
0087 if length(evoked) == 0
0088     fclose(fid);
0089     error(me,'Could not find evoked data');
0090 end
0091 %
0092 %   Identify the aspects
0093 %
0094 naspect = 0;
0095 is_smsh = [];
0096 for k = 1:length(evoked)
0097     sets(k).aspects = fiff_dir_tree_find(evoked(k),FIFF.FIFFB_ASPECT);
0098     sets(k).naspect = length(sets(k).aspects);
0099     if sets(k).naspect > 0
0100         is_smsh = [ is_smsh zeros(1,sets(k).naspect) ];
0101         naspect = naspect + sets(k).naspect;
0102     end
0103     saspects  = fiff_dir_tree_find(evoked(k), FIFF.FIFFB_SMSH_ASPECT);
0104     nsaspects = length(saspects);
0105     if nsaspects > 0
0106         sets(k).naspect = sets(k).naspect + nsaspects;
0107         sets(k).aspects = [ sets(k).aspects saspects ];
0108         is_smsh = [ is_smsh ones(1,sets(k).naspect) ];
0109         naspect = naspect + nsaspects;
0110     end
0111 end
0112 fprintf(1,'\t%d evoked data sets containing a total of %d data aspects in %s\n',length(evoked),naspect,fname);
0113 if setno > naspect || setno < 1
0114     fclose(fid);
0115     error(me,'Data set selector out of range');
0116 end
0117 %
0118 %   Next locate the evoked data set
0119 %
0120 p = 0;
0121 goon = true;
0122 for k = 1:length(evoked)
0123     for a = 1:sets(k).naspect
0124         p = p + 1;
0125         if p == setno
0126             my_evoked = evoked(k);
0127             my_aspect = sets(k).aspects(a);
0128             goon = false;
0129             break;
0130         end
0131     end
0132     if ~goon
0133         break;
0134     end
0135 end
0136 %
0137 %   The desired data should have been found but better to check
0138 %
0139 if ~exist('my_evoked','var') || ~exist('my_aspect','var')
0140     fclose(fid);
0141     error(me,'Desired data set not found');
0142 end
0143 %
0144 %   Now find the data in the evoked block
0145 %
0146 nchan = 0;
0147 sfreq = -1;
0148 q = 0;
0149 for k = 1:my_evoked.nent
0150     kind = my_evoked.dir(k).kind;
0151     pos  = my_evoked.dir(k).pos;
0152     switch kind
0153         case FIFF.FIFF_COMMENT
0154             tag = fiff_read_tag(fid,pos);
0155             comment = tag.data;
0156         case FIFF.FIFF_FIRST_SAMPLE
0157             tag = fiff_read_tag(fid,pos);
0158             first = tag.data;
0159         case FIFF.FIFF_LAST_SAMPLE
0160             tag = fiff_read_tag(fid,pos);
0161             last = tag.data;
0162         case FIFF.FIFF_NCHAN
0163             tag = fiff_read_tag(fid,pos);
0164             nchan = tag.data;
0165         case FIFF.FIFF_SFREQ
0166             tag = fiff_read_tag(fid,pos);
0167             sfreq = tag.data;
0168         case FIFF.FIFF_CH_INFO
0169             q = q+1;
0170             tag = fiff_read_tag(fid,pos);
0171             chs(q) = tag.data;
0172     end
0173 end
0174 if ~exist('comment','var')
0175     comment = 'No comment';
0176 end
0177 %
0178 %   Local channel information?
0179 %
0180 if nchan > 0
0181     if ~exist('chs','var')
0182         fclose(fid);
0183         error(me, ...
0184             'Local channel information was not found when it was expected.');
0185     end
0186     if length(chs) ~= nchan
0187         fclose(fid);
0188         error(me, ...
0189             'Number of channels and number of channel definitions are different');
0190     end
0191     info.chs   = chs;
0192     info.nchan = nchan;
0193     fprintf(1, ...
0194         '\tFound channel information in evoked data. nchan = %d\n',nchan);
0195     if sfreq > 0
0196         info.sfreq = sfreq;
0197     end
0198 end
0199 nsamp = last-first+1;
0200 fprintf(1,'\tFound the data of interest:\n');
0201 fprintf(1,'\t\tt = %10.2f ... %10.2f ms (%s)\n',...
0202     1000*double(first)/info.sfreq,1000*double(last)/info.sfreq,comment);
0203 if ~isempty(info.comps)
0204     fprintf(1,'\t\t%d CTF compensation matrices available\n',length(info.comps));
0205 end
0206 %
0207 % Read the data in the aspect block
0208 %
0209 nepoch = 0;
0210 for k = 1:my_aspect.nent
0211     kind = my_aspect.dir(k).kind;
0212     pos  = my_aspect.dir(k).pos;
0213     switch kind
0214         case FIFF.FIFF_COMMENT
0215             tag = fiff_read_tag(fid,pos);
0216             comment = tag.data;
0217         case FIFF.FIFF_ASPECT_KIND
0218             tag = fiff_read_tag(fid,pos);
0219             aspect_kind = tag.data;
0220         case FIFF.FIFF_NAVE
0221             tag = fiff_read_tag(fid,pos);
0222             nave = tag.data;
0223         case FIFF.FIFF_EPOCH
0224             nepoch = nepoch + 1;
0225             tag = fiff_read_tag(fid,pos);
0226             epoch(nepoch) = tag;
0227     end
0228 end
0229 if ~exist('nave','var')
0230     nave = 1;
0231 end
0232 fprintf(1,'\t\tnave = %d aspect type = %d\n',...
0233     nave,aspect_kind);
0234 if nepoch ~= 1 && nepoch ~= info.nchan
0235     fclose(fid);
0236     error(me,'Number of epoch tags is unreasonable (nepoch = %d nchan = %d)',nepoch,info.nchan);
0237 end
0238 %
0239 if nepoch == 1
0240     %
0241     %   Only one epoch
0242     %
0243     all = epoch(1).data;
0244     %
0245     %   May need a transpose if the number of channels is one
0246     %
0247     if size(all,2) == 1 && info.nchan == 1
0248         all = all';
0249     end
0250 else
0251     %
0252     %   Put the old style epochs together
0253     %
0254     all = epoch(1).data';
0255     for k = 2:nepoch
0256         all = [ all ; epoch(k).data' ];
0257     end
0258 end
0259 if size(all,2) ~= nsamp
0260     fclose(fid);
0261     error(me,'Incorrect number of samples (%d instead of %d)',...
0262         size(all,2),nsamp);
0263 end
0264 %
0265 %   Calibrate
0266 %
0267 for k = 1:info.nchan
0268     cals(k) = info.chs(k).cal;
0269 end
0270 all = diag(cals)*all;
0271 %
0272 %   Put it all together
0273 %
0274 data.info = info;
0275 data.evoked.aspect_kind = aspect_kind;
0276 data.evoked.is_smsh     = is_smsh(setno);
0277 if exist('nave','var')
0278     data.evoked.nave = nave;
0279 else
0280     data.evoked.nave  = 1;
0281 end
0282 data.evoked.first = first;
0283 data.evoked.last  = last;
0284 if exist('comment','var')
0285     data.evoked.comment = comment;
0286 end
0287 %
0288 %   Times for convenience and the actual epoch data
0289 %
0290 data.evoked.times = double(data.evoked.first:1:data.evoked.last)/data.info.sfreq;
0291 data.evoked.epochs = all;
0292 
0293 fclose(fid);
0294 
0295 return;
0296 
0297     function [tag] = find_tag(node,findkind)
0298         
0299         for p = 1:node.nent
0300             kind = node.dir(p).kind;
0301             pos  = node.dir(p).pos;
0302             if kind == findkind
0303                 tag = fiff_read_tag(fid,pos);
0304                 return;
0305             end
0306         end
0307         tag = [];
0308         return
0309     end
0310 
0311 end
0312 
0313