ctf_read_res4.m

读取CTF脑磁图数据的Matlab代码

function [ctf] = ctf_read_res4(folder,SETUP);

% ctf_read_res4 - Read a CTF .res4 file
%
% ctf = ctf_read_res4( [folder], [SETUP])
% 
% This function reads the resource information from a CTF .ds folder.  This
% resource information must be read before reading the data from a .meg4
% file.
%
% INPUTS
%
%   folder: the .ds directory containing the data to be read.  With
%   no input, it will prompt with a gui folder locator.  If only
%   the folder is input, SETUP = 1.
%
%   SETUP: If SETUP = 1, display 'ctf.setup' structure (default)
%          If SETUP = 0, do not display 'ctf.setup' structure
%
% OUTPUTS
%
%   ctf.folder: path of the .ds folder
%
%   ctf.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: 1.3 $ $Date: 2005/02/18 04:53:42 $

% Licence:  GNU GPL, no express or implied warranties
% Modified: 2/2005, Fred Carver
%                   - modified to read new and old format ctf datasets
%                   - modified to read virtual SAM channels
% Modified: 11/2003, Darren.Weber_at_radiology.ucsf.edu
%                    - modified from NIH code readresfile.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ver = '$Revision: 1.3 $';
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('SETUP','var'), SETUP = 1; end

[rootpath,rootname] = fileparts(ctf.folder);
res4file = [ctf.folder,filesep,rootname,'.res4'];
[fid,message] = fopen(res4file,'rb','s');

% 'ieee-be.l64' or 's' - IEEE floating point with big-endian byte ordering
%                        and 64 bit long data type.

if fid < 0,
  error(message);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% READ HEADER

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

% check for the right format
%if strmatch('MEG41RS',ctf.header.res4),
  % OK, we can handle this format
%else
%  msg = sprintf('May not read "%s" format correctly',ctf.header.res4);
%  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,numparam,'double');
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 SETUP,
  fprintf('...data format       : %s\n',ctf.header.res4);
  fprintf('...date collected    : %s\n',ctf.setup.date);
  fprintf('...time collected    : %s\n',ctf.setup.time);
  fprintf('...run name          : %s\n',ctf.setup.run_name);
  fprintf('...run title         : %s\n',ctf.setup.run_title);
  fprintf('...subject           : %s\n',ctf.setup.subject);
  fprintf('...run description   : %s\n',ctf.setup.run_description);
  fprintf('...operator          : %s\n',ctf.setup.operator);
  fprintf('...number of channels: %d\n',ctf.setup.number_channels);
  fprintf('...number of samples : %d per trial\n',ctf.setup.number_samples);
  fprintf('...sample rate       : %g samples/sec\n',ctf.setup.sample_rate);
  fprintf('...number of trials  : %g trials (average of %g)\n',ctf.setup.number_trials,ctf.setup.number_trials_averaged);
  fprintf('...duration          : %g sec total\n',ctf.setup.duration_total);
  fprintf('...duration          : %g sec / trial\n',ctf.setup.duration_trial);
  fprintf('...pretrigger points : %g samples\n',ctf.setup.pretrigger_samples);
  fprintf('...pretrigger msec   : %g msec\n',ctf.setup.pretrigger_msec);
  fprintf('...sensor file name  : %s\n',ctf.setup.sensor_file_name);
  fprintf('...head zeroing      : %s\n',ctf.setup.head_zero);
  fprintf('...number of filters : %d\n',ctf.setup.number_filters);
  for i = 1:ctf.setup.number_filters,
    fprintf(' -Filter # \t%g\n',i)
    fprintf(' -Frequency: \t%g Hz\n',ctf.setup.filters{i}.freq)
    fprintf(' -Class: \t%g\n',ctf.setup.filters{i}.class)
    fprintf(' -Type: \t%g\n',ctf.setup.filters{i}.type)
    if ~isempty(ctf.setup.filters{i}.params)
      fprintf(' -Parameter(s): \t%g\n',ctf.setup.filters{i}.params)
    end
  end
  fprintf('\n');
end






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

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

% 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');
  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

ctf.sensor.type.eeg = 9;
ctf.sensor.type.meg = 5;
ctf.sensor.type.ref = [0 1];
ctf.sensor.type.stim = 11;
ctf.sensor.type.adc = 18;
ctf.sensor.type.vc  = 15;

eegsens = find([ctf.sensor.info.index] == 9);
megsens = find([ctf.sensor.info.index] == 5);
refsens = find([ctf.sensor.info.index] == 1);
refsens = [refsens,find([ctf.sensor.info.index] == 0)];
virtsens = find([ctf.sensor.info.index] == 15);
allsens = [eegsens,megsens,refsens,virtsens];
othersens = setdiff([1:ctf.setup.number_channels],allsens);

ctf.sensor.index = struct(...
  'eeg',{eegsens},...
  'meg',{megsens},...
  'ref',{refsens},...
  'vc',{virtsens},...
  'other',{othersens} ...
  );



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Location coordinates of channels in centimeters, in patient head space

for chan = 1:ctf.setup.number_channels,
  
  switch ctf.sensor.info(chan).index,
    
    case {0,1,5},
      %0=Reference Channels, 
      %1=More Reference Channels, 
      %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

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

return