SYNOPSIS:
supFunSimprovides a comprehensive simulations framework for solving the forward and the inverse problems in electroencephalography (EEG).
supFunSim (provided as supFunSim.org file) contains both ordinary
human language and (discrete) computing machine code that are
available in this document in the form of plain text. This text
tells a story about the reconstruction of neuronal activity
originating from discrete nodes spread among separate brain
regions. This reconstruction is based on the EEG signal.
First neuronal activity is realistically simulated (we also refer to that activity as time series in source-space). The solution to a forward problem (FP) is based on (1) the geometry and conductivity of head compartments and (2) the position of electrodes on the scalp. This solution allows for generation of simulated EEG signals (i.e., time-series in sensor-space that result from the activity of neuronal sources).
Later, sensor-space time-series are used as an input for a variety of spatial filtering methods implemented here. This produces reconstruction (OR estimation) of the original activity of the neuronal sources (i.e., each spatial filter constitutes (OR provides) an inverse problem solution). Next, filtering output is compared with the original source-space time-series in order to estimate reconstruction errors of a given spatial filter.
Some of the filters implemented here are known in the literature and some are novel generalizations of the aforementioned. A number of parameters is available for manipulation in our simulations setup, these are, among others:
- number of considered sources,
- signal to noise ratios for a number of noise types (described below).
We also provide a number of reconstruction error evaluation methods (that are relevant for the present neuroscience research).
Herein we adapt number of solutions provided in toolboxes developed by others, namely:
Also, please note that creation of this document is an ongoing process as we intend to expand it (therefore we are not seeking for it to be complete in any near future).
This document was written using GNU Emacs in plain text exercising advantages of org-mode and (org-)Babel. This Emacs mode allows for the application of a literate approach to programming. Briefly speaking, the source code blocks are interspersed with ordinary human language blocks that provide explanations and some insights into governing dynamics and intrinsic mechanics of the code.
org-mode allows for easy export of whole documents to HTML,
\LaTeX, PDF and other formats. Therefore, various versions of
this document may also be available, but the plain text version is
favorable for studying the simulations setup (as it provides the
most interactive encounter with the code). org-mode also allows
for latex formatted equations, e.g.,
$ei π + 1 = 0$
GNU Emacs (with org-mode and Babel enabled) provides an
environment in which this document can be read and at the same
time the code snippets can be executed (evaluated) within the same
environment and the code execution results can be incorporated in
this document as a plain text, figures or in another form
(provided that, apart from GNU Emacs, the given machine can
execute MATLAB or GNU Octave scripts).
In order to execute all of the simulations steps in a single
run (or multiple runs) please proceed to the section BATCH (or
LOOP respectively). Those sections extensively exert on the
tangling, a feature of the org-mode that allows for source
code blocks content to be written to separate files that can be
later executed (e.g., using MATLAB).
For this simulations files that are tangled from this
org-mode document are stored in the same directory as the
supFunSim.org file and have names that follow the pattern
tg_*.
Simulations are hosted on GitHub.com:
We have tested supFunSim using:
- Lenovo ThinkPad W540:
- Intel(R) Core(TM) i7-4910MQ CPU @ 2.90GHz,
- 32768 MB RAM 1600 MHz DDR3,
- GNU with Linux operating system:
- SUSE Linux Enterprise Desktop 12 SP2,
- x86_64-suse-linux-gnu,
- kernel release 4.4.21-81-default,
- software:
- GNU Emacs 25.1.2 (GTK+ Version 3.20.9),
- MathWorks MATLAB R2017a (9.2.0.538062).
Copyright (c) 2014-2017 Jan Nikadon, Tomasz Piotrowski, Dania Gutiérrez
supFunSim contains both ordinary human language text and discrete
computing machine code (both) provided in this document. Ordinary
human language text is copyrighted by authors. Digital computing
machine code is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
supFunSim is provided “AS IS”, WITHOUT WARRANTY of ANY KIND,
EXPRESS or IMPLIED, INCLUDING but NOT LIMITED to the WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE and
NONINFRINGEMENT. In NO event shall the authors or copyright holders
be liable for ANY CLAIM, DAMAGES or OTHER LIABILITY, WHETHER in an
ACTION of CONTRACT, TORT or OTHERWISE, arising from, out of or in
connection with the software or the use or other dealings in the
software.
See the GNU General Public License for more details. You should
have received a copy of the GNU General Public License along with
supFunSim. If not, see http://www.gnu.org/licenses/
The following code block prepares MATLAB workspace, adds paths to the necessary toolboxes and ./fun directory containing some additional functions. It also loads initial data:
- head compartments geometry (cortex)
- mesh containing candidates for the location of deep sources (based on thalami)
- mesh containing candidates for the location of deep sources (based on icosahedron642)
- cortex geometry with (anatomical) ROI parcellation
- volume conduction model (head-model)
- geometry of electrode positions
Additionally we also provide:
- pre-computed leadfields for the above mentioned source candidates
Contents of the above mentioned files can be checked with:
whos('-file','./mat/sel_atl.mat')
whos('-file','./mat/sel_ele.mat')
whos('-file','./mat/sel_geo_deep_icosahedron642.mat')
whos('-file','./mat/sel_geo_deep_thalami.mat')
whos('-file','./mat/sel_mri00.mat')
whos('-file','./mat/sel_msh.mat')
whos('-file','./mat/sel_src.mat')
whos('-file','./mat/sel_vol.mat')
Load each file separately (mostly for testing and didactic purpose).
load('./mat/sel_atl.mat')
load('./mat/sel_ele.mat')
load('./mat/sel_geo_deep_icosahedron642.mat')
load('./mat/sel_geo_deep_thalami.mat')
load('./mat/sel_mri00.mat')
load('./mat/sel_msh.mat')
load('./mat/sel_src.mat')
load('./mat/sel_vol.mat')
NB: For the sake of backward compatibility and to facilitate bug tracking, it also provides initial settings for the simulations. Pleas remember that ut supra those settings (most likely) need not to be modified. During batch execution of simulations initial settings are overwritten by values defined in sections BATCH and LOOP (the very last sections of this document.)
NB: In the code we make extensive use of a number of acronyms and abbreviations, these include:
- SrcActiv
- Activity of interest (biological)
- IntNoise
- Interference (biological) noise
- BcgNoise
- Background (biological) noise
- MesNoise
- Measurement noise
We assume that this file is located in (or symbolically linked to)
~/supFunSim on machines running GNU with Linux (or any similar
decent operating system). Similarly the necessary toolboxes are
expected to be found in sub-directories of the ~/toolboxes
directory (or symlinked there). Otherwise the following code might
need to be adjusted (e.g., by changing value of the TMP_TBX_PATH
variable etc.) Some additional necessary functions are located in
the ./fun directory.
Simulation parameters are controlled using SETUP variable. The
geometrical arrangement and number of cortical sources in ROIs can
be controlled using SRCS field (SETUP.SRCS), which is a 3
column <int> array.
- Rows of
SETUP.SRCSreppresent consequent ROIs. - Cols of
SETUP.SRCSrepresent signal types, namely:SrcActiv,IntNoiseandBcgNoise,
respectively.
- Integers represent a number of sources in the given ROI for the given signal type.
Variable SETUP also contains field DEEP (SETUP.DEEP) which
defines number of signals in the brain center (around thalami) that
belong to particular signal type (same as for cortical signal types
these include: SrcActiv, IntNoise and BcgNoise).
% Tidy up workspace and change working directory
clc; close all;clearvars('-except','fPath');
if exist('fPath'), cd(fPath); else, try, cd('~/supFunSim/'); catch, warningMessage = 'Problem encoutered while trying to change working directory to ''~/supFunSim/''.'; end; end; disp(['CYBERCRAFT:: pwd is: ',pwd])
% Add path to additional functions
addpath([pwd,'/fun']);
% Add path to toolboxes
if 0, disp('COMMENT: if having problems with toolboxes consider wise use of:'); restoredefaultpath; end
clearvars TMP_TOOLB_PATH; if exist('~/toolboxes/','dir') == 7, TMP_TOOLB_PATH = '~/toolboxes/'; else, warningMessage = sprintf('CYBERCRAFT:: Warning:: toolboxes directory was not found (use ''~/toolboxes/'')'); end; if exist('TMP_TOOLB_PATH','var'), disp(['CYBERCRAFT:: looking for toolboxes in: ',TMP_TOOLB_PATH]), end;
addpath([TMP_TOOLB_PATH,'/arfit']);
addpath([TMP_TOOLB_PATH,'/mvarica']);
addpath([TMP_TOOLB_PATH,'spm12/']);
addpath([TMP_TOOLB_PATH,'spm12/toolbox/aal/']);
addpath([TMP_TOOLB_PATH,'eeglab/']);
addpath([TMP_TOOLB_PATH,'fieldtrip/']); ft_defaults;
if exist([TMP_TOOLB_PATH,'fieldtrip/privatePublic/',filesep]), addpath([TMP_TOOLB_PATH,'fieldtrip/privatePublic/',filesep]); else, copyfile([TMP_TOOLB_PATH,'fieldtrip/private/',filesep], [TMP_TOOLB_PATH,'fieldtrip/privatePublic/',filesep]); addpath([TMP_TOOLB_PATH,'fieldtrip/privatePublic/',filesep]); end;
if exist([TMP_TOOLB_PATH,'fieldtrip/forward/privatePublic/',filesep]), addpath([TMP_TOOLB_PATH,'fieldtrip/forward/privatePublic/',filesep]); else, copyfile([TMP_TOOLB_PATH,'fieldtrip/forward/private/',filesep],[TMP_TOOLB_PATH,'fieldtrip/forward/privatePublic/',filesep]); addpath([TMP_TOOLB_PATH,'fieldtrip/forward/privatePublic/',filesep]); end;
addpath([TMP_TOOLB_PATH,'fieldtrip/utilities']);
if exist([TMP_TOOLB_PATH,'fieldtrip/utilities/privatePublic/',filesep]), addpath([TMP_TOOLB_PATH,'fieldtrip/utilities/privatePublic/',filesep]); else, copyfile([TMP_TOOLB_PATH,'fieldtrip/utilities/private/',filesep],[TMP_TOOLB_PATH,'fieldtrip/utilities/privatePublic/',filesep]); addpath([TMP_TOOLB_PATH,'fieldtrip/utilities/privatePublic/',filesep]); end;
% Initialize FieldTrip
ft_defaults;
% Load head geometry, electrode positions and ROIs data.
load('./mat/sel_msh.mat'); % head compartments geometry (cortex)
load('./mat/sel_geo_deep_thalami.mat'); % mesh containing candidates for lacation of deep sources (based on thalami)
load('./mat/sel_geo_deep_icosahedron642.mat'); % mesh containing candidates for lacation of deep sources (based on icosahedron642)
load('./mat/sel_atl.mat'); % cortex geometry with (anatomical) ROI parcellation
load('./mat/sel_vol.mat'); % volume conduction model (head-model)
load('./mat/sel_ele.mat'); % geometry of electrode positions
load('./mat/sel_src.mat'); % all cx leadfields
clearvars ii jj kk mm nn tmp* SETUP
% Simulations main setup
SETUP.rROI = logical(0); % random (1) or predefined (0) ROIs
SETUP.rPNT = logical(0); % random (1) or predefined (0) candidate points for source locations: if 0,
% number of sources as in SETUP.SRCS(1,1) will be fixed and in close locations
SETUP.SRCS = []; % Cortical sources (avoid placing more than 10 sources in single ROI)
SETUP.SRCS = [ SETUP.SRCS; 8 0 3 ];
SETUP.SRCS = [ SETUP.SRCS; 2 0 0 ];
SETUP.SRCS = [ SETUP.SRCS; 2 1 0 ];
SETUP.SRCS = [ SETUP.SRCS; 0 4 0 ];
SETUP.SRCS = [ SETUP.SRCS; 0 4 0 ];
SETUP.SRCS = [ SETUP.SRCS; 0 4 0 ];
SETUP.SRCS = [ SETUP.SRCS; 0 4 3 ];
SETUP.SRCS = [ SETUP.SRCS; 0 0 3 ];
SETUP.SRCS = [ SETUP.SRCS; 0 0 3 ];
SETUP.SRCS = [ SETUP.SRCS; 0 0 3 ];
SETUP.DEEP = [ 2 1 6 ]; % deep sources
SETUP.ERPs = 5; % Add ERPs (timelocked activity)
SETUP.rROI = 0; % random (1) or predefined (0) ROIs
SETUP.ELEC = size(sel_ele.elecpos,1); % number of electrodes
SETUP.n00 = 500; % number of time samples per trial
SETUP.K00 = 1; % number of independent realizations of signal and noise based on generated MVAR model
% note: covariance matrix R of observed signal and noise covariance matrix N are estimated from samples originating
% from all realizations of signal and noise
SETUP.P00 = 6; % order of the MVAR model used to generate time-courses for signal of interest
SETUP.FRAC = 0.20; % proportion of ones to zeros in off-diagonal elements of the MVAR coefficients masking array
SETUP.STAB = 0.99; % MVAR stability limit for MVAR eigenvalues (less than 1.0 results in more stable model producing more stationary signals)
SETUP.RNG = [0,2.8]; % range for pseudo-random sampling of eigenvalues for MVAR coefficients range
SETUP.ITER = 5e5; % iterations limit for MVAR pseudo-random sampling and stability verification
SETUP.PDC_RES = [0:0.01:0.5]; % resolution vector for normalized PDC estimation
SETUP.TELL = 1; % provide additional comments during code execution ("tell me more")
SETUP.PLOT = 1; % plot figures during the intermediate stages
SETUP.SCRN = get(0,'MonitorPositions'); % get screens positions
SETUP.DISP = SETUP.SCRN(end,:); % force figures to be displayed on (3dr) screenscreen
SETUP.FIXED_SEED = 0; % Settings for seed selection
if SETUP.FIXED_SEED, SETUP.SEED = rng(1964);else,SETUP.SEED = rng(round(1e3*randn()^2*sum(clock)));end
SETUP.RANK_EIG = sum(SETUP.SRCS(:,1)); % rank of EIG-LCMV filter: set to number of active sources
SETUP.fltREMOVE = 1; % to keep (0) or remove (1) selected filters
SETUP.SHOWori = 1; % to show (1) or do not show (0) Original and Dummy signals on Figures
SETUP.IntLfgRANK = round(0.3*sum(SETUP.SRCS(:,2))); % rank of patch-constrained reduced-rank leadfield
SETUP.supSwitch = 'loc'; % 'rec': run reconstruction of sources activity, 'loc': find active sources
NB: The following code is not executed during simulations, it is left here only to facilitate bug tracking and for explanatory/illustrative purpose.
The above applies to all Checkup sections (usually marked with
org-todo-keyword: !opt).
whos
SETUP
chkSim___tg_s01_PRE___001(SETUP);
sel_atl
sel_ele
sel_ele.label
sel_msh
{sel_msh.bnd.inf}'
sel_vol
SETUP
function chkSim___tg_s01_PRE___000_PARSE_SETUP(SETUP,sel_atl)
% Simulations setup parser
disp('Checking consistency of ''SETUP'' definitions')
% check if number of ROIs is smaller than number of cortex parcels
tmp_MaxROIs = size(sel_atl.Atlas(sel_atl.atl).Scouts,2);
tmp_ActROIs = size(SETUP.SRCS,1);
if 0
disp(SETUP.SRCS)
tmp_MaxROIs
tmp_ActROIs
end
if tmp_ActROIs > tmp_MaxROIs,
tmp_msg = ['please avoid putting more than ', num2str(tmp_MaxROIs) ' source ROIs (decrease the number of rows in the ''SETUP.SRCS''.)'];
error(tmp_msg);
end
% check if no roi contains more sources than allowed (allowed number
% of sources is the number of nodes in the ROI that contains the
% least nodes)
tmp_MaxSrc = []; for ii = 1:tmp_MaxROIs, tmp_MaxSrc = [tmp_MaxSrc; size(sel_atl.Atlas(sel_atl.atl).Scouts(ii).Vertices,2)]; end;
tmp_MaxSrc = min(tmp_MaxSrc);
[tmp_ActSrcs, tmp_ActColMaxIdx] = max(sum(SETUP.SRCS,2));
if 0
SETUP.SRCS
tmp_MaxSrc
tmp_ActSrcs
tmp_ActColMaxIdx
end
if tmp_ActSrcs > min(tmp_MaxSrc),
tmp_msg = ['please avoid putting more than ', num2str(tmp_MaxSrc) ' sources in any of the ROIs (decrease sum of the column number ', num2str(tmp_ActColMaxIdx), ' in the ''SETUP.SRCS''.)'];
error(tmp_msg);
end
% Check if a number of ERPs is not greater than the number of signals
% of interest.
tmp_MaxERPs = sum(SETUP.SRCS(:,1));
tmp_ActERPs = SETUP.ERPs;
if 0
SETUP.SRCS
tmp_MaxERPs
tmp_ActERPs
end
if tmp_ActERPs > tmp_MaxERPs,
tmp_msg = ['please avoid putting more than ', num2str(tmp_MaxERPs) ' ERPs (this is a number of sources of activity of interest, please change the sum in the 1st column of ''SETUP.SRCS'' or the ''SETUP.ERPs'' value.)'];
error(tmp_msg);
end
end
function chkSim___tg_s01_PRE___001(SETUP)
fprintf(' \n');
fprintf('CYBERCRAFT:: Number of signal types per CORTEX ROI:\n\n');
disp(array2table([[1:size(SETUP.SRCS,1)]',SETUP.SRCS,sum(SETUP.SRCS,2)],'VariableNames',{'ROI','SrcActiv','IntNoise','BcgNoise','TOTAL'}))
fprintf(' \n');
fprintf('CYBERCRAFT:: TOTAL number of signals for CORTEX:\n\n');
disp(array2table([size(SETUP.SRCS,1);sum(SETUP.SRCS(:,1));sum(SETUP.SRCS(:,2));sum(SETUP.SRCS(:,3));sum(sum(SETUP.SRCS(:,1:3)))]','VariableNames',{'ROI','SrcActiv','IntNoise','BcgNoise','TOTAL'}));
fprintf(' \n');
fprintf('CYBERCRAFT:: Number of signals for DEEP sources:\n\n');
disp(array2table([size(SETUP.SRCS,1)+1;sum(SETUP.DEEP(1,1));sum(SETUP.DEEP(1,2));sum(SETUP.DEEP(1,3));sum(SETUP.DEEP(1,1:3))]','VariableNames',{'ROI','SrcActiv','IntNoise','BcgNoise','TOTAL'}));
fprintf(' \n');
fprintf('CYBERCRAFT:: TOTAL number of signals for CORTEX and DEEP sources:\n\n');
disp(array2table([NaN;sum(SETUP.SRCS(:,1))+sum(SETUP.DEEP(1,1));sum(SETUP.SRCS(:,2))+sum(SETUP.DEEP(1,2));sum(SETUP.SRCS(:,3))+sum(SETUP.DEEP(1,3));sum(sum(SETUP.SRCS(:,1:3)))+sum(SETUP.DEEP(1,1:3))]','VariableNames',{'ROI','SrcActiv','IntNoise','BcgNoise','TOTAL'}));
fprintf(' \n');
disp(['CYBERCRAFT:: MesNoise: ',num2str(SETUP.ELEC)])
end
This section provides further simulation settings.
NB: For the sake of backward compatibility and to facilitate bug tracking, it also provides initial settings for the simulations. Pleas remember that ut supra those settings (most likely) need not to be modified because for execution of simulations in batch they will be overwritten by values defined in sections BATCH and LOOP, the two sections that can be found very close to the end of this document.
SETUP.CUBE = 20; % perturbation of the leadfields based on the shift of source position within a cube of given edge length (centered at the original leadfields positions), in [mm]
SETUP.CONE = pi/32; % perturbation of the leadfields based on the rotation of source orientation (azimuth TH, elevation PHI), in [rad]
SETUP.H_Src_pert = 0; % use original (0) or perturbed (1) leadfield for signal reconstruction
SETUP.H_Int_pert = 0; % use original (0) or perturbed (1) leadfield for nulling constrains
SETUP.SINR = 5; % signal to interference noise power ratio expressed in dB (both measured on electrode level)
SETUP.SBNR = 10; % signal to biological noise power ratio expressed in dB (both measured on electrode level)
SETUP.SMNR = 15; % signal to measurment noise power ratio expressed in dB (both measured on electrode level)
SETUP.WhtNoiseAddFlg = 1; % white noise admixture in biological noise interference noise (FLAG)
SETUP.WhtNoiseAddSNR = 3; % SNR of BcgNoise and WhiNo (dB)
SETUP.SigPre = 0; SETUP.IntPre = 0; SETUP.BcgPre = 1; SETUP.MesPre = 1; % final signal components for pre-interval (use zero or one for signal, interference noise, biological noise, measurement noise)
SETUP.SigPst = 1; SETUP.IntPst = 1; SETUP.BcgPst = 1; SETUP.MesPst = 1; % final signal components for post-interval (as above)
% For localization, the default is to consider random locations of ROI with some fixed candidate points such that
% SETUP.SRCS(1,1) of them are in close locations; additionally, no interfering sources are considered
% in the localization model. You may comment out the 'if' section below to experiment with other settings.
if(SETUP.supSwitch == 'loc')
SETUP.rROI = logical(1); % random (1) or predefined (0) ROIs
SETUP.rPNT = logical(0); % random (1) or predefined (0) candidate points for source locations: if 0,
% number of sources as in SETUP.SRCS(1,1) will be fixed and in close locations
SETUP.SigPre = 0; SETUP.IntPre = 0; SETUP.BcgPre = 1; SETUP.MesPre = 1; % final signal components for pre-interval (use zero or one for signal, interference noise, biological noise, measurement noise)
SETUP.SigPst = 1; SETUP.IntPst = 0; SETUP.BcgPst = 1; SETUP.MesPst = 1; % final signal components for post-interval (as above)
end
if SETUP.rPNT
disp('CC: using random source locations')
else
disp('CC: using predefined source locations')
end
whos
SETUP
chkSim___tg_s01_PRE___001(SETUP)
chkSim___tg_s02_PRE___001(SETUP)
function chkSim___tg_s02_PRE___001(SETUP)
fprintf('\n');
fprintf('CYBERCRAFT:: Perturbation cube:\n\n');
disp(array2table(SETUP.CUBE,'VariableNames',{'mm'},'RowNames',{'shift distance'}))
fprintf('\n');
fprintf('CYBERCRAFT:: Perturbation cone:\n\n');
disp(array2table([SETUP.CONE,rawRad2Deg(SETUP.CONE)],'VariableNames',{'rad','deg'},'RowNames',{'rotation angle'}));
fprintf('\n');
fprintf('CYBERCRAFT:: Ratios for signal to:\n\n');
disp(array2table([SETUP.SINR;SETUP.SBNR;SETUP.SMNR],'VariableNames',{'SxNR'},'RowNames',{'IntNoise','BcgNoise','MesNoise'}));
fprintf('\n');
fprintf('CYBERCRAFT:: Signal components:\n\n');
disp(array2table([SETUP.SigPre,SETUP.IntPre,SETUP.BcgPre,SETUP.MesPre;SETUP.SigPst,SETUP.IntPst,SETUP.BcgPst,SETUP.MesPst]','VariableNames',{'pre','post'},'RowNames',{'SrcActiv','IntNoise','BcgNoise','MesNoise'}));
fprintf('\n');
end
Code in this section produces timeseries for the signals described in the table below.
| Variable name | Description |
|---|---|
| sim_sig_SrcActiv | Activity of interest (biological) |
| sim_sig_IntNoise | Interference (biological) noise |
| sim_sig_BcgNoise | Background (biological) noise |
| sim_sig_MesNoise | Measurement noise |
| sim_sig_AdjSNRs | All the above signals combined and adjusted for SNR |
In this section multivariate timeseries are generated for:
sim_sig_SrcActiv, sim_sig_IntNoise, sim_sig_BcgNoise.
Timeseries are based on Multivariate Autoregressive (MVAR) model.
See:
- ./toolboxes/arfit/arfit.pdf
- ./toolboxes/arfit/arfit_alg.pdf
- ./toolboxes/mvarica/arfit/arfit.pdf
- ./toolboxes/mvarica/arfit/arfit_alg.pdf
- 2001 - Neumaier,Schneider - Algorithm 808: ARFIT —A Matlab Package for the Estimation of Parameters and Eigenmodes of Multivariate Autoregressive Models.pdf
- https://sccn.ucsd.edu/wiki/Chapter_3._Multivariate_Autoregressive_Modeling
${{v}ν}=w+\underset{l=1}{\overset{p}{\mathop ∑ }}\,{{A}l}{{v}ν-l}+{{ε}ν}$, where
- ${{v}ν}$ are the
$m$ -dimensional state vectors of (stationary)
time series,
-
$p$ is the order of the model, - matrices ${{A}1}\ldots{{A}p} ∈ Rm×{m}$ are the
coefficient matrices of the AR model,
- ${{ε}ν}$ are the
$m$ -dimensional uncorrelated
random vectors with mean zero and covariance matrix $C ∈ Rm×{m}$
- 2001 - Baccala, Sameshima, Partial directed coherence: a new concept in neural structure determination.pdf
- 2014 - SAMESHIMA,BACCALA - Methods in Brain Connectivity Inference through Multivariate Time Series Analysis.pdf
See
- Wiki: EN: Normalized frequency (unit)
-
Wiki: EN: Normalized frequency (unit): Alternative normalizations
- Some programs (such as MATLAB) that design filters with
real-valued coefficients use the Nyquist frequency
(
$\textstyle f_s/2$ ) as the normalization constant. The resultant normalized frequency has units of half-cycles/sample or equivalently cycles per 2 samples. - Sometimes, the unnormalized frequency is represented in units
of radians/second (angular frequency),
and denoted by
$\textstyle ω$ . When$\textstyle ω$ is normalized by the sample-rate (samples/sec), the resulting units are radians/sample. The normalized Nyquist frequency is π radians/sample, and the normalized sample-rate is 2π radians/sample. - The following table shows examples of normalized frequencies
for a 1 kHz signal, a sample rate
$\textstyle f_\mathrm{s}$ = 44.1 kHz, and 3 different choices of normalized units. Also shown is the frequency region containing one cycle of the discrete-time Fourier transform, which is always a periodic function.Units Domain Remarks Computation Value cycles/sample [-½, ½] or [0,1] 1000 / 44100 0.02268 half-cycles/sample [-1,1] or [0,2] MATLAB 1000 / 22050 0.04535 radians/sample [-π,π] or [0,2π] 2 ”π” 1000 / 44100 0.1425
- Some programs (such as MATLAB) that design filters with
real-valued coefficients use the Nyquist frequency
(
- http://www.mathworks.com/help/signal/ug/frequency-response.html#zmw57dd0e1407
- http://www.mathworks.com/help/signal/ref/freqz.html
- https://www.mathworks.com/matlabcentral/newsreader/view_thread/286192
- https://www.mathworks.com/matlabcentral/answers/50575-normalized-frequency-in-analog-filter-design
- https://www.codementor.io/tips/8102473731/understanding-normalized-frequency-in-matlab
- http://dsp.stackexchange.com/questions/17490/how-to-normalize-frequency-in-matlab
- http://dsp.stackexchange.com/questions/15219/convert-normalized-frequency-to-real-frequency-in-ar-model
clearvars sim_sig_SrcActiv tmp*
sim_sig_SrcActiv = cccSim___makeSimSig(SETUP,1,'MVAR based signal for sources activity');
clearvars ii jj kk nn tmp*
if SETUP.ERPs > 0
sim_sig_SrcActiv.sigSRC_pst_orig = sim_sig_SrcActiv.sigSRC_pst;
sim_sig_SrcActiv.sigSRC_pre_orig = sim_sig_SrcActiv.sigSRC_pre;
tmp_LB = -5;
tmp_UB = 7;
tmp_NN = SETUP.n00;
tmp_PP = 1;
for ee = 1:1:SETUP.ERPs
tmp_PP = mod(ee,3)+1;
sim_sig_SrcActiv.ERPs(ee,:) = gauswavf(tmp_LB,tmp_UB,tmp_NN,tmp_PP);
end
sim_sig_SrcActiv.ERPs_pst = transpose(sim_sig_SrcActiv.ERPs);
sim_sig_SrcActiv.ERPs_pst(size(sim_sig_SrcActiv.sigSRC_pst,1),size(sim_sig_SrcActiv.sigSRC_pst,2)) = 0;
sim_sig_SrcActiv.ERPs_pst = cat(3,repmat(sim_sig_SrcActiv.ERPs_pst,[1,1,SETUP.K00]));
sim_sig_SrcActiv.sigSRC_pst = sim_sig_SrcActiv.sigSRC_pst + 20*sim_sig_SrcActiv.ERPs_pst;
end
NB.: Fields:
w01,A01,C01,SBC01,FPE01andth01
(in general *01) are estimated from the sigMVAR_pst signal.
function sim_sig_SrcActiv = cccSim___makeSimSig(SETUP,set_COL,varargin)
sim_sig_SrcActiv.inf = '';
sim_sig_SrcActiv.S00 = [];
sim_sig_SrcActiv.P00 = [];
sim_sig_SrcActiv.M00 = [];
sim_sig_SrcActiv.A00 = zeros([0,1]);
sim_sig_SrcActiv.mdlStabH = [];
sim_sig_SrcActiv.mdlStabM = [];
sim_sig_SrcActiv.w00 = [];
sim_sig_SrcActiv.C00 = [];
sim_sig_SrcActiv.n00 = [];
sim_sig_SrcActiv.K00 = [];
sim_sig_SrcActiv.sigSRC_pre = permute([],[3,2,1]);
sim_sig_SrcActiv.sigSRC_pst = permute([],[3,2,1]);
sim_sig_SrcActiv.w01 = [];
sim_sig_SrcActiv.A01 = zeros([0,1]);
sim_sig_SrcActiv.C01 = [];
sim_sig_SrcActiv.SBC01 = [];
sim_sig_SrcActiv.FPE01 = [];
sim_sig_SrcActiv.th01 = [];
if ~isempty(varargin)
sim_sig_SrcActiv.inf = varargin{end};
end
sim_sig_SrcActiv.S00 = sum(SETUP.SRCS(:,set_COL))+SETUP.DEEP(1,set_COL); % number of signals
if sim_sig_SrcActiv.S00 > 0
sim_sig_SrcActiv.P00 = SETUP.P00; % order of the MVAR model used to generate time-courses
sim_sig_SrcActiv.M00 = cccSim___diagonMask(sim_sig_SrcActiv.S00,SETUP.FRAC); % MVAR coefficients mask
sim_sig_SrcActiv.A00 = cccSim___stableMVAR(sim_sig_SrcActiv.S00,sim_sig_SrcActiv.P00,sim_sig_SrcActiv.M00,SETUP.RNG,SETUP.STAB,SETUP.ITER);
[sim_sig_SrcActiv.mdlStabH,sim_sig_SrcActiv.mdlStabM] = rawIsStableMVAR(sim_sig_SrcActiv.A00,1);
sim_sig_SrcActiv.w00 = zeros(sim_sig_SrcActiv.S00,1); % expected value for time-courses
sim_sig_SrcActiv.C00 = eye(sim_sig_SrcActiv.S00); % covariance of MVAR model's "driving noise"
sim_sig_SrcActiv.n00 = SETUP.n00; % number of time samples
sim_sig_SrcActiv.K00 = SETUP.K00; % number of independent realizations of the driving AR models
sim_sig_SrcActiv.sigSRC_pre = 10*arsim(sim_sig_SrcActiv.w00,sim_sig_SrcActiv.A00,sim_sig_SrcActiv.C00,[sim_sig_SrcActiv.n00,2*sim_sig_SrcActiv.K00],1e3); % 10e-6; % Current density (10nA*mm)
sim_sig_SrcActiv.sigSRC_pst = sim_sig_SrcActiv.sigSRC_pre(:,:,[end/2+1:end]);
sim_sig_SrcActiv.sigSRC_pre = sim_sig_SrcActiv.sigSRC_pre(:,:,[1:end/2]);
if SETUP.TELL,disp('Fitting MVAR model to generated signal...');end;
[sim_sig_SrcActiv.w01,sim_sig_SrcActiv.A01,sim_sig_SrcActiv.C01,sim_sig_SrcActiv.SBC01,sim_sig_SrcActiv.FPE01,sim_sig_SrcActiv.th01] = arfit(sim_sig_SrcActiv.sigSRC_pst,sim_sig_SrcActiv.P00-round(sim_sig_SrcActiv.P00/2),sim_sig_SrcActiv.P00+round(sim_sig_SrcActiv.P00/2));
end
end
function A00 = cccSim___stableMVAR(S00,P00,M00,set_RNG,set_STAB,set_ITER)
% procedure inspired by function stablemvar that is a part of MVARICA Toolbox.
if 0
A00 = cccSim___stableMVAR(S00,P00,M00,set_RNG,set_ITER,set_STAB)
end
tmp_iterNow = 0;
tmp_lambda = Inf;
while any(abs(tmp_lambda)>set_STAB) && tmp_iterNow < set_ITER
tmp_V = orth(rand(S00*P00,S00*P00));
tmp_U = orth(rand(S00*P00,S00*P00));
lambdatmp = set_RNG(1)+(set_RNG(2)-set_RNG(1))*rand(S00*P00,1);
tmp_A00 = tmp_V*diag(lambdatmp)*tmp_U';
A00 = tmp_A00(1:S00,:);
A00 = A00.*repmat(M00,1,P00); % nulling of some coefficients based on mask
tmp_lambda = eig([A00; eye((P00-1)*S00) zeros((P00-1)*S00,S00)]);
tmp_iterNow = tmp_iterNow + 1;
end
if tmp_iterNow >=set_ITER
A0=[];
error('Could not generate stable MVAR model in given number of iterations (set_ITER)');
end
if 0
S00 = sim_sig00.S00
P00 = sim_sig00.P00
w00 = sim_sig00.w00
C00 = sim_sig00.C00
n00 = sim_sig00.n00
K00 = sim_sig00.K00
end
end
function M00 = cccSim___diagonMask(S00,M00_frc)
% cccSim___diagonMask produces masking square array M (SxS) whose all diagonal and some off-diagonal elements are ones. All
% othere elements are zeros. The off-diagonal elements are sampled pseudo-randomly and their number is declared by
% second argument of this function which should be a number between zero and one that describes the proportion of ones
% to zeros among all off-diagonal elements of the resulting array.
%
% Usage:
%
% MSK = cccSim___diagonMask(S,frac)
% Copyright (C) 2000-2015, Jan (CyberCraft) Nikadon (nikadon-AT-SIGN-gmail.com)
%
% This file is free software and a part of CyberCraft Toolkit. You can redistribute it and/or modify it under the
% terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the
% License, or (at your option) any later version.
%
% CyberCraft is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CyberCraft. If not, see <http://www.gnu.org/licenses/>.
if 0
% Example useage:
S00 = 10;
M00_frc = 0.1;
M00 = cccSim___diagonMask(S00,M00_frc)
figure(1);clf;
imagesc(M00);colorbar;
(numel(find(M00))-S00) / (numel(M00)-S00)
imagesc(sim_sig00.M00);colorbar;
end
M00 = eye(S00);
M00_idx = find(~M00);
M00_smp = datasample(M00_idx,round(M00_frc*numel(M00_idx)),'Replace',false);
M00(M00_smp) = deal(1);
end
size(sim_sig_SrcActiv.sigSRC_pst_orig)
size(sim_sig_SrcActiv.sigSRC_pst)
size(sim_sig_SrcActiv.ERPs_pst)
size(sim_sig_SrcActiv.ERPs)
figure, imagesc(sim_sig_SrcActiv.sigSRC_pst_orig(:,:,1))
figure, imagesc(sim_sig_SrcActiv.sigSRC_pst(:,:,1))
figure, imagesc(sim_sig_SrcActiv.ERPs_pst(:,:,1))
ee = 0
ee = ee + 1
tmp_LB = -5;
tmp_UB = 7;
tmp_NN = SETUP.n00;
tmp_PP = mod(ee,3)+1;
tmp_PP
[PSI,X] = gauswavf(tmp_LB,tmp_UB,tmp_NN,tmp_PP);
close all
plot(PSI)
size(cat(2,20*repmat(transpose(sim_sig_SrcActiv.ERPs),[1,1,SETUP.K00]),zeros(SETUP.n00,sum(SETUP.SRCS(:,1))-SETUP.ERPs,SETUP.K00)))
close all
figure
plot(sim_sig_SrcActiv.ERPs(3,:))
figure
plot(sim_sig_SrcActiv.sigSRC_pst(:,3,1))
figure
plot(mean(sim_sig_SrcActiv.sigSRC_pst(:,3,:),3))
whos
chkSim___tg_s01_PRE___001(SETUP)
sim_sig_SrcActiv
tmp_arrC = [rawSize(sim_sig_SrcActiv.sigSRC_pre,[1,2,3]); rawSize(sim_sig_SrcActiv.sigSRC_pst,[1,2,3])];
tmp_varN = {'Timesamples','Signals','Realizations'};
tmp_rowN = {'sim_sig_SrcActiv.sigSRC_pre','sim_sig_SrcActiv.sigSRC_pst'};
disp(' ');disp(array2table(tmp_arrC,'VariableNames',tmp_varN,'RowNames',tmp_rowN))
tmp_arrC = [size(sim_sig_SrcActiv.A00),size(sim_sig_SrcActiv.A00,2)/size(sim_sig_SrcActiv.A00,1);size(sim_sig_SrcActiv.A01),size(sim_sig_SrcActiv.A01,2)/size(sim_sig_SrcActiv.A01,1)];
tmp_varN = {'Signals','Signals_cdot_ModOrd','ModOrd'};
tmp_rowN = {'sim_sig_SrcActiv.A00','sim_sig_SrcActiv.A01'};
disp(' ');disp(array2table(tmp_arrC,'VariableNames',tmp_varN,'RowNames',tmp_rowN))
~isequal(sim_sig_SrcActiv.sigSRC_pre,sim_sig_SrcActiv.sigSRC_pst)
~isequal(sim_sig_SrcActiv.A00,sim_sig_SrcActiv.A01)
Matrix produced by cccSim___diagonMask (mask for nulling of
some MVAR coefficients).
figure(101);clf;set(gcf,'Position',SETUP.DISP);imagesc(sim_sig_SrcActiv.M00);colormap(rawHotColdColorMap(10));colorbar;
set(gcf,'color','w');
Plotted are:
A00the original coefficient matrix used for timeseries generation andA01coefficient matrix obtained by fitting to the generated timeseries.
figure(105);clf;set(gcf,'Position',SETUP.DISP);subplot(2,1,1);cccSim___rawDispA00(sim_sig_SrcActiv.A00);title('A00');subplot(2,1,2);cccSim___rawDispA00(sim_sig_SrcActiv.A01);title('A01');
set(gcf,'color','w');
figure(120);clf;set(gcf,'Position',SETUP.DISP);cccSim___rawPlotA00(sim_sig_SrcActiv.A00);
set(gcf,'color','w');
figure(121);clf;set(gcf,'Position',SETUP.DISP);cccSim___rawPlotA00(sim_sig_SrcActiv.A01);
set(gcf,'color','w');
function cccSim___rawDispA00(A00)
rawImgSC(A00,8);
tmp_base = 255;
tmp_colorMapMat = rawHotColdColorMap(tmp_base);
caxis([-max(abs(min(min(A00))),abs(max(max(A00)))) max(abs(min(min(A00))),abs(max(max(A00))))]);
colormap(tmp_colorMapMat);
colorbar;
hold on;
tmp_stem = 0.5+size(A00,1):size(A00,1):size(A00,2)-0.5;
tmp_vals = size(A00,1)+0.5*ones(size(tmp_stem));
stem(tmp_stem,tmp_vals,'Color','k','LineWidth',0.5,'Marker', 'none');
set(gca,'XTick',[0:size(A00,1):size(A00,2)])
end
function cccSim___rawPlotA00(A00)
rawPlotPDC(abs(PDC(A00,SETUP.PDC_RES)),SETUP.PDC_RES);
end
In general this signal is intended to be highly correlated (or anti-correlated) with the signal of interest. Also some white noise can be added.
clearvars sim_sig_IntNoise tmp*
sim_sig_IntNoise = sim_sig_SrcActiv;
sim_sig_IntNoise.inf = 'MVAR based signal for interference (biological) noise activity sources';
tmp_IntNoiseCol = 2;
if SETUP.TELL,disp('Getting base for IntNoise...');end;
sim_sig_IntNoise.sigSRC_pre = -sim_sig_IntNoise.sigSRC_pre;
sim_sig_IntNoise.sigSRC_pst = -sim_sig_IntNoise.sigSRC_pst;
[sim_sig_IntNoise.w01,sim_sig_IntNoise.A01,sim_sig_IntNoise.C01,sim_sig_IntNoise.SBC01,sim_sig_IntNoise.FPE01,sim_sig_IntNoise.th01] = deal([]);
if SETUP.TELL,disp('Populating signals for IntNoise...');end;
sim_sig_IntNoise.sigSRC_pre = repmat(sim_sig_IntNoise.sigSRC_pre,1,ceil( (sum(SETUP.SRCS(:,tmp_IntNoiseCol))+SETUP.DEEP(1,tmp_IntNoiseCol))/size(sim_sig_IntNoise.sigSRC_pre,2)),1);
sim_sig_IntNoise.sigSRC_pst = repmat(sim_sig_IntNoise.sigSRC_pst,1,ceil( (sum(SETUP.SRCS(:,tmp_IntNoiseCol))+SETUP.DEEP(1,tmp_IntNoiseCol))/size(sim_sig_IntNoise.sigSRC_pst,2)),1);
if SETUP.TELL,disp('Reducing dimensionality of IntNoise...');end;
sim_sig_IntNoise.sigSRC_pre = sim_sig_IntNoise.sigSRC_pre(:,1:(sum(SETUP.SRCS(:,tmp_IntNoiseCol))+SETUP.DEEP(1,tmp_IntNoiseCol)),:);
sim_sig_IntNoise.sigSRC_pst = sim_sig_IntNoise.sigSRC_pst(:,1:(sum(SETUP.SRCS(:,tmp_IntNoiseCol))+SETUP.DEEP(1,tmp_IntNoiseCol)),:);
if SETUP.TELL,disp('Backing-up IntNoise before white noise admixture...');end;
sim_sig_IntNoise.sigSRC_pre_b4admix = sim_sig_IntNoise.sigSRC_pre;
sim_sig_IntNoise.sigSRC_pst_b4admix = sim_sig_IntNoise.sigSRC_pst;
if SETUP.WhtNoiseAddFlg
if SETUP.TELL,disp('Adding some white noise to the biological noise activity...');end;
sim_sig_IntNoise.noiseAdmix_pre = randn(size(sim_sig_IntNoise.sigSRC_pre));
sim_sig_IntNoise.noiseAdmix_pst = randn(size(sim_sig_IntNoise.sigSRC_pst));
for kk = 1:SETUP.K00
sim_sig_IntNoise.noiseAdmix_pre(:,:,kk) = rawAdjTotSNRdB(sim_sig_IntNoise.sigSRC_pre(:,:,kk),sim_sig_IntNoise.noiseAdmix_pre(:,:,kk),SETUP.WhtNoiseAddSNR);
sim_sig_IntNoise.noiseAdmix_pst(:,:,kk) = rawAdjTotSNRdB(sim_sig_IntNoise.sigSRC_pst(:,:,kk),sim_sig_IntNoise.noiseAdmix_pst(:,:,kk),SETUP.WhtNoiseAddSNR);
end
sim_sig_IntNoise.sigSRC_pre = sim_sig_IntNoise.sigSRC_pre+sim_sig_IntNoise.noiseAdmix_pre;
sim_sig_IntNoise.sigSRC_pst = sim_sig_IntNoise.sigSRC_pst+sim_sig_IntNoise.noiseAdmix_pst;
end
if SETUP.TELL,disp('Fitting MVAR model to generated signal...');end;
[sim_sig_IntNoise.w01,sim_sig_IntNoise.A01,sim_sig_IntNoise.C01,sim_sig_IntNoise.SBC01,sim_sig_IntNoise.FPE01,sim_sig_IntNoise.th01] = arfit(sim_sig_IntNoise.sigSRC_pst,sim_sig_IntNoise.P00-round(sim_sig_IntNoise.P00/2),sim_sig_IntNoise.P00+round(sim_sig_IntNoise.P00/2));
if 0
sim_sig_IntNoise = rmfield(sim_sig_IntNoise,{'sigSRC_pre_b4admix','sigSRC_pst_b4admix','noiseAdmix_pre','noiseAdmix_pst'})
end
clearvars ii jj kk nn tmp*
whos
sim_sig_SrcActiv
sim_sig_IntNoise
chkSim___tg_s01_PRE___001(SETUP)
kk = 1;
jj = 1:min(size(sim_sig_SrcActiv.sigSRC_pre,2),size(sim_sig_IntNoise.sigSRC_pre,2));
tmp_arrC = [rawSNRdB(sim_sig_IntNoise.sigSRC_pre_b4admix(:,:,kk),sim_sig_IntNoise.noiseAdmix_pre(:,:,kk));rawSNRdB(sim_sig_IntNoise.sigSRC_pst_b4admix(:,:,kk),sim_sig_IntNoise.noiseAdmix_pst(:,:,kk))];
tmp_varN = {'SNR'};
tmp_rowN = {'sim_sig_IntNoise.sigSRC_pre_b4admix','sim_sig_IntNoise.sigSRC_pst_b4admix'};
disp(' ');disp(array2table(tmp_arrC,'VariableNames',tmp_varN,'RowNames',tmp_rowN))
tmp_arrC = corrcoef(sim_sig_SrcActiv.sigSRC_pre(:,jj,kk),sim_sig_IntNoise.sigSRC_pre(:,jj,kk));
tmp_varN = {'sim_sig_SrcActiv_a','sim_sig_IntNoise_a'};
tmp_rowN = tmp_varN;
disp(' ');disp(array2table(tmp_arrC,'VariableNames',tmp_varN,'RowNames',tmp_rowN))
tmp_arrC = corrcoef(sim_sig_SrcActiv.sigSRC_pst(:,jj,kk),sim_sig_IntNoise.sigSRC_pst(:,jj,kk));
tmp_varN = {'sim_sig_SrcActiv_b','sim_sig_IntNoise_b'};
tmp_rowN = tmp_varN;
disp(' ');disp(array2table(tmp_arrC,'VariableNames',tmp_varN,'RowNames',tmp_rowN))
tmp_arrC = corrcoef(sim_sig_SrcActiv.sigSRC_pre(:,jj,kk),sim_sig_IntNoise.sigSRC_pst(:,jj,kk));
tmp_varN = {'sim_sig_SrcActiv_a','sim_sig_IntNoise_b'};
tmp_rowN = tmp_varN;
disp(' ');disp(array2table(tmp_arrC,'VariableNames',tmp_varN,'RowNames',tmp_rowN));disp('CYBERCRAFT:: This should be negligible')
tmp_arrC = corrcoef(sim_sig_SrcActiv.sigSRC_pst(:,jj,kk),sim_sig_IntNoise.sigSRC_pre(:,jj,kk));
tmp_varN = {'sim_sig_SrcActiv_b','sim_sig_IntNoise_a'};
tmp_rowN = tmp_varN;
disp(' ');disp(array2table(tmp_arrC,'VariableNames',tmp_varN,'RowNames',tmp_rowN));disp('CYBERCRAFT:: This should be negligible')
clearvars ii jj kk nn tmp*
whos
sim_sig_SrcActiv
sim_sig_IntNoise
disp(~isequal(sim_sig_SrcActiv.sigSRC_pre,sim_sig_SrcActiv.sigSRC_pst))
disp(~isequal(sim_sig_SrcActiv.sigSRC_pre,sim_sig_IntNoise.sigSRC_pre))
disp(~isequal(sim_sig_SrcActiv.sigSRC_pre,sim_sig_IntNoise.sigSRC_pst))
tmp_arrC = [rawSize(sim_sig_SrcActiv.sigSRC_pre,[1,2,3]); rawSize(sim_sig_SrcActiv.sigSRC_pst,[1,2,3]); rawSize(sim_sig_IntNoise.sigSRC_pre,[1,2,3]); rawSize(sim_sig_IntNoise.sigSRC_pst,[1,2,3])];
tmp_varN = {'Timesamples','Signals','Realizations'};
tmp_rowN = {'sim_sig_SrcActiv.sigSRC_pre','sim_sig_SrcActiv.sigSRC_pst','sim_sig_IntNoise.sigSRC_pre','sim_sig_IntNoise.sigSRC_pst'};
disp(' ');disp(array2table(tmp_arrC,'VariableNames',tmp_varN,'RowNames',tmp_rowN));
tmp_arrC = [size(sim_sig_SrcActiv.A00),size(sim_sig_SrcActiv.A00,2)/size(sim_sig_SrcActiv.A00,1);size(sim_sig_SrcActiv.A01),size(sim_sig_SrcActiv.A01,2)/size(sim_sig_SrcActiv.A01,1);size(sim_sig_IntNoise.A00),size(sim_sig_IntNoise.A00,2)/size(sim_sig_IntNoise.A00,1);size(sim_sig_IntNoise.A01),size(sim_sig_IntNoise.A01,2)/size(sim_sig_IntNoise.A01,1)];
tmp_varN = {'Signals','Signals_cdot_ModOrd','ModOrd'};
tmp_rowN = {'sim_sig_SrcActiv.A00','sim_sig_SrcActiv.A01','sim_sig_IntNoise.A00','sim_sig_IntNoise.A01'};
disp(' ');disp(array2table(tmp_arrC,'VariableNames',tmp_varN,'RowNames',tmp_rowN));
clearvars ii jj kk nn tmp*
close all
figure(55);clf;set(gcf, 'Position', SETUP.DISP);
subplot(4,1,1);cccSim___rawDispA00(sim_sig_SrcActiv.A00);title('sim\_sig00.A00');
subplot(4,1,2);cccSim___rawDispA00(sim_sig_SrcActiv.A01);title('sim\_sig00.A01');
subplot(4,1,3);cccSim___rawDispA00(sim_sig_IntNoise.A00);title('sim\_sig30.A00');
subplot(4,1,4);cccSim___rawDispA00(sim_sig_IntNoise.A01);title('sim\_sig30.A01');
set(gcf,'color','w');
close all
figure(100);clf;set(gcf, 'Position', SETUP.DISP);cccSim___rawPlotA00(sim_sig_SrcActiv.A00);
set(gcf,'color','w');
figure(101);clf;set(gcf, 'Position', SETUP.DISP);cccSim___rawPlotA00(sim_sig_SrcActiv.A01);
set(gcf,'color','w');
figure(130);clf;set(gcf, 'Position', SETUP.DISP);cccSim___rawPlotA00(sim_sig_IntNoise.A00);
set(gcf,'color','w');
figure(131);clf;set(gcf, 'Position', SETUP.DISP);cccSim___rawPlotA00(sim_sig_IntNoise.A01);
set(gcf,'color','w');
clearvars sim_sig_BcgNoise tmp*
sim_sig_BcgNoise = cccSim___makeSimSig(SETUP,3,'MVAR based signal for background (biological) noise activity sources');
clearvars ii jj kk nn tmp*
whos
sim_sig_SrcActiv
sim_sig_IntNoise
sim_sig_BcgNoise
chkSim___tg_s01_PRE___001(SETUP)
% Check size of signals generated
tmp_arrC = [rawSize(sim_sig_SrcActiv.sigSRC_pre,[1,2,3]); rawSize(sim_sig_SrcActiv.sigSRC_pst,[1,2,3]); rawSize(sim_sig_IntNoise.sigSRC_pre,[1,2,3]); rawSize(sim_sig_IntNoise.sigSRC_pst,[1,2,3]); rawSize(sim_sig_BcgNoise.sigSRC_pre,[1,2,3]); rawSize(sim_sig_BcgNoise.sigSRC_pst,[1,2,3])];
tmp_varN = {'Timesamples','Signals','Realizations'};
tmp_rowN = {'sim_sig_SrcActiv.sigSRC_pre','sim_sig_SrcActiv.sigSRC_pst','sim_sig_IntNoise.sigSRC_pre','sim_sig_IntNoise.sigSRC_pst','sim_sig_BcgNoise.sigSRC_pre','sim_sig_BcgNoise.sigSRC_pst'};
disp(' ');disp(array2table(tmp_arrC,'VariableNames',tmp_varN,'RowNames',tmp_rowN));
jj = 1
kk = 1
tmp_arrC = corrcoef(sim_sig_SrcActiv.sigSRC_pre(:,jj,kk),sim_sig_BcgNoise.sigSRC_pre(:,jj,kk));
tmp_varN = {'sim_sig_SrcActiv_a','sim_sig_BcgNoise_a'};
tmp_rowN = tmp_varN;
disp(' ');disp(array2table(tmp_arrC,'VariableNames',tmp_varN,'RowNames',tmp_rowN));disp('CYBERCRAFT:: This should be negligible')
close all
figure(105);clf;set(gcf,'Position',SETUP.DISP);subplot(2,1,1);cccSim___rawDispA00(sim_sig_BcgNoise.A00);title('A00');subplot(2,1,2);cccSim___rawDispA00(sim_sig_BcgNoise.A01);title('A01');
set(gcf,'color','w');
figure(100);clf;set(gcf,'Position',SETUP.DISP);cccSim___rawPlotA00(sim_sig_SrcActiv.A01);
set(gcf,'color','w');
figure(130);clf;set(gcf,'Position',SETUP.DISP);cccSim___rawPlotA00(sim_sig_IntNoise.A01);
set(gcf,'color','w');
figure(140);clf;set(gcf,'Position',SETUP.DISP);cccSim___rawPlotA00(sim_sig_BcgNoise.A01);
set(gcf,'color','w');
close all
figure(55);clf;set(gcf, 'Position', SETUP.DISP);
subplot(6,1,1);cccSim___rawDispA00(sim_sig_SrcActiv.A00);title('sim\_sig\_SrcActiv.A00')
subplot(6,1,2);cccSim___rawDispA00(sim_sig_SrcActiv.A01);title('sim\_sig\_SrcActiv.A01')
subplot(6,1,3);cccSim___rawDispA00(sim_sig_IntNoise.A00);title('sim\_sig\_IntNoise.A00')
subplot(6,1,4);cccSim___rawDispA00(sim_sig_IntNoise.A01);title('sim\_sig\_IntNoise.A01')
subplot(6,1,5);cccSim___rawDispA00(sim_sig_BcgNoise.A00);title('sim\_sig\_BcgNoise.A00')
subplot(6,1,6);cccSim___rawDispA00(sim_sig_BcgNoise.A01);title('sim\_sig\_BcgNoise.A01')
clearvars sim_sig_MesNoise tmp*
sim_sig_MesNoise.inf = 'RANDN based signal for measurment noise';
sim_sig_MesNoise.sigSNS_pre = randn([SETUP.n00,size(sel_ele.chanpos,1),SETUP.K00]);
sim_sig_MesNoise.sigSNS_pst = randn([SETUP.n00,size(sel_ele.chanpos,1),SETUP.K00]);
whos
sim_sig_SrcActiv
sim_sig_IntNoise
sim_sig_BcgNoise
sim_sig_MesNoise
tmp_arrC = [rawSize(sim_sig_SrcActiv.sigSRC_pre,[1,2,3]); rawSize(sim_sig_SrcActiv.sigSRC_pst,[1,2,3]); rawSize(sim_sig_IntNoise.sigSRC_pre,[1,2,3]); rawSize(sim_sig_IntNoise.sigSRC_pst,[1,2,3]); rawSize(sim_sig_BcgNoise.sigSRC_pre,[1,2,3]); rawSize(sim_sig_BcgNoise.sigSRC_pst,[1,2,3]); rawSize(sim_sig_MesNoise.sigSNS_pre,[1,2,3]); rawSize(sim_sig_MesNoise.sigSNS_pst,[1,2,3])];
tmp_varN = {'Timesamples','Signals','Realizations'};
tmp_rowN = {'sim_sig_SrcActiv.sigSRC_pre','sim_sig_SrcActiv.sigSRC_pst','sim_sig_IntNoise.sigSRC_pre','sim_sig_IntNoise.sigSRC_pst','sim_sig_BcgNoise.sigSRC_pre','sim_sig_BcgNoise.sigSRC_pst','sim_sig_BcgNoise.sigSNS_pre','sim_sig_BcgNoise.sigSNS_pst'};
disp(' ');disp(array2table(tmp_arrC,'VariableNames',tmp_varN,'RowNames',tmp_rowN));disp('NB: measurement noise should be already in electrode space not in the source space.')
chkSim___tg_s01_PRE___001(SETUP)
chkSim___tg_s02_PRE___001(SETUP)
close all
figure('Color',[1,1,1]);clf;
sel_geo_deep_icosahedron642
trisurf( ...
sel_geo_deep_icosahedron642.tri, ...
sel_geo_deep_icosahedron642.pnt(:,1), ...
sel_geo_deep_icosahedron642.pnt(:,2), ...
sel_geo_deep_icosahedron642.pnt(:,3), ...
'facealpha',0.2, ...
'facecolor','m', ...
'edgecolor','m', ...
'edgealpha',0.4);
hold on
ccrender([-150,150],'finish','matte')
close all
figure('Color',[1,1,1]);clf;
sel_geo_deep_thalami
trisurf( ...
sel_geo_deep_thalami.tri, ...
sel_geo_deep_thalami.pnt(:,1), ...
sel_geo_deep_thalami.pnt(:,2), ...
sel_geo_deep_thalami.pnt(:,3), ...
'facealpha',0.05, ...
'facecolor','m', ...
'edgecolor','m', ...
'edgealpha',0.1);
hold on
ccrender([-50,50],'finish','matte')
ft_plot_mesh( ...
sel_atl, ...
'facecolor',[0.9 0.9 0.9], ...
'facealpha',0.0, ...
'edgecolor',[0.7 0.7 0.7], ...
'edgealpha',0.2);
hold on
ccrender([-160,160],'finish','matte')
view(0,90)
view(90,0)
view(180,0)
close all
figure('Color',[1,1,1]);clf;
ft_plot_mesh( ...
sel_msh.bnd(1), ...
'facecolor','m', ...
'facealpha',0.1, ...
'edgecolor','m', ...
'edgealpha',0.05);
ccrender([-160,160],'finish','matte')
ft_plot_mesh( ...
sel_msh.bnd(2), ...
'facecolor','c', ...
'facealpha',0.1, ...
'edgecolor','c', ...
'edgealpha',0.1);
ccrender([-160,160],'finish','matte')
ft_plot_mesh( ...
sel_msh.bnd(3), ...
'facecolor','skin', ...
'facealpha',0.1, ...
'edgecolor','red', ...
'edgealpha',0.2);
ccrender([-160,160],'finish','matte')
stem3( sel_ele.elecpos(:,1), ...
sel_ele.elecpos(:,2), ...
sel_ele.elecpos(:,3), ...
'.m','filled','LineStyle','none','MarkerSize',12);
% for ii=1:size(tmp_elecpos,1),sel_elec00.label{ii,1} = sprintf('Ch%02d',ii);end
text( sel_ele.elecpos(:,1) * 1.2, ...
sel_ele.elecpos(:,2) * 1.2, ...
sel_ele.elecpos(:,3) * 1.2, ...
sel_ele.label,'Color',[0.5 0.3 0.5],'FontSize',10);
view(0,90)
view(90,0)
view(180,0)
view(225,60)
clearvars sim_geo_cort ii jj kk nn tmp*
sim_geo_cort.numROIs = size(sel_atl.Atlas(sel_atl.atl).Scouts,2);
if SETUP.rROI
disp('CC: using random ROI locations')
sim_geo_cort.lstROIs = randsample([1:sim_geo_cort.numROIs],size(SETUP.SRCS,1));
else
disp('CC: using predefined ROI locations')
sim_geo_cort.lstROIs = [31 30 34 29 56 51 84 53 83 58 57];
sim_geo_cort.lstROIs = sim_geo_cort.lstROIs(1:size(SETUP.SRCS,1));
end
if SETUP.rPNT
disp('CC: using random source locations')
else
disp('CC: using predefined source locations')
sim_geo_cort.lstROIs(1) = 31;
end
SETUP.rROI
SETUP.rPNT
sim_geo_cort
sim_geo_cort.lstROIs
chkSim___tg_s01_PRE___001(SETUP)
chkSim___tg_s02_Geometry___01_RandSamp(sim_geo_cort,sel_atl)
function chkSim___tg_s02_Geometry___01_RandSamp(sim_geo_cort,sel_atl)
fprintf('\n');
disp(cell2table([num2cell(sim_geo_cort.lstROIs)',{sel_atl.Atlas(sel_atl.atl).Scouts(sim_geo_cort.lstROIs).Label}'],'VariableNames',{'ROI_number','ROI_name'}));
fprintf('\n');
end
Select all vertices and all triangles for each ROI.
sim_geo_cort.indROIs = [];
for ii = 1:length(sim_geo_cort.lstROIs)
tmp_roi = sim_geo_cort.lstROIs(ii);
sim_geo_cort.indROIs(ii).pntNum00 = sel_atl.Atlas(sel_atl.atl).Scouts(tmp_roi).Vertices';
sim_geo_cort.indROIs(ii).triNum00 = find(all(ismember(sel_atl.tri,sim_geo_cort.indROIs(ii).pntNum00),2));
end
clearvars ii jj kk nn tmp*
sim_geo_cort
sim_geo_cort.indROIs
ii = 1
sim_geo_cort.indROIs(ii)
sim_geo_cort.indROIs(ii).pntNum00
sim_geo_cort.indROIs(ii).triNum00
function chkSim___tg_s02_Geometry___02_Indices___Plot_000(sel_atl)
ft_plot_mesh( ...
sel_atl, ...
'facecolor',[0.95 0.95 0.95], ...
'facealpha',0.1, ...
'edgecolor',[0.85 0.85 0.85], ...
'edgealpha',0.2);
hold on
ccrender([-160,160],'finish','matte')
end
function chkSim___tg_s02_Geometry___02_Indices___Plot_001(sel_atl,sim_geo_cort)
hold on
tmp_colorBar = lines(length(sim_geo_cort.lstROIs));
for ii = 1:length(sim_geo_cort.lstROIs)
trisurf( ...
sel_atl.tri(sim_geo_cort.indROIs(ii).triNum00,:),...
sel_atl.pnt(:,1)-0.1*sel_atl.vn1(:,1), ...
sel_atl.pnt(:,2)-0.1*sel_atl.vn1(:,2), ...
sel_atl.pnt(:,3)-0.1*sel_atl.vn1(:,3), ...
'facealpha',0.2, ...
'facecolor',tmp_colorBar(ii,:), ...
'edgecolor',tmp_colorBar(ii,:), ...
'edgealpha',0.4);
hold on
end
ccrender([-160,160],'finish','matte')
% caxis([1 length(sim_geo_cort.lstROIs)]);
% colormap(tmp_colorBar)
% colorbar
end
function chkSim___tg_s02_Geometry___02_Indices___Plot_002(sim_geo_cort)
hold on
ii = 1
qh = quiver3(...
sim_geo_cort.pos_orig{ii}(:,1),sim_geo_cort.pos_orig{ii}(:,2),sim_geo_cort.pos_orig{ii}(:,3),...
sim_geo_cort.ori_orig{ii}(:,1),sim_geo_cort.ori_orig{ii}(:,2),sim_geo_cort.ori_orig{ii}(:,3),...
1,...
'color','k');
set(qh,'linewidth',1);
qh = quiver3(...
sim_geo_cort.pos_pert{ii}(:,1),sim_geo_cort.pos_pert{ii}(:,2),sim_geo_cort.pos_pert{ii}(:,3),...
sim_geo_cort.ori_pert{ii}(:,1),sim_geo_cort.ori_pert{ii}(:,2),sim_geo_cort.ori_pert{ii}(:,3),...
1,':',...
'color','k');
set(qh,'linewidth',1);
ii = 2
qh = quiver3(...
sim_geo_cort.pos_orig{ii}(:,1),sim_geo_cort.pos_orig{ii}(:,2),sim_geo_cort.pos_orig{ii}(:,3),...
sim_geo_cort.ori_orig{ii}(:,1),sim_geo_cort.ori_orig{ii}(:,2),sim_geo_cort.ori_orig{ii}(:,3),...
1,...
'color','r');
set(qh,'linewidth',1);
qh = quiver3(...
sim_geo_cort.pos_pert{ii}(:,1),sim_geo_cort.pos_pert{ii}(:,2),sim_geo_cort.pos_pert{ii}(:,3),...
sim_geo_cort.ori_pert{ii}(:,1),sim_geo_cort.ori_pert{ii}(:,2),sim_geo_cort.ori_pert{ii}(:,3),...
1,':',...
'color','r');
set(qh,'linewidth',1);
ii = 3
qh = quiver3(...
sim_geo_cort.pos_orig{ii}(:,1),sim_geo_cort.pos_orig{ii}(:,2),sim_geo_cort.pos_orig{ii}(:,3),...
sim_geo_cort.ori_orig{ii}(:,1),sim_geo_cort.ori_orig{ii}(:,2),sim_geo_cort.ori_orig{ii}(:,3),...
1,...
'color','b');
set(qh,'linewidth',1);
qh = quiver3(...
sim_geo_cort.pos_pert{ii}(:,1),sim_geo_cort.pos_pert{ii}(:,2),sim_geo_cort.pos_pert{ii}(:,3),...
sim_geo_cort.ori_pert{ii}(:,1),sim_geo_cort.ori_pert{ii}(:,2),sim_geo_cort.ori_pert{ii}(:,3),...
1,':',...
'color','b');
set(qh,'linewidth',1);
ccrender([-160,160],'finish','matte')
end
Figure giving ultimate situation view.
close all
chkSim___tg_s01_PRE___001(SETUP)
chkSim___tg_s02_Geometry___01_RandSamp(sim_geo_cort,sel_atl)
chkSim___tg_s02_Geometry___02_Indices___Plot_000(sel_atl)
chkSim___tg_s02_Geometry___02_Indices___Plot_001(sel_atl,sim_geo_cort)
view(-115,55)
view(-115,-55)
view(0,60)
view(45,60)
view(90,60)
view(135,60)
view(180,60)
view(225,60)
view(270,60)
view(315,60)
view(315,-20)
view(135,45)
sim_geo_cort.distrROIs = SETUP.SRCS;
sim_geo_cort.distrROIsSum = sum(SETUP.SRCS,2);
sim_geo_cort.bulkSRC = {};
for ii = 1:length(sim_geo_cort.distrROIsSum)
sim_geo_cort.bulkSRC{ii,1} = randsample(sim_geo_cort.indROIs(ii).pntNum00,sim_geo_cort.distrROIsSum(ii));
end
sim_geo_cort.splitSRC = {};
for ii = 1:length(sim_geo_cort.distrROIsSum)
sim_geo_cort.splitSRC(ii,:) = mat2cell(sim_geo_cort.bulkSRC{ii},sim_geo_cort.distrROIs(ii,:),1)';
end
if 0
sim_geo_cort.splitSRC(1,1)
sim_geo_cort.splitSRC{1,1}
end
if SETUP.rPNT
disp('CC: using random source locations')
else
disp('CC: using predefined source locations')
sim_geo_cort.splitSRC{1,1} = [5441;5506;5481;5283;5144;5063;5823;5987;6154;6065;6166;6346;6087;6446;6367;6726;6859;6829;6613;6698;6966;7162;7227;6949;6996;7365;7247;7480];
sim_geo_cort.splitSRC{1,1} = sim_geo_cort.splitSRC{1,1}(1:SETUP.SRCS(1,1));
if 0
sim_geo_cort.splitSRC{1,1}
end
end
sim_geo_cort.mergeSRC = {};
for ii = 1:3
sim_geo_cort.mergeSRC{ii} = cat(1,sim_geo_cort.splitSRC{:,ii});
end
sim_geo_cort.pos_orig = {};
sim_geo_cort.ori_orig = {};
sim_geo_cort.pos_pert = {};
sim_geo_cort.ori_pert = {};
for ii = 1:3
sim_geo_cort.pos_orig{ii} = sel_atl.pnt(sim_geo_cort.mergeSRC{ii},:);
sim_geo_cort.ori_orig{ii} = sel_atl.vn1(sim_geo_cort.mergeSRC{ii},:);
sim_geo_cort.pos_pert{ii} = sim_geo_cort.pos_orig{ii};
sim_geo_cort.ori_pert{ii} = sim_geo_cort.ori_orig{ii};
end
clearvars ii jj kk nn tmp*
SETUP.SRCS
SETUP.DEEP
sim_geo_cort
sim_geo_cort.distrROIs
sim_geo_cort.distrROIsSum
sim_geo_cort.bulkSRC
sim_geo_cort.splitSRC
sim_geo_cort.distrROIs
sim_geo_cort.mergeSRC
sim_geo_cort.mergeSRC{1}
chkSim___tg_s02_Geometry___01_RandSamp(sim_geo_cort,sel_atl)
chkSim___tg_s01_PRE___001(SETUP)
sim_geo_cort.mergeSRC
sim_geo_cort.pos_orig
sim_geo_cort.ori_orig
sim_geo_cort.pos_pert
sim_geo_cort.ori_pert
clearvars sel_geo_deep ii jj kk ll nn tmp*
sim_geo_deep = sel_geo_deep_icosahedron642;
sim_geo_deep = sel_geo_deep_thalami;
sim_geo_deep.bulkSRC = randsample(1:size(sim_geo_deep.pnt,1),sum(SETUP.DEEP))';
[sim_geo_deep.mergeSRC{1},sim_geo_deep.mergeSRC{2},sim_geo_deep.mergeSRC{3}] = deal([]);
sim_geo_deep.mergeSRC{1,1} = sim_geo_deep.bulkSRC(1:SETUP.DEEP(1));
sim_geo_deep.mergeSRC{1,2} = sim_geo_deep.bulkSRC(1+SETUP.DEEP(1):SETUP.DEEP(1)+SETUP.DEEP(2));
sim_geo_deep.mergeSRC{1,3} = sim_geo_deep.bulkSRC(1+SETUP.DEEP(1)+SETUP.DEEP(2):SETUP.DEEP(1)+SETUP.DEEP(2)+SETUP.DEEP(3));
sim_geo_deep.pos_orig = {};
sim_geo_deep.ori_orig = {};
sim_geo_deep.pos_pert = {};
sim_geo_deep.ori_pert = {};
for ii = 1:3
sim_geo_deep.pos_orig{ii} = sim_geo_deep.pnt(sim_geo_deep.mergeSRC{ii},:);
sim_geo_deep.ori_orig{ii} = sim_geo_deep.vn1(sim_geo_deep.mergeSRC{ii},:);
sim_geo_deep.pos_pert{ii} = sim_geo_deep.pos_orig{ii};
sim_geo_deep.ori_pert{ii} = sim_geo_deep.ori_orig{ii};
end
clearvars ii jj kk nn tmp*
sim_geo_deep
trisurf( ...
sim_geo_deep.tri, ...
sim_geo_deep.pnt(:,1), ...
sim_geo_deep.pnt(:,2), ...
sim_geo_deep.pnt(:,3), ...
'facealpha',0.2, ...
'facecolor','m', ...
'edgecolor','m', ...
'edgealpha',0.4);
ccrender([-150,150],'finish','matte')
whos sel_*
whos sim_*
whos sim_geo*
chkSim___tg_s01_PRE___001(SETUP)
disp(SETUP.DEEP)
[sim_geo_deep.bulkSRC,[sim_geo_deep.mergeSRC{1}; sim_geo_deep.mergeSRC{2};sim_geo_deep.mergeSRC{3}]]
sim_geo_deep
sim_geo_cort
sim_geo = sim_geo_cort;
for ii = 1:3,
sim_geo.pos_orig{ii} = [sim_geo.pos_orig{ii};sim_geo_deep.pos_orig{ii}];
sim_geo.ori_orig{ii} = [sim_geo.ori_orig{ii};sim_geo_deep.ori_orig{ii}];
sim_geo.pos_pert{ii} = [sim_geo.pos_pert{ii};sim_geo_deep.pos_pert{ii}];
sim_geo.ori_pert{ii} = [sim_geo.ori_pert{ii};sim_geo_deep.ori_pert{ii}];
end;
whos sim* sel*
sim_geo
sim_geo_deep
sim_geo_cort
To facilitate bug tracking and for explanatory/illustrative purpose it might be plausible to modify the perturbation parameters here.
SETUP.CUBE
SETUP.CUBE = 30
SETUP.CUBE
SETUP.CONE
SETUP.CONE = pi/16
SETUP.CONE
for ii = 1:3
sim_geo.ori_pert{ii} = normr(rawSphToCart(rawCartToSph(sim_geo.ori_orig{ii}) + [SETUP.CONE*0.5*(2*rand(size(sim_geo.ori_pert{ii},1),2)-1),zeros(size(sim_geo.ori_pert{ii},1),1)]));
end
tmp_count = 0;
for ii = 1:3
for jj = 1:size(sim_geo.pos_pert{ii},1)
tmp_srcInside = false;
while ~tmp_srcInside
tmp_count = tmp_count+1;
sim_geo.pos_pert{ii}(jj,:) = sim_geo.pos_orig{ii}(jj,:)+SETUP.CUBE*0.5*(2*rand([1,3])-1);
tmp_srcInside = bounding_mesh(sim_geo.pos_pert{ii}(jj,:),sel_msh.bnd(1).pnt,sel_msh.bnd(1).tri);
end
end
end
if SETUP.TELL
disp(['CYBERCRAFT:: Perturbation took ',tmp_count,' iterations'])
end
clearvars ii jj kk nn tmp*
NOTE 1: The following procedure is applied to the orientation vectors of activity sources: in order to achieve their random perturbation:
- Original orientations are vertex vormal vectors of the triangulation mesh;
- Coordinates are transformed to spherical;
SETUP.CONEis the max perturbation (angle of max rotation) considered- Multiplied here by random number ∈ [-1,1];
- Azimuth [TH] and elevation [PHI] are modified by adding the above random numbers;
- Radius is not changed (the added matrix has last column padded with zeros);
- Coordinates are transformed back to the cartesian system.
ii = 1
SETUP.CUBE
min(sim_geo.pos_orig{ii}-sim_geo.pos_pert{ii})
max(sim_geo.pos_orig{ii}-sim_geo.pos_pert{ii})
SETUP.CONE
min(rawCartToSph(sim_geo.ori_orig{ii})-rawCartToSph(sim_geo.ori_pert{ii}))
max(rawCartToSph(sim_geo.ori_orig{ii})-rawCartToSph(sim_geo.ori_pert{ii}))
sim_geo
whos
close all
% figure(1)
figure('Color',[1 1 1]);
clf
chkSim___tg_s01_PRE___001(SETUP)
chkSim___tg_s02_Geometry___01_RandSamp(sim_geo,sel_atl)
% mesh for ROIs on cortex
chkSim___tg_s02_Geometry___02_Indices___Plot_001(sel_atl,sim_geo)
% mesh for deep sources ROI
chkSim___tg_s02_Geometry___02_Indices___Plot_003(sim_geo_deep)
% sources
chkSim___tg_s02_Geometry___02_Indices___Plot_002(sim_geo)
legend('SrcActiv\_orig','SrcActiv\_pert','IntNoise\_orig','IntNoise\_pert','BcgNoise\_orig','BcgNoise\_pert')
% chkSim___tg_s02_Geometry___02_Indices___Plot_000(sel_atl)
view(0,90)
view(90,0)
view(180,0)
view(-115,55)
view(-115,-55)
view(0,60)
view(45,60)
view(90,60)
view(135,60)
view(180,60)
view(225,60)
view(270,60)
view(315,60)
view(315,-20)
view(135,45)
This one works for deep sources!!!
function chkSim___tg_s02_Geometry___02_Indices___Plot_003(sim_geo_deep)
hold on
trisurf( ...
sim_geo_deep.tri, ...
sim_geo_deep.pnt(:,1), ...
sim_geo_deep.pnt(:,2), ...
sim_geo_deep.pnt(:,3), ...
'facealpha',0.1, ...
'facecolor','m', ...
'edgecolor','m', ...
'edgealpha',0.2);
ccrender([-150,150],'finish','matte')
end
Please note that if meshes for brain, skull and scalp change
the headmodel also needs to be recalculated
Leadfields are named with reference to signals.
| Variable name | Description |
|---|---|
| sim_lfg_SrcActiv_orig | Original activity of interest (biological) |
| sim_lfg_IntNoise_orig | Original interference (biological) noise |
| sim_lfg_BcgNoise_orig | Original background (biological) noise |
| sim_lfg_SrcActiv_pert | Perturbed activity of interest (biological) |
| sim_lfg_IntNoise_pert | Perturbed interference (biological) noise |
| sim_lfg_BcgNoise_pert | Perturbed background (biological) noise |
% WARNING:: suppress warning about MATLAB version and CFG tracking
if 0
TMP_S = warning('off','all');
end
tmp_cfg = [];
tmp_cfg.grid.unit = 'mm';
tmp_cfg.reducerank = 'no';
tmp_cfg.normalize = 'no';
tmp_cfg.elec = sel_ele;
tmp_cfg.vol = sel_vol;
ii = 1;
tmp_cfg.grid.pos = sim_geo.pos_orig{ii};
tmp_cfg.grid.mom = transpose(sim_geo.ori_orig{ii});
sim_lfg_SrcActiv_orig.lfg = ft_prepare_leadfield(tmp_cfg);
sim_lfg_SrcActiv_orig.LFG = cat(2,sim_lfg_SrcActiv_orig.lfg.leadfield{:});
tmp_cfg.grid.pos = sim_geo.pos_pert{ii};
tmp_cfg.grid.mom = transpose(sim_geo.ori_pert{ii});
sim_lfg_SrcActiv_pert.lfg = ft_prepare_leadfield(tmp_cfg);
sim_lfg_SrcActiv_pert.LFG = cat(2,sim_lfg_SrcActiv_pert.lfg.leadfield{:});
ii = 2;
tmp_cfg.grid.pos = sim_geo.pos_orig{ii};
tmp_cfg.grid.mom = transpose(sim_geo.pos_orig{ii});
sim_lfg_IntNoise_orig.lfg = ft_prepare_leadfield(tmp_cfg);
sim_lfg_IntNoise_orig.LFG = cat(2,sim_lfg_IntNoise_orig.lfg.leadfield{:});
tmp_cfg.grid.pos = sim_geo.pos_pert{ii};
tmp_cfg.grid.mom = transpose(sim_geo.pos_pert{ii});
sim_lfg_IntNoise_pert.lfg = ft_prepare_leadfield(tmp_cfg);
sim_lfg_IntNoise_pert.LFG = cat(2,sim_lfg_IntNoise_pert.lfg.leadfield{:});
ii = 3;
tmp_cfg.grid.pos = sim_geo.pos_orig{ii};
tmp_cfg.grid.mom = transpose(sim_geo.pos_orig{ii});
sim_lfg_BcgNoise_orig.lfg = ft_prepare_leadfield(tmp_cfg);
sim_lfg_BcgNoise_orig.LFG = cat(2,sim_lfg_BcgNoise_orig.lfg.leadfield{:});
tmp_cfg.grid.pos = sim_geo.pos_pert{ii};
tmp_cfg.grid.mom = transpose(sim_geo.pos_pert{ii});
sim_lfg_BcgNoise_pert.lfg = ft_prepare_leadfield(tmp_cfg);
sim_lfg_BcgNoise_pert.LFG = cat(2,sim_lfg_BcgNoise_pert.lfg.leadfield{:});
clearvars ii jj kk nn tmp*
clearvars ii jj kk nn tmp*
for kk = 1:SETUP.K00
sim_sig_SrcActiv.sigSNS_pre(:,:,kk) = (sim_lfg_SrcActiv_orig.LFG*sim_sig_SrcActiv.sigSRC_pre(:,:,kk)')';
sim_sig_SrcActiv.sigSNS_pst(:,:,kk) = (sim_lfg_SrcActiv_orig.LFG*sim_sig_SrcActiv.sigSRC_pst(:,:,kk)')';
sim_sig_IntNoise.sigSNS_pre(:,:,kk) = (sim_lfg_IntNoise_orig.LFG*sim_sig_IntNoise.sigSRC_pre(:,:,kk)')';
sim_sig_IntNoise.sigSNS_pst(:,:,kk) = (sim_lfg_IntNoise_orig.LFG*sim_sig_IntNoise.sigSRC_pst(:,:,kk)')';
sim_sig_BcgNoise.sigSNS_pre(:,:,kk) = (sim_lfg_BcgNoise_orig.LFG*sim_sig_BcgNoise.sigSRC_pre(:,:,kk)')';
sim_sig_BcgNoise.sigSNS_pst(:,:,kk) = (sim_lfg_BcgNoise_orig.LFG*sim_sig_BcgNoise.sigSRC_pst(:,:,kk)')';
end
clearvars ii jj kk nn tmp*
We introduced desired SNR by adjusting power for:
- biological interference,
- background biological activity,
- measurement noise.
In order to allow for evaluation of:
- reconstruction accuracy, and
- functional dependencies estimation (PDC, DTF, dDTF etc.)
we did not change the power of the sources signal!
clearvars ii jj kk nn tmp*
sim_sig_AdjSNRs.SrcActivPre = sim_sig_SrcActiv.sigSNS_pre;
sim_sig_AdjSNRs.SrcActivPst = sim_sig_SrcActiv.sigSNS_pst;
sim_sig_AdjSNRs.IntNoisePre = sim_sig_IntNoise.sigSNS_pre;
sim_sig_AdjSNRs.IntNoisePst = sim_sig_IntNoise.sigSNS_pst;
sim_sig_AdjSNRs.BcgNoisePre = sim_sig_BcgNoise.sigSNS_pre;
sim_sig_AdjSNRs.BcgNoisePst = sim_sig_BcgNoise.sigSNS_pst;
sim_sig_AdjSNRs.MesNoisePre = sim_sig_MesNoise.sigSNS_pre;
sim_sig_AdjSNRs.MesNoisePst = sim_sig_MesNoise.sigSNS_pst;
for kk = 1:SETUP.K00
sim_sig_AdjSNRs.IntNoisePre(:,:,kk) = rawAdjTotSNRdB(sim_sig_AdjSNRs.SrcActivPre(:,:,kk),sim_sig_AdjSNRs.IntNoisePre(:,:,kk),SETUP.SINR);
sim_sig_AdjSNRs.IntNoisePst(:,:,kk) = rawAdjTotSNRdB(sim_sig_AdjSNRs.SrcActivPst(:,:,kk),sim_sig_AdjSNRs.IntNoisePst(:,:,kk),SETUP.SINR);
sim_sig_AdjSNRs.BcgNoisePre(:,:,kk) = rawAdjTotSNRdB(sim_sig_AdjSNRs.SrcActivPre(:,:,kk),sim_sig_AdjSNRs.BcgNoisePre(:,:,kk),SETUP.SBNR);
sim_sig_AdjSNRs.BcgNoisePst(:,:,kk) = rawAdjTotSNRdB(sim_sig_AdjSNRs.SrcActivPst(:,:,kk),sim_sig_AdjSNRs.BcgNoisePst(:,:,kk),SETUP.SBNR);
sim_sig_AdjSNRs.MesNoisePre(:,:,kk) = rawAdjTotSNRdB(sim_sig_AdjSNRs.SrcActivPre(:,:,kk),sim_sig_AdjSNRs.MesNoisePre(:,:,kk),SETUP.SMNR);
sim_sig_AdjSNRs.MesNoisePst(:,:,kk) = rawAdjTotSNRdB(sim_sig_AdjSNRs.SrcActivPst(:,:,kk),sim_sig_AdjSNRs.MesNoisePst(:,:,kk),SETUP.SMNR);
end
clearvars ii jj kk nn tmp*
whos
SETUP
sim_sig_SrcActiv % signal
sim_sig_IntNoise % interference
sim_sig_BcgNoise % background
sim_sig_MesNoise % measurment
sim_sig_AdjSNRs % merged signals (bulk) with SNR adjusted
kk = 1
2*pow2db(rawNrm(sim_sig_AdjSNRs.SrcActivPre(:,:,kk)))
2*pow2db(rawNrm(sim_sig_AdjSNRs.SrcActivPst(:,:,kk)))
rawSNRdB(sim_sig_AdjSNRs.SrcActivPre(:,:,kk),sim_sig_AdjSNRs.IntNoisePre(:,:,kk))
rawSNRdB(sim_sig_AdjSNRs.SrcActivPst(:,:,kk),sim_sig_AdjSNRs.IntNoisePst(:,:,kk))
rawSNRdB(sim_sig_AdjSNRs.SrcActivPre(:,:,kk),sim_sig_AdjSNRs.BcgNoisePre(:,:,kk))
rawSNRdB(sim_sig_AdjSNRs.SrcActivPst(:,:,kk),sim_sig_AdjSNRs.BcgNoisePst(:,:,kk))
rawSNRdB(sim_sig_AdjSNRs.SrcActivPre(:,:,kk),sim_sig_AdjSNRs.MesNoisePre(:,:,kk))
rawSNRdB(sim_sig_AdjSNRs.SrcActivPst(:,:,kk),sim_sig_AdjSNRs.MesNoisePst(:,:,kk))
% sim_sig_SrcActiv.sigSNS_pre should remain unchanged after SNR adjustments
isequal(sim_sig_AdjSNRs.SrcActivPre,sim_sig_SrcActiv.sigSNS_pre)
% sim_sig_SrcActiv.sigSNS_pre should remain unchanged after SNR adjustments
isequal(sim_sig_AdjSNRs.SrcActivPst,sim_sig_SrcActiv.sigSNS_pst)
whos
clearvars ii jj kk nn tmp*
% clearvars ii jj kk nn tmp*
% Note: in this and subsequent subsections we aim at following closely notation from the paper.
if SETUP.H_Src_pert
H_Src = sim_lfg_SrcActiv_pert.LFG;
else
H_Src = sim_lfg_SrcActiv_orig.LFG;
end
if SETUP.H_Int_pert
H_Int = sim_lfg_IntNoise_pert.LFG;
else
H_Int = sim_lfg_IntNoise_orig.LFG;
end
% reduced-rank leadfield for patch constraints
[U_H_Int Si_H_Int V_H_Int] = svd(H_Int);
for ii = SETUP.IntLfgRANK+1:size(H_Int,2)
Si_H_Int(ii,ii)=0;
end
H_Int = U_H_Int*Si_H_Int*V_H_Int';
y_Pre = SETUP.SigPre*sim_sig_AdjSNRs.SrcActivPre ...
+ SETUP.IntPre*sim_sig_AdjSNRs.IntNoisePre ...
+ SETUP.BcgPre*sim_sig_AdjSNRs.BcgNoisePre ...
+ SETUP.MesPre*sim_sig_AdjSNRs.MesNoisePre;
y_Pst = SETUP.SigPst*sim_sig_AdjSNRs.SrcActivPst ...
+ SETUP.IntPst*sim_sig_AdjSNRs.IntNoisePst ...
+ SETUP.BcgPst*sim_sig_AdjSNRs.BcgNoisePst ...
+ SETUP.MesPst*sim_sig_AdjSNRs.MesNoisePst;
% background activity (as in Pre segment) in Pst segment
y_PstNoise = SETUP.SigPre*sim_sig_AdjSNRs.SrcActivPst ...
+ SETUP.IntPre*sim_sig_AdjSNRs.IntNoisePst ...
+ SETUP.BcgPre*sim_sig_AdjSNRs.BcgNoisePst ...
+ SETUP.MesPre*sim_sig_AdjSNRs.MesNoisePst;
N = cov(reshape(permute(y_Pre,[1 3 2]),[],size(y_Pre,2),1));
% N = cov(reshape(permute(y_PstNoise,[1 3 2]),[],size(y_PstNoise,2),1));
R = cov(reshape(permute(y_Pst,[1 3 2]),[],size(y_Pst,2),1));
G = H_Src'*pinv(N)*H_Src;
S = H_Src'*pinv(R)*H_Src;
G_SrcInt = [H_Src H_Int]'*pinv(N)*[H_Src H_Int];
S_SrcInt = [H_Src H_Int]'*pinv(R)*[H_Src H_Int];
C_SrcInt = pinv(S_SrcInt) - pinv(G_SrcInt);
% Estimated C for both regular and nulling filters
C_NL = C_SrcInt(1:size(H_Src,2),1:size(H_Src,2));
C = pinv(S)-pinv(G);
clear U_H_Int Si_H_Int V_H_Int G_SrcInt S_SrcInt;
whos
size(H_Src)
size(H_Int)
size(y_Pre)
size(y_Pst)
size(N)
size(R)
size(G)
size(S)
size(C)
% mean values of signals
y_Pst_RSH = reshape(permute(y_Pst,[1 3 2]),[],size(y_Pst,2),1);
max(abs(mean(y_Pst_RSH)))
clear y_Pst_RSH;
SrcPst = reshape(permute(sim_sig_AdjSNRs.SrcActivPst,[1 3 2]),[],size(sim_sig_AdjSNRs.SrcActivPst,2),1);
max(abs(mean(SrcPst)))
clear SrcPst;
IntPst = reshape(permute(sim_sig_AdjSNRs.IntNoisePst,[1 3 2]),[],size(sim_sig_AdjSNRs.IntNoisePst,2),1);
max(abs(mean(IntPst)))
clear IntPst;
BcgPst = reshape(permute(sim_sig_AdjSNRs.BcgNoisePst,[1 3 2]),[],size(sim_sig_AdjSNRs.BcgNoisePst,2),1);
max(abs(mean(BcgPst)))
clear BcgPst;
BcgPre = reshape(permute(sim_sig_AdjSNRs.BcgNoisePre,[1 3 2]),[],size(sim_sig_AdjSNRs.BcgNoisePre,2),1);
% note: max and min values vector across electrodes
max(abs(BcgPre))'
min(abs(BcgPre))'
clear BcgPre;
MesPst = reshape(permute(sim_sig_AdjSNRs.MesNoisePst,[1 3 2]),[],size(sim_sig_AdjSNRs.MesNoisePst,2),1);
max(abs(mean(MesPst)))
clear MesPst;
% correlation between signals check and # of highly correlated electrodes for Src, Bcg, Mes signals
SrcPst = reshape(permute(sim_sig_AdjSNRs.SrcActivPst,[1 3 2]),[],size(sim_sig_AdjSNRs.SrcActivPst,2),1);
BcgPst = reshape(permute(sim_sig_AdjSNRs.BcgNoisePst,[1 3 2]),[],size(sim_sig_AdjSNRs.BcgNoisePst,2),1);
Corr_SrcPst = corrcoef(SrcPst);
for ii = 1:length(Corr_SrcPst)
[a k]=find(Corr_SrcPst(:,ii)>.9);
x(ii)=sum(k);
end
x
x=[];
Corr_BcgPst = corrcoef(BcgPst);
for ii = 1:length(Corr_BcgPst)
[a k]=find(Corr_BcgPst(:,ii)>.9);
x(ii)=sum(k);
end
x
x=[];
Corr_SrcBcgPst = corrcoef([SrcPst BcgPst]);
Cross_SrcBcgPst = Corr_SrcBcgPst(length(Corr_SrcPst)+1:end,1:length(Corr_SrcPst));
max([max(Cross_SrcBcgPst) abs(min(Cross_SrcBcgPst))])
clear Corr_SrcPst Corr_BcgPst Corr_SrcBcgPst Cross_SrcBcgPst x;
MesPst = reshape(permute(sim_sig_AdjSNRs.MesNoisePst,[1 3 2]),[],size(sim_sig_AdjSNRs.MesNoisePst,2),1);
Corr_SrcPst = corrcoef(SrcPst);
for ii = 1:length(Corr_SrcPst)
[a k]=find(Corr_SrcPst(:,ii)>.9);
x(ii)=sum(k);
end
x
x=[];
Corr_MesPst = corrcoef(MesPst);
for ii = 1:length(Corr_MesPst)
[a k]=find(Corr_MesPst(:,ii)>.9);
x(ii)=sum(k);
end
x
x=[];
Corr_SrcMesPst = corrcoef([SrcPst MesPst]);
Cross_SrcMesPst = Corr_SrcMesPst(length(Corr_SrcPst)+1:end,1:length(Corr_SrcPst));
max([max(Cross_SrcMesPst) abs(min(Cross_SrcMesPst))])
clear Corr_MesPst Corr_SrcMesPst Cross_SrcMesPst x;
Corr_BcgMesPst = corrcoef([BcgPst MesPst]);
Cross_BcgMesPst = Corr_BcgMesPst(length(Corr_SrcPst)+1:end,1:length(Corr_SrcPst));
max([max(Cross_BcgMesPst) abs(min(Cross_BcgMesPst))])
clear BcgPst MesPst Corr_BcgMesPst Cross_BcgMesPst;
IntPst = reshape(permute(sim_sig_AdjSNRs.IntNoisePst,[1 3 2]),[],size(sim_sig_AdjSNRs.IntNoisePst,2),1);
Corr_SrcIntPst = corrcoef([SrcPst IntPst]);
Cross_SrcIntPst = Corr_SrcIntPst(length(Corr_SrcPst)+1:end,1:length(Corr_SrcPst));
max([max(Cross_SrcIntPst) abs(min(Cross_SrcIntPst))])
clear ii a k SrcPst Corr_SrcPst IntPst Corr_SrcIntPst Cross_SrcIntPst;
rec_opt
whos
% NL filter as proposed in Hui & Leahy 2010
H_SrcInt = [H_Src H_Int];
P_NL = [eye(size(H_Src,2)) zeros(size(H_Src,2),size(H_Int,2))];
% EIG-LCMV filter
[UR, SiR, ~] = svd(R);
P_EIG = UR(:,1:SETUP.RANK_EIG)*UR(:,1:SETUP.RANK_EIG)';
clear UR SiR;
% sMVP filters
H_Src_R = pinv(sqrtm(R))*H_Src;
H_Src_N = pinv(sqrtm(N))*H_Src;
% sMVP_MSE filter
% note: S and G use only H_Src, as we want to recover only activity of interest, not interference
% note: pinv(H_Src_R)*pinv(H_Src_R)' = pinv(S)
K_MSE = pinv(S)-2*C;
[U_K_MSE, W_K_MSE] = eig(K_MSE);
[~, p_K_MSE] = sort(diag(W_K_MSE));
% note: for selecting optimal rank, Theobald's Th. can be used without evaluation of the cost function as is done here
U_K_MSE = U_K_MSE(:,p_K_MSE);
sMVP_MSE_ranks = zeros(1,size(H_Src,2));
for jj = 1:size(H_Src,2)
sMVP_MSE_ranks(jj) = trace(U_K_MSE(:,1:jj)*U_K_MSE(:,1:jj)'*K_MSE);
end
[~, sMVP_MSE_rank_opt] = min(sMVP_MSE_ranks);
P_sMVP_MSE_opt = U_K_MSE(:,1:sMVP_MSE_rank_opt)*U_K_MSE(:,1:sMVP_MSE_rank_opt)';
rec_opt.ranks.sMVP_MSE = sMVP_MSE_rank_opt;
clear K_MSE U_K_MSE W_K_MSE p_K_MSE sMVP_MSE_ranks sMVP_MSE_rank_opt;
% sMVP_R filter
% note: S and G use only H_Src, as we want to recover only activity of interest, not interference
% note: pinv(H_Src_R)*pinv(H_Src_R)' = pinv(S)
K_MSE = pinv(S)-2*C;
K_R = K_MSE+2*C;
[U_K_R, W_K_R] = eig(K_R);
[~, p_K_R] = sort(diag(W_K_R));
U_K_R = U_K_R(:,p_K_R);
sMVP_R_ranks = zeros(1,size(H_Src,2));
for jj = 1:size(H_Src,2)
sMVP_R_ranks(jj) = trace(U_K_R(:,1:jj)*U_K_R(:,1:jj)'*K_MSE);
end
[~, sMVP_R_rank_opt] = min(sMVP_R_ranks);
P_sMVP_R_opt = U_K_R(:,1:sMVP_R_rank_opt)*U_K_R(:,1:sMVP_R_rank_opt)';
rec_opt.ranks.sMVP_R = sMVP_R_rank_opt;
clear K_MSE K_R U_K_R W_K_R p_K_R sMVP_R_ranks;
% sMVP_N filter
% note: S and G use only H_Src, as we want to recover only activity of interest, not interference
% note: pinv(H_Src_N)*pinv(H_Src_N)' = pinv(G)
K_MSE_N = pinv(G)-C;
K_N = K_MSE_N+C;
[U_K_N, W_K_N] = eig(K_N);
[~, p_K_N] = sort(diag(W_K_N));
U_K_N = U_K_N(:,p_K_N);
sMVP_N_ranks = zeros(1,size(H_Src,2));
for jj = 1:size(H_Src,2)
sMVP_N_ranks(jj) = trace(U_K_N(:,1:jj)*U_K_N(:,1:jj)'*K_MSE_N);
end
[~, sMVP_N_rank_opt] = min(sMVP_N_ranks);
P_sMVP_N_opt = U_K_N(:,1:sMVP_N_rank_opt)*U_K_N(:,1:sMVP_N_rank_opt)';
rec_opt.ranks.sMVP_N = sMVP_N_rank_opt;
clear K_MSE_N K_N U_K_N W_K_N p_K_N sMVP_N_ranks;
% sMVP filters
clear H_Src_R H_Src_N;
% sMVP-NL filters
H_Src_R = pinv(sqrtm(R))*H_Src;
H_Src_N = pinv(sqrtm(N))*H_Src;
H_Int_R = pinv(sqrtm(R))*H_Int;
H_Int_N = pinv(sqrtm(N))*H_Int;
P_sMVP_NL_R = eye(size(H_Int_R,1))-H_Int_R*pinv(H_Int_R);
P_sMVP_NL_N = eye(size(H_Int_N,1))-H_Int_N*pinv(H_Int_N);
% sMVP_NL_MSE filter
K_NL_MSE = pinv(P_sMVP_NL_R*H_Src_R)*P_sMVP_NL_R*(pinv(P_sMVP_NL_R*H_Src_R))'-2*C_NL;
[U_K_NL_MSE, W_K_NL_MSE] = eig(K_NL_MSE);
[~, p_K_NL_MSE] = sort(diag(W_K_NL_MSE));
% note: for selecting optimal rank, Theobald's Th. can be used without evaluation of the cost function as is done here
U_K_NL_MSE = U_K_NL_MSE(:,p_K_NL_MSE);
sMVP_NL_MSE_ranks = zeros(1,size(H_Src,2));
for jj = 1:size(H_Src,2)
sMVP_NL_MSE_ranks(jj) = trace(U_K_NL_MSE(:,1:jj)*U_K_NL_MSE(:,1:jj)'*K_NL_MSE);
end
[~, sMVP_NL_MSE_rank_opt] = min(sMVP_NL_MSE_ranks);
P_sMVP_NL_MSE_opt = U_K_NL_MSE(:,1:sMVP_NL_MSE_rank_opt)*U_K_NL_MSE(:,1:sMVP_NL_MSE_rank_opt)';
rec_opt.ranks.sMVP_NL_MSE = sMVP_NL_MSE_rank_opt;
clear K_NL_MSE U_K_NL_MSE W_K_NL_MSE p_K_NL_MSE sMVP_NL_MSE_ranks sMVP_NL_MSE_rank_opt;
% sMVP_NL_R filter
K_NL_MSE = pinv(P_sMVP_NL_R*H_Src_R)*P_sMVP_NL_R*(pinv(P_sMVP_NL_R*H_Src_R))'-2*C_NL;
%K_NL_R = pinv(P_sMVP_NL_R*H_Src_R)*P_sMVP_NL_R*(pinv(P_sMVP_NL_R*H_Src_R))';
K_NL_R = K_NL_MSE+2*C_NL;
[U_K_NL_R, W_K_NL_R] = eig(K_NL_R);
[~, p_K_NL_R] = sort(diag(W_K_NL_R));
U_K_NL_R = U_K_NL_R(:,p_K_NL_R);
sMVP_NL_R_ranks = zeros(1,size(H_Src,2));
for jj = 1:size(H_Src,2)
sMVP_NL_R_ranks(jj) = trace(U_K_NL_R(:,1:jj)*U_K_NL_R(:,1:jj)'*K_NL_MSE);
end
[~, sMVP_NL_R_rank_opt] = min(sMVP_NL_R_ranks);
P_sMVP_NL_R_opt = U_K_NL_R(:,1:sMVP_NL_R_rank_opt)*U_K_NL_R(:,1:sMVP_NL_R_rank_opt)';
rec_opt.ranks.sMVP_NL_R = sMVP_NL_R_rank_opt;
clear K_NL_MSE K_NL_R U_K_NL_R W_K_NL_R p_K_NL_R sMVP_NL_R_ranks sMVP_NL_R_rank_opt;
% sMVP_NL_N filter
K_NL_MSE_N = pinv(P_sMVP_NL_N*H_Src_N)*P_sMVP_NL_N*(pinv(P_sMVP_NL_N*H_Src_N))'-C_NL;
%K_NL_N = pinv(P_sMVP_NL_N*H_Src_N)*P_sMVP_NL_N*(pinv(P_sMVP_NL_N*H_Src_N))';
K_NL_N = K_NL_MSE_N+C_NL;
[U_K_NL_N, W_K_NL_N] = eig(K_NL_N);
[~, p_K_NL_N] = sort(diag(W_K_NL_N));
U_K_NL_N = U_K_NL_N(:,p_K_NL_N);
sMVP_NL_N_ranks = zeros(1,size(H_Src,2));
for jj = 1:size(H_Src,2)
sMVP_NL_N_ranks(jj) = trace(U_K_NL_N(:,1:jj)*U_K_NL_N(:,1:jj)'*K_NL_MSE_N);
end
[~, sMVP_NL_N_rank_opt] = min(sMVP_NL_N_ranks);
% note the simplified notation for P_N_opt
P_sMVP_NL_N_opt = U_K_NL_N(:,1:sMVP_NL_N_rank_opt)*U_K_NL_N(:,1:sMVP_NL_N_rank_opt)';
rec_opt.ranks.sMVP_NL_N = sMVP_NL_N_rank_opt;
clear K_NL_MSE_N K_NL_N U_K_NL_N W_K_NL_N p_K_NL_N sMVP_NL_N_ranks sMVP_NL_N_rank_opt;
% sMVP-NL filters
clear H_Int_R H_Int_N;
% note: S and G use only H_Src, as we want to recover only activity of interest, not interference
rec_flt.LCMV_R = pinv(S)*H_Src'*pinv(R);
rec_flt.LCMV_N = pinv(G)*H_Src'*pinv(N);
% As proposed in Hui & Leahy 2010
rec_flt.NL = P_NL*pinv(H_SrcInt'*pinv(R)*H_SrcInt)*H_SrcInt'*pinv(R);
rec_flt.MMSE = C*H_Src'*pinv(R);
rec_flt.MMSE_INT = C_SrcInt(1:size(H_Src,2),:)*H_SrcInt'*pinv(R);
rec_flt.ZF = pinv(H_Src);
rec_flt.RANDN = randn(size(H_Src'));
rec_flt.ZEROS = zeros(size(H_Src'));
rec_flt.EIG_LCMV_R = rec_flt.LCMV_R*P_EIG;
rec_flt.EIG_LCMV_N = rec_flt.LCMV_N*P_EIG;
% not applicable in our case as it assumes access to covariance matrix of interference + background noise + measurement noise, whereas in our case N does not contain interference
rec_flt.sMVP_MSE = P_sMVP_MSE_opt*rec_flt.LCMV_R;
rec_flt.sMVP_R = P_sMVP_R_opt*rec_flt.LCMV_R;
rec_flt.sMVP_N = P_sMVP_N_opt*rec_flt.LCMV_N;
rec_flt.sMVP_NL_MSE = P_sMVP_NL_MSE_opt*pinv(P_sMVP_NL_R*H_Src_R)*P_sMVP_NL_R*pinv(sqrtm(R));
rec_flt.sMVP_NL_R = P_sMVP_NL_R_opt*pinv(P_sMVP_NL_R*H_Src_R)*P_sMVP_NL_R*pinv(sqrtm(R));
rec_flt.sMVP_NL_N = P_sMVP_NL_N_opt*pinv(P_sMVP_NL_N*H_Src_N)*P_sMVP_NL_N*pinv(sqrtm(N));
if(SETUP.fltREMOVE)
rec_flt = rmfield(rec_flt,'RANDN');
rec_flt = rmfield(rec_flt,'ZEROS');
end
rawFixStrJoin
clearvars R N C C_NL S G P_* H_* jj;
whos
rec_flt
rec_opt.ranks
Reconstruct signal in the sources using previously prepared spatial filters.
clearvars rec_sig ii jj kk nn tmp*
rec_sig.Original = sim_sig_SrcActiv.sigSRC_pst;
rec_sig.Dummy = sim_sig_SrcActiv.sigSRC_pre;
tmp_fltFields = fieldnames(rec_flt);
for nn = 1:length(tmp_fltFields),
for kk = 1:SETUP.K00,
rec_sig.(tmp_fltFields{nn})(:,:,kk) = (rec_flt.(tmp_fltFields{nn}) * y_Pst(:,:,kk)')';
end,
end
clearvars ii jj kk nn tmp*
whos
rec_flt
rec_opt
rec_sig
whos
Calculate MVAR model coefficients and PDC.
NB: Keep in mind that P00 is restricted to single value, i.e. the actual order of the MVAR model used to generate time-courses for signal of interest (it is not being estimated here).
clearvars ii jj kk nn rec_funDep_A00 rec_funDep_PDC
rec_funDep_A00.Original = sim_sig_SrcActiv.A00;
[~,rec_funDep_A00.Dummy,~,~,~,~] = arfit(sim_sig_SrcActiv.sigSRC_pre,sim_sig_SrcActiv.P00,sim_sig_SrcActiv.P00);
tmp_fltFields = fieldnames(rec_flt);
for nn = 1:length(tmp_fltFields),[~,rec_funDep_A00.(tmp_fltFields{nn}),~,~,~,~] = arfit(rec_sig.(tmp_fltFields{nn}),SETUP.P00,SETUP.P00);end
tmp_allFields = fieldnames(rec_funDep_A00);
for nn = 1:length(tmp_allFields),rec_funDep_PDC.(tmp_allFields{nn}) = abs(PDC(rec_funDep_A00.(tmp_allFields{nn}),SETUP.PDC_RES));end
whos
whos rec*
rec_flt
rec_opt
rec_opt.ranks
rec_sig
rec_funDep_A00
rec_funDep_PDC
% pre-computed leadfields for source candidates
load('./mat/sel_src.mat');
structure = sel_src.lfg;
% lf is the array of cells containing all leafields within ROI
lf = structure.leadfield;
% position of sources
pos_Src = sim_lfg_SrcActiv_orig.lfg.pos;
[~, pos_Src_ind]=ismember(pos_Src,structure.pos,'rows');
% Results will be stored in structs with fields:
% leadfield, source positions, localization error, rank selected
% 1 = MAI(2011), 2 = MAI(2013), 3 = MPZ(2011), 4 = MPZ(2013)
% 5 = MAI_RR_I, 6 = MAI_RR_C, 7 = MPZ_RR_I, 8 = MPZ_RR_C
RES = repmat(struct('H',[],'sources',[]),8,1);
% to be used for the proposed heuristic rank-selection criterion
% leading eigenvalues of R*pinv(N) are equal to those of C*G0+1
x = sort(eig(R*pinv(N)),'descend');
x = x(1:sum(SETUP.SRCS(:,1)));
% percentage of variance explained
y = x./sum(x);
% percent of explained variance
rs_opt = min(find(cumsum(y)>0.8));
clear x y;
for i = 1:sum(SETUP.SRCS(:,1))
tic
i
% for REAL WORLD APPLICABLE activity indices, rank is fixed
RES_tmp = repmat({zeros(length(lf),1)}, 1, 8);
rs = min([i rs_opt]);
for k = 1:8
for j = 1:length(lf)
h_curr = lf{j};
% we use H obtained for a given index at previous iteration with respect to i,
% starting with null matrix
h = [RES(k).H h_curr];
G = h'*pinv(N)*h;
S = h'*pinv(R)*h;
T = h'*pinv(R)*N*pinv(R)*h;
[U D ~] = svd(S);
MAI_eigenv = sort(real(eig(G*pinv(S))),'descend');
MPZ_eigenv = sort(real(eig(S*pinv(T))),'descend');
if ismember(k,[6 8]) & i>rs_opt
% Csq is rank-deficient if h = lf{j} matches previously found sources
Csq = real(sqrtm(pinv(S)-pinv(G)));
h2 = h*Csq;
G2 = h2'*pinv(N)*h2;
S2 = h2'*pinv(R)*h2;
T2 = h2'*pinv(R)*N*pinv(R)*h2;
[U2 D2 ~] = svd(S2);
end
% REAL-WORLD APPLICABLE family of indices
switch k
case 1
% MAI(2011) index
RES_tmp{k}(j) = sum(MAI_eigenv);
case 2
% MAI(2013) index
RES_tmp{k}(j) = sum(MAI_eigenv(1:rs));
case 3
% MPZ(2011) index
RES_tmp{k}(j) = sum(MPZ_eigenv);
case 4
% MPZ(2013) index
RES_tmp{k}(j) = sum(MPZ_eigenv(1:rs));
case 5
% MAI_RR_I index
P = U(:,1:rs)*U(:,1:rs)';
RES_tmp{k}(j) = trace(G*pinv(S)*P);
clear P;
case 6
% MAI_RR_C index
if i<=rs_opt
RES_tmp{k}(j) = RES_tmp{k-1}(j);
elseif i>rs_opt % in this case rs=rs_opt
P2 = U2(:,1:rs_opt)*U2(:,1:rs_opt)';
RES_tmp{k}(j) = trace(G2*pinv(S2)*P2);
clear P2;
end
case 7
% MPZ_RR_I index
P = U(:,1:rs)*U(:,1:rs)';
RES_tmp{k}(j) = trace(S*pinv(T)*P);
clear P;
case 8
% MPZ_RR_C index
if i<=rs_opt
RES_tmp{k}(j) = RES_tmp{k-1}(j);
elseif i>rs_opt % in this case rs=rs_opt
P2 = U2(:,1:rs_opt)*U2(:,1:rs_opt)';
RES_tmp{k}(j) = trace(S2*pinv(T2)*P2);
clear P2;
end
end
% explicit clear to emphasize that all matrices
% are created anew at each iteration for each j and k
clear h_curr h G S T Csq h2 G2 S2 T2 U D U2 D2 P P2;
% below: end of j = 1:length(lf) loop
end
[~, w]=max(RES_tmp{k});
[dis dis_ind] = pdist2(pos_Src,structure.pos(w,:),'chebychev','Smallest',1);
% store results
RES(k).H = [RES(k).H lf{w}];
RES(k).sources(i).candidate_number = w;
RES(k).sources(i).rank_opt = rs;
RES(k).sources(i).source_discovered = dis_ind;
RES(k).sources(i).distance = dis;
% below: end of k = 1:8 loop
end
% below: end of i = 1:SETUP.RANK_EIG loop
toc
end
clear structure lf pos_Src;Replace original signal with sine and the reconstructed signals with waves that produce easy to predict error values…
if 0
clearvars rec* ii jj kk nn tmp*
tmp_arg = linspace(0,8*pi,SETUP.n00)';
size(tmp_arg)
rec_sig.Original = repmat(sin(tmp_arg),[1,sum(SETUP.SRCS(:,1)),SETUP.K00]);
rec_sig.Dummy = repmat(cos(tmp_arg),[1,sum(SETUP.SRCS(:,1)),SETUP.K00]);
rec_sig.Test00 = rec_sig.Original
rec_sig.Test01 = -rec_sig.Original
rec_sig.Test02 = rec_sig.Original+0.2*randn(size(rec_sig.Original))
rec_sig.Test03 = -rec_sig.Original+0.2*randn(size(rec_sig.Original))
rec_sig.Test04 = zeros(size(rec_sig.Original))
end
clearvars ii jj kk nn tmp*
whos
whos rec*
rec_sig
clearvars ii jj kk nn tmp*
tmp_allFields = fieldnames(rec_sig);
for nn = 1:length(tmp_allFields)
for ii = 1:SETUP.n00
for jj = 1:SETUP.K00
rec_sig.(tmp_allFields{nn})(ii,:,jj)=rec_sig.(tmp_allFields{nn})(ii,:,jj)/norm(rec_sig.(tmp_allFields{nn})(ii,:,jj));
end
end
end
clearvars ii jj kk nn tmp* rec_sigAmp_ErrEuclid
tmp_allFields = fieldnames(rec_sig);
for nn = 1:length(tmp_allFields)
tmp_error = zeros(SETUP.n00,SETUP.K00);
for ii = 1:SETUP.n00
for jj = 1:SETUP.K00
tmp_error(ii,jj)=norm(rec_sig.(tmp_allFields{1})(ii,:,jj)-rec_sig.(tmp_allFields{nn})(ii,:,jj))^2;
end
end
rec_sigAmp_ErrEuclid.(tmp_allFields{nn})=mean(mean(tmp_error));
endtmp_allFields = fieldnames(rec_sig);
disp(tmp_allFields{1})
isequal(rec_sig.(tmp_allFields{1}),sim_sig_SrcActiv.sigSRC_pst)
norm(rec_sig.(tmp_allFields{1})(:)-sim_sig_SrcActiv.sigSRC_pst(:))
size(rec_sig.(tmp_allFields{1}))
size(sim_sig_SrcActiv.sigSRC_pre)
size(sim_sig_AdjSNRs.SrcActivPre)
clearvars ii jj kk nn tmp* rec_sigAmp_ErrCorrCf
tmp_allFields = fieldnames(rec_sig);
for nn = 1:length(tmp_allFields)
tmp_error = zeros(sum(SETUP.SRCS(:,1)),SETUP.K00);
for kk = 1:sum(SETUP.SRCS(:,1))
for jj = 1:SETUP.K00
tmp_x=corrcoef(rec_sig.(tmp_allFields{1})(:,kk,jj),rec_sig.(tmp_allFields{nn})(:,kk,jj));
tmp_error(kk,jj)=tmp_x(1,end);
end
end
rec_sigAmp_ErrCorrCf.(tmp_allFields{nn}) = mean(mean(tmp_error));
end
clearvars ii jj kk nn tmp* rec_funDep_A00_ErrMonday rec_funDep_A00_ErrEuclid
tmp_allFields = fieldnames(rec_sig);
for nn = 1:length(tmp_allFields)
rec_funDep_A00_ErrEuclid.(tmp_allFields{nn}) = (norm(rec_funDep_A00.(tmp_allFields{1}) - rec_funDep_A00.(tmp_allFields{nn}),'fro')^2);
end
clearvars ii jj kk nn tmp* rec_funDep_A00_ErrCorrCf
tmp_allFields = fieldnames(rec_sig);
for nn = 1:length(tmp_allFields)
tmp_error = zeros(sum(SETUP.SRCS(:,1)),1);
for kk = 1:sum(SETUP.SRCS(:,1))
tmp_x=corrcoef(rec_funDep_A00.(tmp_allFields{1})(kk,:),rec_funDep_A00.(tmp_allFields{nn})(kk,:));
tmp_error(kk)=tmp_x(1,end);
end
rec_funDep_A00_ErrCorrCf.(tmp_allFields{nn}) = mean(tmp_error);
end
clearvars ii jj kk nn tmp* rec_funDep_PDC_err rec_funDep_PDC_ErrEuclid
tmp_allFields = fieldnames(rec_sig);
for nn = 1:length(tmp_allFields)
tmp_error = zeros(sum(SETUP.SRCS(:,1)),1);
for ii = 1:sum(SETUP.SRCS(:,1))
tmp_error(ii) = norm(squeeze(rec_funDep_PDC.(tmp_allFields{1})(ii,:,:)) - squeeze(rec_funDep_PDC.(tmp_allFields{nn})(ii,:,:)),'fro')^2;
end
rec_funDep_PDC_ErrEuclid.(tmp_allFields{nn}) = mean(tmp_error);
end
clearvars ii jj kk nn tmp* rec_funDep_PDC_ErrCorrCf
tmp_allFields = fieldnames(rec_sig);
for nn = 1:length(tmp_allFields)
tmp_error = zeros(sum(SETUP.SRCS(:,1)),1);
for ii = 1:sum(SETUP.SRCS(:,1))
tmp_x=corrcoef(reshape(squeeze(rec_funDep_PDC.(tmp_allFields{1})(ii,:,:)),[],1),reshape(squeeze(rec_funDep_PDC.(tmp_allFields{nn})(ii,:,:)),[],1));
tmp_error(ii)=tmp_x(1,end);
end
rec_funDep_PDC_ErrCorrCf.(tmp_allFields{nn}) = mean(tmp_error);
end
clearvars ii jj kk nn tmp* rec*vec
rec_sigAmp_ErrEuclid_vec = cellfun(@(fn) rec_sigAmp_ErrEuclid.(fn), fieldnames(rec_sigAmp_ErrEuclid), 'UniformOutput', false);
rec_sigAmp_ErrEuclid_vec = squeeze(vertcat(rec_sigAmp_ErrEuclid_vec{:}));
rec_sigAmp_ErrCorrCf_vec = cellfun(@(fn) rec_sigAmp_ErrCorrCf.(fn), fieldnames(rec_sigAmp_ErrCorrCf), 'UniformOutput', false);
rec_sigAmp_ErrCorrCf_vec = squeeze(vertcat(rec_sigAmp_ErrCorrCf_vec{:}));
rec_funDep_A00_ErrEuclid_vec = cellfun(@(fn) rec_funDep_A00_ErrEuclid.(fn), fieldnames(rec_funDep_A00_ErrEuclid), 'UniformOutput', false);
rec_funDep_A00_ErrEuclid_vec = squeeze(vertcat(rec_funDep_A00_ErrEuclid_vec{:}));
rec_funDep_A00_ErrCorrCf_vec = cellfun(@(fn) rec_funDep_A00_ErrCorrCf.(fn), fieldnames(rec_funDep_A00_ErrCorrCf), 'UniformOutput', false);
rec_funDep_A00_ErrCorrCf_vec = squeeze(vertcat(rec_funDep_A00_ErrCorrCf_vec{:}));
rec_funDep_PDC_ErrEuclid_vec = cellfun(@(fn) rec_funDep_PDC_ErrEuclid.(fn), fieldnames(rec_funDep_PDC_ErrEuclid), 'UniformOutput', false);
rec_funDep_PDC_ErrEuclid_vec = squeeze(vertcat(rec_funDep_PDC_ErrEuclid_vec{:}));
rec_funDep_PDC_ErrCorrCf_vec = cellfun(@(fn) rec_funDep_PDC_ErrCorrCf.(fn), fieldnames(rec_funDep_PDC_ErrCorrCf), 'UniformOutput', false);
rec_funDep_PDC_ErrCorrCf_vec = squeeze(vertcat(rec_funDep_PDC_ErrCorrCf_vec{:}));
clearvars ii jj kk nn tmp* rec_res
if ~isequal(fieldnames(rec_sigAmp_ErrEuclid),fieldnames(rec_sigAmp_ErrCorrCf))
error('check fieldnames')
elseif ~isequal(fieldnames(rec_sigAmp_ErrEuclid),fieldnames(rec_funDep_A00_ErrEuclid))
error('check fieldnames')
elseif ~isequal(fieldnames(rec_sigAmp_ErrEuclid),fieldnames(rec_funDep_A00_ErrCorrCf))
error('check fieldnames')
elseif ~isequal(fieldnames(rec_sigAmp_ErrEuclid),fieldnames(rec_funDep_PDC_ErrEuclid))
error('check fieldnames')
elseif ~isequal(fieldnames(rec_sigAmp_ErrEuclid),fieldnames(rec_funDep_PDC_ErrCorrCf))
error('check fieldnames')
else
rec_res.table_arrC = [ rec_sigAmp_ErrEuclid_vec, rec_sigAmp_ErrCorrCf_vec, rec_funDep_A00_ErrEuclid_vec, rec_funDep_A00_ErrCorrCf_vec, rec_funDep_PDC_ErrEuclid_vec, rec_funDep_PDC_ErrCorrCf_vec];
rec_res.table_varN = {'rec_sigAmp_ErrEuclid_vec','rec_sigAmp_ErrCorrCf_vec','rec_funDep_A00_ErrEuclid_vec','rec_funDep_A00_ErrCorrCf_vec','rec_funDep_PDC_ErrEuclid_vec','rec_funDep_PDC_ErrCorrCf_vec'};
rec_res.table_rowN = fieldnames(rec_sigAmp_ErrEuclid);
rec_res.table = array2table(rec_res.table_arrC,'VariableNames',rec_res.table_varN,'RowNames',rec_res.table_rowN);
end
rec_res.table_arrC = [ rec_sigAmp_ErrEuclid_vec, rec_sigAmp_ErrCorrCf_vec ];
rec_res.table_varN = {'rec_sigAmp_ErrEuclid_vec','rec_sigAmp_ErrCorrCf_vec'};
rec_res.table_rowN = fieldnames(rec_sigAmp_ErrEuclid);
rec_res.table = array2table(rec_res.table_arrC,'VariableNames',rec_res.table_varN,'RowNames',rec_res.table_rowN);
whos rec*
whos rec*vec
rec_res
sortrows(rec_res.table,'rec_sigAmp_ErrEuclid_vec')
sortrows(rec_res.table,'rec_sigAmp_ErrCorrCf_vec','descend')
sortrows(rec_res.table,'rec_funDep_A00_ErrEuclid_vec')
sortrows(rec_res.table,'rec_funDep_A00_ErrCorrCf_vec','descend')
sortrows(rec_res.table,'rec_funDep_PDC_ErrEuclid_vec')
sortrows(rec_res.table,'rec_funDep_PDC_ErrCorrCf_vec','descend')
whos rec*
whos rec*vec
rec_sigAmp_ErrEuclid
rec_sigAmp_ErrCorrCf
rec_funDep_A00_ErrEuclid
rec_funDep_A00_ErrMonday
rec_funDep_PDC_ErrEuclid
rec_funDep_PDC_ErrMonday
rec_sig
rec_sigAmp_ErrEuclid
[r_opt2,r_opt3,size(vSRC_pst,2)]
TMP_COV = cov(sim_sig_SrcActiv.sigSRC_pst(:,:,1));
TMP_COR = corrcoef(sim_sig_SrcActiv.sigSRC_pst(:,:,1));
figure(1);clf;
rawImgSC(TMP_COV);
tmp_colorMapMat = rawHotColdColorMap(2048);
caxis([-max(abs(min(min(TMP_COV))),abs(max(max(TMP_COV)))) max(abs(min(min(TMP_COV))),abs(max(max(TMP_COV))))]);
colormap(tmp_colorMapMat);
colorbar;
figure(2);clf;
rawImgSC(TMP_COR);
tmp_colorMapMat = rawHotColdColorMap(2048);
caxis([-max(abs(min(min(TMP_COR))),abs(max(max(TMP_COR)))) max(abs(min(min(TMP_COR))),abs(max(max(TMP_COR))))]);
colormap(tmp_colorMapMat);
colorbar;
sort(eig(K))
figure(501);clf;set(gcf, 'Position', SETUP.DISP);rawImgSC(sort(eig(K)));colorbar;
rank(Hp)
svd(Hp)
svd(S1)
R
rec_sig
rec_sigAmp_ErrEuclid
[r_opt2,r_opt3,size(vSRC_pst,2)]
rec_funDep_A00
rec_funDep_A00_ErrEuclid
rec_funDep_PDC
rec_funDep_PDC_ErrEuclid
rec_funDep_PDC_ErrMonday
figure(100000)
set(gcf, 'Position', SETUP.DISP);
subplot(2,1,1)
plot(rec_funDep_PDC_ErrEuclid);colorbar;
subplot(2,1,2)
plot(rec_funDep_PDC_ErrMonday);colorbar;
figure(30000);clf;
set(gcf, 'Position', SETUP.DISP);
for nn = 1:length(tmp_allFields)
subplot(length(tmp_allFields),1,nn)
imagesc(rec_funDep_A00.(tmp_allFields{nn}));
tmp_colorMapMat = jet(255);
tmp_colorMapMat(128,:) = [1 1 1];
caxis([-max(abs(min(min(rec_funDep_A00.(tmp_allFields{nn})))),abs(max(max(rec_funDep_A00.(tmp_allFields{nn}))))) max(abs(min(min(rec_funDep_A00.(tmp_allFields{nn})))),abs(max(max(rec_funDep_A00.(tmp_allFields{nn})))))]);
colormap(tmp_colorMapMat);
colorbar;
hold on;
tmp_stem = 0.5+sim_sig_SrcActiv.S00:sim_sig_SrcActiv.S00:sim_sig_SrcActiv.S00*(sim_sig_SrcActiv.P00-1)+0.5;
tmp_vals = sim_sig_SrcActiv.S00+0.5*ones(size(tmp_stem));
stem(tmp_stem,tmp_vals,'Color','k','LineWidth',0.5,'Marker', 'none');
set(gca,'XTick',[1:sim_sig_SrcActiv.S00*sim_sig_SrcActiv.P00])
title(tmp_allFields{nn})
end
rec_sigAmp_ErrEuclid
[r_opt2,r_opt3,size(vSRC_pst,2)]
rec_funDep_PDC_ErrEuclid
rec_funDep_PDC_ErrMonday
for nn = 1:length(tmp_allFields)
figure(40000+nn);clf;
set(gcf, 'Position', SETUP.DISP);
rawPlotPDC(rec_funDep_PDC.(tmp_allFields{nn}),SETUP.PDC_RES)
title(tmp_allFields{nn})
end
Single run for tangled matlab-scripts.
The below code block can be used to test simulations single run.
NB: before running the code below remember to tangle code blocks after any modifications !!!
% if exist('fPath'), cd(fPath); else, try, cd('~/supFunSim/'); catch, warningMessage = 'Problem encoutered while trying to change working directory to ''~/supFunSim/''.'; end; end; disp(['CYBERCRAFT:: pwd is: ',pwd])
run('./tg_s07_BATCH___01_Prelude.m'); whos
SETUP
chkSim___tg_s01_PRE___001(SETUP);
chkSim___tg_s01_PRE___000_PARSE_SETUP(SETUP,sel_atl);
run('./tg_s07_BATCH___02_Signals.m'); whos
run('./tg_s07_BATCH___03_Leadfields.m'); whos
% if BcgNoise (biological background noise) is to be white Gaussian,
% then please uncomment the following two lines.
% sim_sig_BcgNoise.sigSRC_pre = randn(size(sim_sig_BcgNoise.sigSRC_pre));
% sim_sig_BcgNoise.sigSRC_pst = randn(size(sim_sig_BcgNoise.sigSRC_pst));
% if one wishes the biological noise on sensors to be white Gaussian
% sim_sig_BcgNoise.sigSNS_pre = randn(size(sim_sig_BcgNoise.sigSNS_pre));
% sim_sig_BcgNoise.sigSNS_pst = randn(size(sim_sig_BcgNoise.sigSNS_pst));
% rawImgSC(squareform(pdist(sim_lfg_SrcActiv_orig.lfg.pos,'euclidean')))
run('./tg_s07_BATCH___04_Preparations.m'); whos
switch SETUP.supSwitch
case {'rec'}
% Reconstruction
disp('Reconstruction')
run('./tg_s07_BATCH___05_Filters.m'); whos
% Check the results
sortrows(rec_res.table,'rec_sigAmp_ErrEuclid_vec')
sortrows(rec_res.table,'rec_sigAmp_ErrCorrCf_vec','descend')
sortrows(rec_res.table,'rec_funDep_A00_ErrEuclid_vec')
sortrows(rec_res.table,'rec_funDep_A00_ErrCorrCf_vec','descend')
sortrows(rec_res.table,'rec_funDep_PDC_ErrEuclid_vec')
sortrows(rec_res.table,'rec_funDep_PDC_ErrCorrCf_vec','descend')
rec_opt.ranks
case {'loc'}
% Localization
disp('Localization')
run('./tg_s07_BATCH___05_Localizers.m'); whos
otherwise
warning('Unexpected supSwitch value. No action taken, only data generated.')
end
After execution of tg_s07_BATCH___02_Signals.m the generated
timeseries can be swapped for simple waveforms (used for to
facilitate bug tracking and for explanatory/illustrative
purpose).
% (just trying simple signals over here)
tmp_arg = linspace(0,8*pi,SETUP.n00)';
size(tmp_arg)
sim_sig_SrcActiv.sigSRC_pre = repmat(sin(tmp_arg), [1,sum(SETUP.SRCS(:,1)),SETUP.K00]);
sim_sig_SrcActiv.sigSRC_pst = repmat(sin(tmp_arg), [1,sum(SETUP.SRCS(:,1)),SETUP.K00]);
sim_sig_IntNoise.sigSRC_pre = repmat(sin(tmp_arg*100),[1,sum(SETUP.SRCS(:,2)),SETUP.K00]);
sim_sig_IntNoise.sigSRC_pst = repmat(sin(tmp_arg*100),[1,sum(SETUP.SRCS(:,2)),SETUP.K00]);
sim_sig_BcgNoise.sigSRC_pre = repmat(sin(tmp_arg*10), [1,sum(SETUP.SRCS(:,3)),SETUP.K00]);
sim_sig_BcgNoise.sigSRC_pst = repmat(sin(tmp_arg*10), [1,sum(SETUP.SRCS(:,3)),SETUP.K00]);
size(sim_sig_SrcActiv.sigSRC_pre)
size(sim_sig_SrcActiv.sigSRC_pst)
size(sim_sig_BcgNoise.sigSRC_pre)
size(sim_sig_BcgNoise.sigSRC_pst)
close all
jj = 1
kk = 1
figure('Name','sim_sig_SrcActiv.sigSRC_pre'); plot(sim_sig_SrcActiv.sigSRC_pre(:,jj,kk))
figure('Name','sim_sig_BcgNoise.sigSRC_pre'); plot(sim_sig_BcgNoise.sigSRC_pre(:,jj,kk))
figure('Name','sim_sig_SrcActiv.sigSNS_pre'); plot(sim_sig_SrcActiv.sigSNS_pre(:,jj,kk))
figure('Name','sim_sig_SrcActiv.sigSNS_pst'); plot(sim_sig_SrcActiv.sigSNS_pst(:,jj,kk))
figure('Name','sim_sig_BcgNoise.sigSRC_pre'); plot(sim_sig_BcgNoise.sigSRC_pre(:,jj,kk))
figure('Name','sim_sig_BcgNoise.sigSRC_pst'); plot(sim_sig_BcgNoise.sigSRC_pst(:,jj,kk))
figure('Name','sim_sig_BcgNoise.sigSNS_pre'); plot(sim_sig_BcgNoise.sigSNS_pre(:,jj,kk))
figure('Name','sim_sig_BcgNoise.sigSNS_pst'); plot(sim_sig_BcgNoise.sigSNS_pst(:,jj,kk))
size(y_Pre)
size(y_Pst)
figure('Name','y_Pre'); plot(y_Pre(:,jj,kk))
figure('Name','y_Pst'); plot(y_Pst(:,jj,kk))
% (just trying simple signals over here)
sim_sig_BcgNoise.sigSRC_pre = randn(size(sim_sig_BcgNoise.sigSRC_pre))
sim_sig_BcgNoise.sigSRC_pst = sim_sig_BcgNoise.sigSRC_pre
% (just trying simple signals over here)
sim_sig_BcgNoise.sigSRC_pst = sim_sig_BcgNoise.sigSRC_pre
The following settings will overwritte simulation default settings. The devault settings are optimized for bug tracking and for explanatory/illustrative purpose. Please modify simulation settings only in this section.
% Tidy up and change working directory.
clc; close all;clearvars('-except','fPath');
if exist('fPath'),cd(fPath);else,try,cd('~/cc_overkill/git/supFunSim');catch,cd('~/supFunSim');end;end;
clear all;
close all;
clc;
run('./tg_s00_Prelude___01_Simulations_main_setup.m')
run('./tg_s00_Prelude___02_SNR_Adjustments.m')
% Rows of SETUP.SRCS reppresent ROIs.
% Cols of SETUP.SRCS represent SrcActiv, IntNoise and BcgNoise, respectively.
SETUP.rROI = logical(1); % random (1) or predefined (0) ROIs
SETUP.rPNT = logical(1); % random (1) or predefined (0) candidate points for source locations: if 0,
% number of sources as in SETUP.SRCS(1,1) will be fixed and in close locations
SETUP.SRCS = []; % Cortical sources (avoid placing more than 10 sources in single ROI)
SETUP.SRCS = [ SETUP.SRCS; 3 3 0 ];
SETUP.SRCS = [ SETUP.SRCS; 3 3 0 ];
SETUP.SRCS = [ SETUP.SRCS; 3 3 0 ];
SETUP.SRCS = [ SETUP.SRCS; 4 3 0 ];
SETUP.SRCS = [ SETUP.SRCS; 0 3 0 ];
SETUP.SRCS = [ SETUP.SRCS; 0 3 0 ];
SETUP.SRCS = [ SETUP.SRCS; 0 3 0 ];
SETUP.SRCS = [ SETUP.SRCS; 0 3 0 ];
SETUP.SRCS = [ SETUP.SRCS; 0 3 3 ];
SETUP.SRCS = [ SETUP.SRCS; 0 0 4 ];
SETUP.DEEP = [ 0 0 20 ]; % deep sources
SETUP.ERPs = 0; % Add ERPs (timelocked activity)
% Basic setting s for sources and noise.
SETUP.n00 = 500; % number of time samples per trial
SETUP.K00 = 1; % number of independent realizations of signal and noise based on generated MVAR model
% note: covariance matrix R of observed signal and noise covariance matrix N are
% estimated from samples originating from all realizations of signal and noise
SETUP.FIXED_SEED = 1; % Settings for seed selection
if SETUP.FIXED_SEED, SETUP.SEED = rng(1964);else,SETUP.SEED = rng(round(1e3*randn()^2*sum(clock)));end
SETUP.RANK_EIG = sum(SETUP.SRCS(:,1)); % rank of EIG-LCMV filter: set to number of active sources
SETUP.fltREMOVE = 1; % to keep (0) or remove (1) selected filters
SETUP.SHOWori = 0; % to show (1) or do not show (0) Original and Dummy signals on Figures
SETUP.IntLfgRANK = round(0.3*sum(SETUP.SRCS(:,2))); % rank of patch-constrained reduced-rank leadfield
SETUP.supSwitch = 'rec'; % 'rec': run reconstruction of sources activity, 'loc': find active sources
SETUP.CUBE = 20; % perturbation of the leadfields based on the shift of source position within a cube of given edge length (centered at the original leadfields positions), in [mm]
SETUP.CONE = pi/32; % perturbation of the leadfields based on the rotation of source orientation (azimuth TH, elevation PHI), in [rad]
SETUP.H_Src_pert = 1; % use original (0) or perturbed (1) leadfield for signal reconstruction
SETUP.H_Int_pert = 0; % use original (0) or perturbed (1) leadfield for nulling constrains
SETUP.SINR = 0; % signal to interference noise power ratio expressed in dB (both measured on electrode level)
SETUP.SBNR = 0; % signal to biological noise power ratio expressed in dB (both measured on electrode level)
SETUP.SMNR = 0; % signal to measurment noise power ratio expressed in dB (both measured on electrode level)
SETUP.WhtNoiseAddFlg = 1; % white noise admixture in biological noise interference noise (FLAG)
SETUP.WhtNoiseAddSNR = 3; % SNR of BcgNoise and WhiNo (dB)
SETUP.SigPre = 0; SETUP.IntPre = 0; SETUP.BcgPre = 1; SETUP.MesPre = 1; % final signal components for pre-interval (use zero or one for signal, interference noise, biological noise, measurement noise)
SETUP.SigPst = 1; SETUP.IntPst = 1; SETUP.BcgPst = 1; SETUP.MesPst = 1; % final signal components for post-interval (as above)
% For localization, the default is to consider random locations of ROI with some fixed candidate points such that
% SETUP.SRCS(1,1) of them are in close locations; additionally, no interfering sources are considered
% in the localization model. You may comment out the 'if' section below to experiment with other settings.
if(SETUP.supSwitch == 'loc')
SETUP.rROI = logical(1); % random (1) or predefined (0) ROIs
SETUP.rPNT = logical(0); % random (1) or predefined (0) candidate points for source locations: if 0,
% number of sources as in SETUP.SRCS(1,1) will be fixed and in close locations
SETUP.SigPre = 0; SETUP.IntPre = 0; SETUP.BcgPre = 1; SETUP.MesPre = 1; % final signal components for pre-interval (use zero or one for signal, interference noise, biological noise, measurement noise)
SETUP.SigPst = 1; SETUP.IntPst = 0; SETUP.BcgPst = 1; SETUP.MesPst = 1; % final signal components for post-interval (as above)
end
if SETUP.rPNT
disp('CC: using random source locations')
else
disp('CC: using predefined source locations')
end
disp('CYBERCRAFT:: Generation of timeseries for bioelectrical activity of interest')
run('./tg_s01_Timeseries___01_SrcActiv.m')
disp('CYBERCRAFT:: Generation of timeseries for bioelectrical interference noise')
run('./tg_s01_Timeseries___02_IntNoise.m')
disp('CYBERCRAFT:: Generation of timeseries for bioelectrical background noise')
run('./tg_s01_Timeseries___03_BcgNoise.m')
disp('CYBERCRAFT:: Measurement noise generation')
run('./tg_s01_Timeseries___04_MesNoise.m')
run('./tg_s02_Geometry___01_RandSamp.m')
run('./tg_s02_Geometry___02_Indices.m')
run('./tg_s02_Geometry___03_Coordinates.m')
run('./tg_s02_Geometry___04_DeepSrc.m')
run('./tg_s02_Geometry___05_Perturbation.m')
run('./tg_s02_Geometry___06_lfg.m')
run('./tg_s03_Forward_modeling.m')
run('./tg_s04_Preparation___00_SNRs_Adjustment.m')
run('./tg_s04_Preparation___01_Signal_Covariances.m')
run('./tg_s05_Spatial_filtering___01_Preparation.m')
run('./tg_s05_Spatial_filtering___02_Definitions.m')
run('./tg_s05_Spatial_filtering___03_Execution.m')
run('./tg_s06_Error_evaluation.m')
run('./tg_s05_Localization___00_Preparations.m')
run('./tg_s05_Localization___01_MAI_and_MPZ_Locallizers.m')
Multiple runs for tangled matlab-scripts.
- Simulation settings can be modified in:
- {s00} Prelude
- Prelude (with optional settings overwrite), above
- the code block below
if exist('fPath'), cd(fPath); else, try, cd('~/supFunSim/'); catch, warningMessage = 'Problem encoutered while trying to change working directory to ''~/supFunSim/''.'; end; end; disp(['CYBERCRAFT:: pwd is: ',pwd])
clearvars ii jj kk nn tmp* LOOP* Loop*;
run('./tg_s07_BATCH___01_Prelude.m');
SETUP
% Get some details for simulation output filename (*.mat file)
LOOP.DATE = datestr(now,'yyyymmdd_HHMMSS');
LOOP.NAME = tempname; [~, LOOP.NAME] = fileparts(LOOP.NAME); % simulation unique name
% Number of simulation runs for each SNRs combination
LOOP.totSimCount = 1000;
LOOP.rngSimCount = 1:LOOP.totSimCount;
% Range of SNRs
LOOP.rngMesNoise = 10;
LOOP.rngBcgNoise = 0;
LOOP.rngIntNoise = [0:5:10];
LOOP
SETUP
run('./tg_s08_LOOP___01_Main_LOOP.m'); whos
switch SETUP.supSwitch
case {'rec'}
run('./tg_s08_Loop_Reconstruction_Plot___02_.m')
case {'loc'}
run('./tg_s08_Loop_Localization_Plot___02_.m')
otherwise
warning('Unexpected supSwitch value. No action taken, only data generated.')
endLOOP
LOOP.table_varN'
LOOP.table_rowN
size(LOOP.table_arrC)
% avg and std arrays #filters x #error_measures x #SXNR_levels
LOOP.table_arrC_avgSimCount = squeeze(mean(LOOP.table_arrC,3))
LOOP.table_arrC_stdSimCount = squeeze(std(LOOP.table_arrC,[],3))
size(LOOP.table_arrC_avgSimCount)
% SAVE THE RESULTS
mySave = ['./res___',LOOP.DATE,'___',LOOP.NAME,'___',datestr(now,'yyyymmdd_HHMMSS')];
save(mySave)
LOOP
close all
for ii = 1:length(LOOP.table_varN)
LOOP.table_varN'
tmp_pickVar = ii;
disp(LOOP.table_varN(tmp_pickVar))
TMP_res.table_arrC = squeeze(LOOP.table_arrC_avgSimCount(:,tmp_pickVar,:))
TMP_res.table_varN = strrep(strrep(strsplit(num2str(LOOP.(LOOP.NamNoise),'%+d ')),'-','n'),'+','p')
TMP_res.table_rowN = LOOP.table_rowN
TMP_res.table = array2table(TMP_res.table_arrC,'VariableNames',TMP_res.table_varN,'RowNames',TMP_res.table_rowN);
sortrows(TMP_res.table,1)
figure('Name',[LOOP.table_varN{tmp_pickVar},' for ',LOOP.NamNoise,' in ',num2str(LOOP.(LOOP.NamNoise),'%+d ')])
bar(TMP_res.table_arrC)
set(gca, 'XTick', 1:length(TMP_res.table_rowN), 'XTickLabel', strrep(TMP_res.table_rowN,'_','\_'));
h1 = legend(strsplit(num2str(LOOP.(LOOP.NamNoise),'%+d ')))
myTitle = strrep([LOOP.table_varN{tmp_pickVar},' for ',LOOP.NamNoise,' in ',num2str(LOOP.(LOOP.NamNoise),'%+d ')],'_','\_');
% title(myTitle)
myName = ['./fig/',LOOP.DATE,'___',LOOP.NAME,'___fig_',num2str(ii),'___',strrep(strrep(myTitle,'\',''),' ','_')];
v1 = get(h1,'title')
% fix me: it does not have to be SINR
set(v1,'string','SINR')
set(gcf, 'PaperUnits', 'centimeters')
set(gcf, 'PaperOrientation', 'portrait')
set(gcf, 'PaperSize', [20, 20])
set(gcf, 'Position', [100, 100, 800, 750])
legend('Location','northeastoutside')
xtickangle(45)
savefig(gcf,myName,'compact')
savefig(gcf,myName)
saveas(gcf,myName,'pdf')
end
LOOP
% USED FOR SAVING FIGURES
myName = ['./fig___',LOOP.DATE,'___',LOOP.NAME,'___',datestr(now,'yyyymmdd_HHMMSS')];
% use LOOP.RES.ITERATION{x}{y}(z) to access results struct for
% z-th activity index (z=1,2,...,8 with z=1 for MAI(2011) etc., as defined in !run MAI and MPZ Locallizers)
% y-th level of SNR of type considered (interference, biological, measurement)
% x-th iteration of the experiment
% XRANK is an array of size |x| by |y|, representing ranks selected for reduced-rank indices
XRANK=zeros(LOOP.totSimCount,length(LOOP.(LOOP.NamNoise)));
% number of simulation runs
for i=1:LOOP.totSimCount
% number of levels of SNR
for j=1:length(LOOP.(LOOP.NamNoise))
XRANK(i,j)=max([LOOP.RES.ITERATION{i}{j}(1).sources.rank_opt]);
end
end
% XRES is an array of size |x| by |y| by |z| by |#sources|, representing errors for each source
% found using iterative procedure
XRES=zeros(LOOP.totSimCount,length(LOOP.(LOOP.NamNoise)),8,sum(SETUP.SRCS(:,1)));
% number of simulation runs
for i=1:LOOP.totSimCount
% number of levels of SNR
for j=1:length(LOOP.(LOOP.NamNoise))
% number of activity indices
for k=1:8
XRES(i,j,k,:)=[LOOP.RES.ITERATION{i}{j}(k).sources.distance]';
end
end
end
% RRES is same like XRES, with the trailing few sources considered
RRES=XRES(:,:,:,end-1:end);
% |x| by |y| by |z| array, with last dimension representing mean of errors across sources discovered
XRES_s = mean(XRES,4);
RRES_s = mean(RRES,4);
% |y| by |z| matrix, with entries representing averaged mean error and its std across experiment repetitions
XRES_sa = squeeze(mean(XRES_s,1));
XRES_sstd = squeeze(std(XRES_s,1));
RRES_sa = squeeze(mean(RRES_s,1));
RRES_sstd = squeeze(std(RRES_s,1));
% change order of and remove some activity indices as:
% 1 = MAI(2011), 2 = MAI(2013), 3 = MAI_RR_I
% 4 = MPZ(2011), 5 = MPZ(2013), 6 = MPZ_RR_I
XRES_s = XRES_s(:,:,[1 2 5 3 4 7]);
RRES_s = RRES_s(:,:,[1 2 5 3 4 7]);
XRES_sa = XRES_sa(:,[1 2 5 3 4 7]);
XRES_sstd = XRES_sstd(:,[1 2 5 3 4 7]);
RRES_sa = RRES_sa(:,[1 2 5 3 4 7]);
RRES_sstd = RRES_sstd(:,[1 2 5 3 4 7]);
% Mann-Whitney U-test to check whether proposed activity indices achieve
% significantly lower average localization error
% XRES_U and RRES_U are |y| by 4 (# of comparisons) storing resulting p-values of performed Mann-Whitney U test
XRES_U=zeros(length(LOOP.(LOOP.NamNoise)),4);
RRES_U=zeros(length(LOOP.(LOOP.NamNoise)),4);
for j=1:length(LOOP.(LOOP.NamNoise))
XRES_U(j,1)=ranksum(XRES_s(:,j,1),XRES_s(:,j,3),'tail','right');
XRES_U(j,2)=ranksum(XRES_s(:,j,2),XRES_s(:,j,3),'tail','right');
XRES_U(j,3)=ranksum(XRES_s(:,j,4),XRES_s(:,j,6),'tail','right');
XRES_U(j,4)=ranksum(XRES_s(:,j,5),XRES_s(:,j,6),'tail','right');
end
for j=1:length(LOOP.(LOOP.NamNoise))
RRES_U(j,1)=ranksum(RRES_s(:,j,1),RRES_s(:,j,3),'tail','right');
RRES_U(j,2)=ranksum(RRES_s(:,j,2),RRES_s(:,j,3),'tail','right');
RRES_U(j,3)=ranksum(RRES_s(:,j,4),RRES_s(:,j,6),'tail','right');
RRES_U(j,4)=ranksum(RRES_s(:,j,5),RRES_s(:,j,6),'tail','right');
end
% pad with zero vectors to increase spacing between bar groups
XRES_sa = [XRES_sa; zeros(2,6)];
XRES_sstd = [XRES_sstd; zeros(2,6)];
RRES_sa = [RRES_sa; zeros(2,6)];
RRES_sstd = [RRES_sstd; zeros(2,6)];
clear gca;
figure
name = {'MAI';'MAI_{ext}';'MAI_{RR-I}';'MPZ';'MPZ_{ext}';'MPZ_{RR-I}'};
bar(XRES_sa',1.2);
hold on
bar(XRES_sstd',1.2,'FaceColor',[0 0 0])
hold off
f=get(gca,'Children');
f=flipud(f(3:end));
lgd=strsplit(num2str(LOOP.(LOOP.NamNoise),'%d '));
% need to skip items in legend
for i=1:length(lgd)
legendInfo{i} = ['SNR = ', lgd{i}, ' dB'];
end
l=legend(f,legendInfo);
l.Location='northwest';
set(gca,'XTickLabel',name);
ylabel('Average localization error [mm]')
% ax = findobj(gcf,'type','axes');
print(gcf,[fullfile(pwd,filesep), strcat(myName,'_XRES_')],'-depsc')
figure
bar(RRES_sa',1.2);
hold on
bar(RRES_sstd',1.2,'FaceColor',[0 0 0])
hold off
f=get(gca,'Children');
f=flipud(f(3:end));
l=legend(f,legendInfo);
l.Location='northwest';
set(gca,'XTickLabel',name);
ylabel('Average localization error [mm]')
print(gcf,[fullfile(pwd,filesep), strcat(myName,'_RRES_')],'-depsc')
clear gca;
figure
plot(XRANK,'-o','MarkerIndices',1:LOOP.totSimCount,'LineWidth',2);
f=get(gca,'Children');
f=flipud(f);
l=legend(f,legendInfo);
l.Location='northwest';
axis([1 LOOP.totSimCount 1 sum(SETUP.SRCS(:,1))])
xlabel('Simulation run')
xticks([1:LOOP.totSimCount])
ylabel('Rank selected')
yticks([1:sum(SETUP.SRCS(:,1))])
print(gcf,[fullfile(pwd,filesep), strcat(myName,'_XRANK_')],'-depsc')
[LOOP.c01RR,LOOP.c02MM,LOOP.c03BB,LOOP.c04II] = ndgrid(LOOP.rngSimCount,LOOP.rngMesNoise,LOOP.rngBcgNoise,LOOP.rngIntNoise);
LOOP.coords = [LOOP.c01RR(:),LOOP.c02MM(:),LOOP.c03BB(:),LOOP.c04II(:)];
2 == length(size(squeeze(double.empty([[length(LOOP.rngMesNoise),length(LOOP.rngBcgNoise),length(LOOP.rngIntNoise)],0]))))
if 2 == length(size(squeeze(double.empty([[length(LOOP.rngMesNoise),length(LOOP.rngBcgNoise),length(LOOP.rngIntNoise)],0]))))
if max([length(LOOP.rngMesNoise),length(LOOP.rngBcgNoise),length(LOOP.rngIntNoise)]) == min([length(LOOP.rngMesNoise),length(LOOP.rngBcgNoise),length(LOOP.rngIntNoise)]),
LOOP.IdxNoise = 0;
else,
[~,LOOP.IdxNoise] = max([length(LOOP.rngMesNoise),length(LOOP.rngBcgNoise),length(LOOP.rngIntNoise)]);
end;
if LOOP.IdxNoise == 1,
LOOP.NamNoise = 'rngMesNoise';
elseif LOOP.IdxNoise == 2,
LOOP.NamNoise = 'rngBcgNoise';
elseif LOOP.IdxNoise == 3,
LOOP.NamNoise = 'rngIntNoise';
else,
LOOP.NamNoise = 'AnyNoise';
end;
else,
error('CYBERCRAFT:: ensure that only one noise type is vector and the remaining are scalars');
end;
LOOP.progress = []; LOOP.table_arrC = []; LOOP.table_varN = {}; LOOP.table_rowN = {};
for LoopSimCount = 1:LOOP.totSimCount
LoopIterCount = [];
run('./tg_s07_BATCH___02_Signals.m');
run('./tg_s07_BATCH___03_Leadfields.m');
for LoopIntNoise = 1:length(LOOP.rngIntNoise)
for LoopBcgNoise = 1:length(LOOP.rngBcgNoise)
for LoopMesNoise = 1:length(LOOP.rngMesNoise)
% used by localization
LoopIterCount = [LoopIterCount;LoopMesNoise,LoopBcgNoise,LoopIntNoise];
% used by reconstruction
LOOP.progress = [LOOP.progress;LoopSimCount,LOOP.rngMesNoise(LoopMesNoise),LOOP.rngBcgNoise(LoopBcgNoise),LOOP.rngIntNoise(LoopIntNoise)];
fprintf('\n');
disp('============================')
disp(['CurrIter: ',num2str(size(LOOP.progress,1)),' of ',num2str(length(LOOP.coords))])
disp('----------------------------')
disp(['SimCount: ',num2str(LoopSimCount),' of ',num2str(LOOP.totSimCount)])
disp(['MesNoise: ',num2str(LoopMesNoise),' of ',num2str(length(LOOP.rngMesNoise)),' (',num2str(LOOP.rngMesNoise(LoopMesNoise)),' dB in ' ,'[ ',num2str(LOOP.rngMesNoise),' ])'])
disp(['BcgNoise: ',num2str(LoopBcgNoise),' of ',num2str(length(LOOP.rngBcgNoise)),' (',num2str(LOOP.rngBcgNoise(LoopBcgNoise)),' dB in ' ,'[ ',num2str(LOOP.rngBcgNoise),' ])'])
disp(['IntNoise: ',num2str(LoopIntNoise),' of ',num2str(length(LOOP.rngIntNoise)),' (',num2str(LOOP.rngIntNoise(LoopIntNoise)),' dB in ' ,'[ ',num2str(LOOP.rngIntNoise),' ])'])
disp('----------------------------')
disp(['SETUP.n00: ',num2str(SETUP.n00)])
disp(['SETUP.K00: ',num2str(SETUP.K00)])
% chkSim___tg_s01_PRE___001(SETUP)
disp('============================')
SETUP.SINR = LOOP.rngIntNoise(LoopIntNoise);
SETUP.SBNR = LOOP.rngBcgNoise(LoopBcgNoise);
SETUP.SMNR = LOOP.rngMesNoise(LoopMesNoise);
run('./tg_s07_BATCH___04_Preparations.m');
switch SETUP.supSwitch
case {'rec'}
% Reconstruction
disp('Reconstruction')
run('./tg_s07_BATCH___05_Filters.m');
% removes Original and Dummy signals from Figures
if(SETUP.SHOWori==0)
rec_res.table_rowN = rec_res.table_rowN(3:end);
rec_res.table_arrC = rec_res.table_arrC(3:end,:);
end
LOOP.table_rowN{size(LOOP.progress,1)} = rec_res.table_rowN;
LOOP.table_varN{size(LOOP.progress,1)} = rec_res.table_varN;
LOOP.table_arrC(:,:,LoopSimCount,LoopMesNoise,LoopBcgNoise,LoopIntNoise) = rec_res.table_arrC;
case {'loc'}
% Localization
disp('Localization')
run('tg_s07_BATCH___05_Localizers.m');
% use LOOP.RES.ITERATION{x}{y}(z) to access results struct for
% z-th activity index (z=1,2,...,8 with z=1 for MAI(2011) etc., as defined in !run MAI and MPZ Locallizers)
% y-th level of SNR of type considered (interference, biological, measurement)
% x-th iteration of the experiment
LOOP.RES.ITERATION{LoopSimCount}{size(LoopIterCount,1)} = RES;
% temporary helper code
% save(['RES__','ITERATION_',num2str(LoopSimCount),'__SINR_',num2str(SETUP.SINR),'__SBNR_',num2str(SETUP.SBNR),'__SMNR_',num2str(SETUP.SMNR)],'RES')
otherwise
warning('Unexpected supSwitch value. No action taken, only data generated.')
end
end
end
end
end
switch SETUP.supSwitch
case {'rec'}
if isequal(LOOP.table_rowN{:}), LOOP.table_rowN = LOOP.table_rowN{1}; else, error('CYBERCRAFT:: row names are inconsistent'); end;
if isequal(LOOP.table_varN{:}), LOOP.table_varN = LOOP.table_varN{1}; else, error('CYBERCRAFT:: var names are inconsistent'); end;
case {'loc'}
otherwise
warning('Unexpected supSwitch value. No action taken, only data generated.')
end
function y = rawTimeSeries3d2d(x)
% y = rawTimeSeries3d2d(x)
% NB.: this function asumes that trials are along third dimention.
y = reshape(permute(x,[1,3,2]),[],size(x,2),1);
end
Convert angles from radians to degrees.
function angDeg = rawRad2Deg(angRad)
% Convert angs from radians to degrees
angDeg = (180/pi) * angRad;
Check whether the MVAR model is stable.
function [H,varargout] = rawIsStableMVAR(A00,varargin)
if isempty(varargin)
STAB = 1;
disp(['Checking if eigenvalues are inside unit circle.']);
elseif length(varargin) == 1
STAB = varargin{1};
disp(['Checking if eigenvalues are inside circle with radius of ',num2str(STAB) ,'.']);
else
error('Too many input arguments');
end
S00 = size(A00,1); % dimension of state vectors
P00 = size(A00,2)/S00; % order of process
tmp_lambda = eig([A00; eye((P00-1)*S00) zeros((P00-1)*S00,S00)]);
H = ~any(abs(tmp_lambda)>STAB);
varargout{1} = max(abs(tmp_lambda));
end
Get size of an array (plain).
function [s] = rawSize(x,varargin)
if isempty(varargin)
s = size(x);
else
s = [];
for ii = 1:length(varargin{1})
s = [s,size(x,ii)];
end
end
end
Pretty table made from 2d array.
function [h1,h2] = rawImgSC(mat,varargin)
h1 = imagesc(mat);
[x,y] = ndgrid(1:size(mat,2),1:size(mat,1));
mat = mat';
if ~isempty(varargin)
fSize = varargin{1};
else
fSize = 18;
end
h2 = text(x(:),y(:),strtrim(cellstr(num2str(mat(:),'%.2f'))),'HorizontalAlignment','center','FontSize',fSize);
if 0
mat = [reshape(-12:12,5,5);reshape(-12:12,5,5);reshape(-12:12,5,5)]
figure(17);clf;rawImgSC(mat);colorbar;
end
end
Pretty color-map for MATLAB figures.
function [y] = rawHotColdColorMap(x)
% [y] = rawHotColdColorMap(x)
% Suggested rationalization: base/2
if 0
x = 255
y = rawHotColdColorMap(x)
end
y = [linspace(0,1,x)',linspace(0,1,x)',linspace(1,1,x)'];
y = [y;fliplr(flipud(y))];
% imagesc(permute(y,[1,3,2]))
end
Plot PDC values.
function rawPlotPDC(pdcData,accf)
tmp_row = size(pdcData,1);
tmp_col = size(pdcData,2);
tmp_cnt = 0;
for i = 1:tmp_row
for j = 1:tmp_col
tmp_cnt = tmp_cnt+1;
subplot(tmp_row,tmp_col,tmp_cnt)
area(accf,squeeze(pdcData(i,j,:)))
axis tight
ylim([0,1])
end
end
end
Adjust the SNR by changing the amplitude of the noise time-series.
function [y] = rawAdjTotSNRdB(x01,x02,newSNR)
y = ((x02 / rawNrm(x02)) * rawNrm(x01)) / (db2pow(0.5*newSNR));
end
Compute norm of the random vector.
function [y] = rawNrm(x)
% Compute norm of the random vector
y = sqrt(sum(rawMom(x,2),2));
end
Compute raw moment of the k-th order for the random vector.
function [M] = rawMom(x01,k)
% Compute raw moment of the k-th order for the random vector
M = mean(x01.^k, 1);
end
Compute raw moment of the k-th order for the random vector.
function [y] = rawSNRdB(x1,x2)
% Compute SNR in dB for two random vectors (signal and noise)
y = 20.*log10(norm(x1,'fro')./norm(x2,'fro'));
end
Prettify 3d MATLAB figures.
function [] = ccrender(varargin)
%%
% CCRENDER --- figure rendering settings, usage:
%
% ccrender( axisLimits, finish, origin, originLimits )
%
% - axisLimits [optional]:
% limits display to specified range defined using either:
% - vector of length 2 (all axes have same limits)
% [ xyzLim_1, xyzLim_2 ],
% - vector of length 6 (each axis has its own limits)
% [ xLim_1, xLim_2, yLim_1, yLim_2, zLim_1, zLim_2 ];
% - finish [optional]:
% changes rendering scheme;
% - origin [optional]:
% true/false;
% - originLim [optional]:
% a number.
%
% Questions and comments: [email protected]
%
% TODO:
% - add "try" envelopes for major blocks
%
%
%%
p = inputParser;
defAxisLimits = [];
% chkAxisLimits = @(x) validateattributes(x,{'numeric'},{'vector','numel',2});
chkAxisLimits = @(x) isnumeric(x) && isvector(x) && ( length(x) == 2 || length(x) == 6 );
addOptional(p,'axisLimits',defAxisLimits,chkAxisLimits)
defOriginLimits = [];
chkOriginLimits = @(x) isnumeric(x) && isvector(x) && ( length(x) == 2 || length(x) == 6 );
addOptional(p,'originLimits',defOriginLimits,chkOriginLimits)
defEqualize = false;
chkEqualize = @(x) validateattributes(x,{'logical'},{});
addParameter(p,'equalize',defEqualize,chkEqualize)
defOrigin = true;
chkOrigin = @(x) validateattributes(x,{'logical'},{});
addParameter(p,'origin',defOrigin,chkOrigin)
defLabels = true;
chkLabels = @(x) validateattributes(x,{'logical'},{});
addParameter(p,'labels',defLabels,chkLabels)
defFinish = 'matte';
vldFinish = {'glossy','matte'};
chkFinish = @(x) any(validatestring(x,vldFinish));
addParameter(p,'finish',defFinish,chkFinish)
s = dbstack;
p.KeepUnmatched = true;
parse(p,varargin{:})
tmp_gcf = gcf
ccdisp('INIT: ccrender')
ccdisp(['Updating figure: ',num2str(tmp_gcf.Number),'...'])
if ~isempty(p.UsingDefaults)
ccdisp([ s(1).name,': Using defaults:'])
disp(p.UsingDefaults')
end
if ~isempty(p.Results)
ccdisp([ s(1).name,': Results:'])
disp(p.Results')
end
if ~isempty(fieldnames(p.Unmatched))
ccdisp([ s(1).name,': Extra (unmatched) inputs:'])
disp(p.Unmatched')
end
%%
if length(p.Results.axisLimits) == 2
if p.Results.axisLimits(1) >= p.Results.axisLimits(2)
help ccrender;
error('CYBERCRAFT: axis limits [xyzLim_1,xyzLim_2] must be increasing.');
else
axis([ p.Results.axisLimits(1) p.Results.axisLimits(2) p.Results.axisLimits(1) p.Results.axisLimits(2) p.Results.axisLimits(1) p.Results.axisLimits(2) ]);
end
elseif length(p.Results.axisLimits) == 6
if p.Results.axisLimits(1) >= p.Results.axisLimits(2) || p.Results.axisLimits(3) >= p.Results.axisLimits(4) || p.Results.axisLimits(5) >= p.Results.axisLimits(6)
help ccrender;
error('CYBERCRAFT: axis limits [xLim_1,xLim_2,yLim_1,yLim_2,zLim_1,zLim_2] must be increasing for each axis.');
else
axis([ p.Results.axisLimits(1) p.Results.axisLimits(2) p.Results.axisLimits(3) p.Results.axisLimits(4) p.Results.axisLimits(5) p.Results.axisLimits(6) ]);
end
end
%{
xlim([-15,15])
ylim([-15,15])
zlim([-15,15])
axis square
axis equal
axis auto
%}
if p.Results.equalize == true
xMin = min(xlim);
xMax = max(xlim);
yMin = min(ylim);
yMax = max(ylim);
zMin = min(zlim);
zMax = max(zlim);
xyzMin = min( [ xMin yMin zMin ] );
xyzMax = max( [ xMax yMax zMax ] );
axis([ xyzMin xyzMax xyzMin xyzMax xyzMin xyzMax ]);
end
if p.Results.labels == true
xlabel('X-axis')
ylabel('Y-axis')
zlabel('Z-axis')
end
set(gca,'DataAspectRatio', [1 1 1])
set(gca,'PlotBoxAspectRatio',[1 1 1])
camproj('perspective'); % perspective | ortographic
%{
set(gca,'Projection','ortographic') % problem here?
%}
set(gca,'CameraViewAngle',10)
az = 60;
el = 30;
view(az, el);
if strcmp(p.Results.finish,'glossy')
lighting gouraud;
material shiny;
camlight;
%{
shading interp;
light;
lighting phong;
%}
end
axis on;
box on;
%%
if p.Results.origin
if length(p.Results.originLimits) == 2
if p.Results.originLimits(1) > 0
help ccrender;
error('CYBERCRAFT: origin axis low limit must be negative number.')
elseif p.Results.originLimits(2) < 0
help ccrender;
error('CYBERCRAFT: origin axis higher limit must be positive number.')
else
xMin = p.Results.originLimits(1);
xMax = p.Results.originLimits(2);
yMin = p.Results.originLimits(1);
yMax = p.Results.originLimits(2);
zMin = p.Results.originLimits(1);
zMax = p.Results.originLimits(2);
end
elseif length(p.Results.originLimits) == 6
if p.Results.originLimits(1) > 0 || p.Results.originLimits(3) > 0 || p.Results.originLimits(5) > 0
help ccrender;
error('CYBERCRAFT: origin axis low limit must be negative number.')
elseif p.Results.originLimits(2) < 0 || p.Results.originLimits(4) < 0 || p.Results.originLimits(6) < 0
help ccrender;
error('CYBERCRAFT: origin axis higher limit must be positive number.')
else
xMin = p.Results.originLimits(1);
xMax = p.Results.originLimits(2);
yMin = p.Results.originLimits(3);
yMax = p.Results.originLimits(4);
zMin = p.Results.originLimits(5);
zMax = p.Results.originLimits(6);
end
else
xMin = min(xlim);
xMax = max(xlim);
yMin = min(ylim);
yMax = max(ylim);
zMin = min(zlim);
zMax = max(zlim);
end
hold on
try
ccfgFigH = evalin( 'base', 'ccfgFigH' );
catch
ccfgFigH = [];
end
try
tmp_fields = fieldnames(ccfgFigH(tmp_gcf.Number).origin);
for ii = 1:length(tmp_fields)
try
delete(ccfgFigH(tmp_gcf.Number).origin.(tmp_fields{ii}));
catch tmp_catched_2
fprintf('\n\nPSEUDO-WARNING[2]: %s\n\n',tmp_catched_2.message);
end
end
catch tmp_catched_1
fprintf('\n\nPSEUDO-WARNING [1]: %s\n\n',tmp_catched_1.message);
end
ccfgFigH(tmp_gcf.Number).origin.xp = quiver3( 0, 0, 0, 1.1*xMax, 0, 0, 'color', [1 0 0], 'LineWidth', 1.0, 'LineStyle', '-' ); hold on;
ccfgFigH(tmp_gcf.Number).origin.yp = quiver3( 0, 0, 0, 0, 1.1*yMax, 0, 'color', [0 1 0], 'LineWidth', 1.0, 'LineStyle', '-' ); hold on;
ccfgFigH(tmp_gcf.Number).origin.zp = quiver3( 0, 0, 0, 0, 0, 1.1*zMax, 'color', [0 0 1], 'LineWidth', 1.0, 'LineStyle', '-' ); hold on;
ccfgFigH(tmp_gcf.Number).origin.xn = quiver3( 0, 0, 0, xMin, 0, 0, 'color', [1 0 0], 'LineWidth', 1.0, 'LineStyle', '-.', 'ShowArrowHead', 'off' ); hold on;
ccfgFigH(tmp_gcf.Number).origin.yn = quiver3( 0, 0, 0, 0, yMin, 0, 'color', [0 1 0], 'LineWidth', 1.0, 'LineStyle', '-.', 'ShowArrowHead', 'off' ); hold on;
ccfgFigH(tmp_gcf.Number).origin.zn = quiver3( 0, 0, 0, 0, 0, zMin, 'color', [0 0 1], 'LineWidth', 1.0, 'LineStyle', '-.', 'ShowArrowHead', 'off' ); hold on;
ccfgFigH(tmp_gcf.Number).origin.xt = text( 1.2*xMax, 0, 0, 'OX', 'color', [1 0 0] );
ccfgFigH(tmp_gcf.Number).origin.yt = text( 0, 1.2*yMax, 0, 'OY', 'color', [0 1 0] );
ccfgFigH(tmp_gcf.Number).origin.zt = text( 0, 0, 1.2*zMax, 'OZ', 'color', [0 0 1] );
assignin( 'base', 'ccfgFigH', ccfgFigH );
end
%%
xLimits = xlim;
yLimits = ylim;
zLimits = zlim;
ccdisp( [ 'xLimits: [', num2str(xLimits(1)), ', ', num2str(xLimits(2)), ']' ] )
ccdisp( [ 'yLimits: [', num2str(yLimits(1)), ', ', num2str(yLimits(2)), ']' ] )
ccdisp( [ 'zLimits: [', num2str(zLimits(1)), ', ', num2str(zLimits(2)), ']' ] )
grid on;
rotate3d on;
hold on;
ccdisp(['Figure: ',num2str(tmp_gcf.Number),' has been updated!'])
ccdisp('DONE: ccrender')
return
Just throw some text to MATLAB console.
function [ ] = ccdisp( str )
% CCDISP display as CYBERCRAFT
%
% Use
%
% ccdisp( str )
%
% add here "try + catch" OR "ifexists"
disp( ['CYBERCRAFT: ', str ] );
return
Convert cartesian to spherical coordinates (row-wise).
function [y] = rawCartToSph(x)
% Convert cartesian to spherical coordinates (row-wise)
if ~isempty(x)
for ii = 1:size(x,1)
[y(ii,1),y(ii,2),y(ii,3)] = cart2sph(x(ii,1),x(ii,2),x(ii,3));
end
else
y = x;
end
end
Convert spherical to cartesian coordinates (row-wise).
function [y] = rawSphToCart(x)
% Convert spherical to cartesian coordinates (row-wise)
if ~isempty(x)
for ii = 1:size(x,1)
[y(ii,1),y(ii,2),y(ii,3)] = sph2cart(x(ii,1),x(ii,2),x(ii,3));
end
else
y=x;
end
end
Some EEGLab plugins (e.g., cleanline) contain function strjoin
which overwrites the original MATLAB function. The following
function overshadows the any external strjoin functions with the
original one.
function [] = rawFixStrJoin()
disp('=== Before ===')
which('strjoin','-all')
addpath([matlabroot,'/toolbox/matlab/strfun/'])
disp('=== After ===')
which('strjoin','-all')
end
rm ./tg_*
rm ./fun/*
- Principal sub-matrices in MVAR model are reconstructed quite robustly from IntNoise signal (that may have lowered dimensionality) when compared with SigIn
- The above also aplies to PDC
(setq org-image-actual-width '(800))
(setq org-todo-keywords '(
(sequence "TODO(t)" "|" "DONE(d)" "SRTD" "NRML")
(sequence "NEW!" "FIX!" "WTF!" "DUE!" "CHCK" "REDO" "TOSS" "TEST" "OPEN" "DOIN" "|" "CNCL" "NVRM")
(sequence "SOON" "LATE" "NEXT" "|" "OPTN" "WAIT" "HOLD")
(sequence "!run" "!std" "!opt" "!fun" "!fig" "|")
(sequence "CYCL" "|")
))
As per:
- cite:g2010—Hui-et-al—Identifying-true-cortical-interactions-in-MEG-using-the-nulling-beamformer
/--------------------------\ /--------------------------\
| sel_atl.mat | | ROIs |
|--------------------------| |--------------------------|
| sel_ele.mat * | *
|--------------------------| |--------------------------|
| sel_ele.mat | | |
|--------------------------| |--------------------------|
| sel_geo_deep_thalami.mat | | |
|--------------------------| |--------------------------|
| sel_mri00.mat | | |
|--------------------------| |--------------------------|
| sel_msh.mat | | |
|--------------------------| |--------------------------|
| sel_vol.mat | | |
\--------------------------/ \--------------------------/
/----------------------------\
| |
| Brain triangulation |
| |
+----------------------------+ /----------------------------\
| | | |
| Skull triangulation *------------------>* Volume conduction model |
| | | |
+----------------------------+ +----------*-----------------+
| | |
| Scalp Triangulation | |
| | |
\----------------------------/ |
|
/----------------------------------------\ |
| | |
| Cortex parcellation | |
| | |
\----------------*-----------------------/ |
| |
| |
V V
/----------------------------------------\ /----------*--------------------\
| | | |
| ROIs random sampling *------>* Leadfield computation |
| | | |
\-----------*----------------------------/ \-------------------------------/
| |
| /------------------------\ |
| | | |
+-->* Volume cond. model | |
| | V
+------------------------+ /----------*--------------------\
| | | |
| Grid (ROI based) *<------* Volume lookup (ROI) |
| | | |
+------------------------+ \-------------------------------/
| |
| Leadfields |
| |
\------------------------/