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

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