function [F,Channel,imegsens,ieegsens,iothersens,irefsens,grad_order_no,no_trials,filter,Time,RunTitle] = ds2brainstorm(ds_directory,VERBOSE,READRES,CHANNELS, TIME, NO_TRIALS, DCOffset);
%DS2BRAINSTORM - Convert a DS CTF dataset into BrainStorm format
% function [F,Channel,imegsens,ieegsens,iothersens,irefsens,grad_order_no,no_trials,filter,Time,RunTitle] = ds2brainstorm(ds_directory,VERBOSE,READRES,CHANNELS,TIME, NO_TRIALS, DCOffset);
% [F,Channel,imegsens,ieegsens,iothersens,irefsens,grad_order_no,no_trials,filter,Time, RunTitle] 
%                      ... = ds2brainstorm(ds_directory,VERBOSE,READRES,CHANNELS,TIME, NO_TRIALS, DCOffset);
% Reads CTF Systems Inc. Data and Resource file formats (4.1 MEG Data Format)
%
% INPUTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DS_DIRECTORY: full path to the .ds directory containing the orginal data set
%
% VERBOSE = 1 : Toggle VERBOSE mode ON
%               -> Activate screen display
%               -> Save channel information in channel_file.txt
%               -> Save coefficient information in coef_file.txt
%
% VERBOSE = 0 : Toggle VERBOSE mode OFF
%
%---> The rest of the input parameters are optionnal (more specific to data block extraction) 
%
% READRES = 0: Do not read the original .res4 file assuming it has been read previously. 
%              Read  _res4.mat instead. If _res4.mat does not exist, automatic switch to READRES = 1.   
% READRES = 1: Read the original .res4 file and writes out a file called **_res4.mat in the .ds folder, 
%              which is a shortened version of the CTF resources file .res4.
%              MEG/EEG data are also read from the .meg4 file ONLY IF other input arguments are specified 
% READRES = 2: Same as READRES = 1, without the long reading of coefficient information
%              but reads ONLY a short version of data information and returns them in F as a structure.
%              MEG/EEG data are nor read from the .meg4 file.
%              This is useful for display in importdata_gui for instance.
% READRES = 3: Do not read .meg4 file (useful when data need to be kept in native file format when importing into BrainStorm)
%              F is then a string indicating the name of the corresponding .ds folder.
%
%
% Set READRES to 0 when you want to extract another block of data after you have already read one block with READRES = 1.
%
% If READRES is left blank, the whole data are read - ATTENTION when data set is large - possible 'out of memory' issues
%
% If the optional fields below are left blank, ds2brainstorm only reads the res4 file and generates a *_res4.mat file in the current .ds folder 
% with all useful information regarding EEG, MEG channels etc. for subsequent data extraction.
%
% CHANNELS: the indices of the channels to be extracted as in IEEGSENS, IMEGSENS or IREFSENS
% If you don't know what are the indices of these channels, run ds2brainstorm with READRES = 1 first.
%
% TIME: the time window to be extracted in seconds within a trial; ie a subset of Time, a vector saved in **_res4.mat (see above)
%
% NO_TRIALS: a vector containing the number(s) of the trial(s) to be extracted.
%
% DCOffset: A flag indicating whether DC offset needs to be removed from every MEG channel.
% Acceptable values:
%                   0 - DC offset is not removed from MEG channels (DEFAULT)
%                   1 - DC offset is removed from MEG channels based on the entire trial length 
%                   2 - DC offset is removed from MEG channels based on pretigger time period 
%
%  Example; ds2brainstorm('tipouic.ds',1,0,[10:30],[-.5 3],2)
% .. extracts in VERBOSE mode from tipouic.ds, assuming tipouic_res4.mat was already created during a previous call to ds2brainstorm with READRES = 1
% channels 10 to 30 (could be MEG or other -> check IEEGSENS and IMEGSENS first), between -0.5sec and 3sec, in trial 2. 
% 
% %     
% OUTPUTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% CHANNEL - MEG,EEG, REFERENCE and OTHER channel information 
% CHANNEL is a cell array of structures following the BrainStorm data format.
% F - Data Matrix: One sensor per row; one time sample per column / alternatively a structure containing all data information
%    (resources contained in .res4 file).
% IMEGSENS - Indices of the MEG data rows in F
% IEEGSENS - Indices of the EEG data rows in F 
% IOTHERSENS - Indices of the OTHER data rows (not including reference coils) in F
% IREFSENS - Indices of the REFERENCE data rows in F
% GRAD_ORDER_NO - Order of the gradient correction to be taken into account in the forward model only
% NO_TRIALS - Number of trials in the .ds folder.
% FILTER - An array of structures; each structure containing information related to each filter
% TIME - A vector containing the time instants at every data sample
% RUNTITLE - A character array labelling the current run.
% 
% Please refer to CTF Systems Inc. Technical Note #1 for details on DataSet file formats

%<autobegin> ---------------------- 09-Jul-2004 22:17:03 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Data Processing
%
% Alphabetical list of external functions (non-Matlab):
%   toolbox\bst_message_window.m
%   toolbox\findclosest.m
%   toolbox\good_channel.m
%
% Subfunctions in this file, in order of occurrence in file:
%   PrintCoilPos(channel_file, CoilRec_ext,nCoils) ;
%   save_sensor_locs(Channel,SensorLocFile)
%
% At Check-in: $Author: Mosher $  $Revision: 31 $  $Date: 7/09/04 8:42p $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 09-Jul-2004
% 
% Principal Investigators and Developers:
% ** Richard M. Leahy, PhD, Signal & Image Processing Institute,
%    University of Southern California, Los Angeles, CA
% ** John C. Mosher, PhD, Biophysics Group,
%    Los Alamos National Laboratory, Los Alamos, NM
% ** Sylvain Baillet, PhD, Cognitive Neuroscience & Brain Imaging Laboratory,
%    CNRS, Hopital de la Salpetriere, Paris, France
% 
% See BrainStorm website at http://neuroimage.usc.edu for further information.
% 
% Copyright (c) 2004 BrainStorm by the University of Southern California
% This software distributed  under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPL
% license can be found at http://www.gnu.org/copyleft/gpl.html .
% 
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%<autoend> ------------------------ 09-Jul-2004 22:17:03 -----------------------


% /---Script Authors-------------------------------------\
% |                                                      |
% |  *** Sylvain Baillet, Ph.D.                          |
% |  Cognitive Neuroscience & Brain Imaging Laboratory   |
% |  CNRS UPR640 - LENA                                  | 
% |  Hopital de la Salpetriere, Paris, France            |
% |  sylvain.baillet@chups.jussieu.fr                    |
% |                                                      |
% |  *** John C. Mosher, Ph.D.                           |
% |  Biophysics Group                                    |
% |                                                      |
% |  -> with contributions from                          |
% |       Line Garnero, Ph.D. and Antoine Ducorps, Ph.D. |
% |       MEG Center &                                   |
% |       Cognitive Neuroscience                         | 
% |       & Brain Imaging Laboratory - CNRS              |
% |       La Salpetriere Hospital, Paris - France        |
% \------------------------------------------------------/

