ctf_read_meg4.m

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

function [ctf] = ctf_read_meg4(folder,ctf,CHAN,TIME,TRIALS);

% ctf_read_meg4 - read meg4 format data from a CTF .ds folder
%
% [ctf] = ctf_read_meg4([folder],[ctf],[CHAN],[TIME],[TRIALS]);
%
% INPUTS
%
% folder - the directory of the .ds data set to read.  By
% default, a gui prompts for the folder.
%
% ctf - a struct with setup, sensor and data fields.  If the setup field is
% missing or empty, this function calls ctf_read_res4.
%
% CHAN - a integer array of channel numbers to read.
% eg, [30:35] reads channels 30 to 35.  Also
% If CHAN = 'eeg', read only eeg channels/sensorIndices
% If CHAN = 'meg', read only meg channels/sensorIndices
% If CHAN = 'ref', read only reference channels/sensorIndices
% IF CHAN = 'vc', read only SAM virtual channels/sensorIndices
% If CHAN = 'other', read only the other channels/sensorIndices
%
% TIME - eg. [0 5] - the desired time interval to read, in sec.
%        If TIME = 'all', all data is read (the default)
%
% TRIALS - If TRIALS = n, the nth trial will be read.
%          If TRIALS = [3,5,8], reads trials 3,5, and 8 such that
%          ctf.data{1} = data for trial 3,
%          ctf.data{2} = data for trial 5, and 
%          ctf.data{3} = data for trial 8.
%          If TRIALS = [3:7], reads trials 3 to 7
%          If TRIALS = 'all', reads all data (the default)
%
% OUTPUTS
% ctf.data                  - cell array of all the data, eg. ctf.data{x}
% ctf.sensor.names          - cell array of sensor names
% ctf.sensor.location       - array of sensor locations for plotting
% ctf.sensor.orientation    - array of sensor orientations
% ctf.header                - used for writing new meg4 file
%
%
%      <>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> %
%      <                                                       > %  
%      <                      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: 11/2003, Darren.Weber_at_radiology.ucsf.edu
%                    - modified from NIH code simply to allocate data into
%                    one large struct (ctf)
%                    - modified channel selection section at the end so
%                    that it doesn't try to get orientation information for
%                    EEG channels
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

if ~isfield(ctf,'setup'),
  ctf = ctf_read_res4(ctf.folder);
end

if ~exist('CHAN','var'),   CHAN   = 'all'; end
if ~exist('TIME','var'),   TIME   = 'all'; end
if ~exist('TRIALS','var'), TRIALS = 'all'; end

if isempty(CHAN),   CHAN   = 'all'; end
if isempty(TIME),   TIME   = 'all'; end
if isempty(TRIALS), TRIALS = 'all'; end

smallCHAN = 0;

switch num2str(CHAN),
  case 'meg',
    CHAN = ctf.sensor.index.meg;
  case 'ref',
    CHAN = ctf.sensor.index.ref;
  case 'eeg',
    CHAN = ctf.sensor.index.eeg;
  case 'other',
    CHAN = ctf.sensor.index.other;
  case 'all',
    CHAN = [1:ctf.setup.number_channels];
  otherwise
    % assume the input is a range of channel numbers
end

switch num2str(TIME),
  case 'all',
    TIME = ctf.setup.time_sec;
  otherwise
    % assume the input is a range of times
end

switch num2str(TRIALS),
  case 'all',
    TRIALS = 1:ctf.setup.number_trials;
  otherwise
    % assume the input is a range of trials
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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


% open the data file
[path,rootname] = fileparts(ctf.folder);
[fid,message] = fopen([ctf.folder,filesep,rootname,'.meg4'],'rb','s');
if fid < 0, error(message); end

% Read the header (data format identifier)
ctf.header.meg4 = char(fread(fid,8,'char'))';

% check the format
%if strmatch('MEG41CP',ctf.header.meg4),
  % OK, we can handle this format
%else
%  msg = sprintf('May not read "%s" format correctly',ctf.header.meg4);
%  warning(msg);
%end


% Get sensor indices
megIndex =   ctf.sensor.index.meg;
refIndex =   ctf.sensor.index.ref;
eegIndex =   ctf.sensor.index.eeg;
vcIndex =    ctf.sensor.index.vc;
otherIndex = ctf.sensor.index.other;

% Setup gains and offsets
channel_gain = zeros(1,ctf.setup.number_channels);
channel_gain(megIndex)   = [ctf.sensor.info(megIndex).proper_gain] .* [ctf.sensor.info(megIndex).q_gain];
channel_gain(refIndex)   = [ctf.sensor.info(refIndex).proper_gain] .* [ctf.sensor.info(refIndex).q_gain];
channel_gain(eegIndex)   = [ctf.sensor.info(eegIndex).q_gain];
channel_gain(vcIndex)   = [ctf.sensor.info(vcIndex).q_gain];
channel_gain(otherIndex) = [ctf.sensor.info(otherIndex).q_gain];



