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);  endelse  if ~exist('ctf','var'),    ctf = ctf_folder(folder);  else    ctf = ctf_folder(folder,ctf);  endendif ~isfield(ctf,'setup'),  ctf = ctf_read_res4(ctf.folder);endif ~exist('CHAN','var'),   CHAN   = 'all'; endif ~exist('TIME','var'),   TIME   = 'all'; endif ~exist('TRIALS','var'), TRIALS = 'all'; endif isempty(CHAN),   CHAN   = 'all'; endif isempty(TIME),   TIME   = 'all'; endif isempty(TRIALS), TRIALS = 'all'; endsmallCHAN = 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 numbersendswitch num2str(TIME),  case 'all',    TIME = ctf.setup.time_sec;  otherwise    % assume the input is a range of timesendswitch num2str(TRIALS),  case 'all',    TRIALS = 1:ctf.setup.number_trials;  otherwise    % assume the input is a range of trialsend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%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 indicesmegIndex =   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 offsetschannel_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 durationduration = 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);endsamples = 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 datactf.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_res4fprintf('...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        endend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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)');endt = toc; fprintf('...done (%6.2f sec)\n\n',t);return