Skip to content

Commit f6f7d1d

Browse files
tiborauertiborauer
committed
NEW - MEEG: fEI
1 parent 9318912 commit f6f7d1d

5 files changed

Lines changed: 228 additions & 1 deletion

File tree

‎aa_modules/aamod_meeg_fei.m‎

Lines changed: 62 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,62 @@
1+
function [aap, resp] = aamod_meeg_fei(aap,task,subj,sess)
2+
3+
resp='';
4+
5+
switch task
6+
case 'report'
7+
% for fn = cellstr(spm_select('FPList',aas_getsesspath(aap,subj,sess),'^diagnostic_.*jpg$'))'
8+
% aap = aas_report_add(aap,subj,'<table><tr><td>');
9+
% aap=aas_report_addimage(aap,subj,fn{1});
10+
% aap = aas_report_add(aap,subj,'</td></tr></table>');
11+
% end
12+
case 'doit'
13+
infname = cellstr(aas_getfiles_bystream(aap,'meeg_session',[subj sess],'meeg'));
14+
15+
[~, EL] = aas_cache_get(aap,'eeglab');
16+
EL.load;
17+
18+
indfnEEG = strcmp(spm_file(infname,'ext'),'set');
19+
EEG = pop_loadset('filepath',spm_file(infname{indfnEEG},'path'),'filename',spm_file(infname{indfnEEG},'filename'));
20+
EL.unload;
21+
22+
data = cell2mat(arrayfun(@(ch) abs(hilbert(EEG.data(ch,:)))', 1:EEG.nbchan, 'UniformOutput',false));
23+
24+
%% DFA
25+
cfg = aas_getsetting(aap,'dfa');
26+
if ~isempty(cfg.windowsize) && ~isempty(cfg.windowoverlap)
27+
expDFA = calculateDFA(data, cfg.windowsize*EEG.srate, cfg.windowoverlap);
28+
else
29+
expDFA = NaN(EEG.nbchan,1);
30+
end
31+
32+
%% fEI
33+
cfg = aas_getsetting(aap,'fei');
34+
fEI = calculateFEI(data, cfg.windowsize*EEG.srate, cfg.windowoverlap);
35+
36+
%% Save
37+
[~, FT] = aas_cache_get(aap,'fieldtrip');
38+
FT.load;
39+
40+
data = keepfields(ft_read_header(infname{indfnEEG}),{'elec','label'});
41+
data.dimord = 'chan';
42+
43+
FT.unload;
44+
45+
% dfa
46+
if ~all(isnan(expDFA))
47+
dfa = data;
48+
dfa.expdfa = expDFA;
49+
save(spm_file(infname{indfnEEG},'prefix','dfa_'),'dfa');
50+
aas_desc_outputs(aap,'meeg_session',[subj sess],'dfa',spm_file(infname{indfnEEG},'prefix','dfa_'));
51+
end
52+
53+
% fei
54+
fei = data;
55+
fei.fei = fEI;
56+
save(spm_file(infname{indfnEEG},'prefix','fei_'),'fei');
57+
aas_desc_outputs(aap,'meeg_session',[subj sess],'fei',spm_file(infname{indfnEEG},'prefix','fei_'));
58+
59+
case 'checkrequirements'
60+
if ~aas_cache_get(aap,'eeglab'), aas_log(aap,true,'EEGLAB is not found'); end
61+
if ~aas_cache_get(aap,'fieldtrip'), aas_log(aap,true,'FieldTrip is not found'); end
62+
end

‎aa_modules/aamod_meeg_fei.xml‎

Lines changed: 31 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,31 @@
1+
<?xml version="1.0" encoding="utf-8"?>
2+
<aap>
3+
<tasklist>
4+
<currenttask domain='meeg_session' desc='Measurement of excitation-inhibition ratio' modality='MEEG'>
5+
6+
<qsub>
7+
<timeBase>0.15</timeBase>
8+
<memoryBase>0.5</memoryBase>
9+
</qsub>
10+
11+
<dfa desc='DFA calculation parameters. If any if not specified, DFA is not calculated'>
12+
<windowsize desc='Array of calculation window size in seconds, [1xN]'></windowsize>
13+
<windowoverlap desc='Fraction of overlap between calculation windows'></windowoverlap>
14+
</dfa>
15+
<fei>
16+
<windowsize desc='Size of calculation window in seconds'></windowsize>
17+
<windowoverlap desc='Fraction of overlap between calculation windows'></windowoverlap>
18+
</fei>
19+
20+
<inputstreams>
21+
<stream ismodified='0'>meeg</stream>
22+
</inputstreams>
23+
24+
<outputstreams>
25+
<stream>dfa</stream>
26+
<stream>fei</stream>
27+
</outputstreams>
28+
29+
</currenttask>
30+
</tasklist>
31+
</aap>

‎external/README.md‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,8 @@ Automatic Analysis relies on several externally-developed packages and functions
1919
## Functions
2020
| Name | Version | Last Updated | URL | Maintainer |
2121
| ------------------------ | -------------------------------------------------- | ----------------- | ---------------------------------------------------------------------------------------------------------------------------------------------------------- | ----------------------------- |
22+
| calculateDFA.m | | 16 April 2026 | https://doi.org/10.6084/m9.figshare.12284099 | Tibor Auer |
23+
| calculateFEI.m | | 16 April 2026 | https://github.com/rhardstone/fEI | Tibor Auer |
2224
| catstruct.m | (mod) | 22 July 2015 | http://www.mathworks.com/matlabcentral/fileexchange/7842-catstruct | Tibor Auer |
2325
| colorspace.m | | 01 August 2015 | http://www.mathworks.com/matlabcentral/fileexchange/28790-colorspace-transformations | Tibor Auer |
2426
| cor2mni.m | | 09 Jun 2016 | http://www.alivelearn.net/?m=201006 | Tibor Auer |
@@ -28,4 +30,4 @@ Automatic Analysis relies on several externally-developed packages and functions
2830
| imagescnan.m | 2.1 | 16 September 2015 | http://www.mathworks.com/matlabcentral/fileexchange/20516-imagescnan-m-v2-1--aug-2009- | Tibor Auer |
2931
| resize_img.m | | 21 March 2017 | http://www0.cs.ucl.ac.uk/staff/gridgway/vbm | Tibor Auer |
3032
| rotateticklabel.m | (mod) | | http://www.mathworks.com/matlabcentral/fileexchange/29800-scenes-objects-classification-toolbox | Tibor Auer |
31-
| sscanfitem.m | | 22 July 2015 | https://surfer.nmr.mgh.harvard.edu/fswiki/DownloadAndInstall | Tibor Auer |
33+
| sscanfitem.m | | 22 July 2015 | https://surfer.nmr.mgh.harvard.edu/fswiki/DownloadAndInstall | Tibor Auer |

‎external/calculateDFA.m‎

Lines changed: 70 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,70 @@
1+
function [DFAExponent, meanDF, windowSizes] = calculateDFA(Signal, windowSizes, windowOverlap)
2+
% Originally created by Richard Hardstone (2020), rhardstone@gmail.com
3+
4+
% This program is free software: you can redistribute it and/or modify
5+
% it under the terms of the GNU General Public License as published by
6+
% the Free Software Foundation, either version 3 of the License, or
7+
% (at your option) any later version.
8+
%
9+
% This program is distributed in the hope that it will be useful,
10+
% but WITHOUT ANY WARRANTY; without even the implied warranty of
11+
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12+
% GNU General Public License for more details.
13+
%
14+
% You should have received a copy of the GNU General Public License
15+
% along with this program. If not, see <http://www.gnu.org/licenses/>.
16+
17+
% Calculates DFA (on a set of window sizes) for signal
18+
% Steps refer to description of DFA algorithm in Figure 1 of:
19+
% Detrended fluctuation analysis: a scale-free view on neuronal oscillations
20+
% Hardstone, Richard, Simon-Shlomo Poil, Giuseppina Schiavone, Rick Jansen, Vadim V. Nikulin, Huibert D. Mansvelder, and Klaus Linkenkaer-Hansen
21+
% Frontiers in physiology 3 (2012): 450.
22+
23+
24+
% Signal dimensions (numSamples,numChannels)
25+
% windowSize in samples
26+
% windowOverlap is fraction of overlap between windows (0-1)
27+
% Signal dimensions (numSamples,numChannels)
28+
% windowSize in samples
29+
% windowOverlap in samples
30+
31+
32+
lengthSignal = size(Signal,1);
33+
numChannels = size(Signal,2);
34+
meanDF = zeros(numChannels,length(windowSizes));
35+
DFAExponent = zeros(numChannels,1);
36+
windowSizes = windowSizes(:);
37+
38+
for i_channel = 1:size(Signal,2)
39+
for i_windowSize = 1:length(windowSizes)
40+
windowOffset = floor(windowSizes(i_windowSize) * (1-windowOverlap));
41+
allWindowIndex = createWindowIndices(lengthSignal, windowSizes(i_windowSize), windowOffset);
42+
originalAmplitude = Signal(:,i_channel);
43+
signalProfile = cumsum(originalAmplitude - mean(originalAmplitude)); %Step A->B
44+
xSignal = signalProfile(allWindowIndex); %Step B->C
45+
dSignal = detrend(xSignal'); %Step C
46+
w_detrendedFluctuations = std(dSignal,1); %std(dSignal,1) means normalized by N instead of N-1
47+
meanDF(i_channel,i_windowSize) = mean(w_detrendedFluctuations);
48+
end
49+
50+
X = [ones(length(windowSizes),1) log10(windowSizes)];
51+
Y = log10(meanDF(i_channel,:))';
52+
if length(windowSizes) > 1
53+
regressOutput = regress(Y,X);
54+
DFAExponent(i_channel) = regressOutput(2,1); %Step D
55+
else
56+
DFAExponent(i_channel) = nan;
57+
end
58+
end
59+
60+
function allWindowIndex = createWindowIndices(lengthSignal, lengthWindow, windowOffset)
61+
%Gets indices for a set of windows of size (lengthWindow) with a set overlap between them
62+
63+
windowStarts = (1:windowOffset:lengthSignal-lengthWindow+1)-1;
64+
numWindows = length(windowStarts);
65+
66+
oneWindowIndex = 1:lengthWindow;
67+
allWindowIndex = repmat(oneWindowIndex,[numWindows,1]);
68+
69+
allWindowIndex = allWindowIndex + repmat(windowStarts',[1 lengthWindow]);
70+

‎external/calculateFEI.m‎

Lines changed: 62 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,62 @@
1+
function [EI, wAmp, wDNF] = calculateFEI(Signal, windowSize, windowOverlap)
2+
% Originally created by Richard Hardstone (2020), rhardstone@gmail.com
3+
% Please note that commercial use of this algorithm is protected by Patent claim (PCT/NL2019/050167) .Method of determining brain activity.; with priority date 16 March 2018
4+
% This code is licensed under creative commons license CC-BY-NC https://creativecommons.org/licenses/by-nc/4.0/legalcode.txt
5+
6+
%%
7+
% Calculates fEI (on a set window size) for signal
8+
% Steps refer to description of fEI algorithm in Figure 2D of paper:
9+
% Measurement of excitation inhibition ratio in autism spectrum disorder using critical brain dynamics
10+
% Scientific Reports (2020)
11+
% Hilgo Bruining*, Richard Hardstone*, Erika L. Juarez-Martinez*, Jan Sprengers*, Arthur-Ervin Avramiea, Sonja Simpraga, Simon J. Houtman, Simon-Shlomo Poil5,
12+
% Eva Dallares, Satu Palva, Bob Oranje, J. Matias Palva, Huibert D. Mansvelder & Klaus Linkenkaer-Hansen
13+
% (*Joint First Author)
14+
%
15+
%
16+
17+
%%Input
18+
% Signal: amplitude envelope with dimensions (numSamples,numChannels)
19+
% windowSize: in samples
20+
% windowOverlap: is fraction of overlap between windows (0-1)
21+
22+
%% Example:
23+
% [EI, wAmp, wDNF] = calculateFEI(randn(100000,1), 5000, 0.8)
24+
% Will calculate fEI on the single channel white noise signal with length 100000
25+
% using a window size of 5000 samples, and 80% overlap between each window
26+
27+
lengthSignal = size(Signal,1);
28+
numChannels = size(Signal,2);
29+
30+
windowOffset = floor(windowSize * (1-windowOverlap));
31+
allWindowIndex = createWindowIndices(lengthSignal, windowSize, windowOffset);
32+
numWindows = size(allWindowIndex,1);
33+
34+
EI = zeros(numChannels,1);
35+
wAmp = zeros(numChannels,numWindows);
36+
wDNF = zeros(numChannels,numWindows);
37+
38+
for i_channel = 1:numChannels
39+
originalAmplitude = Signal(:,i_channel);
40+
signalProfile = cumsum(originalAmplitude - mean(originalAmplitude)); % Step ii -> iii
41+
w_originalAmplitude = mean(originalAmplitude(allWindowIndex),2); % Calculate mean amplitude for each window
42+
xAmp = repmat(w_originalAmplitude,[1 windowSize]);
43+
xSignal = signalProfile(allWindowIndex); % Arrange Signals into windows
44+
xSignal = (xSignal ./ xAmp)'; % Step iii -> iv
45+
dSignal = detrend(xSignal); % Step iv -> v
46+
w_detrendedNormalizedFluctuations = std(dSignal,1)'; % Step v -> vi
47+
EI(i_channel) = 1-corr(w_detrendedNormalizedFluctuations,w_originalAmplitude); % Step vi -> vii
48+
wAmp(i_channel,:) = w_originalAmplitude;
49+
wDNF(i_channel,:) = w_detrendedNormalizedFluctuations;
50+
end
51+
52+
53+
function allWindowIndex = createWindowIndices(lengthSignal, lengthWindow, windowOffset)
54+
%Gets indices for a set of windows of size (lengthWindow) with a set overlap between them
55+
56+
windowStarts = (1:windowOffset:lengthSignal-lengthWindow+1)-1;
57+
numWindows = length(windowStarts);
58+
59+
oneWindowIndex = 1:lengthWindow;
60+
allWindowIndex = repmat(oneWindowIndex,[numWindows,1]);
61+
62+
allWindowIndex = allWindowIndex + repmat(windowStarts',[1 lengthWindow]);

0 commit comments

Comments
 (0)