ctf_read_res4.m

function [ctf] = ctf_read_res4(folder,VERBOSE,COEFS);

% ctf_read_res4 - Read a CTF .res4 file
%
% ctf = ctf_read_res4( [folder], [verbose], [coefs])
% 
% This function reads the resource information from a CTF .ds folder.  This
% resource information must be read before reading the .meg4 data file.
% All input arguments are optional.
%
% INPUTS
%
%   folder: the .ds directory containing the data to be read.  With
%           no input, it will prompt with a gui folder locator.
%
%   verbose: If verbose = 1, display 'ctf.setup' structure (default)
%            If verbose = 0, do not display 'ctf.setup' structure
%
%   coefs: an option to read the sensory coefficients, which give the
%          weights for calculation of synthetic 2nd or 3rd order
%          gradiometers.
%          If coefs = 1, read the sensor coefficients
%          If coefs = 0, do not read the sensor coefficients (default)
% 
% OUTPUTS
%
%   ctf.folder - path of the .ds folder
%
%   ctf.res4.file - data file path/name
%   ctf.res4.header - data format header
%
%   ctf.setup - a header structure consisting of date, time, run name, run
%   title, subject, run description, operator, number of channels, number
%   of samples, sample rate, number of trials, duration, pretrigger_samples,
%   sensor filename, head zeroing, and number of filters.
%   
%   ctf.sensor.index - a sensor structure consisting of EEG sensors, MEG
%   sensors, reference sensors, and other sensors.
%   
%   ctf.sensor.info - a structure with gain and offset information consisting
%   of proper gain, Q gain, io gain, io offset, and index.
%
%      <>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
%      <                                                      >
%      <                    DISCLAIMER:                       >
%      <                                                      >
%      < THIS PROGRAM IS INTENDED FOR RESEARCH PURPOSES ONLY. >
%      < THIS PROGRAM IS IN NO WAY INTENDED FOR CLINICAL OR   >
%      <                    OFFICIAL USE.                     >
%      <                                                      >
%      <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>
%

% $Revision: 253 $ $Date: 2004/08/19 03:17:10 $

% Copyright (C) 2003  Darren L. Weber
% 
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
% 
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% 
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

% Modified: 11/2003, Darren.Weber_at_radiology.ucsf.edu
%                    - modified from NIH code readresfile.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ver = '$Revision: 253 $';
fprintf('\nCTF_READ_RES4 [v %s]\n',ver(11:15)); tic;

if ~exist('folder','var'),
  ctf = ctf_folder;
else
  ctf = ctf_folder(folder);
end

if ~exist('VERBOSE','var'), VERBOSE = true; end
if ~exist('COEFS','var'), COEFS = false; end

[folderPath,folderName,folderExt] = fileparts(ctf.folder);
ctf.res4.file = findres4file(ctf.folder);

%----------------------------------------------------------------
% open the data file

[fid,message] = fopen(ctf.res4.file,'rb','ieee-be.l64');
if fid < 0, error('cannot open .res4 file'); end


%-------------------------------------------------------------
% READ HEADER

fseek(fid,0,-1);
ctf.res4.header = char(fread(fid,8,'char'))';

% check for the right format
if strmatch('MEG41RS',ctf.res4.header),
  % OK, we can handle this format
else
  msg = sprintf('This function is designed to read MEG41RS format.\nIt may not read "%s" format correctly',ctf.res4.header);
  warning(msg);
end


%-------------------------------------------------------------
% READ SETUP


%---DATE/TIME
fseek(fid, 778,-1);
ctf.setup.time = char(fread(fid,255,'char'))';
fseek(fid,1033,-1);
ctf.setup.date = char(fread(fid,255,'char'))';

