/trunk/examples/toolboxes/fieldtrip/external/ctf/writeCTFds.m
MATLAB | 1616 lines | 1126 code | 163 blank | 327 comment | 175 complexity | b75e119c345beec5a724138aae97b89d MD5 | raw file
Possible License(s): BSD-3-Clause, AGPL-1.0, LGPL-2.0, GPL-3.0, GPL-2.0, LGPL-2.1, LGPL-3.0, BSD-2-Clause
Large files files are truncated, but you can click here to view the full file
- function ds=writeCTFds(datasetname,ds,data,unit);
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % %
- % This program creates datasets that can be analyzed by CTF software. %
- % %
- % Datasets created by this program MUST NOT BE USED FOR CLINICAL APPLICATIONS. %
- % %
- % Please do not redistribute it without permission from VSM MedTech Ltd. %
- % %
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % Author : Harold Wilson
- % Version 1.3 5 October 2007 Spelling errors in some variables corrected.
- % Version 1.2 24 April 2007 Modified to write both MEG and fMEG .hc files.
- % writeCTFds.m Version 1.1 Prepares a CTF-format data set. MATLAB to CTF conversion.
- % Adds record of filter changes to hist and newds files.
- % Adds no-clinical-use messages.
- % Creates multiple meg4 files when the total size of the data array
- % >536870910 elements.
- % Operation :
- % 1. User reads a data set using readCTFds and getTrial2.
- % size(data)= [SAMPLES, CHANNELS, TRIALS].
- % 2. After working on data in MATLAB, user adjusts ds structure to reflect changes.
- % (writeCTFds will adjust the number of channels, channel names and trial structure.)
- % 3. This program then creates a new CTF data set.
- % datasetname must include the complete path description.
- % If a data set with name datasetname already exists, writeCTFds will issue an error message.
- % The new directory contains files for each field of structure ds. If the field is
- % missing no file is created. If the field is empty, an empty file is created.
- % Files default.* are not created by writeCTFds.
- % 4. The following fields of structure ds.res4 are modified based on the size of array data :
- % no_samples
- % no_channels
- % no_trials.
- % If size(ds.res4.chanNames,1)<no_channels, additional channel names are created
- % as required. The additional channels are called MXT%%%. It is assumed that
- % there will be <1000 new channels.
- % Inputs : datasetname : Output dataset including path. Extension '.ds' is optional.
- % ds : Structure produced by readCTFds.
- % data : MEG data array. size(data)=[no_samples no_channels no_trials]
- % Array data may be single or double.
- % unit : Determines the unit of the SQUID and EEG signals:
- % If unit is missing or unit==[], the unit is set to 'fT'.
- % 'ft' or 'fT' : Convert to fT (MEG), uV (EEG)
- % 't' or 'T' : Convert to T (MEG), V (EEG)
- % 'phi0' : Convert to phi0 (MEG), uV (EEG)
- % 'int': Read plain integers from *.meg4-file
- % Outputs : - ds : The ds structure of the output data set.
- % - datasetout : the name of the output data set.
- % - A data set. The .hz and .hz2 subdirectories are not included.
- % Function calls
- % Included in this listing:
- % - check_senres: Does simple checks on the fields of the senres table.
- % - writeHc: Creates the new .hc file
- % - checkMrk: Checks structure ds.mrk to see if it is valid marker set.
- % - writeEEG: Creates the new .EEG file.
- % - writeBadSegments: Creates the new bad.segments file.
- % - writeClassFile: Creates the new ClassFile.cls file.
- % - writeVirtualChannels: Creates the new VirtualChannels file.
- % - updateDescriptors: Adds non-clinical-use and creation software messages to
- % infods, res4, newds and hist fields of ds.
- % - updateHLC: Adds head-coil movement information to infods.
- % - updateDateTime : Resets dattime fields of res4 and infods.
- % - updateBandwidth: Resets bandwidth of newds and infods. Adds res4 filter
- % description to ds.hist.
- % - getArrayField : Extracts one field of a structure array to make it easier to
- % manipulate.
- % - writeMarkerFile: Creates the new MarkerFile.mrk file.
- % - writeCPersist: Creates the new .acq and .infods files.
- % Other calls:
- % - writeRes4: Writes the new .res4 file.
- % Output files are opened with 'w' permission and 'ieee-be' machine format in order
- % to be compatible with the Linux acquisition and analysis software. Do not open files
- % with 'wt' permission because this will add an extra char(13) byte at the end of each
- % line of text.
- persistent printWarning bandwidthMessage
- delim=filesep;
- if nargin==0 & nargout==0 % Print a version number
- fprintf(['\twriteCTFds: Version 1.3 5 October 2007 ',...
- 'Creates v4.1 and v4.2 CTF data sets.\n',...
- '\tCall: ds=writeCTFds(datasetname,ds,data,unit);\n',...
- '\t\tdatasetname = Name of the new dataset including the path.\n',...
- '\t\tds = Structure describing the new dataset (ds.hc must be v1.2 format).\n',...
- '\t\tdata = data that will be written to the new dataset .meg4 file.\n',...
- '\t\tunit = physical units of the data.\n\n']);
- return
- end
- % Allowed 8-byte headers for res4 and meg4 files.
- res4_headers=strvcat(['MEG41RS',char(0)],['MEG42RS',char(0)]);
- meg4_headers=strvcat(['MEG41CP',char(0)],['MEG42CP',char(0)]);
- maxMEG4Size=2^31; % Maximum MEG$ file in bytes. (Limit set by Linux software)
- MAX_COILS=8; % Parameter that determines the size of the ds.res.senres structure.
- lenSensorName=32; % Channel names must be 32 characters
- printDefaultBandwidthMessage=0; % Print a message about default bandwidth?
- default_flp=0.25; % Default is flp=0.25*ds.res4.sample_rate
- clinicalUseMessage='NOT FOR CLINICAL USE';
- creatorSoftware='writeCTFds'; % Added to .infods file
- meg4ChunkSize=2^20; % Write new .meg4 file in chunks of 4*meg4ChunkSize bytes.
- DATASET_HZ_UNKNOWN=round(2^31-1); % Peculiar requirement of DataEditor as of 23 Oct. 2006
- if ~exist('clinicalUseMessage');
- clinicalUseMessage=char([]);
- end
- % Check inputs
- if nargin<3
- fprintf(['\nwriteCTFds: Must supply inputs datasetname,ds,data. ',...
- 'Only %d input arguments are present.\n\n'],nargin);
- ds=-1; % Force an error in the calling program.
- return
- end
- % Check input argument unit. Convert unit to lower case.
- if exist('unit')~=1
- unit='ft'; % default
- elseif isempty(unit)
- unit='ft'; % default
- elseif ischar(unit)
- unit=lower(unit);
- if ~strcmp(unit,'int') & ~strcmp(unit,'ft') & ~strcmp(unit,'t') & ~strcmp(unit,'phi0')
- fprintf(['\nwriteCTFds : unit=%s Not a valid option. Must be ',...
- '''fT'', ''T'', ''phi0'' or ''int''\n\n'],unit);
- ds=-1; % Force an error in the calling program.
- return
- end
- end
- % Check argument type
- if ~isstruct(ds) | ~isnumeric(data) | ~ischar(unit) | ~ischar(datasetname)
- fprintf('\nwriteCTFds: Some of the inputs are the wrong type.\n');
- whos datasetname ds data unit;
- ds=-1;
- return
- elseif ~isfield(ds,'res4') | ~isfield(ds,'meg4')
- fprintf('\nwriteCTFds: Fields res4 and meg4 must be present in structure ds.\n\n');
- ds % List the fields of structure ds.
- ds=-1; % Force an error in the calling program.
- return
- end
- % Refuse to write a data set with balanced reference gradiometers.
- balancedGref=0;
- for k=1:ds.res4.no_channels
- balancedGref=(ds.res4.senres(k).sensorTypeIndex==1 & ds.res4.senres(k).grad_order_no~=0);
- end
- if balancedGref
- fprintf('\nwriteCTFds: ds.res4.senres indicates balanced reference gradiometers.\n\n');
- ds=-1; % Force an error in the calling program.
- return
- end
- clear k balancedGref;
- % Separate datasetname into a path and the baseName
- datasetname=deblank(datasetname);
- ksep=max([0 findstr(datasetname,delim)]);
- baseName=datasetname((ksep+1):length(datasetname));
- path=datasetname(1:ksep); % String path is terminated by the file delimiter (or path=[]).
- % Remove the last .ds from baseName.
- kdot=max(findstr(baseName,'.ds'));
- if kdot==(length(baseName)-2)
- baseName=baseName(1:(max(kdot)-1));
- else
- datasetname=[datasetname,'.ds'];
- end
- clear ksep kdot;
- % Save the name already in structure ds, and change to the new datset name.
- if isfield(ds,'baseName')
- olddatasetname=[ds.baseName,'.ds'];
- if isfield(ds.path)
- olddatasetname=[ds.path,olddatasetname];
- end
- else
- olddatasetname=' None ';
- end
- ds.path=path;
- ds.baseName=baseName;
- % Does the dataset already exist?
- if exist(datasetname)==7
- fprintf('\nwriteCTFds: Dataset %s already exists. Use a different name.\n\n',...
- datasetname);
- ds=-1; % Force an error in the calling program.
- return
- end
- if size(ds.res4.chanNames,2)~=lenSensorName
- fprintf(['\nwriteCTFds : size(ds.res4.chanNames)=[%d %d] ? Must ',...
- 'have %d-character channel names.\n\n'],size(ds.res4.chanNames),lenSensorName);
- ds=-1;
- return
- end
- % Check that the channel names have a sensor-file identification extensions.
- % If it is missing, print a warning message.
- % Sensor type indices : SQUIDs (0:7), ADCs (10), DACs(14), Clock (17), HLC (13,28,29)
- % See Document CTF MEG File Formats (PN 900-0088), RES4 File Format/
- for index=[0:7 10 13 14 17 28 29]
- for k=find([ds.res4.senres.sensorTypeIndex]==index);
- if isempty(strfind(ds.res4.chanNames(k,:),'-'))
- fprintf(['writeCTFds: Channel %3d %s No sensor-file identification.',...
- ' (''-xxxx'' appended to channel name).\n',...
- '\t\tSome CTF software may not work with these channel names.\n'],...
- k,deblank(ds.res4.chanNames(k,:)));
- break;
- end
- end
- if isempty(strfind(ds.res4.chanNames(k,:),'-'));break;end
- end
- clear index k chanName;
- % Update the data description in the ds.res4 structure.
- [nSample, nChan, trials]=size(data);
- % Update ds.res4 fields to match the size of array data.
- ds.res4.no_trials=trials;
- ds.res4.no_channels=nChan;
- ds.res4.no_samples=nSample;
- ds.res4.epoch_time=nSample*trials/ds.res4.sample_rate;
- % Check if channels have been added or removed from array data.
- [no_chanNames len_chanNames]=size(ds.res4.chanNames);
- if no_chanNames<nChan
- % Assume that new channel are at the end of the data set. Add a fake extension.
- for kx=1:(nChan-no_chanNames)
- ds.res4.chanNames=...
- strvcat(ds.res4.chanNames,['MXT' num2str(kx,'%3.3d') '-0001' char(0)]);
- end
- fprintf('\tAdded %d SQUID channels to the end of ds.res4.chanNames table.\n',...
- nChan-no_chanNames);
- elseif no_chanNames>nChan
- fprintf(['\nlength(chanNames)=%d, but only %d channels of data. ',...
- 'writeCTFds cannot tell which channel names to remove.\n\n'],no_chanNames,nChan);
- ds=-1;
- return
- end
- clear no_chanNames len_chanNames;
- % The senres table may have been changed, especially if channels are added or removed.
- % Check structure senres, print error messages, and possibly fix the errors.
- [ds.res4.senres,status]=check_senres(ds.res4.senres,ds.res4.no_channels);
- if status<0;ds=-1;return;end
- clear status;
- % Check that ds.res4.numcoef is the size of structure array ds.res4.scrr.
- if isfield(ds.res4,'scrr')
- if ~isequal(size(ds.res4.scrr,2),ds.res4.numcoef)
- fprintf('Error in ds.res4: ds.res4.numcoef=%d, but size(ds.res4.scrr)=[',...
- ds.res4.numcoef);
- fprintf(' %d',size(ds.res4.scrr));fprintf('] ?\n');
- return
- end
- elseif ds.res4.numcoef~=0
- fprintf(['Error in ds.res4: ds.res4.numcoef=%d,',...
- ' but scrr is not a field of ds.res4\n'],ds.res4.numcoef);
- return
- end
- % Before converting data to integers, save HLC channels for motion analysis in function
- % make_new_infods. Pass HLCdata to function make_new_infods
- if isempty(strmatch('HLC',ds.res4.chanNames))
- HLCdata=[];
- else
- % Make a list of head-coil channels
- coil=0;
- HLClist=[];
- while ~isempty(strmatch(['HLC00',int2str(coil+1)],ds.res4.chanNames))
- coil=coil+1;
- for k=1:3
- HLClist=[HLClist strmatch(['HLC00',int2str(coil),int2str(k)],ds.res4.chanNames)];
- end
- end
- HLCdata=reshape(double(data(:,HLClist,:)),ds.res4.no_samples,3,coil,ds.res4.no_trials);
- clear coil k HLClist;
- end
- % Convert data to integers because CTF data sets are stored as raw numbers and
- % not as physical qunatities. The res4 file contains the calibrations for
- % converting back to physical units. Array data may be single precision, so
- % convert to double before doing any adjustments to the data.
- % Already checked that unit is valid.
- if strcmp(unit,'int')
- data=reshape(data,nSample,nChan*trials);
- for k=1:nChan*trials
- if strcmp(class(data),'single')
- data(:,k)=single(round(double(data(:,k))));
- else
- data(:,k)=round(double(data(:,k)));
- end
- end
- data=reshape(data,nSample,nChan,trials);
- clear k;
- else
- for chan=1:nChan % Convert EEGs from uV to V, SQUIDs from fT to T
- SQUIDtype=any(ds.res4.senres(chan).sensorTypeIndex==[0:7]);
- EEGtype=any(ds.res4.senres(chan).sensorTypeIndex==[8 9]);
- if EEGtype & (strcmp(unit,'ft') | strtcmp(unit,'phi0'))
- alphaG=1e-6;
- elseif SQUIDtype & strcmp(unit,'ft')
- alphaG=1e-15;
- elseif SQUIDtype & strcmp(unit,'phi0')
- alphaG=1./(ds.res4.senres(chan).properGain*ds.res4.senres(chan).ioGain);
- else
- alphaG=1;
- end
- % Convert from physical units to integers using the gains in the senres table.
- for kt=1:trials
- buff=round(double(data(:,chan,kt))*(alphaG*...
- (ds.res4.senres(chan).properGain*ds.res4.senres(chan).qGain*ds.res4.senres(chan).ioGain)));
- if strcmp(class(data),'single')
- data(:,chan,kt)=single(buff);
- else
- data(:,chan,kt)=buff;
- end
- end
- end
- clear chan alphaG SQUIDtype EEGtype buff kt;
- end
- % Create the output dataset
- [status,msg]=mkdir(path,[baseName '.ds']);
- if status==0 | ~isempty(msg)
- fprintf('\nwriteCTFds: Failed to create directory.\n');
- fprintf(' [status,msg]=mkdir(%s)\n',datasetname);
- fprintf(' returns status=%d, msg=%s',status,msg);fprintf('\n\n');
- ds=-1;
- return
- end
- clear msg status;
- % Write the data file (meg4 file). Check ds.meg4.header and make sure that
- % the output .meg4 file has an acceptable header.
- headerMessage=1;
- if ~isfield(ds.meg4,'header'); % If ds.meg4.header is missing, add it.
- nChar=length(deblank(meg4_headers(1,:)));
- ds.meg4.header=[meg4_headers(1,1:min(7,nChar)) char(zeros(1,7-nChar))];
- elseif isempty(strmatch(ds.meg4.header(1:7),meg4_headers(:,1:7),'exact'))
- ds.meg4.header=meg4_headers(1,1:7);
- else
- headerMessage=0;
- end
- if headerMessage;
- fprintf('writeCTFds: Set ds.meg4.header=%s\n',ds.meg4.header);
- end
- clear headerMessage;
- if isempty(printWarning)
- fprintf(['\nwriteCTFds: The data you are writing have been processed by software not\n',...
- '\tmanufactured by VSM MedTech Ltd. and that has not received marketing clearance\n',...
- '\tfor clinical applications. These data should not be later employed for clinical\n',...
- '\tand/or diagnostic purposes.\n\n']);
- printWarning=1;
- end
- % Write the meg4 file(s). If there are more than maxMEG4Size-8 bytes, then additional meg4
- % files will be created.
- % Convert data to a 1-D array
- ndata=prod(size(data));
- data=reshape(data,ndata,1);
- ptsPerTrial=nSample*nChan;
- maxPtsPerFile=ptsPerTrial*floor((maxMEG4Size-8)/(4*ptsPerTrial));
- pt=0; % Last point written to the output file(s).
- while pt<ndata
- endPt=pt+min(ndata-pt,maxPtsPerFile);
- if pt==0
- meg4Ext='.meg4';
- else
- meg4Ext=['.',int2str(floor(pt/maxPtsPerFile)),'_meg4'];
- end
- fidMeg4=fopen([path,baseName,'.ds\',baseName,meg4Ext],'w','ieee-be');
- fwrite(fidMeg4,[ds.meg4.header(1:7),char(0)],'uint8');
- while pt<endPt
- pt1=min(pt+meg4ChunkSize,endPt); % Convert to double in case data is
- fwrite(fidMeg4,double(data((pt+1):pt1)),'int32'); % is single and write in short
- pt=pt1; % pieces.
- end
- fclose(fidMeg4);
- end
- % Update the .meg4 part of structure ds.
- ds.meg4.fileSize=4*ndata+8*(1+floor(ndata/maxPtsPerFile));
- clear data pt pt1 ndata fidMeg4 ptsPerTrial maxPtsPerFile meg4Ext;
- % Add dataset names to .hist
- if ~isfield(ds,'hist');ds.hist=char([]);end
- ds.hist=[ds.hist char(10) char(10) datestr(now) ' :' char(10) ...
- ' Read into MATLAB as data set ' olddatasetname char(10) ...
- ' Rewritten by writeCTFds as data set ' datasetname char(10)];
- % If infods doesn't exist or is empty create it.
- if ~isfield(ds,'infods');ds.infods=[];end
- if isempty(ds.infods)
- ds.infods=make_dummy_infods(isfield(ds,'hc'),~isempty(HLCdata),ds.res4.sample_rate);
- end
- % Update text fields of res4,infods, newds and hist.
- ds=updateDescriptors(ds,clinicalUseMessage,creatorSoftware);
- % Add HLC data to infods
- ds.infods=updateHLC(ds.infods,HLCdata);
- % Analyze structure array ds.res4.filters to make text info for .hist file and
- % bandwidth parameters for .newds file.
- fhp=0; % High-pass cutoff assuming no filters
- flp=default_flp*ds.res4.sample_rate; % Assumed lowpass cutoff. SHOULD THIS BE CHANGED?
- if isempty(bandwidthMessage) & printDefaultBandwidthMessage
- fprintf('writeCTFds: Lowpass filter set to flp=%0.2f*sample_rate\n',default_flp);
- bandwidthMessage=1;
- end
- ds=updateBandwidth(ds,fhp,flp);
- % Update date/time fields of infods and res4
- ds=updateDateTime(ds);
- % Create the .res4 file in the output dataset.
- ds.res4=writeRes4([path,baseName,'.ds\',baseName,'.res4'],ds.res4,MAX_COILS);
- if ds.res4.numcoef<0
- fprintf('\nwriteCTFds: writeRes4 returned ds.res4.numcoef=%d (<0??)\n\n',...
- ds.res4.numcoef);
- % Kill the output dataset.
- rmdir([path,baseName,'.ds'],'s');
- ds=-1;
- return
- end
- % Create .hist file
- histfile=[path,baseName,'.ds\',baseName,'.hist'];
- fid=fopen(histfile,'w');
- fwrite(fid,ds.hist,'char');
- fclose(fid);
- % New .newds file
- if isfield(ds,'newds')
- fid=fopen([path,baseName,'.ds\',baseName,'.newds'],'w');
- fwrite(fid,ds.newds,'char');
- fclose(fid);
- end
- % New infods file.
- if isfield(ds,'infods')
- writeCPersist([path,baseName,'.ds',delim,baseName,'.infods'],ds.infods);
- end
- % new hc file
- if isfield(ds,'hc')
- ds.hc=writeHc([path,baseName,'.ds',delim,baseName,'.hc'],ds.hc,HLCdata(:,:,1));
- end
- clear HLCdata;
- % MarkerFile.mrk
- if checkMrk(ds)
- writeMarkerFile([path,baseName,'.ds',delim,'MarkerFile.mrk'],ds.mrk);
- end
- % .EEG
- if isfield(ds,'EEG')
- writeEEG([path,baseName,'.ds',delim,baseName,'.EEG'],ds.EEG);
- end
- % .acq
- if isfield(ds,'acq');
- % Check that ds.acq has the correct fields
- if isfield(ds.acq,'name') & isfield(ds.acq,'type') & isfield(ds.acq,'data')
- acqFilename=[path,baseName,'.ds',delim,baseName,'.acq'];
- writeCPersist(acqFilename,ds.acq);
- end
- end
- % bad.segments file
- if isfield(ds,'badSegments')
- writeBadSegments([path,baseName,'.ds',delim,'bad.segments'],ds.badSegments,...
- ds.res4.no_trials,ds.res4.no_samples/ds.res4.sample_rate);
- end
- % BadChannels
- if isfield(ds,'BadChannels');
- if ischar(ds.BadChannels) & ~isempty(ds.BadChannels)
- fid=fopen([path,baseName,'.ds',delim,'BadChannels'],'w','ieee-be');
- for k=1:size(ds.BadChannels,1)
- fprintf(fid,'%s\n',deblank(ds.BadChannels(k,:)));
- end
- fclose(fid);
- end
- end
- % ClassFile
- if check_cls(ds)
- writeClassFile([path,baseName,'.ds',delim,'ClassFile.cls'],ds.TrialClass);
- end
- % VirtualChannels
- if isfield(ds,'Virtual')
- writeVirtualChannels([path,baseName,'.ds',delim,'VirtualChannels'],ds.Virtual);
- end
- % processing.cfg
- if isfield(ds,'processing');
- if ischar(ds.processing)
- fid=fopen([datasetname,delim,'processing.cfg'],'w','ieee-be');
- fwrite(fid,ds.processing,'char');
- fclose(fid);
- end
- end
- % Update the data set path and name
- ds.path=path;
- ds.baseName=baseName;
- return
- % *************** End of function writeCTFds ********************************************
- % **************************************************************************************
- % **************************************************************************************
- % *************** Function check_senres **************************
- function [senres,status]=check_senres(senres,numChan);
- % A user may have augmented the senres table, so check that all the fields have the
- % correct size. This will cause errors in programs that try to compute
- % sensor response using the geometry in the senres table.
- % Does "sanity" checks on the senres table. If there are obviously incorrect entries, it
- % tries to fix them, and prints a message. If the senres table does not specify coil
- % positions, orientations or areas, set them to zero, but give them the correct array
- % size.
- newChannelType=4; % Create the fake sensors as MEG magnetometers. This way offsets can be removed.
- status=-1;
- % Does the senres table have the correct no. of channels?
- no_senres=length(senres);
- if no_senres<numChan
- % Add channels. Assume that they are gradiometers.
- ioGain=1;
- qGain=2^20;
- gain=0.3;
- properGain=1e15/(qGain*ioGain*gain); % sensor gain in phi0/fT
- for kx=(no_senres+1):numChan
- senres(kx)=struct(...
- 'sensorTypeIndex',newChannelType,'originalRunNum',0,'coilShape',0,...
- 'properGain',properGain,'qGain',qGain,'ioGain',ioGain,...
- 'ioOffset',0,'numCoils',1,...
- 'grad_order_no',0,...
- 'gain',gain,...
- 'pos0',zeros(3,1),'ori0',zeros(3,1),...
- 'area',0.1,'numturns',1,...
- 'pos',[0 0 21]','ori',zeros(3,1));
- end
- fprintf(['\tAdded %d SQUID channels to senres table. Nominal gain of each = ',...
- '%8.4f fT/step\n'],numChan-no_senres,gain);
- clear kx gain qGain ioGain properGain;
- elseif no_senres>numChan
- % Channels have been removed from the data set, but can't tell which elements of the
- % senres table to remove.
- fprintf(['length(senres)=%d, but only %d channels of data. writeCTFds can''t',...
- ' tell which channels to remove from senres table.\n'],no_senres,numChan);
- return
- end
- no_senres=length(senres);
- % Previous version of check_senres ensures that several fields had size [1 1].
- % Seems unnecessary, so removed it.
- for k=1:no_senres
- % Check the fields that define pickup loops. Force the field defining the pickup
- % loops to have the correct size. It is not clear why the EEG channels and
- % the type=13,28,29 HLC channels need numCoils=1, and area>0.
- if any(senres(k).sensorTypeIndex==[0:7]) % SQUID channels
- correct_numCoils=rem(senres(k).sensorTypeIndex,4)+1;
- elseif any(senres(k).sensorTypeIndex==[13 28 29]) % HLC channels
- correct_numCoils=1;
- elseif any(senres(k).sensorTypeIndex==[8 9]) % EEG channels
- correct_numCoils=1;
- else
- correct_numCoils=0;
- end
- if senres(k).numCoils~=correct_numCoils & any(senres(k).sensorTypeIndex==[0:7])
- fprintf('writeCTFds_test: senres(%d).sensorTypeIndex=%d but numCoils=%d??\n',...
- k,senres(k).sensorTypeIndex,senres(k).numCoils);
- fprintf(' Set numCoils=%d\n',correct_numCoils);
- senres(k).numCoils=correct_numCoils;
- end
- numCoils=senres(k).numCoils;
- if size(senres(k).pos0)~=[3 numCoils];pos0=zeros(3,numCoils);end
- if size(senres(k).ori0)~=[3 numCoils];ori0=zeros(3,numCoils);end
- if size(senres(k).pos)~=[3 numCoils];pos=zeros(3,numCoils);end
- if size(senres(k).ori)~=[3 numCoils];ori=zeros(3,numCoils);end
- if size(senres(k).numturns)~=[1 numCoils];numturns=zeros(1,numCoils);end
- if size(senres(k).area)~=[1 numCoils];area=zeros(3,numCoils);end
- end
- status=1;
- return
- % ************* End of function check_senres*********************************************
- % **************************************************************************************
- % **************************************************************************************
- % ************* function make_dummy_infods*********************************************
- function Tag=make_dummy_infods(exist_hc,exist_HLC,sample_rate);
- % If the user does not supply ds.infos, this code makes a dummy version. It has all of
- % the tags that Acq created in file Richard_SEF_20060606_04.infods.
- % ***********************************************************************************
- % ***********************************************************************************
- DATASET_HZ_UNKNOWN=round(2^31-1); % Peculiar requirement of DataEditor as of 23 Oct. 2006
- fprintf('writeCTDds (make_dummy_infods): Creating a dummy infods file with all the tags.\n');
- Tag(1)=struct('name','WS1_','type',0,'data',[]);
- Tag(length(Tag)+1)=struct('name','_PATIENT_INFO','type',2,'data',[]);
- Tag(length(Tag)+1)=struct('name','WS1_','type',0,'data',[]);
- Tag(length(Tag)+1)=struct('name','_PATIENT_UID','type',10,...
- 'data','2.16.124.113000.000000.00000000000000.000000000.00000000.0011');
- Tag(length(Tag)+1)=struct('name','_PATIENT_NAME_FIRST','type',10,'data','');
- Tag(length(Tag)+1)=struct('name','_PATIENT_NAME_MIDDLE','type',10,'data','');
- Tag(length(Tag)+1)=struct('name','_PATIENT_NAME_LAST','type',10,'data','');
- Tag(length(Tag)+1)=struct('name','_PATIENT_ID','type',10,'data','x');
- Tag(length(Tag)+1)=struct('name','_PATIENT_BIRTHDATE','type',10,'data','19500101000000');
- Tag(length(Tag)+1)=struct('name','_PATIENT_SEX','type',5,'data',2);
- Tag(length(Tag)+1)=struct('name','_PATIENT_PACS_NAME','type',10,'data','');
- Tag(length(Tag)+1)=struct('name','_PATIENT_PACS_UID','type',10,'data','');
- Tag(length(Tag)+1)=struct('name','_PATIENT_INSTITUTE','type',10,'data','');
- Tag(length(Tag)+1)=struct('name','EndOfParameters','type',-1,'data',[]);
- Tag(length(Tag)+1)=struct('name','_PROCEDURE_INFO','type',2,'data',[]);
- Tag(length(Tag)+1)=struct('name','WS1_','type',0,'data',[]);
- Tag(length(Tag)+1)=struct('name','_PROCEDURE_VERSION','type',5,'data',1);
- Tag(length(Tag)+1)=struct('name','_PROCEDURE_UID','type',10,...
- 'data','2.16.124.113000.000000.00000000000000.000000000.00000000.0041');
- Tag(length(Tag)+1)=struct('name','_PROCEDURE_ACCESSIONNUMBER','type',10,'data','0');
- Tag(length(Tag)+1)=struct('name','_PROCEDURE_TITLE','type',10,'data','0');
- Tag(length(Tag)+1)=struct('name','_PROCEDURE_SITE','type',10,'data','');
- Tag(length(Tag)+1)=struct('name','_PROCEDURE_STATUS','type',5,'data',1);
- Tag(length(Tag)+1)=struct('name','_PROCEDURE_TYPE','type',5,'data',2); % Research type
- Tag(length(Tag)+1)=struct('name','_PROCEDURE_STARTEDDATETIME','type',10,...
- 'data','20060606164306');
- Tag(length(Tag)+1)=struct('name','_PROCEDURE_CLOSEDDATETIME','type',10,...
- 'data','19000100000000');
- Tag(length(Tag)+1)=struct('name','_PROCEDURE_COMMENTS','type',10,'data','');
- Tag(length(Tag)+1)=struct('name','_PROCEDURE_LOCATION','type',10,'data','');
- Tag(length(Tag)+1)=struct('name','_PROCEDURE_ISINDB','type',5,'data',0);
- Tag(length(Tag)+1)=struct('name','EndOfParameters','type',-1,'data',[]);
- Tag(length(Tag)+1)=struct('name','_DATASET_INFO','type',2,'data',[]);
- Tag(length(Tag)+1)=struct('name','WS1_','type',0,'data',[]);
- Tag(length(Tag)+1)=struct('name','_DATASET_VERSION','type',5,'data',2);
- Tag(length(Tag)+1)=struct('name','_DATASET_UID','type',10,...
- 'data','2.16.124.113000.000000.00000000000000.000000000.00000000.0042');
- Tag(length(Tag)+1)=struct('name','_DATASET_PATIENTUID','type',10,...
- 'data','2.16.124.113000.000000.00000000000000.000000000.00000000.0011');
- Tag(length(Tag)+1)=struct('name','_DATASET_PROCEDUREUID','type',10,...
- 'data','2.16.124.113000.000000.00000000000000.000000000.00000000.0041');
- Tag(length(Tag)+1)=struct('name','_DATASET_STATUS','type',10,'data','');
- Tag(length(Tag)+1)=struct('name','_DATASET_RPFILE','type',10,'data','default.rp');
- Tag(length(Tag)+1)=struct('name','_DATASET_PROCSTEPTITLE','type',10,'data','');
- Tag(length(Tag)+1)=struct('name','_DATASET_PROCSTEPPROTOCOL','type',10,'data','');
- Tag(length(Tag)+1)=struct('name','_DATASET_PROCSTEPDESCRIPTION','type',10,...
- 'data','');
- Tag(length(Tag)+1)=struct('name','_DATASET_COLLECTIONDATETIME','type',10,...
- 'data','Unknown');
- Tag(length(Tag)+1)=struct('name','_DATASET_COLLECTIONSOFTWARE','type',10,...
- 'data','Acq ');
- Tag(length(Tag)+1)=struct('name','_DATASET_CREATORDATETIME','type',10,...
- 'data',sprintf('%d',floor(clock)));
- Tag(length(Tag)+1)=struct('name','_DATASET_CREATORSOFTWARE','type',10,...
- 'data','Acq ');
- Tag(length(Tag)+1)=struct('name','_DATASET_KEYWORDS','type',10,'data','');
- Tag(length(Tag)+1)=struct('name','_DATASET_COMMENTS','type',10,...
- 'data','Dummy infods.');
- Tag(length(Tag)+1)=struct('name','_DATASET_OPERATORNAME','type',10,'data','');
- Tag(length(Tag)+1)=struct('name','_DATASET_LASTMODIFIEDDATETIME','type',10,...
- 'data',sprintf('%d',floor(clock)));
- if exist_hc
- nominalPositions=0; % Measured
- else
- nominalPositions=1; % Nominal
- end
- Tag(length(Tag)+1)=struct('name','_DATASET_NOMINALHCPOSITIONS','type',5,...
- 'data',nominalPositions);
- Tag(length(Tag)+1)=struct('name','_DATASET_COEFSFILENAME','type',10,...
- 'data','ds.res4.scrr');
- Tag(length(Tag)+1)=struct('name','_DATASET_SENSORSFILENAME','type',10,...
- 'data','ds.res4.senres');
- %Tag(length(Tag)+1)=struct('name','_DATASET_COEFSFILENAME','type',10,...
- % 'data','/opt/ctf-5.1/hardware/M015/M015_1609.coef');
- %Tag(length(Tag)+1)=struct('name','_DATASET_SENSORSFILENAME','type',10,...
- % 'data','/opt/ctf-5.1/hardware/M015/M015_1609.sens');
- Tag(length(Tag)+1)=struct('name','_DATASET_SYSTEM','type',10,'data','DSQ-2010');
- Tag(length(Tag)+1)=struct('name','_DATASET_SYSTEMTYPE','type',10,'data','Untitled');
- Tag(length(Tag)+1)=struct('name','_DATASET_LOWERBANDWIDTH','type',4,'data',0);
- Tag(length(Tag)+1)=struct('name','_DATASET_UPPERBANDWIDTH','type',4,'data',...
- round(0.25*sample_rate));
- Tag(length(Tag)+1)=struct('name','_DATASET_ISINDB','type',5,'data',0);
- if exist_HLC
- HZ_MODE=5;
- elseif exist_hc
- HZ_MODE=1;
- else
- HZ_MODE=DATASET_HZ_UNKNOWN;
- end
- Tag(length(Tag)+1)=struct('name','_DATASET_HZ_MODE','type',5,'data',HZ_MODE);
- Tag(length(Tag)+1)=struct('name','_DATASET_MOTIONTOLERANCE','type',4,'data',0.005);
- Tag(length(Tag)+1)=struct('name','_DATASET_MAXHEADMOTION','type',4,'data',0.005);
- Tag(length(Tag)+1)=struct('name','_DATASET_MAXHEADMOTIONTRIAL','type',7,'data',0);
- Tag(length(Tag)+1)=struct('name','_DATASET_MAXHEADMOTIONCOIL','type',10,'data','1');
- Tag(length(Tag)+1)=struct('name','EndOfParameters','type',-1,'data',[]);
- Tag(length(Tag)+1)=struct('name','EndOfParameters','type',-1,'data',[]);
- return
- % ************* End of function make_dummy_infods*********************************************
- % **************************************************************************************
- % **************************************************************************************
- % ************* Function writeHc*********************************************
- function hc=writeHc(hcFileName,hc,HLCdata);
- % Modified for v1.2 ds.hc structures. ds.hc.names has the exact, complete names of the
- % head coils. The coordinates relative to the subject may be ds.hc.head(MEG) OR
- % ds.hc.abdomen (fMEG).
- % Creates a .hc file in a CTF dataset
- % Inputs: hcFileName : Complete name of .hc file including path, basename and .hc ext.
- % hc : structure with nasion, left, right head coil positions in dewar and
- % CTF head coordinates.
- % HLCdata : head coil locations at the start of 1st trial. unit=m
- % coordinates=dewar. Used only if structure hc is empty and it is MEG
- % data (not fMEG data).
- % Output : hc : hc structure. hc=struct([]) on failure. hc is set to the coil positions
- % in array HLCdata if hc=[] on entry.
- % Check inputs
- if exist('hc')~=1
- hc=struct([]);
- elseif ~isstruct(hc)
- hc=struct([]);
- end
- hcfields=fieldnames(hc);
- if exist('HLCdata')~=1
- HLCdata=[];
- elseif ~isnumeric(HLCdata) | size(HLCdata,1)~=3 | size(HLCdata,2)<3
- HLCdata=[];
- end
- % Check hc Filename
- if exist('hcFileName')~=1;hcFileName=char([]);end
- if ~ischar(hcFileName) | isempty(hcFileName)
- fprintf('writeCTFds (writeHc): Called writeHc with bad file name.\n');
- hcFileName
- return
- end
- % Both hc and HLCdata bad?
- if length(hcfields)~=4 & isempty(HLCdata)
- fprintf('writeCTFds (writeHc): Called writeHc with bad hc and bad HLCdata.\n');
- hc=struct([])
- return
- elseif length(hcfields)~=4
- rstandard=8/sqrt(2)*[1 1 0;-1 1 0;1 -1 0]';
- rstandard(3,:)=-27;
- rdewar=100*HLCdata(:,1:3);
-
- % Convert from dewar coordinates to CTF head coordinates.
- originCTF=0.5*(hc.dewar(:,2)+hc.dewar(:,3));
- % Unit vectors for the CTF coordinates
- uCTF(:,1)=hc.dewar(:,1)-originCTF;
- uCTF(:,3)=cross(uxCTF,hc.dewar(:,2)-hc.dewar(:,3));
- uCTF(:,2)=cross(uCTF(:,3),uCTF(:,1));
- uCTF=uCTF./(ones(3,1)*sqrt(sum(uCTF.^2,1)));
- rCTF=uCTF'*(rdewar-originCTF*ones(1,3))
- hc=struct('names',strvcat('nasion','left ear','right ear'),...
- 'standard',rstandard,'dewar',rdewar,'head',rCTF);
- clear originCTF uCTF rstandard rdewar rCTF;
- end
- % Character strings for generating the .hc text file
- % Should never have both 'head' and 'abdomen' fields.
- labelword=strvcat('standard','measured','measured');
- printField=[strmatch('standard',hcfields) strmatch('dewar',hcfields) ...
- strmatch('head',hcfields) strmatch('abdomen',hcfields)];
- if ~strmatch('names',hcfields) | length(printField)~=3
- fprintf(['writeCTFds (writeHc): Structure hc does not have all of the required fields.\n',...
- ' No .hc file will appear in the output dataset.\n']);
- hc;
- hc=struct([]);
- return
- end
- relative=strvcat('dewar','dewar',hcfields{printField(3)});
- coilcoord=strvcat('standard','dewar',hcfields{printField(3)});
- comp='xyz';
- coilname=hc.names;
- fid=fopen(hcFileName,'w','ieee-be');
- for k=1:size(coilcoord,1)
- rcoil=getfield(hc,coilcoord(k,:));
- for coil=1:size(hc.names,1)
- clName=deblank(hc.names(coil,:));
- fwrite(fid,[labelword(k,:) ' ' clName ' coil position relative to ',...
- deblank(relative(k,:)) ' (cm):' char(10)],'char');
- for m=1:3
- fwrite(fid,[char(9) comp(m) ' = ' num2str(rcoil(m,coil),'%7.5f') char(10)],'char');
- end
- end
- end
- fclose(fid);
- status=0;
- return
- % ************* End of function writeHc*********************************************
- % **************************************************************************************
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %%%%%%%%%%Function checkMrk %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- function MarkerFileOK=checkMrk(ds);
- % Examines structure ds to see if a sensible MarkerFile can be created from ds.mrk.
- % Output: MarkerFileOK=1 : ds.mrk looks OK
- % MarkerFileOK=0 : ds.mrk cannot be a set of valid markers for these data.
- MarkerFileOK=isfield(ds,'mrk');
- if MarkerFileOK
- MarkerFileOK=~isempty(ds.mrk);
- end
- if MarkerFileOK
- % Are the markers appropriate?
- minMarkerTrial=[];
- minMarkerTime=[];
- maxMarkerTrial=[];
- maxMarkerTime=[];
- for k=1:length(ds.mrk)
- maxMarkerTrial=max([maxMarkerTrial max(ds.mrk(k).trial)]);
- maxMarkerTime=max([maxMarkerTime max(ds.mrk(k).time)]);
- minMarkerTrial=min([minMarkerTrial min(ds.mrk(k).trial)]);
- minMarkerTime=min([minMarkerTime min(ds.mrk(k).time)]);
- end
- if isempty(maxMarkerTrial) | isempty(maxMarkerTime)
- MarkerFileOK=0; % Do not create MarkerFile.mrk if all of the marker classes are empty.
- else
- MarkerFileOK=(maxMarkerTrial<=ds.res4.no_trials & minMarkerTrial>=1 & ...
- maxMarkerTime<=(ds.res4.no_samples/ds.res4.sample_rate) & ...
- minMarkerTime>=(-ds.res4.preTrigPts/ds.res4.sample_rate));
- if ~MarkerFileOK
- fprintf(['\nwriteCTFds (checkMrk): ds.mrk cannot possibly be a set of markers ',...
- 'for array(data).\n']);
- fprintf([' minMarkerTrial=%d (must be >=1) ',...
- 'maxMarkerTrial=%d (must be <=%d)\n'],...
- minMarkerTrial,maxMarkerTrial,ds.res4.no_trials);
- fprintf([' minMarkerTime=%7.4f (must be >=%7.4f) ',...
- 'maxMarkerTrial=%7.4f (must be <=%7.4f )\n'],...
- minMarkerTime,-ds.res4.preTrigPts/ds.res4.sample_rate,...
- maxMarkerTime,ds.res4.no_samples/ds.res4.sample_rate);
- fprintf(' MarkerFile.mrk will not be created.\n\n');
- end
- end
- end
- return
- %%%%%%%%%%%%%% end of checkMrk %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %%%%%%%%% Function writeEEG %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- function writeEEG(eegFileName,EEG);
- % Reads the EEG file of a dtaset and stores the infoemation in strucure array EEG.
- % EEG(k).chanName = channel name in the dataset (EEGmmm where mmm=channel number)
- % EEG(k).name = channel name given by the user (e.g. Fp4)
- % EEG(k).pos = electrode position in cm in CTF head coordinates
- % Check inputs
- if exist('eegFileName')~=1;eegFileName=char([]);end
- if ~ischar(eegFileName) | isempty(eegFileName)
- fprintf('writeCTFds (writeEEG): Called writeEEG with bad file name.\n');
- eegFileName
- EEG=struct([]);
- end
- if exist('EEG')~=1
- EEG=struct([]);
- elseif ~isstruct(EEG)
- EEG=struct([]);
- end
- if isempty(EEG);return;end
- fid=fopen(eegFileName,'w','ieee-be');
- if fid<0
- fprintf('writeCTFds (writeEEG): Could not open file %s\n',eegFileName);
- return
- end
- nEEG=length(EEG);
- for k=1:nEEG
- fprintf(fid,'%d\t%s\t%7.5f\t%7.5f\t%7.5f\n',EEG(k).chanNum,EEG(k).name,EEG(k).pos);
- end
- fclose(fid);
- return
- %%%%%%%%% End of writeEEG %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %%%%%%%%% Function writeBadSegments %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- function writeBadSegments(badSegmentsFile,badSegments,nTrial,tTrial);
- % Creates a bad.segements file in a CTF data set from the information in structure
- % badSegments which is created by read_badSegments, or by the user.
- % If structure badSegments is empty, then the file is not created.
- % badSegments structure:
- % badSegments.trial = List of trial numbers
- % badSegments.StartTime = List of bad segment start times (relative to trial).
- % badSegments.EndTime = List of bad segment end times.
- % Check badSegmentsFile
- if exist('badSegmentsFile')~=1;badSegmentsFile=char([]);end
- if isempty(badSegmentsFile) | ~ischar(badSegmentsFile)
- fprintf('writeCTFds(writeBadSegments): Bad file name.\n');
- badSegmentsFile
- return
- end
- % Check that structure badSegments is defined correctly
- if exist('badSegments')~=1 | exist('nTrial')~=1
- return
- elseif ~isstruct(badSegments) | isempty(badSegments)
- return
- elseif ~isfield(badSegments,'trial') | ~isfield(badSegments,'StartTime') | ...
- ~isfield(badSegments,'EndTime')
- return
- elseif isempty(badSegments.trial) | isempty(badSegments.StartTime) | ...
- isempty(badSegments.EndTime)
- return
- elseif ~isequal(size(badSegments.trial),size(badSegments.StartTime),...
- size(badSegments.EndTime))
- fprintf(['\nwriteCTFds (writeBadSegments): ',...
- 'The fields of structure badSegments do not all have the same size.\n']);
- return
- elseif any(badSegments.trial>nTrial) | any(badSegments.trial<1) | ...
- any(badSegments.EndTime>tTrial)
- fprintf(['\nwriteCTFds (writeBadSegments): ds.badSegments cannot possibly describe ',...
- 'bad segments for these data.\n',...
- '\tmin(badSegments.trial)=%d max(badSegments.trial)=%d ',...
- 'max(badSegments.EndTime)=%0.4f s\n\t\tDataset: nTrial=%d tTrial=%0.4f s\n'],...
- min(badSegments.trial),max(badSegments.trial),max(badSegments.EndTime),...
- nTrial,tTrial);
- fprintf('\t\tbad.segments file will not be created.\n\n');
- return
- end
- % Convert all fields to simple vectors
- nSeg=prod(size(badSegments.trial));
- trial=reshape(badSegments.trial,1,nSeg);
- StartTime=reshape(badSegments.StartTime,1,nSeg);
- EndTime=reshape(badSegments.EndTime,1,nSeg);
- fid=fopen(badSegmentsFile,'w','ieee-be');
- if fid<0
- fprintf('writeCTFds (writeBadSegments): Could not open file %s\n',badSegmentsFile);
- return
- end
- % Extra tabs are inserted to reproduce the format of bad.segments files produced by
- % DataEditor (5.3.0-experimental-linux-20060918).
- fprintf(fid,'%0.6g\t\t%0.6g\t\t%0.6g\t\t\n',[trial;StartTime;EndTime]);
- fclose(fid);
- return
- %%%%%%%%% End of writeBadSegments %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %%%%%%%%%%Function check_cls %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- function ClassFileOK=check_cls(ds);
- % Examines structure ds to see if a sensible ClassFile can be created from ds.TrialClass.
- ClassFileOK=isfield(ds,'TrialClass');
- if ClassFileOK
- ClassFileOK=~isempty(ds.TrialClass);
- end
- if ClassFileOK
- % Are the class trials appropriate?
- minClassTrial=[];
- maxClassTrial=[];
- for k=1:length(ds.TrialClass)
- maxClassTrial=max([maxClassTrial max(ds.TrialClass(k).trial)]);
- minClassTrial=min([minClassTrial min(ds.TrialClass(k).trial)]);
- end
- % Create ClassFile.cls even when the trail classes are empty.
- if ~isempty(maxClassTrial)
- ClassFileOK=(maxClassTrial<=ds.res4.no_trials & minClassTrial>=1);
- if ~ClassFileOK
- fprintf(['\nwriteCTFds (check_cls): ds.TrialClass cannot possibly be a set of ',...
- 'trial classes for array(data).\n minClassTrial=%d (must be >=1) ',...
- 'maxClassTrial=%d (must be <=%d)\n'],...
- minClassTrial,maxClassTrial,ds.res4.no_trials);
- fprintf(' ClassFile.cls will not be created.\n');
- end
- end
- end
- return
- %%%%%%%%%%%%%% end of check_cls %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %%%%%%%%%% Function writeClassFile %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- function writeClassFile(ClassFile,TrialClass);
- % Write the ClassFile of a CTF data set.
- % The CLassFile allows a user to store a list of trial classifications in a data set.
- % The ClassFile format is defined in document CTF MEG File Formats, PN900-0088.
- % This format is rigid.
- % Inputs :
- % ClassFile : marker file including the full path and extension .mrk.
- % TrialClass : Structure creted by read_ClassFile.
- % Output : ClassFile.cls.
- % Check input TrialClass.
- if exist('TrialClass')~=1;TrialClass=[];end
- if isempty(TrialClass) | ~isstruct(TrialClass)
- fprintf('writeCTFds (writeClassFile): TrialClass is empty or is not a structure.\n');
- TrialClass
- return
- end
- % Check ClassFile
- if exist('ClassFile')~=1;ClassFile=char([]);end
- if isempty(ClassFile) | ~ischar(ClassFile)
- fprintf('writeCTFds (writeClassFile): Bad file name.\n');
- ClassFile
- end
- fid=fopen(ClassFile,'w','ieee-be');
- if fid<0
- fprintf('writeCTFds (writeClassFile): Could not open file %s\n',ClassFile);
- return
- end
- nClass=length(TrialClass);
- % Generate datasetname from ClassFIle.
- ksep=max([0 strfind(ClassFile,filesep)]);
- datasetname=ClassFile(1:ksep-1);
- if isempty(datasetname);datasetname=cd;end
- fprintf(fid,'PATH OF DATASET:\n%s\n\n\n',datasetname);
- fprintf(fid,'NUMBER OF CLASSES:\n%d\n\n\n',nClass);
- for k=1:nClass
- if k==1 % Add sign character to make output match the output of Acq.
- sgn=char([]); % There should be no real significance to this.
- else % Why does DataEditor places the + sign only on ClassID 2,3,...?
- sgn='+';
- end
- No_of_Trials=prod(size(TrialClass(k).trial));
- fprintf(fid,'CLASSGROUPID:\n%s%d\nNAME:\n%s\nCOMMENT:\n%s\n',...
- sgn,TrialClass(k).ClassGroupId,TrialClass(k).Name,TrialClass(k).Comment);
- fprintf(fid,'COLOR:\n%s\nEDITABLE:\n%s\nCLASSID:\n%s%d\nNUMBER OF TRIALS:\n%d\n',...
- TrialClass(k).Color,TrialClass(k).Editable,sgn,TrialClass(k).ClassId,No_of_Trials);
- fprintf(fid,'LIST OF TRIALS:\nTRIAL NUMBER\n');
- % Subtract one so trial numbering starts at 0 in ClassFile.cls
- fprintf(fid,'%20d\n',reshape(TrialClass(k).trial,1,No_of_Trials)-1);
- fprintf(fid,'\n\n');
- end
- fclose(fid);
- return
- %%%%%%%%%%%%%% end of writeClassFile %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %%%%%%%%%% Function writeVirtualChannels %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- function writeVirtualChannels(VirtualChannelsFile,Virtual);
- % Writes a VirtualChannels file using the information in structure Virtual.
- % S\tructure Virtual is prepared by read_VirtualChannels
- % Check VirtualChannelsFile
- if exist('VirtualChannelsFile')~=1;VirtualChannelsFile=char([]);end
- if isempty(VirtualChannelsFile) | ~ischar(VirtualChannelsFile)
- fprintf('write_VirtualChannelsFile: Bad file name.\n');
- VirtualChannelsFile
- return
- end
- % Check structure array Virtual
- if exist('Virtual')~=1
- fprintf('writeVirtualChannels: Must specify structure Virtual.\n');
- return
- elseif isempty(Virtual) | ~isstruct(Virtual)
- return
- elseif ~isfield(Virtual,'Name') | ~isfield(Virtual,'wt');
- return
- elseif isempty(Virtual(1).Name)
- return
- end
- fid=fopen(VirtualChannelsFile,'w','ieee-be');
- if fid<0
- fprintf('writeCTFds (writeVirtualChannels): Could not open file %s\n',VirtualChannelsFile);
- return
- end
- fprintf(fid,'//Virtual channel configuration\n\n');
- for k=1:length(Virtual)
- fprintf(fid,'VirtualChannel\n{\n\tName:\t%s\n\tUnit:\t%s\n',...
- Virtual(k).Name,Virtual(k).Unit);
- for q=1:size(Virtual(k).wt)
- % Floating format chosen to match VirtualChanels file creatd by
- % DataEditor (5.3.0-experimental-linux-20060918).
- fprintf(fid,'\tRef:\t%s,%0.6g\n',deblank(Virtual(k).chan(q,:)),Virtual(k).wt(q));
- end
- fprintf(fid,'}\n\n');
- end
- fclose(fid);
- return
- %%%%%%%%%%%%%% end of writeVirtualChannels %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %%%%%%%%%%%%%% Function updateDescriptors %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- function ds=updateDescriptors(ds,comment,creatorSoftware,HLCcntrl);
- % Makes sure that certain tags are in structure infods, and makes sure that information
- % is transfered to infods and newds from res4. Updates several fields of res4 with
- % clinical use message (=comment) and name of creator program.
- % Inputs: ds : ds structure produced by readCTFds.
- % comment : Character string that is added to infods tags listed in addCommentTag
- % and res4 fields listed in addCommentField
- % creatorSoftware : Character string indicating that the data set was
- % created by writeCTFds. Added to infods tags listed in addCreatorTag and
- % appName field of res4.
- % HLCcntrl: If ds.infods is missing or empty, HLCcntrl determines the
- % _DATA_HZ_MODE tage of infods. if noit present or empty, HLCcntrl=0.
- % Creates a dummy infods structure if necessary using function make_dummy_infods.
- % Changes infods and newds to match information in
- % ds.res4.nf_run_title
- % ds.res4.collect_descriptor
- % ds.res4.nf_subject_id
- % ds.res4.nf_operator
- % Adds comment (clinical use message) and creator software name to infods, newds and hist files.
- if ~isfield(ds,'infods');ds.infods=[];end
- % infods tags that will display the comment.
- addCommentTag=strvcat('_PATIENT_ID','_PATIENT_INSTITUTE','_PROCEDURE_COMMENTS',...
- '_DATASET_COMMENTS','_DATASET_STATUS','_DATASET_PROCSTEPTITLE');
- % res4 text fields that will display the comment
- addCommentField=strvcat('appName','dataOrigin','dataDescription',...
- 'nf_run_title','nf_subject_id','run_description');
- % res4 text string lengths (from document "CTF MEG File Formats", PN900-0088)
- addCommentLength=[256 256 256 256 32 -1]'; % -1 indicates variable length
- % infods tags that will display the creator software.
- addCreatorTag=strvcat('_PROCEDURE_COMMENTS','_DATASET_STATUS',...
- '_DATASET_COLLECTIONSOFTWARE','_DATASET_CREATORSOFTWARE');
- % res4 text fields that will display the creator software.
- addCreatorField=strvcat('appName','run_description');
- addCreatorLength=[256 -1]';
- % List of infods tags that will be set to ds.res4.fields (not necessarily fields listed in
- % addCommentField or addCreatorField
- addRes4Tag=strvcat('_PATIENT_ID','_DATASET_PROCSTEPTITLE',...
- '_DATASET_PROCSTEPDESCRIPTION','_DATASET_OPERATORNAME');
- % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % Add the Tags to infods. In most cases they will already be there, but just to be sure.
- addTag=strvcat(addRes4Tag,addCommentTag,addCreatorTag);
- tagClass=strvcat('_PATIENT','_PROCEDURE','_DATASET');
- % List of all the tag names as a character array
- tagName=getArrayField(ds.infods,'name')';
- if length(ds.infods)==1;tagName=tagName';end
- tagType=getArrayField(ds.infods,'type');
- tagPtr=1; % If a class is missing, inject the class starting at tagPtr+1.
- for k=1:size(tagClass,1)
- addIndex=strmatch(deblank(tagClass(k,:)),addTag)';
- % Are there any tags to be added in this class?
- if isempty(addIndex);continue;end
- % List of infods tags in the tagClass, but excluding the CPerist type (which marks the
- % start of a class.
- infodsIndex=strmatch(deblank(tagClass(k,:)),tagName);
- if isempty(infodsIndex) % Create a new class of tags.
- if strcmp(deblank(tagName(tagPtr+1,:)),'EndOfParameters') & k>1;
- tagPtr=tagPtr+1;
- end
- nTag=length(ds.infods);
- tagName((tagPtr+4):(nTag+3),:)=tagName((tagPtr+1):nTag,:);
- …
Large files files are truncated, but you can click here to view the full file