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