[Master Index]
[Index for Toolbox]
multidipole_simulation_eeg_data_create
(Toolbox/multidipole_simulation_eeg_data_create.m in BrainStorm 2.0 (Alpha))
Function Synopsis
[Data,G] = montreal_data_create(SNR);
Help Text
MONTREAL_DATA_CREATE - Create simulated EEG data using EGI array
function [Data,G] = montreal_data_create(SNR);
function Data = montreal_data_create(SNR);
Cross-Reference Information
This function calls
Listing of function C:\BrainStorm_2001\Toolbox\multidipole_simulation_eeg_data_create.m
function [Data,G] = montreal_data_create(SNR);
%MONTREAL_DATA_CREATE - Create simulated EEG data using EGI array
% function [Data,G] = montreal_data_create(SNR);
% function Data = montreal_data_create(SNR);
%<autobegin> ---------------------- 02-Nov-2004 16:14:27 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Help
%
% Alphabetical list of external functions (non-Matlab):
% toolbox\blk_diag.m
% toolbox\bst_headmodeler.m
% toolbox\find_brainstorm_files.m
% toolbox\save_fieldnames.m
%
% Group : Preference data and their calls in this file:
% BRAINSTORM = getpref('BrainStorm','brainstormHomeDir');
%
% At Check-in: $Author: $ $Revision: $ $Date: $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 02-Nov-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> ------------------------ 02-Nov-2004 16:14:27 -----------------------
% ----------------------------- Script History ---------------------------------
% JCM 27-May-2004 Creation
% ----------------------------- Script History ---------------------------------
if ~exist('SNR','var') | isempty(SNR),
SNR = 100; % signal to noise ratio
end
PRESTIM = 100; % number of samples for prestim
% where is the home dir
BRAINSTORM = getpref('BrainStorm','brainstormHomeDir');
% setup the directory locations
SUBJECT = fullfile(BRAINSTORM,'Phantom','montreal_eeg_subject');
DATA = fullfile(BRAINSTORM,'Phantom','montreal_eeg_data');
% load the Image file
cd(SUBJECT);
% load the BST MRI file
Image = load('EGI_256_warped_subjectimage');
DipoleLocs = Image.Landmarks.MRImmXYZ; % the preset dipole locations in MRI mm.
NumSources = size(DipoleLocs,2); % number of sources
% convert to subject coordinate system
DipoleLocs = Image.SCS.R * DipoleLocs + Image.SCS.T(:,ones(1,NumSources));
% convert to MKS (meters)
DipoleLocs = DipoleLocs / 1000;
% so each column of DipoeLocs is a desired dipole location
% Load the channel information
cd(DATA)
Channel = load('EGI_256_channel');
Channel = Channel.Channel; % a 256 channel helmet hair net array, average reference
% 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 4.5/100]');
% setup the head model for EEG
% if it doesn't already exist
[DataDir,DataPopup,Leader] = find_brainstorm_files('headmodel',DATA);
if isempty(DataPopup),
Options = bst_headmodeler; % get the options structure
% from the headmodel Gui fit to this head
%Radii = [8.1 8.6 9.3]';
%Conductivity = [0.33 0.004 2]';
%Center = [0.7 -0.1 3.8]';
%[mu_berg,lam_berg] = berg(Radii, Conductivity);
% create the gain matrix for dipoles (order -1).
Options.ChannelFile = fullfile(DATA,'EGI_256_channel.mat');
Options.Scalp.FileName = fullfile(SUBJECT,'EGI_256_warped_4layer_tess.mat');
Options.Scalp.iGrid = 3;
Options.Cortex.FileName = fullfile(SUBJECT,'EGI_256_warped_4layer_tess.mat');
Options.Cortex.iGrid = 4;
Options.ChannelType = 'EEG';
Options.ImageGridFile = 'default';
Options.HeadModelFile = 'default';
Options.VolumeSourceGridSpacing = 1.0; % cm
cd(DATA)
[G,ReturnOptions] = bst_headmodeler(fullfile(DATA,'EGI_256_brainstormstudy.mat'),Options); % generate the head model
else
G = fullfile(DATA,'EGI_256_headmodelSurfGrid_CD.mat');
end
HeadModel = load(G);
% CHEAT, possible flaw in the saving convention of bst_headmodeler, JCM fix:
HeadModel.GridLoc = strrep(HeadModel.GridLoc,fullfile(BRAINSTORM,'Phantom'),'');
save_fieldnames(HeadModel,G);
G= feval(HeadModel.Function{1},DipoleLocs,Channel,HeadModel.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
switch 'upwards'
case 'strong'
for i = 1:NumSources,
[ignore,ignore,temp] = svd(G(:,[-2:0]+i*3),0);
DipoleOrient(:,i) = temp(:,1); % strongest direction.
end
case 'upwards'
for i = 1:NumSources,
DipoleOrient(:,i) = [sin(pi/6) 0 cos(pi/6)]'; % strongest direction.
end
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(PRESTIM + length(Mother) + (NumSources)*PadLength,NumSources); % one column per source
for i = 1:NumSources,
ndx = PRESTIM + (i)*PadLength + [1:length(Mother)]; % next source
S(ndx,i) = (-1)^i*Mother;
end
% reset these waveforms to be synchronous
S(:,[3 4]) = S([[-9:0]+end [1:(end-10)]],[1 2]);
% 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
%Vprestim = randn(size(B,1),100);
%Vprestim = sqrt(PowerV/(Vprestim(:)'*Vprestim(:)))*Vprestim; % rescaled
F = B + V; % the noisy data
%F = [Vprestim F];
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(F,2)-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('EGI_256_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('EGI_256_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);
FullName = strrep(FullName,'_data_','_true_');
save(FullName,'B', 'S', 'G', 'DipoleOrient', 'DipoleLocs')
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