trial_size  = 4 * ctf.setup.number_channels * ctf.setup.number_samples;

small_trial = 4 * (min(CHAN)-1) * ctf.setup.number_samples;

large_trial = 4 * (ctf.setup.number_channels - max(CHAN)) * ctf.setup.number_samples; 

% check the duration

duration = TIME(end) - TIME(1);

if duration > ctf.setup.duration_trial,
  fprintf('...TIME input too large for trial\n')
  duration = ctf.setup.duration_trial;
  fprintf('...setting TIME= %g seconds (ctf.setup.duration_trial)',ctf.setup.duration_trial);
end

samples = round((duration) * ctf.setup.sample_rate) + 1;
intime  = round((TIME(1) - ctf.setup.time_sec(1)) * ctf.setup.sample_rate) + 1;

% s_rt = 1/(ctf.setup.time_sec(2) - ctf.setup.time_sec(1));
% samples = round((duration)*s_rt)+1;
% intime = round((TIME(1)-ctf.setup.time_sec(1))*s_rt)+1;



% Read trial data

ctf.setup.number_trials = length(TRIALS);
ctf.setup.number_channels = length(CHAN);
channels = length([min(CHAN):max(CHAN)]);
ctf.data = cell(ctf.setup.number_trials,1);

tr = 0;
for trial = TRIALS,
  
  tr = tr + 1;
  
  ctf.data{tr} = zeros(samples,channels);
  
  if tr > 1,
    backspaces = '\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b';
    fprintf([backspaces,'...reading %4d of %4d trials\n'],tr,ctf.setup.number_trials);
  else
    fprintf('...reading %4d of %4d trials\n',tr,ctf.setup.number_trials);
  end
  
  % Define data file bytes to read
  if trial == TRIALS(1);
    % 1st trial
    bytes = [ (trial-1) * trial_size ] + small_trial + [ 4 * (intime-1) ];
    fseek(fid,bytes,0);
  else
    bytes = [ (trial - TRIALS(tr-1) -1) * trial_size ] + small_trial + large_trial;
    fseek(fid,bytes,0);
  end
  
  %Read data, channels in columns, data samples in rows
  %fseek(fid,ftell(fid) + 4,-1);
  
  ctf.data{tr} = fread(fid,[samples channels],[num2str(samples),'*int32=>int32'],4*(ctf.setup.number_samples-samples));
  
  % [DLW] Why is this CHAN - min(CHAN)+1?  The next line only uses CHAN.
  ctf.data{tr} = ctf.data{tr}(:,CHAN - min(CHAN)+1);
  
  ctf.data{tr} = double(ctf.data{tr}) * diag(1./channel_gain(CHAN));
  
end


% assign sensor locations and orientations for selected channels, this
% section will simplify the data allocated by ctf_read_res4

fprintf('...sorting %d sensors\n',ctf.setup.number_channels);

ctf.sensor.location = zeros(3,ctf.setup.number_channels);
ctf.sensor.orientation = zeros(3,ctf.setup.number_channels);

j=0;
for channel = CHAN,
  
  % not sure this next line makes sense [DLW]
  j=j+1;
  
  % All channels have a label
  ind = strfind(ctf.sensor.info(channel).label,'-');
  if isempty(ind),
    ctf.sensor.label{1,j} = ctf.sensor.info(channel).label;
  else
    ctf.sensor.label{1,j} = ctf.sensor.info(channel).label(1:ind-1);
  end
  
  % All channels have a location
  
  % EEG channels do not have any orientation
  
  switch ctf.sensor.info(channel).index,
    
    case {ctf.sensor.type.meg, ctf.sensor.type.ref},
      %0=Reference Channels, 
      %1=More Reference Channels, 
      %5=MEG Channels
      
      % MEG channels are radial gradiometers, so they have an inner (1) and
      % an outer (2) location - it might be better to take the average of
      % their locations
      if ~isempty(ctf.sensor.info(channel).location),
        ctf.sensor.location(:,j) = ctf.sensor.info(channel).location(:,1);
      end
      
      if ~isempty(ctf.sensor.info(channel).orientation),
        ctf.sensor.orientation(:,j) = ctf.sensor.info(channel).orientation(:,1);
      end
      
    case ctf.sensor.type.eeg,
      %9=EEG Channels
      
      if ~isempty(ctf.sensor.info(channel).location),
        ctf.sensor.location(:,j) = ctf.sensor.info(channel).location(:,1);
      end
      
  end
end







%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NEED TO CHECK ctf.setup parameters here, to adjust for any changes
% required by the CHAN, TIME, TRIALS inputs
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ctf.setup.number_channels ~= size(ctf.data,2),
  %error('setup.number_channels ~= size(ctf.data,2)');
end





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

return