[Master Index] [Index for Toolbox]

# montreal_data_create

## (Toolbox/montreal_data_create.m in BrainStorm 2.0 (Alpha))

### Function Synopsis

    Data = montreal_data_create(SNR);

### Help Text

    MONTREAL_DATA_CREATE - Read landmarks from the Montreal Phantom and create simulated data
     function Data = montreal_data_create(SNR);

### Cross-Reference Information

This function calls

    blk_diag                           C:\BrainStorm_2001\Toolbox\blk_diag.m
    os_meg                             C:\BrainStorm_2001\Toolbox\os_meg.m
    save_fieldnames                    C:\BrainStorm_2001\Toolbox\save_fieldnames.m

### Listing of function C:\BrainStorm_2001\Toolbox\montreal_data_create.m

    function Data = montreal_data_create(SNR);
    %MONTREAL_DATA_CREATE - Read landmarks from the Montreal Phantom and create simulated data
    % function Data = montreal_data_create(SNR);

    %<autobegin> ---------------------- 14-Jun-2004 17:10:33 -----------------------
    % --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
    %
    % CATEGORY: Help
    %
    % Alphabetical list of external functions (non-Matlab):
    %   toolbox\blk_diag.m
    %   toolbox\os_meg.m
    %   toolbox\save_fieldnames.m
    %
    % Group : Preference data and their calls in this file:
    %   BRAINSTORM = getpref('BrainStorm','brainstormHomeDir');
    %
    % At Check-in: $Author: Mosher $  $Revision: 4 $  $Date: 6/14/04 3:38p $
    %
    % This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 14-Jun-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> ------------------------ 14-Jun-2004 17:10:33 -----------------------

    % ----------------------------- Script History ---------------------------------
    % JCM 27-May-2004  Creation
    % ----------------------------- Script History ---------------------------------

    if ~exist('SNR','var') | isempty(SNR),
       SNR = 100; % signal to noise ratio
    end

    % where is the home dir
    BRAINSTORM = getpref('BrainStorm','brainstormHomeDir');

    % setup the directory locations
    SUBJECT = fullfile(BRAINSTORM,'Phantom','montreal_subject');
    DATA = fullfile(BRAINSTORM,'Phantom','montreal_data');

    % load the Image file

    cd(SUBJECT);

    % load the BST MRI file
    Image = load('montreal_subjectimage');

    DipoleLocs = Image.Landmarks.MRImmXYZ; % the preset dipole locations in MRI mm.

    % convert to subject coordinate system

    DipoleLocs = Image.SCS.R * DipoleLocs + Image.SCS.T(:,[1 1 1]);

    % convert to MKS (meters)

    DipoleLocs = DipoleLocs / 1000;

    % so each column of DipoeLocs is a desired dipole location

    NumSources = size(DipoleLocs,2); % number of sources

    % Load the channel information

    cd(DATA)
    Channel = load('montreal_channel');
    Channel = Channel.Channel; % a 151 channel helmet gradiometer array

    % We will run the single sphere model. A good location for the sphere center in
    % this phantom is [0 0 6] cm, i.e. 60 mm up from the SCS origin.

    NumChannels = length(Channel); % number channels
    [Param(1:NumChannels).Center] = deal([0 0 6/100]');

    % create the gain matrix for dipoles (order -1).

    G = os_meg(DipoleLocs,Channel,Param,-1);

    % Each three columns represents one dipole location.
    % Since this is a single spherical model, set the orientation to a non-radial
    % direction.

    DipoleOrient = zeros(3,NumSources); % one column for each dipole

    for i = 1:NumSources,
       [ignore,ignore,temp] = svd(G(:,[-2:0]+i*3),0);
       DipoleOrient(:,i) = temp(:,1); % strongest direction.
    end

    % now create time series for each source.

    Mother = hann(101) * 10e-9; % basic hump shape, max of 10 nA-m
    PadLength = 30; % number of sample to slide
    % each additional source slides a PadLength
    S = zeros(length(Mother) + (NumSources-1)*PadLength,NumSources); % one column per source

    for i = 1:NumSources,
       ndx = (i-1)*PadLength + [1:length(Mother)]; % next source
       S(ndx,i) = Mother;
    end

    % now create the noiseless forward data

    A = G*blk_diag(DipoleOrient,1); % multiply each unconstrained dipole by its orientation

    B = A*S'; % now multiple each constrained dipole by its time series.

    B = B*1e15; % convert to fT for convenience

    PowerB = B(:)'*B(:); % sum squared power of the noiseless data

    V = randn(size(B)); % white random noise, unity elements

    % now scale the noise to the desired SNR

    PowerV = PowerB/SNR; % desired power

    V = sqrt(PowerV/(V(:)'*V(:))) * V; % rescaled

    F = B + V; % the noisy data

    F = F/1e15; % back to MKS units

    % now save out the data in BST format, see help_data_data

    Data = struct('ChannelFlag',[],'Comment',[],'Device',[],'System',[],...
       'F',[],'NoiseCov',[],'Projector',[],'SourceCov',[],'Time',[]);

    Data.ChannelFlag = ones(1,NumChannels); % all good
    Data.Comment = 'Simulated phantom data';
    Data.Device = 'CTF'; % simulated machine on which acquired
    Data.System = 'CTF'; % coordinate system used
    Data.F = F; % the spatio-temporal data matrix
    Data.Time = [0:(size(S,1)-1)] * 1e-3; % one millisecond sampling

    % now save to a unique filename
    CLOCK = clock; % get the current time
    TimeHash = CLOCK([-1:0]+end); % minutes, seconds
    Filename = sprintf('montreal_data_%02.0f%02.0f.mat',TimeHash); % unique name
    FullName = fullfile(DATA,Filename);

    while exist(FullName,'file'), % name already exists, rare occurrence
       % subtract another minute
       TimeHash(1) = mod(TimeHash(1) - 1,60);
       Filename = sprintf('montreal_data_%02.0f%02.0f.mat',TimeHash); % unique name
       FullName = fullfile(DATA,Filename);
    end

    save_fieldnames(Data,FullName); % saves out the fields of Data into the file
    % can load directly as Data = load(FullName);

Produced by color_mat2html, a customized BrainStorm 2.0 (Alpha) version of mat2html on Tue Nov 2 16:17:59 2004
Cross-Directory links are: ON