% Script History ----------------------------------------------------------------------------------------
%
% 11/28/00 - SB   : Updated to read the 3rd-order gradient information
% 03/09/01 - SB   : Data block extraction
% 15/11/01 - SB   : Extraction of the ADC Input channels in iothersens, combined with the STIM channel
% 19/11/01 - SB   : Save iothersens indices in the _res4.mat file
% JCM 06-Jun-2002 : made waitbar update less often, autocomments
% SB 26-Jul-2002  : Added the READRES = 2 condition to read basic file information for display importdata_gui.
% ..............    Removed creation of channel_file.txt. 
% SB 27-Jul-2002  : Minor alterations in the READRES conditions - see help header. 
% SB 28-Jul-2002  : When CHANNELS is emptu, all channels, including reference channels are read.
% ..............    Changed name of sensor location result file (for visualization in MRITool or 3DViewer)
% SB 10-Sep-2002  : Minor layout editing.
% SB 07-Nov-2002  : Display all progress information into bst_message_window
% SB 15-Nov-2002  : EEG channel .loc field is now only 3x1 (was 3x2, second column filled with zeros).
%                   Added DCOffset input argurment to remove DC Offset from MEG channels.
% SB 20-Nov-2002  : Assign "MEG REF" type to MEG reference channels for fast retrieval with GOOD_CHANNEL
%                   Removed .Gcoef and .i*sens fields form Channel structure.
% SB 23-Dec-2002  : Time now returns only the time window on which data has been extracted (not the whole original time span of data) 
% SB 21-Jan-2003  : Fixed 1-sample bug for total Time length
% SB 23-Jan-2003  : Added READRES = 3 option.
% SB 06-Feb-2003  : Fixed bug when accessing data with multiple trials specificied in NO_TRIALS
%                   From the second NO_TRIALS entry on, data was always read from beguinning of trial 
%                   not from the first entry specified in TIME. 
% JCM 10-Dec-2003 : Added additional blank SensorNames to account for type "17" in the 275 channel arrays. Moved around
%                   redundant definitions of SensorNames. Thanks to Fred Carver NIMH for finding this.
% JCM 11-Dec-2003 : Updating SensorNames based on email from CTF, plus correcting 0 indexing scheme
% SB  09-Jul-2004 : Issuing a message when time clipping occurs (requested time window is outside original trial time window)
% --------------------------------------------------------------------------------------------------------

% Check input arguments

if nargin == 2
   READRES = 1;
   DCOffset = 0;
elseif nargin==5 % No Marker file and marker range are defined
   MARKER_FILE = 0;
   MARKER_RANGE = [];
   DCOffset = 0;
elseif nargin==3
   % Read only file information.
   % Check that READRES has acceptable values - if not, the whole file could be read
   % which could end-up in loading very large files into Matlab format 
   % and cause memory crash.
   if READRES ~= 1 | READRES ~= 2
      READRES = 2; % Read data information only.
   end
   MARKER_FILE = 0;
   MARKER_RANGE = [];
   CHANNELS = [];
   TIME = [];
   NO_TRIALS = [];   
   DCOffset = 0;
end