%---NUMBER OF SAMPLES
fseek(fid,1288,-1);
ctf.setup.number_samples = (fread(fid,1,'int32')');

%---NUMBER OF CHANNELS
fseek(fid,1292,-1);
ctf.setup.number_channels = (fread(fid,1,'int16')');

%---SAMPLE RATE
fseek(fid,1296,-1);
ctf.setup.sample_rate = fread(fid,1,'double')';
ctf.setup.sample_msec = 1000 / ctf.setup.sample_rate;
ctf.setup.sample_sec  =    1 / ctf.setup.sample_rate;

%---NUMBER OF TRIALS
fseek(fid,1312,-1);
ctf.setup.number_trials = fread(fid,1,'int16')';
fseek(fid, 776,-1);
ctf.setup.number_trials_averaged = fread(fid,1,'int16');

%---DURATION
fseek(fid,1304,-1);
ctf.setup.duration_total = fread(fid,1,'double');
ctf.setup.duration_trial = ctf.setup.duration_total / ctf.setup.number_trials;

%---PRE_TRIG POINTS
fseek(fid,1316,-1);
ctf.setup.pretrigger_samples = fread(fid,1,'int32');
ctf.setup.pretrigger_msec = (ctf.setup.pretrigger_samples / ctf.setup.sample_rate) * 1000;

%---HEAD ZEROING
fseek(fid,1348,-1);
h_zero = fread(fid,1,'int32')';
no_yes = {'no','yes'};
ctf.setup.head_zero = no_yes{h_zero+1};

%---RUN NAME
fseek(fid,1360,-1);
ctf.setup.run_name = char(fread(fid,32,'char'))';

%---RUN TITLE
fseek(fid,1392,-1);
ctf.setup.run_title = char(fread(fid,256,'char'))';

%---SUBJECT
fseek(fid,1712,-1);
ctf.setup.subject = char(fread(fid,32,'char'))';

%---OPERATOR
fseek(fid,1744,-1);
ctf.setup.operator = char(fread(fid,32,'char'))';

%---SENSOR FILE NAME
fseek(fid,1776,-1);
ctf.setup.sensor_file_name = char(fread(fid,60,'char'))';

%---RUN DESCRIPTION & FILTERS
fseek(fid,1836,-1);
run_size = fread(fid,1,'int32');
fseek(fid,1844,-1);
ctf.setup.run_description = char(fread(fid,run_size,'char'));

ctf.setup.number_filters = fread(fid,1,'int16');

for i = 1:ctf.setup.number_filters,
  ctf.setup.filters(i).freq     = fread(fid,1,'double');
  ctf.setup.filters(i).class    = fread(fid,1,'int32');
  ctf.setup.filters(i).type     = fread(fid,1,'int32');
  ctf.setup.filters(i).numparam = fread(fid,1,'int16');
  ctf.setup.filters(i).params   = fread(fid,ctf.setup.filters(i).numparam,'double');
end

if(COEFS)
    b = 1846 + run_size;
    if(ctf.setup.number_filters == 0)
        np = 0;
    else
        np = ctf.setup.filters.numparam;
        warning('3rd gradient + hardware filter parameters not fully tested! let''s see what happens... :)');
    end
    nf = ctf.setup.number_filters;
    f = ( nf * 18 ) + ( np * 8 );
    offset = b + f + ctf.setup.number_channels * 1360;
end


%-------------------------------------------------------------
% CREATE TIME ARRAYS

% the time arrays must be based on increments of the sample_msec
ctf.setup.time_msec = [0:ctf.setup.number_samples - 1]' * ctf.setup.sample_msec;
ctf.setup.time_msec = ctf.setup.time_msec - ctf.setup.pretrigger_msec;

% adjust the sample point closest to zero so that it is zero, if it
% is reasonably close to zero, say within 3 decimal places for msec timing
zero_index = find(abs(ctf.setup.time_msec) == min(abs(ctf.setup.time_msec)));
zero_value = ctf.setup.time_msec(zero_index);
if (-0.0001 < zero_value) & (zero_value < 0.0001),
  ctf.setup.time_msec(zero_index) = 0;
end

ctf.setup.start_msec = ctf.setup.time_msec(1);
ctf.setup.end_msec   = ctf.setup.time_msec(end);

ctf.setup.time_sec  = ctf.setup.time_msec / 1000;
ctf.setup.start_sec = ctf.setup.time_sec(1);
ctf.setup.end_sec   = ctf.setup.time_sec(end);


%-------------------------------------------------------------
% PRINT SETUP
if VERBOSE,
  ctf_print_setup(ctf);
end





%-------------------------------------------------------------
% READ SENSOR INFORMATION

ctf.sensor.info = struct(...
  'proper_gain',[],...
  'q_gain',[],...
  'io_gain',[],...
  'io_offset',[],...
  'index',[],...
  'extra',[],...
  'label',[],...
  'grad_order_no',[]);

% read channel names (not trivial!)
for chan = 1:ctf.setup.number_channels,
  temp = fread(fid,32,'char');
  temp(temp>127) = 0;
  temp(temp<0) = 0;
  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