meg4read

(Toolbox/meg4read.m in BrainStorm 2.0 (Alpha))

Function Synopsis

F = meg4read(dsfolder,Trials,TimeIndices,CHANNELS,VERBOSE)

Help Text

MEG4READ - READ CTF MEG4 file format
 function F = meg4read(dsfolder,Trials,TimeIndices,CHANNELS,VERBOSE)
 reads directly the meg4 file in a ds folder once the res4 was read using ds2brainstorm
 useful when extracting data by block when using native file format and when dat awas not translated in
 BraInStorm data format (mainly for disk space saving purposes)
 
 dsfolder : character string containing the name of the original ds folder
 Trials : integer specifying the trial to extract from the data file
 TimeIndices : a vector specifying the indices of the time samples to extract (default is extract entire trial length)
 CHANNELS : indices of channels to extract
 VERBOSE : turn on/off verbose mode (default is off : 0 )

 F : a matrix containing the data.

Cross-Reference Information

bst_message_window                 C:\BrainStorm_2001\Toolbox\bst_message_window.m
ds2brainstorm                      C:\BrainStorm_2001\Toolbox\ds2brainstorm.m
offset                             C:\BrainStorm_2001\Toolbox\offset.m

bst_sourceimaging                  C:\BrainStorm_2001\Toolbox\bst_sourceimaging.m
dataplot_cb                        C:\BrainStorm_2001\Toolbox\dataplot_cb.m
lcmvbf                             C:\BrainStorm_2001\Toolbox\lcmvbf.m
minnorm                            C:\BrainStorm_2001\Toolbox\minnorm.m

Listing of function C:\BrainStorm_2001\Toolbox\meg4read.m

function F = meg4read(dsfolder,Trials,TimeIndices,CHANNELS,VERBOSE)
%MEG4READ - READ CTF MEG4 file format
% function F = meg4read(dsfolder,Trials,TimeIndices,CHANNELS,VERBOSE)
% reads directly the meg4 file in a ds folder once the res4 was read using ds2brainstorm
% useful when extracting data by block when using native file format and when dat awas not translated in
% BraInStorm data format (mainly for disk space saving purposes)
% 
% dsfolder : character string containing the name of the original ds folder
% Trials : integer specifying the trial to extract from the data file
% TimeIndices : a vector specifying the indices of the time samples to extract (default is extract entire trial length)
% CHANNELS : indices of channels to extract
% VERBOSE : turn on/off verbose mode (default is off : 0 )
%
% F : a matrix containing the data.

%<autobegin> ---------------------- 26-May-2004 11:30:56 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Utility - General
%
% Alphabetical list of external functions (non-Matlab):
%   toolbox\bst_message_window.m
%   toolbox\ds2brainstorm.m
%
% At Check-in: $Author: Mosher $  $Revision: 14 $  $Date: 5/26/04 9:59a $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 24-May-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> ------------------------ 26-May-2004 11:30:56 -----------------------


DCOffset = 0; % Remove DCoffset based on entire trial length

% Load Res4 information
[PATHSTR,NAME,EXT,VERSN] = fileparts(dsfolder) ;

if ~isempty(PATHSTR)
    cd(PATHSTR)
end

try 
    [path,rootname] = fileparts(dsfolder);
    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
    
    try 
        load(res4_mat)
    catch
        [f,Channel,imegsens,ieegsens,iothersens,irefsens,grad_order_no,no_trials,filter,Time,RunTitle] = ds2brainstorm(dsfolder,VERBOSE,3,[],[],[],1);
        res4_mat = fullfile('BrainStorm',[rootname,'_res4.mat']);
        load(res4_mat)
    end
    
    nref=length(irefsens);
catch
    errordlg(...
        {sprintf('Please make sure that resource file %s is available',[NAME,'_res4.mat']),...
            sprintf('in %s.', dsfolder)},...
        'res4 file missing'...
        )
    return
end

if nargin == 3
    VERBOSE = 0;
    CHANNELS = [1:length(Channel)];
elseif nargin == 4
    CHANNELS = [1:length(Channel)];
end


meg4file = ([NAME,'.meg4']);

[meg,message] = fopen(meg4file,'rb','s'); % Big-endian byte ordering
if meg < 0
    errordlg(message,mfilename)
    return
end

% no_samples = round((TIME(end)-TIME(1))*gSetUp.sample_rate) +1; % Number of time samples in a trial

% Extraction time range (sample indices)
if nargin == 2 % Not defined: take entire trial length
    TimeIndices = 1:length(Time);
else
    if length(TimeIndices) > length(Time)
        errordlg(...
            sprintf('Only %d time samples are available in each trial. Please adjust time extraction window',length(Time)));
        return
    end
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.

no_samples = length(TimeIndices);
no_trials = 1:gSetUp.no_trials;

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

F = zeros(no_channels,length(Trials)*no_samples);

for k = 1:length(Trials)
    itrial = Trials(k); %0
    %itrial = itrial+1;
    if VERBOSE
        bst_message_window({...
                sprintf('ds2brainstorm -> Trial %d / %d',itrial,length(no_trials))...
            })
    end
    
    Ftmp = zeros(no_channels,no_samples);
    
    if itrial == no_trials(1) % Read first trial
        fseek(meg,(itrial-1)*ByteSizeOfTrial + FirstChannelSkip + (TimeIndices(1)-1)*4 ,-1);
    else % just shift from the size of a trial
        fseek(meg,LastChannelSkip + (itrial-1)*ByteSizeOfTrial + FirstChannelSkip + (TimeIndices(1)-1)*4 ,-1);
    end
    
    Ftmp = fread(meg,[no_samples no_channels],[num2str(no_samples),'*int32=>int32'], samples_skip*4)';
    %F = F(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
    Ftmp = diag(1./gain_chan(CHANNELS))*double(Ftmp);
    
    % 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
            [tmp ia ib] = intersect(CHANNELS,imegsens);
            
            if ~isempty(ia)
                Ftmp(ia,:) = Ftmp(ia,:) - repmat(mean(Ftmp(ia,:)')',1,size(Ftmp,2));
            end
            
        else % Based on pretrigger
            TimeNeg = intersect(find(Time(TimeIndices(1):TimeIndices(end)) >= TIME(1)),find(Time(TimeIndices(1):TimeIndices(end))<0)); % Find pretrigger time points
            if ~isempty(TimeNeg)
                [tmp ia ib] = intersect(CHANNELS,imegsens);
                if ~isempty(ia)
                    Ftmp(ia,:) = Ftmp(ia,:) - repmat(mean(Ftmp(ia,TimeNeg)')',1,size(Ftmp,2));
                end
                
            else
                % Do nothing
            end
        end
        
        if VERBOSE
            bst_message_window({...
                    'ds2brainstorm -> Removing DC offset - DONE'...
                })
        end
        
    otherwise
        % Do nothing
    end
    
    F(:,(itrial-1)*no_samples + TimeIndices) = Ftmp;
    
end

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

fclose(meg);