% SensorNames ={'Ref Magnetometer','Ref Gradiometer' ,''  , '' , '' ,'MEG Sensor','' , '', '','EEG Sensor','ADC Input','Stimulation input',''}; % CHEAT - added the 13th '' ofor testing purposes
% 10 December 2003 Fred Carver at NIMH found that 275 channel arrays had a type 17. Added extra blanks to account for.
% 11 December 2003, based on email directly from Dough McKenzie at CTF:
% Comments on types, email of 10 December 2003:
%The new types are
%
%   SAM Sensor       -- Synthetic electrode channel from SAM
%   Virtual Channel  -- Linear combination of sensor channels
%   System Clock     -- 32 bit counter incremented at 12,000
%                       ticks per second for Omega 2000 system
%                       ( new systems like NIMH )
%                       and incremented at 12,500 ticks per
%                       second for older systems.
%
%   Video Time       -- A 32 bit encoded value that can be recorded from
%                       a SONY Video time output from a video player or 
%                       recorder
%
%                       The format is 'hhmmssff'  where each byte contains
%                            hh    hours in decimal ( 0--
%                            mm    minutes          ( 0 ~ 59 )
%                            ss    seconds          ( 0 ~ 59 )
%                            ff    frames           ( 0 ~ 29 )
%
%   ADC Input Voltage -- Used if the ADC channels are measuring
%                        voltage. == 11 if measuring current
%
%NOTE: There will be some additional channels types defined
%      mid next year.
SensorNames ={...
      'Ref Magnetometer',... % Sensor Type Index of 0
      'Ref Gradiometer' ,... % Index of 1
      ''  ,...               % 2
      '' ,...                % 3
      '' ,...                % 4
      'MEG Sensor',...       % 5
      '' ,...                % 6
      '',...                 % 7
      '',...                 % 8
      'EEG Sensor',...       % 9
      'ADC Input Current',...% 10 ADC Input Current (Amps)
      'Stimulation input',...% 11
      'Video Time',...       % 12
      '',...                 % 13
      '',...                 % 14
      'SAM Sensor',...       % 15
      'Virtual Channel',...  % 16
      'System Clock',...     % 17 System Time Ref
      'ADC Input Voltage',...% 18 ADC Input Voltage (Volts)
   };


%% Checking files
cd(ds_directory)
[path,rootname] = fileparts(ds_directory);
meg4file = [rootname,'.meg4'];
rec4file = [rootname,'.res4'];
if isdir('BrainStorm') % Was a BrainStorm directory created in .ds folder ? Yes if ds2brainstorm called from importdata
   res4_mat = fullfile('BrainStorm',[rootname,'_res4.mat']);
else
   res4_mat = [rootname,'_res4.mat'];
end

if ~exist(rec4file,'file') | ~exist(meg4file,'file') 
   errordlg([rec4file ' or ' meg4file ' missing'])
   return
end

% If .mat version of the .res4 does not exist while we supposed it did - read it.
if ~exist(res4_mat,'file') & READRES == 0  
   % Force reading of original .res4 
   READRES = 1;
end

if VERBOSE
   bst_message_window(['Working in...',ds_directory ])
end


%*********************************************    

if READRES > 0 % Read the .res4 file     
   
   %% Define constants
   MAX_COILS = 8;
   MAX_BALANCING = 50;
   SENSOR_LABEL = 31;
   MAX_AVERAGE_BINS = 8;
   
   %***** nfSetUp **** Data Structure
   gSetUp = struct('no_samples',[],'no_channels',[],'sample_rate',[],...
      'epoch_time',[],'no_trials',[],'preTrigPts',[],'no_trials_done',[],'no_trials_display',[],...
      'save_trials',[],'primaryTrigger',[],'secondaryTrigger',[],'triggerPolarityMask',[],...
      'trigger_mode',[],'accept_reject_Flag',[],'run_time_display',[],'zero_Head_Flag',[],...
      'artifact_mode',[]);
   
   nfSetUp = struct('nf_run_name','','nf_run_title','','nf_instruments','','nf_collectdescriptor','',...
      'nf_subject_id','','nf_operator','','nf_sensorFileName','');
   
   GenRes4 = struct('appName','','dataOrigin','','dataDescription','',...
      'no_trials_avgd',[],'data_time',[],'data_date',[],'gSetUp',gSetUp,'nfSetUp',nfSetUp,'rdlen',[]);
   
   %*********************************************
   
   
   %% Reading Resources------------------------------------------------------------------------
   [rec,message] = fopen(rec4file,'rb','s'); % Big-endian byte ordering
   if rec < 0
      errordlg(message)
      return
   end
   
   % Read HEADER
   header = fread(rec,8,'char')';
   if VERBOSE
      bst_message_window({...
            sprintf('ds2brainstorm -> %s file format',char(header))...
         })
   end
   
   
   % Read nfSetUp
   GenRes4.appName = char(fread(rec,256,'char')');
   GenRes4.dataOrigin = char(fread(rec,256,'char')');
   GenRes4.dataDescription = char(fread(rec,256,'char')');
   
   GenRes4.no_trials_avgd = (fread(rec,1,'int16')');
   GenRes4.data_time = char(fread(rec,255,'char')');
   GenRes4.data_date = char(fread(rec,255,'char')');
   
   gSetUp.no_samples = (fread(rec,1,'int32')');
   gSetUp.no_channels = (fread(rec,1,'int16')');
   
   fseek(rec,ceil(ftell(rec)/8)*8,-1);
   gSetUp.sample_rate = (fread(rec,1,'double')');
   fseek(rec,ceil(ftell(rec)/8)*8,-1);
   gSetUp.epoch_time = (fread(rec,1,'double')');
   gSetUp.no_trials =  (fread(rec,1,'int16')');
   fseek(rec,ceil(ftell(rec)/4)*4,-1);
   gSetUp.preTrigPts =  (fread(rec,1,'int32')');
   gSetUp.no_trials_done = (fread(rec,1,'int16')');
   gSetUp.no_trials_bst_message_windowlay = (fread(rec,1,'int16')');
   fseek(rec,ceil(ftell(rec)/4)*4,-1);
   gSetUp.save_trials = (fread(rec,1,'int32')');
   gSetUp.primaryTrigger = char(fread(rec,1,'uchar')');
   gSetUp.secondaryTrigger = char(fread(rec,MAX_AVERAGE_BINS,'uchar')');
   gSetUp.triggerPolarityMask = char(fread(rec,1,'uchar')');
   
   gSetUp.trigger_mode = (fread(rec,1,'int16')');
   fseek(rec,ceil(ftell(rec)/4)*4,-1);
   gSetUp.accept_reject_Flag = (fread(rec,1,'int32')');
   gSetUp.run_time_bst_message_windowlay = (fread(rec,1,'int16')');
   fseek(rec,ceil(ftell(rec)/4)*4,-1);
   gSetUp.zero_Head_Flag = (fread(rec,1,'int32')');
   fseek(rec,ceil(ftell(rec)/4)*4,-1);
   gSetUp.artifact_mode = (fread(rec,1,'int32')');
   gSetUp.padding = (fread(rec,1,'int32')');
   
   
   nfSetUp.nf_run_name = char(fread(rec,32,'char')');
   nfSetUp.nf_run_title = char(fread(rec,256,'char')');
   
   RunTitle = nfSetUp.nf_run_title ;
   
   nfSetUp.nf_instruments = char(fread(rec,32,'char')');
   nfSetUp.nf_collect_descriptor = char(fread(rec,32,'char')');
   nfSetUp.nf_subject_id = char(fread(rec,32,'char')');
   nfSetUp.nf_operator = char(fread(rec,32,'char')');
   tmp = fread(rec,60,'char')';
   tmp(tmp<0) = 0; % prevent out of range character conversion warning
   nfSetUp.nf_sensorFileName = tmp; clear tmp
   
   fseek(rec,ceil(ftell(rec)/4)*4,-1);
   nfSetUp.rdlen = fread(rec,1,'int32')';
   
   GenRes4.gSetUp = gSetUp;
   GenRes4.nfSetUp = nfSetUp;
   
   
   % Time span of data record
   Time =  linspace(-GenRes4.gSetUp.preTrigPts/GenRes4.gSetUp.sample_rate,...
      GenRes4.gSetUp.epoch_time/GenRes4.gSetUp.no_trials-(GenRes4.gSetUp.preTrigPts+1)/GenRes4.gSetUp.sample_rate,...
      GenRes4.gSetUp.no_samples);
   %     
   %----------------------------------------------------------------------------------------
   
   if VERBOSE
      ON_OFF = {'Off','On'};
      bst_message_window({...
            sprintf('ds2brainstorm -> Header information'),...
            sprintf('                 Collected %s starting at %s', GenRes4.data_date, GenRes4.data_time),...
            sprintf('                 Run Name: %s', nfSetUp.nf_run_name),...
            sprintf('                 Run Title: %s', nfSetUp.nf_run_title),...
            sprintf('                 Col Desc: %s', nfSetUp.nf_collect_descriptor),...
            sprintf('                 Run Desc: %s', GenRes4.dataDescription),...
            sprintf('                 Operator: %s', nfSetUp.nf_operator),...
            sprintf('                 Subject : %s', nfSetUp.nf_subject_id),...
            sprintf('                 Channels: %d',gSetUp.no_channels),...
            sprintf('                 Samples : %d per trial', gSetUp.no_samples),...
            sprintf('                 Rate    : %g samples/sec', gSetUp.sample_rate),...
            sprintf('                 Trials  : %d (average of %d)', gSetUp.no_trials, GenRes4.no_trials_avgd),...
            sprintf('                 Duration: %g seconds/trial', gSetUp.epoch_time/gSetUp.no_trials),...
            sprintf('                 Pre-trig: %g samples', gSetUp.preTrigPts),...
            sprintf('                 Sensor file name : %s', nfSetUp.nf_sensorFileName),...
            sprintf('                 Head zeroing: %s', ON_OFF{gSetUp.zero_Head_Flag+1}),...
            sprintf('_________________')...
         })
   end
   
   
   %-------------------------------------------------------------------------------
   
   %----------------------------------READ FILTERS---------------------------------
   
   fseek(rec,1844,-1);
   
   % Run Description
   rundescript = char(fread(rec,nfSetUp.rdlen,'char'));
   if VERBOSE
      bst_message_window({...
            sprintf('ds2brainstorm -> Run Description: %s', rundescript )...
         })  
   end
   
   classType = {'CLASSERROR','BUTTERWORTH'};
   filtType = {'TYPERROR','LOWPASS','HIGHPASS','NOTCH'};
   
   
   % Number of filters
   no_filters = fread(rec,1,'int16');
   if VERBOSE
      bst_message_window({...
            'ds2brainstorm -> Filter information',...
            sprintf('Number of filters: %d', no_filters)...
         })
   end
   
   % JCM fixes 10/26/00, each of these should be looped through, not read in bulk
   % old
   if(0)
      filter.freq = fread(rec,no_filters,'double');
      filter.fClass = classType{fread(rec,no_filters,'int32')+1};
      filter.fType = filtType{fread(rec,no_filters,'int32')+1};
      filter.numParam = fread(rec,no_filters,'int16');
      
      filter.params = [];
      for fi = 1:no_filters
         for p = 1: filter(fi).numParam 
            filter(fi).params= fread(rec,filter(fi).numParam,'double');
         end
      end
   else % new
      [filter(1:no_filters)] = struct('freq',[],'fClass',[],'fType',[],'numParam',[],'params',[]);
      for fi = 1:no_filters,
         filter(fi).freq = fread(rec,1,'double');
         filter(fi).fClass = classType{fread(rec,1,'int32')+1};
         filter(fi).fType = filtType{fread(rec,1,'int32')+1};
         filter(fi).numParam = fread(rec,1,'int16');
         filter(fi).params= fread(rec,filter(fi).numParam,'double');
      end
   end
   
   if VERBOSE
      for fi = 1:no_filters
         bst_message_window({...
               sprintf('  Filter - %d',fi),...
               sprintf('         -> Frequency: %g Hz',filter(fi).freq),...
               sprintf('         -> Class: %s',filter(fi).fClass),...
               sprintf('         -> Type: %s',filter(fi).fType),...
               sprintf('         -> Number of parameters: %d',filter(fi).numParam)...
            })
         
         if ~isempty(filter(fi).params)
            bst_message_window({...
                  sprintf('         -> Parameter Value(s): %g',filter(fi).params)...
               })
         end
      end 
      
      bst_message_window({...
            sprintf('ds2brainstorm -> Reading Filter Information - DONE'),...
            sprintf('ds2brainstorm -> Reading Channel Information. . .')...
         })
   end
   
   % Channel Names
   for chan = 1:gSetUp.no_channels 
      channel_name{chan} = fread(rec,32,'char')'; 
      tmp = channel_name{chan};
      tmp(tmp>127) = 0; 
      tmp(tmp<0) = 0;
      channel_name{chan} = strtok(tmp,char(0));
      ChannelName{chan} = char(strtok(channel_name{chan},'-'));
   end
   
   % Sensor Resources
   CoilType = {'CIRCULAR','SQUARE','???'};
   
   for chan = 1:gSetUp.no_channels 
      oldtell = ftell(rec);
      SensorRes(chan).sensorTypeIndex = fread(rec,1,'int16');
      SensorRes(chan).originalRunNum = fread(rec,1,'int16');
      
      id = fread(rec,1,'int32')+1;
      
      if isempty(id)
         id = -1;
      end
      
      if id > 3 | id <0 
         id = 3;
      end
      
      SensorRes(chan).coilShape = CoilType{id};
      
      SensorRes(chan).properGain = fread(rec,1,'double');
      SensorRes(chan).qGain = fread(rec,1,'double');
      SensorRes(chan).ioGain= fread(rec,1,'double');
      SensorRes(chan).ioOffset = fread(rec,1,'double');
      SensorRes(chan).numCoils = fread(rec,1,'int16');
      SensorRes(chan).grad_order_no = fread(rec,1,'int16');
      
      padding = fread(rec,1,'int32'); 
      
      for coil = 1:MAX_COILS
         SensorRes(chan).coilTbl(coil).position.x = fread(rec,1,'double');
         SensorRes(chan).coilTbl(coil).position.y = fread(rec,1,'double');
         SensorRes(chan).coilTbl(coil).position.z = fread(rec,1,'double');
         SensorRes(chan).coilTbl(coil).position.junk = fread(rec,1,'double');
         SensorRes(chan).coilTbl(coil).orient.x = fread(rec,1,'double');
         SensorRes(chan).coilTbl(coil).orient.y = fread(rec,1,'double');
         SensorRes(chan).coilTbl(coil).orient.z = fread(rec,1,'double');
         SensorRes(chan).coilTbl(coil).orient.junk = fread(rec,1,'double');
         SensorRes(chan).coilTbl(coil).numturns = fread(rec,1,'int16');
         padding = fread(rec,1,'int32');
         padding = fread(rec,1,'int16');
         SensorRes(chan).coilTbl(coil).area = fread(rec,1,'double');
      end
      
      for coil = 1:MAX_COILS
         SensorRes(chan).HdcoilTbl(coil).position.x = fread(rec,1,'double');
         SensorRes(chan).HdcoilTbl(coil).position.y = fread(rec,1,'double');
         SensorRes(chan).HdcoilTbl(coil).position.z = fread(rec,1,'double');
         SensorRes(chan).HdcoilTbl(coil).position.junk = fread(rec,1,'double');
         SensorRes(chan).HdcoilTbl(coil).orient.x = fread(rec,1,'double');
         SensorRes(chan).HdcoilTbl(coil).orient.y = fread(rec,1,'double');
         SensorRes(chan).HdcoilTbl(coil).orient.z = fread(rec,1,'double');
         SensorRes(chan).HdcoilTbl(coil).orient.junk = fread(rec,1,'double');
         SensorRes(chan).HdcoilTbl(coil).numturns = fread(rec,1,'int16');
         padding = fread(rec,1,'int32'); 
         padding = fread(rec,1,'int16'); 
         SensorRes(chan).HdcoilTbl(coil).area = fread(rec,1,'double');
      end
      
   end
   
   if VERBOSE
      bst_message_window({...
            'ds2brainstorm -> Reading Channel Information - DONE',...
            'ds2brainstorm -> Converting Channel Information. . .'...
         })
      
   end
   
   nchan= length(channel_name); % Number of channels
   
   % Explode data according to channel type
   imegsens = find([SensorRes.sensorTypeIndex] == 5); % Indices of MEG sensors
   ieegsens = find([SensorRes.sensorTypeIndex] == 9); % Indices of EEG sensors
   iothersens = [find([SensorRes.sensorTypeIndex] == 11),find([SensorRes.sensorTypeIndex] == 10)]; % Indices of OTHER sensors ('Stimulation input' and 'ADC')
   irefsens = find([SensorRes.sensorTypeIndex] == 0); % Reference Channels
   irefsens = [irefsens,find([SensorRes.sensorTypeIndex] == 1)]; % Reference Channels
   
   if READRES == 2 % Enough data information - send it out of here
      % Initialize output arguments 
      [F,Channel,grad_order_no,no_trials,Time, RunTitle] = deal([]);
      F = GenRes4;
      F.filter = filter;
      F.imegsens = imegsens;
      F.ieegsens = ieegsens;
      F.irefsens = irefsens;
      F.iothersens = iothersens;
      F.ChannelNames = char(ChannelName);
      F.grad_order_no = [SensorRes(:).grad_order_no];        
      F.nchan = length([imegsens,ieegsens,irefsens,iothersens]);
      F.Time = Time;
      return
   end
   
   ieeg = 0; imeg = 0; 
   iother = 0; 
   iother2 = 0; 
   eegID = [];
   otherID = [];
   other2ID = [];
   megID = [];
   
   iref = 0;
   refID = [];
   
   chaneeg = [];% Index of EEG channels 
   chanmeg = [];chanother = [];
   
   [Channel(1:nchan)] = deal(struct('Loc',[],'Orient',[],'Comment',[],'Weight',[],'Type',[],'Name',''));
   
   for chan = 1:nchan
      switch SensorRes(chan).sensorTypeIndex
      case {0,1} %'References coils   
         iref = iref+1;
         MAX_COILS = 2; % Take only into account the 2 first coils for every sensor 
         tmp = [SensorRes(chan).HdcoilTbl(:).position];              
         tmpx  = [tmp(1:MAX_COILS).x]/100; % In meters
         tmpy  = [tmp(1:MAX_COILS).y]/100;
         tmpz  = [tmp(1:MAX_COILS).z]/100;
         
         Channel(chan).Loc = ...
            [tmpx;tmpy;tmpz];
         
         tmp = [SensorRes(chan).HdcoilTbl(:).orient];
         tmpx  = [tmp(1:MAX_COILS).x];
         tmpy  = [tmp(1:MAX_COILS).y];
         tmpz  = [tmp(1:MAX_COILS).z];
         
         Channel(chan).Orient = ...
            [tmpx;tmpy;tmpz];
         
         
         % Reorient along the outward pointing normal
         loc =  Channel(chan).Loc';
         ps1=sum((loc(1,:).*Channel(chan).Orient(:,1)')')';
         ps2=sum((loc(2,:).*Channel(chan).Orient(:,2)')')';
         Channel(chan).Orient(:,1) = (sign(ps1)*ones(1,3).*Channel(chan).Orient(:,1)')';
         Channel(chan).Orient(:,2) = (sign(ps2)*ones(1,3).*Channel(chan).Orient(:,2)')';
         
         Channel(chan).Comment = ''; % SB 20-Nov-2002
         Channel(chan).Name = [SensorNames{SensorRes(chan).sensorTypeIndex+1},' ',int2str(chan)];% SB 20-Nov-2002
         Channel(chan).Type = 'MEG REF'; % SB 20-Nov-2002
         
         % Channel(chan).Type = SensorNames{SensorRes(chan).sensorTypeIndex+1} ;
         %             if isempty(deblank(Channel(chan).Type))
         %                 Channel(chan).Type = 'MEG REF';
         %             end
         
         chanmeg = [chanmeg,chan];
         refID(iref) = chan; 
         Channel(chan).Weight = [1 -1] ;  
         
      case 5 % MEG
         imeg = imeg+1;
         MAX_COILS = 2; % Take only into account the 2 first coils for every sensor 
         tmp = [SensorRes(chan).HdcoilTbl(:).position];              
         tmpx  = [tmp(1:MAX_COILS).x]/100; % In meters
         tmpy  = [tmp(1:MAX_COILS).y]/100;
         tmpz  = [tmp(1:MAX_COILS).z]/100;
         
         Channel(chan).Loc = ...
            [tmpx;tmpy;tmpz];
         
         tmp = [SensorRes(chan).HdcoilTbl(:).orient];
         tmpx  = [tmp(1:MAX_COILS).x];
         tmpy  = [tmp(1:MAX_COILS).y];
         tmpz  = [tmp(1:MAX_COILS).z];
         
         Channel(chan).Orient = ...
            [tmpx;tmpy;tmpz];
         
         % Reorient along the outward pointing normal
         loc =  Channel(chan).Loc';
         ps1=sum((loc(1,:).*Channel(chan).Orient(:,1)')')';
         ps2=sum((loc(2,:).*Channel(chan).Orient(:,2)')')';
         Channel(chan).Orient(:,1) = (sign(ps1)*ones(1,3).*Channel(chan).Orient(:,1)')';
         Channel(chan).Orient(:,2) = (sign(ps2)*ones(1,3).*Channel(chan).Orient(:,2)')';
         
         Channel(chan).Comment = ' ';%[SensorNames{SensorRes(chan).sensorTypeIndex+1},' ',int2str(chan)];
         Channel(chan).Type = 'MEG';%SensorNames{SensorRes(chan).sensorTypeIndex+1} ;
         
         chanmeg = [chanmeg,chan];
         megID(imeg) = chan; 
         Channel(chan).Weight = [1 -1] ;  
         Channel(chan).Name = ChannelName{chan}; %char(strtok(channel_name{chan},'-'));
         
      case 9 % EEG
         ieeg = ieeg + 1;
         MAX_COILS = 1; % Take only into account the first sensor location for EEG (second pseudo-coil location is [0 0 0] as writtent in .ds files)
         tmp = [SensorRes(chan).HdcoilTbl(:).position]; 
         tmpx  = [tmp(1:MAX_COILS).x]/100;% In meters
         tmpy  = [tmp(1:MAX_COILS).y]/100;
         tmpz  = [tmp(1:MAX_COILS).z]/100;
         
         Channel(chan).Loc = ...
            [tmpx;tmpy;tmpz];
         
         Channel(chan).Orient = [];
         
         Channel(chan).Type = 'EEG';%SensorNames{SensorRes(chan).sensorTypeIndex+1} ;
         Channel(chan).Comment = ' ';%[SensorNames{SensorRes(chan).sensorTypeIndex+1},' ',int2str(chan)];
         Channel(chan).Weight = [];   
         Channel(chan).Name = ChannelName{chan}; %char(strtok(channel_name{chan},'-'));%int2str(ieeg) ;  
         chaneeg = [chaneeg,chan];
         eegID(ieeg) = chan; 
         
      case 11 % STIM
         iother = iother+1;
         Channel(chan).Loc = [];
         Channel(chan).Orient = [];
         Channel(chan).Type = SensorNames{SensorRes(chan).sensorTypeIndex+1} ;
         if isempty(deblank(Channel(chan).Type ))
            Channel(chan).Type = 'STIM';
         end
         Channel(chan).Comment = [SensorNames{SensorRes(chan).sensorTypeIndex+1},' ',int2str(chan)];
         Channel(chan).Weight = [];   
         chanother = [chanother,chan];
         otherID(iother) = chan;
         
      otherwise
         iother2 = iother2+1;
         Channel(chan).Loc = [];
         Channel(chan).Orient = [];
         Channel(chan).Type = SensorNames{SensorRes(chan).sensorTypeIndex+1} ;
         if isempty(deblank(Channel(chan).Type ))
            Channel(chan).Type = 'OTHER';
         end
         Channel(chan).Comment = [SensorNames{SensorRes(chan).sensorTypeIndex+1},' ',int2str(chan)];
         Channel(chan).Weight = [];   
         chanother = [chanother,chan];
         other2ID(iother) = chan;
         
      end
   end
   
   if VERBOSE
      bst_message_window('ds2brainstorm -> Converting Channel Information  - DONE')
   end
   %-------------------------------------------------------------------------------
   
   
   %----------------------------------READ Coefficients---------------------------------
   
   % Number of coefficient records
   nrec = fread(rec,1,'int16');
   
   channel_name = cellstr(char(channel_name{:}));
   hexadef = {'00000000','47314252','47324252','47334252','47324f49','47334f49'};
   strdef = {'NOGRAD','G1BR','G2BR','G3BR','G2OI','G3OI'};
   
   CoefInfo = cell(length(channel_name),length(strdef)-1);
   
   %SensorCoefResRec
   if VERBOSE
      h = waitbar(0,'Reading Coefficient Information...');
   end
   
   for k = 1:nrec
      SensorCoefResRec(k).sensorName = char(fread(rec,32,'char')');
      SensorCoefResRec(k).coefType = (fread(rec,1,'bit32'));
      padding = fread(rec,1,'int32');
      SensorCoefResRec(k).CoefResRec.num_of_coefs = fread(rec,1,'int16');
      
      SensorCoefResRec(k).CoefResRec.sensor_list = char(fread(rec,[SENSOR_LABEL,MAX_BALANCING],'uchar')');  
      SensorCoefResRec(k).CoefResRec.coefs_list = fread(rec,MAX_BALANCING,'double');
      
      cchannel = strmatch( deblank(SensorCoefResRec(k).sensorName),channel_name);
      CoefType = find(hex2dec(hexadef)==(SensorCoefResRec(k).coefType));
      
      if CoefType > 0
         
         CoefInfo{cchannel,CoefType-1}.num_of_coefs =  SensorCoefResRec(k).CoefResRec.num_of_coefs;
         
         for i=1:CoefInfo{cchannel,CoefType-1}.num_of_coefs 
            CoefInfo{cchannel,CoefType-1}.sensor_list(i) = irefsens(strmatch(strtok(SensorCoefResRec(k).CoefResRec.sensor_list(i,:),char(0)),channel_name(irefsens)));   
            CoefInfo{cchannel,CoefType-1}.coefs(i) = SensorCoefResRec(k).CoefResRec.coefs_list(i);
         end
      end  
      if VERBOSE
         if(~rem(k,floor(nrec/10))), % only occasionally
            waitbar(k/nrec);
         end
      end
      
   end
   
   if VERBOSE
      delete(h)
   end
   
   % Channel Gains
   gain_chan = zeros(size(channel_name,1),1);
   gain_chan(imegsens) = ([SensorRes(imegsens).properGain]'.*[SensorRes(imegsens).qGain]');
   gain_chan(irefsens) = ([SensorRes(irefsens).properGain]'.*[SensorRes(irefsens).qGain]');
   %  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
   
   if VERBOSE
      bst_message_window({...
            'ds2brainstorm -> Coefficient Information',...
            sprintf('                 Number of coefficient records = %d ', nrec)...
         })
      
      %         [coef_file, message] = fopen('coef_file.txt','wt+');
      %         
      %         if coef_file < 0
      %             errordlg(message)
      %             return
      %         end
      %         
      %         fprintf(coef_file,'Number of coefficient records = %d \n', nrec);
      %         
      %         for k = 1:nrec
      %             coefType = dec2hex(SensorCoefResRec(k).coefType);
      %             coefType = strdef{hex2dec(hexadef)==hex2dec(coefType)};
      %             
      %             fprintf(coef_file,'Channel %s type %s, ',SensorCoefResRec(k).sensorName, coefType);
      %             
      %             num_of_coefs = SensorCoefResRec(k).CoefResRec.num_of_coefs;
      %             fprintf(coef_file,'number = %d:\n', num_of_coefs);
      %             if 0
      %                 for j=1:num_of_coefs
      %                     fprintf(coef_file,' %s ', SensorCoefResRec(k).CoefResRec.sensor_list{j});
      %                     fprintf(coef_file,' %g ',SensorCoefResRec(k).CoefResRec.coefs_list(j));
      %                 end
      %                 fprintf(coef_file,'\n');
      %             end
      %         end
      %         
      %         bst_message_window({...
      %                 sprintf('ds2brainstorm -> Coefficient Information saved in coef_file.txt'),...        
      %                 sprintf('ds2brainstorm -> Reading Coefficient Information - DONE'),...
      %                 ' '...
      %             })
      bst_message_window({...
            sprintf('ds2brainstorm -> Reading Coefficient Information - DONE'),...
            ' '...
         })
   end
   
   
   % 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(imegsens),length(min(irefsens):max(irefsens)));
   grad_order_no = [SensorRes(:).grad_order_no];        
   for k = 1:length(imegsens)
      
      % Reference coils for channel k
      if grad_order_no(imegsens(k)) == 0
         %Data is saved as RAW
         %Save 3rd order gradient sensor-list for subsequent correction if requested later by the user
         [refs] = (CoefInfo{imegsens(k),3}.sensor_list);
         Gcoef(k,refs-min(irefsens)+1) = CoefInfo{imegsens(k),3}.coefs ... 
            .* gain_chan(refs)'/gain_chan(imegsens(k)); 
      else
         [refs] = (CoefInfo{imegsens(k),grad_order_no(imegsens(k))}.sensor_list);
         Gcoef(k,refs-min(irefsens)+1) = CoefInfo{imegsens(k),grad_order_no(imegsens(k))}.coefs ... 
            .* gain_chan(refs)'/gain_chan(imegsens(k)); 
      end
      
   end
   
   Channel(imegsens(1)).Comment = Gcoef;
   %     Channel(imegsens(1)).imegsens = imegsens;
   %     Channel(imegsens(1)).ieegsens = imegsens;
   %     Channel(imegsens(1)).irefsens = irefsens;
   %     Channel(imegsens(1)).RefChannel = Channel; % Store MEG reference information for future call to MEG forward routines and proper processing of gradient correction
   
   SensorLocFile = strrep(ds_directory,[fileparts(ds_directory),filesep],'');
   SensorLocFile = strrep(SensorLocFile,'.ds','_SensorLoc_Results.mat');
   save_sensor_locs(Channel,SensorLocFile)
   
   if VERBOSE
      bst_message_window('wrap',...
         sprintf('ds2brainstorm -> Sensor locations can be visualized with the MRI Tool by loading the file: %s',SensorLocFile)...
         )
   end
   
   fclose('all');
   
   % Create the BrainStorm res4 file for subsequent fast access to the data file
   no_channels = length(Channel);
   gain_chan(ieegsens)=1./gain_chan(ieegsens);
   no_trials = [1:gSetUp.no_trials];
   
   save(res4_mat,'gSetUp','meg4file','ieegsens','irefsens', 'iothersens', 'imegsens','Time','no_channels','gain_chan','channel_name','Channel','grad_order_no','filter','RunTitle','no_trials');
   
   if READRES == 1
      READRES = 0; %read the res4.mat file just created and read the data by block.
      if nargin == 2 % Stop there
         return
         %CHANNELS = 1:length(Channel); % With all available channels...
         %TIME = Time; % whole time duration...
         %NO_TRIALS = 1:gSetUp.no_trials; % All trials
      elseif nargin > 2
         % Fill empty arguments
         if isempty(CHANNELS)
            CHANNELS = 1:length(Channel); % Read all available channels
         end
         if isempty(TIME)
            TIME = Time; % Read entire trial duration 
         end
         if isempty(NO_TRIALS)
            NO_TRIALS = 1:gSetUp.no_trials; % Read all trials
         end
      end
   end
   
else % READRES == 0
   if VERBOSE
      bst_message_window({...
            sprintf('ds2brainstorm -> Loading resources : %s',res4_mat),...
         })
   end
   load(res4_mat);
end


%--------------------------- Readind DATA FILE ------------------------------

if VERBOSE
   bst_message_window({...
         sprintf('ds2brainstorm -> Reading Data')...
      })
end

if nargin == 2 & READRES ~= 3 % Read entire block of data
   
   no_trials = gSetUp.no_trials;
   
   [meg,message] = fopen(meg4file,'rb','s'); % Big-endian byte ordering
   if meg < 0
      errordlg(message)
      return
   end
   
   F = cell(gSetUp.no_trials,1); % allocate space for the data array % ASSUME SINGLE TRIAL FOR THE MOMENT - ie enough memory space
   header = char(fread(meg,8,'char')');
   
   nref=length(irefsens);
   
   for trial = 1:gSetUp.no_trials
      
      if VERBOSE
         bst_message_window({...
               sprintf('ds2brainstorm -> Trial %d/%d',trial,gSetUp.no_trials)...
            })
      end
      
      F{trial} = zeros(gSetUp.no_channels,gSetUp.no_samples);
      F{trial} = fread(meg,[gSetUp.no_samples gSetUp.no_channels],'int32')';
      
      if VERBOSE
         bst_message_window({...
               sprintf('ds2brainstorm -> Done')...
            })
      end
      % IMPORTANT NOTICE : Applying Gradient Correction
      % Data are saved in a given nth-order gradient correction
      % Applying gradient correction si only needed for forward model computation
      % or if it is desired to reverse to lower-order gradient correction (see importdata.m for instance).
      
      % Apply channel gains
      F{trial}(ieegsens,:) = diag(1./gain_chan(ieegsens))*F{trial}(ieegsens,:) ; %To get Microvolts->already inverted when saved _res4.mat
      F{trial}(imegsens,:) = diag(1./gain_chan(imegsens))*F{trial}(imegsens,:)  ;
      
      % Removing DC Offset
      switch(DCOffset)
      case {1,2}
         if VERBOSE
            bst_message_window({...
                  'ds2brainstorm -> Removing DC offset'...
               })
         end
         
         if DCOffset == 1 % Based on entire trial length
            F{trial}(imegsens,:) = F{trial}(imegsens,:) - repmat(mean(F{trial}(imegsens,:)')',1,size(F{trial},2));
            F{trial}(irefsens,:) = F{trial}(irefsens,:) - repmat(mean(F{trial}(irefsens,:)')',1,size(F{trial},2));    
         else % Based on pretrigger
            TimeNeg = TIME(TIME<0); % Find pretrigger time points
            if ~isempty(TimeNeg)
               F{trial}(imegsens,:) = F{trial}(imegsens,:) - repmat(mean(F{trial}(imegsens,TimeNeg)')',1,size(F{trial},2));
               F{trial}(irefsens,:) = F{trial}(irefsens,:) - repmat(mean(F{trial}(irefsens,TimeNeg)')',1,size(F{trial},2));
            else
               % Do nothing
            end
         end
         
         if VERBOSE
            bst_message_window({...
                  'ds2brainstorm -> Removing DC offset - DONE'...
               })
         end
         
      otherwise
         % Do nothing
      end
      
      F{trial}(irefsens,:) = diag(1./gain_chan(irefsens))*F{trial}(irefsens,:) ;
      F{trial}(iothersens,:) = diag(1./gain_chan(iothersens))*F{trial}(iothersens,:) ;
      
   end
   
   fclose('all');
   
   %------------------------------------------------------------------------------------------------------------------------------------
   
else % Reads sub-block of data
   
   
   if READRES == 3
      
      F = ds_directory;
      
      if isempty(NO_TRIALS) % argument NO_TRIALS not defined
         no_trials = [1:gSetUp.no_trials];
      else
         no_trials = [NO_TRIALS];
      end
      
   else
      
      %ds2brainstorm(ds_directory,VERBOSE,READRES,CHANNELS,TIME, NO_TRIALS);
      if isempty(NO_TRIALS) % argument NO_TRIALS not defined
         no_trials = [1:gSetUp.no_trials];
      else
         no_trials = [NO_TRIALS];
      end
      
      if isempty(CHANNELS)
         CHANNELS = [1:length(Channel)];
      end
      
      fclose('all');
      [meg,message] = fopen(meg4file,'rb','s'); % Big-endian byte ordering
      if meg < 0
         errordlg(message)
         return
      end
      
      % no_samples = round((TIME(end)-TIME(1))*gSetUp.sample_rate) +1; % Number of time samples in a trial
      
      F = cell(length(no_trials),1); % allocate space for the data array % ASSUME SINGLE TRIAL FOR THE MOMENT - ie enough memory space
      % Adjust time range to the closest time samples;
      if isempty(TIME)
         return
      end
      
      header = char(fread(meg,8,'char')');
      implicit.sample_rate = 1/(Time(2) - Time(1)); % i.e. the one used in the data file given the time begin_time end period.
      
      if ((min(TIME) < Time(1)) |(max(TIME) > Time(end))) & VERBOSE
          bst_message_window('wrap',...
              'ds2brainstorm -> Requested time window outside origial trial time edges; some clipping will occur.')
      end
      
      
      tmp = findclosest([TIME],Time);
      t_in = tmp(1);
      t_out = tmp(2);
      no_samples = length(t_in:t_out);
      
      try 
         Time = Time(t_in:t_out); % SB 23-Dec-2002
      catch 
         error(...
            sprintf('Data time ranges from %3.1f to %3.2f ms. Please adjust time extraction window',Time(1),Time(end)));
      end
      
      diff_trials = diff(no_trials)-1; % Number of trials between each selected trials (useful when skipping a few trials)
      
      ByteSizeOfTrial= no_channels*gSetUp.no_samples*4; % Byte size of a single trial of data (Int32 coding) 
      
      samples_skip = gSetUp.no_samples-no_samples; %Number of time samples to skip per channel
      
      LastChannelSkip = (no_channels - max(CHANNELS))*gSetUp.no_samples*4; % Skip data from last channels
      
      channels = [min(CHANNELS):max(CHANNELS)]; % Block of channels to extract.
      
      FirstChannelSkip = (min(CHANNELS)-1)*gSetUp.no_samples*4; % Skip data from first channels
      no_channels = length(channels);    
      
      itrial = 0;
      
      for trial = no_trials
         itrial = itrial+1;
         if VERBOSE
            bst_message_window({...
                  sprintf('ds2brainstorm -> Trial %d / %d',itrial,length(no_trials))...
               })
         end
         
         F{itrial} = zeros(no_channels,no_samples);
         
         if trial == no_trials(1) % Read first trial
            fseek(meg,(trial-1)*ByteSizeOfTrial + FirstChannelSkip + (t_in-1)*4 ,0);
         else % just shift from the size of a trial
            fseek(meg,LastChannelSkip + diff_trials(itrial-1)*ByteSizeOfTrial + FirstChannelSkip + (t_in-1)*4 ,0);
         end
         
         F{itrial} = fread(meg,[no_samples no_channels],[num2str(no_samples),'*int32=>int32'], samples_skip*4)';
         F{itrial} = F{itrial}(CHANNELS-min(CHANNELS)+1,:);
         
         % IMPORTANT NOTICE : Applying Gradient Correction
         % Data are saved in a given nth-order gradient correction
         % Applying gradient correction is only needed for forward model computation
         % or if it is desired to reverse to lower-order gradient correction (see importdata.m for instance).
         
         % Apply channel gains
         tmp = gain_chan(CHANNELS);
         tmp(tmp == 0) = eps;
         gain_chan(CHANNELS) = tmp; clear tmp
         F{itrial} = diag(1./gain_chan(CHANNELS))*double(F{itrial});
         
         % Removing DC Offset
         if isempty(DCOffset)
            DCOffset = 0;
         end
         switch(DCOffset)
         case {1,2}
            if VERBOSE
               bst_message_window({...
                     'ds2brainstorm -> Removing DC offset'...
                  })
            end
            
            if DCOffset == 1 % Based on entire trial length
               F{itrial}(imegsens,:) = F{itrial}(imegsens,:) - repmat(mean(F{itrial}(imegsens,:)')',1,size(F{itrial},2));
               F{itrial}(irefsens,:) = F{itrial}(irefsens,:) - repmat(mean(F{itrial}(irefsens,:)')',1,size(F{itrial},2));
            else % Based on pretrigger
               TimeNeg = intersect(find(Time(t_in:t_out) >= TIME(1)),find(Time(t_in:t_out)<0)); % Find pretrigger time points
               if ~isempty(TimeNeg)
                  F{itrial}(imegsens,:) = F{itrial}(imegsens,:) - repmat(mean(F{itrial}(imegsens,TimeNeg)')',1,size(F{itrial},2));
                  F{itrial}(irefsens,:) = F{itrial}(irefsens,:) - repmat(mean(F{itrial}(irefsens,TimeNeg)')',1,size(F{itrial},2));
               else
                  % Do nothing
               end
            end
            
            if VERBOSE
               bst_message_window({...
                     'ds2brainstorm -> Removing DC offset - DONE'...
                  })
            end
            
         otherwise
            % Do nothing
         end
         
         
         if VERBOSE
            bst_message_window({...
                  sprintf('ds2brainstorm -> Done')...
               })
         end
         
      end
      
      fclose('all');
   end
   
end


function PrintCoilPos(channel_file, CoilRec_ext,nCoils) ; 

% Subroutine to DS2BRAINSTORM 
% Formatted output of channel locations in a file (channel_file is the pointer to this file)
% as defined in theCoilRec_ext structure, according to the number of coils nCoils for every gradiometer

for coil=1:nCoils
   fprintf(channel_file,'Coil %d:', coil);
   fprintf(channel_file,' pX %g,', CoilRec_ext(coil).position.x);
   fprintf(channel_file,' pY %g,', CoilRec_ext(coil).position.y);
   fprintf(channel_file,' pZ %g ', CoilRec_ext(coil).position.z);
   fprintf(channel_file,'-');
   fprintf(channel_file,' oX %g,', CoilRec_ext(coil).orient.x);
   fprintf(channel_file,' oY %g,', CoilRec_ext(coil).orient.y);
   fprintf(channel_file,' oZ %g ', CoilRec_ext(coil).orient.z);
   fprintf(channel_file,'-');
   fprintf(channel_file,' turns %d,', CoilRec_ext(coil).numturns);
   fprintf(channel_file,' area %g', CoilRec_ext(coil).area);
   
   if coil < nCoils
      fprintf(channel_file,' --- ');
   end
   
   fprintf(channel_file,'\n');
end

return

function  save_sensor_locs(Channel,SensorLocFile)
% Save Channel location in a pseudo-result file for visualization in the MRITool

nchan = length(Channel);
meg = good_channel(Channel,ones(nchan,1),'MEG');
eeg = good_channel(Channel,ones(nchan,1),'EEG');

nchan= length([meg,eeg]);
SourceLoc = cell(1,nchan);

i = 0;
for k = [meg,eeg]
   i = i+1;
   SourceLoc{i} = Channel(k).Loc(:,1);
   SourceOrder(i) = -1;
   Comment = 'Sensor Locations';
end
DataFlag = 'Sensors';
if isdir('BrainStorm')
   save(fullfile('BrainStorm',SensorLocFile),'SourceLoc','SourceOrder','DataFlag','Comment')
else
   save sensor_result SourceLoc SourceOrder DataFlag Comment
end

return
