ctf_read_res4.m

含有多种ICA算法的eeglab工具箱

  temp = strtok(temp,char(0));
  ctf.sensor.info(chan).label = char(temp');
end

for chan = 1:ctf.setup.number_channels,
  
  %ftell(fid);
  
  ctf.sensor.info(chan).index = fread(fid,1,'int16');
  ctf.sensor.info(chan).extra = fread(fid,1,'int16');
  id = fread(fid,1,'int32')+1;
  ctf.sensor.info(chan).proper_gain = fread(fid,1,'double');
  ctf.sensor.info(chan).q_gain = fread(fid,1,'double');
  ctf.sensor.info(chan).io_gain = fread(fid,1,'double');
  ctf.sensor.info(chan).io_offset = fread(fid,1,'double');
  fread(fid,1,'int16');
  ctf.sensor.info(chan).grad_order_no = fread(fid,1,'int16');
  fread(fid,1,'int32');
  
  %fseek(fid,ftell(fid)+6,0);
  
  for pos = 1:8,
    ctf.sensor.info(chan).coil(pos).position.x = fread(fid,1,'double');
    ctf.sensor.info(chan).coil(pos).position.y = fread(fid,1,'double');
    ctf.sensor.info(chan).coil(pos).position.z = fread(fid,1,'double');
    fread(fid,1,'double');
    ctf.sensor.info(chan).coil(pos).orient.x = fread(fid,1,'double');
    ctf.sensor.info(chan).coil(pos).orient.y = fread(fid,1,'double');
    ctf.sensor.info(chan).coil(pos).orient.z = fread(fid,1,'double');
    fread(fid,1,'double');
    fread(fid,1,'int16');
    fread(fid,1,'int32');
    fread(fid,1,'int16');
    fread(fid,1,'double');
    
    %fseek(fid,ftell(fid)+56,0);
    %fseek(fid,ftell(fid)-80,0);
  end
  
  for pos = 1:8,
    ctf.sensor.info(chan).hcoil(pos).position.x = fread(fid,1,'double');
    ctf.sensor.info(chan).hcoil(pos).position.y = fread(fid,1,'double');
    ctf.sensor.info(chan).hcoil(pos).position.z = fread(fid,1,'double');
    fread(fid,1,'double');
    ctf.sensor.info(chan).hcoil(pos).orient.x = fread(fid,1,'double');
    ctf.sensor.info(chan).hcoil(pos).orient.y = fread(fid,1,'double');
    ctf.sensor.info(chan).hcoil(pos).orient.z = fread(fid,1,'double');
    fread(fid,1,'double');
    fread(fid,1,'int16');
    fread(fid,1,'int32');
    fread(fid,1,'int16');
    fread(fid,1,'double');
    
    %fseek(fid,ftell(fid)+56,0);
    %fseek(fid,ftell(fid)+80,0);
  end
  %fseek(fid,ftell(fid)+1288,-1);
end




%-------------------------------------------------------------
% Find channel types and define channel sets, see the
% System Administrators .pdf, 'Channel Sets Configuration'

ctf = ctf_channel_sets(ctf);



%-------------------------------------------------------------
% Channel coordinates, in centimeters, in subject head space

for chan = 1:ctf.setup.number_channels,
  
  switch ctf.sensor.info(chan).index,
    
    case {0,1,5},
      %0=Reference Magnetometers
      %1=Reference Gradiometers
      %5=MEG Channels
      
      coord = [ctf.sensor.info(chan).hcoil(1:2).position];
      ctf.sensor.info(chan).location = [coord.x; coord.y; coord.z];
      
      orient = [ctf.sensor.info(chan).hcoil(1).orient];
      ctf.sensor.info(chan).orientation = [orient.x; orient.y; orient.z];
      
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % This ensures that the orientation of the sensor is away from the
      % center of a sphere.  It uses the sign of the dot product between
      % the orientation vector and the location vector.
      tmp = ctf.sensor.info(chan).orientation' * ctf.sensor.info(chan).location;
      tmp = sign(tmp(1));
      ctf.sensor.info(chan).orientation = tmp * ctf.sensor.info(chan).orientation;
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      
    case 9,
      %EEG Channels    
      
      coord = [ctf.sensor.info(chan).hcoil(1:2).position];
      ctf.sensor.info(chan).location = [coord.x; coord.y; coord.z];
      ctf.sensor.info(chan).orientation = [];
      
  end
end



%%%%%% Coefficient reading appended by SSD %%%%%%

if(COEFS)
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% NUMBER OF COEFFICIENTS, byte offset = b+f+nc*1360, byte size = 1
	
	fseek(fid,offset,-1);
	ctf.res4.numberCoefficients = fread(fid,1,'int16');
	
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% SENSOR COEFFICIENT RECORDS, byte offset = b+f+nc*1360+2, byte size = 1992
	
	if VERBOSE & ctf.res4.numberCoefficients,
		fprintf('...reading %d coefficients\n',ctf.res4.numberCoefficients);
	end
	
	SENSOR_LABEL  = 31;
	MAX_NUM_COEFS = 50;
	MAX_BALANCING = MAX_NUM_COEFS;
	
	hexadef = {'00000000','47314252','47324252','47334252','47324f49','47334f49'};
	strdef = {'NOGRAD','G1BR','G2BR','G3BR','G2OI','G3OI'};
	
	for i = 1:ctf.res4.numberCoefficients,
		
		% read the sensor name (channel label)
		temp = fread(fid,[1,32],'char');
		temp(temp>127) = 0;
		temp(temp<0) = 0;
		temp = strtok(temp,char(0));
		temp = strtok(temp,'-');
		sensorName = char(temp);
		sensorIndex = strmatch( sensorName, {ctf.sensor.info.label} );
		
		% read the coefficient type
		coefType = fread(fid,1,'bit32');
		
		padding = fread(fid,1,'int32'); % not sure why this is needed???
		
		% read the coefficient record
		numberCoefs = fread(fid,1,'int16');
		
		if numberCoefs > MAX_NUM_COEFS,
			msg = sprintf('numberCoefs > MAX_NUM_COEFS\n');
			warning(msg);
		end
      
		sensor_list = char(fread(fid,[SENSOR_LABEL,MAX_BALANCING],'uchar')');
		% clean-up the sensor_list
		sensor_list = sensor_list(1:numberCoefs,:);
		for j=1:numberCoefs,
			temp = strtok(sensor_list(j,:),char(0));
			temp = strtok(temp,'-');
			% check if this sensor is a reference
			refLabels = {ctf.sensor.info(ctf.sensor.index.meg_ref).label};
			refIndex = strmatch(temp,refLabels);
			if refIndex,
				% ensure this one has the same label
				temp = refLabels{refIndex};
			end
       
			new_sensor_list(j) = refIndex;
		end
		sensor_list = ctf.sensor.index.meg_ref(:,new_sensor_list)';
		
		coefs_list = fread(fid,MAX_BALANCING,'double');
		% clean-up the coefs_list
		coefs_list = coefs_list(1:numberCoefs,:)';
		
		% allocate the coefficient parameters into the ctf struct
		ctf.res4.sensorCoef(i).sensorName = sensorName;
		ctf.res4.sensorCoef(i).coefType   = coefType;
		ctf.res4.sensorCoef(i).coefRec.numberCoefs = numberCoefs;
		ctf.res4.sensorCoef(i).coefRec.sensor_list = sensor_list;
		ctf.res4.sensorCoef(i).coefRec.coefs_list  = coefs_list;
			
          
          
		% DLW:
		% This is a brainstorm variable, note the use of coefType
		% Not clear why this is checked and allocated as such
		coefType = find( hex2dec(hexadef) == coefType );
		if coefType,
				CoefInfo{sensorIndex,coefType-1}.numberCoefs = numberCoefs;
				CoefInfo{sensorIndex,coefType-1}.sensor_list = sensor_list;
				CoefInfo{sensorIndex,coefType-1}.coefs = coefs_list;
			end
		end	
		
		
            
        % Channel Gains
        gain_chan = zeros(size(ctf.setup.number_channels,1),1);
        gain_chan(ctf.sensor.index.meg_sens) = ([ctf.sensor.info(ctf.sensor.index.meg_sens).proper_gain]'.*[ctf.sensor.info(ctf.sensor.index.meg_sens).q_gain]');
        gain_chan(ctf.sensor.index.meg_ref) = ([ctf.sensor.info(ctf.sensor.index.meg_ref).proper_gain]'.*[ctf.sensor.info(ctf.sensor.index.meg_ref).q_gain]');
        %  gain_chan(ieegsens) = 1./([SensorRes(ieegsens).qGain]'*1e-6);
		%     gain_chan(ieegsens) = 1./([SensorRes(ieegsens).qGain]');
		%     gain_chan(iothersens) = ([SensorRes(iothersens).qGain]'); % Don't know exactly which gain to apply here
        
       
        
        % Calculus of the matrix for nth-order gradient correction
        % Coefficients for unused reference channels are weigthed by zeros in
        % the correction matrix.
        Gcoef = zeros(length(ctf.sensor.index.meg_sens),length(min(ctf.sensor.index.meg_ref):max(ctf.sensor.index.meg_ref)));
        grad_order_no = 3*ones(306,1);
        for k = 1:length(ctf.sensor.index.meg_sens)
            
            % Reference coils for channel k
            if grad_order_no(ctf.sensor.index.meg_sens(k)) == 0
                %Data is saved as RAW
                %Save 3rd order gradient sensor-list for subsequent correction if requested later by the user
                [refs] = (CoefInfo{ctf.sensor.index.meg_sens(k),3}.sensor_list);
                Gcoef(k,refs-min(ctf.sensor.index.meg_ref)+1) = CoefInfo{ctf.sensor.index.meg_sens(k),3}.coefs ... 
                    .* gain_chan(refs)'/gain_chan(ctf.sensor.index.meg_sens(k)); 
            else
                [refs] = (CoefInfo{ctf.sensor.index.meg_sens(k),grad_order_no(ctf.sensor.index.meg_sens(k))}.sensor_list);
                Gcoef(k,refs-min(ctf.sensor.index.meg_ref)+1) = CoefInfo{ctf.sensor.index.meg_sens(k),grad_order_no(ctf.sensor.index.meg_sens(k))}.coefs ... 
                    .* gain_chan(refs)/gain_chan(ctf.sensor.index.meg_sens(k)); 
            end
            ctf.sensor.info(ctf.sensor.index.meg_sens(k)).Gcoef = Gcoef(k,:);
        end
end   %% end COEF block


fclose(fid);

t = toc; fprintf('...done (%6.2f sec)\n\n',t);

return



% find file name if truncated or with uppercase extension
% added by Arnaud Delorme June 15, 2004
% -------------------------------------------------------
function res4name = findres4file( folder )

res4name = dir([ folder filesep '*.res4' ]);
if isempty(res4name)
    res4name = dir([ folder filesep '*.RES4' ]);
end

if isempty(res4name)
    error('No file with extension .res4 or .RES4 in selected folder');
else
    res4name = [ folder filesep res4name.name ];
end;
return