Attachment 'process_beamformer_lcmv.m'
Download 1 function varargout = process_beamformer_lcmv( varargin )
2 % PROCESS_LCMV_BEAMFORMER: Linearly-constrained minimum variance beamformer
3
4 % @=============================================================================
5 % This software is part of the Brainstorm software:
6 % http://neuroimage.usc.edu/brainstorm
7 %
8 % Copyright (c)2000-2015 University of Southern California & McGill University
9 % This software is distributed under the terms of the GNU General Public License
10 % as published by the Free Software Foundation. Further details on the GPL
11 % license can be found at http://www.gnu.org/copyleft/gpl.html.
12 %
13 % FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
14 % UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
15 % WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
16 % MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
17 % LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
18 %
19 % For more information type "brainstorm license" at command prompt.
20 % =============================================================================@
21 %
22 % Authors: Hui-Ling Chan, Francois Tadel, Sylvain Baillet, 2014
23
24 macro_methodcall;
25 end
26
27
28 %% ===== GET DESCRIPTION =====
29 function sProcess = GetDescription() %#ok<DEFNU>
30 % Description the process
31 sProcess.Comment = 'Beamformer: LCMV (20141022) [Experimental]';
32 sProcess.FileTag = '';
33 sProcess.Category = 'Custom';
34 sProcess.SubGroup = 'Sources';
35 sProcess.Index = 340;
36 % Definition of the input accepted by this process
37 sProcess.InputTypes = {'data'};
38 sProcess.OutputTypes = {'results'};
39 sProcess.nInputs = 1;
40 sProcess.nMinFiles = 1;
41 % Definition of the options
42 sProcess.options.comment.Comment = 'Comment: ';
43 sProcess.options.comment.Type = 'text';
44 sProcess.options.comment.Value = '';
45 % === OPTIONS: Cortial constrained orientation (eig -> singular value, corrected in 20141009)
46 sProcess.options.label_orient.Comment = '<HTML><BR> Sources orientation:';
47 sProcess.options.label_orient.Type = 'label';
48 sProcess.options.oriconstraint.Comment = {'Unconstrained (max singular value)','Unconstrained (trace)', 'Constrained (normal to cortex)'};
49 sProcess.options.oriconstraint.Type = 'radio';
50 sProcess.options.oriconstraint.Value = 1;
51 % === ACTIVE TIME RANGE
52 sProcess.options.label_range.Comment = '<HTML><BR> ';
53 sProcess.options.label_range.Type = 'label';
54 sProcess.options.nai_range.Comment = 'Time range of interest: ';
55 sProcess.options.nai_range.Type = 'poststim';
56 sProcess.options.nai_range.Value = [];
57 % === ACTIVE TIME WINDOW SIZE
58 sProcess.options.active_window_size.Comment = 'Active window size: ';
59 sProcess.options.active_window_size.Type = 'value';
60 sProcess.options.active_window_size.Value = {0, 'ms', 1};
61 % === ACTIVE TEMPORAL RESOLUTION
62 sProcess.options.nai_tresolution.Comment = 'Temporal resolution: ';
63 sProcess.options.nai_tresolution.Type = 'value';
64 sProcess.options.nai_tresolution.Value = {0, 'ms', 1};
65 % === REGULARIZATION
66 sProcess.options.reg.Comment = 'Regularization parameter: ';
67 sProcess.options.reg.Type = 'value';
68 sProcess.options.reg.Value = {1, '%', 4}; % enlarge the default value (20141009)
69 % === SENSOR TYPES
70 sProcess.options.label_sensor.Comment = '<HTML><BR> ';
71 sProcess.options.label_sensor.Type = 'label';
72 sProcess.options.sensortypes.Comment = 'Sensor types or names (empty=all): ';
73 sProcess.options.sensortypes.Type = 'text';
74 sProcess.options.sensortypes.Value = 'MEG';
75 % === SAVE SPATIAL FILTER
76 sProcess.options.savefilter.Comment = 'Save spatial filters';
77 sProcess.options.savefilter.Type = 'checkbox';
78 sProcess.options.savefilter.Value = 0;
79 end
80
81
82 %% ===== FORMAT COMMENT =====
83 function Comment = FormatComment(sProcess) %#ok<DEFNU>
84 Comment = sProcess.Comment;
85 end
86
87
88 %% ===== RUN =====
89 function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
90 % Initialize returned list of files
91 OutputFiles = {};
92 % Get option values
93 NAIRange = sProcess.options.nai_range.Value{1};
94 Reg = sProcess.options.reg.Value{1};
95 SensorTypes = sProcess.options.sensortypes.Value;
96 ActiveWindowSize = sProcess.options.active_window_size.Value{1};
97 NAITResolu = sProcess.options.nai_tresolution.Value{1};
98 Comment = sProcess.options.comment.Value;
99 OrientConstr = sProcess.options.oriconstraint.Value;
100 isSaveFilter = sProcess.options.savefilter.Value;
101
102 % profile on
103
104 % ===== LOAD THE DATA =====
105 bst_progress('start', 'Beamformer: LCMV', 'Loading input data...');
106 % Read the first file in the list, to initialize the loop
107 DataMat = in_bst_data(sInputs(1).FileName);
108 nChannels = size(DataMat.F,1);
109 nTime = size(DataMat.F,2);
110 Time = DataMat.Time;
111 % Lock input time window on the real time vector
112 iNAIRange = panel_time('GetTimeIndices', Time, NAIRange);
113 NAIRange = [Time(iNAIRange(1)), Time(iNAIRange(end))];
114 % Convert active window to a number of samples
115 sFreq = 1 ./ (Time(2) - Time(1));
116 ActiveWindowSmp = ActiveWindowSize / sFreq;
117
118 % ===== CHECK THE TIME WINDOWS =====
119 % Verify the coherence of the inputs
120 if (NAIRange(1) > NAIRange(2))
121 bst_report('Error', sProcess, sInputs, 'The setting of time range of interest is incorrect.');
122 return;
123 end
124 if (NAIRange(1) < Time(1))
125 bst_report('Warning', sProcess, sInputs, 'The start for time range of interest is reset to the first time point of data');
126 NAIRange(1) = Time(1);
127 end
128 if (NAIRange(2) > Time(end))
129 bst_report('Warning', sProcess, sInputs, 'The end for time range of interest is reset to the end point of data');
130 NAIRange(2) = Time(end);
131 end
132 if ((NAIRange(1)+ActiveWindowSize) > NAIRange(2))
133 bst_report('Warning', sProcess, sInputs, 'The active window size is reset to the same as the time range of interest.');
134 ActiveWindowSize = NAIRange(2) - NAIRange(1);
135 end
136 % Temporal resolution=0 : Active window size set to all the resolution of interest
137 if (NAITResolu == 0) || (ActiveWindowSize == 0)
138 NAITResolu = 1;
139 ActiveWindowSize = NAIRange(2) - NAIRange(1);
140 end
141
142 % Get time points of interest
143 NAIPoints = panel_time('GetTimeIndices', Time, [NAIRange(1)+ActiveWindowSize, NAIRange(2)]);
144 if (length(NAIPoints) <= 1)
145 nNAI = 1;
146 else
147 nNAI = length((NAIRange(1)+ActiveWindowSize):NAITResolu:NAIRange(2));
148 end
149
150 NAITimeList= zeros(nNAI,2);
151 for i=1:nNAI
152 NAITimeList(i,:) = NAIRange(1) + NAITResolu*(i-1) + [0 ActiveWindowSize] ;
153 end
154
155 MinVarTime = [NAIRange(1) NAIRange(2)];
156
157 nActiveWindowPoints = length(panel_time('GetTimeIndices', Time, NAITimeList(1,:)));
158 iActiveWindowTime = zeros(nNAI, nActiveWindowPoints);
159 for i = 1:nNAI
160 if (NAITimeList(i,1) < Time(1) || NAITimeList(i,2) > Time(end))
161 bst_report('Error', sProcess, sInputs, 'One active time window is not within the data time range.');
162 return;
163 end
164 single_iActiveWindow = panel_time('GetTimeIndices', Time, NAITimeList(i,:));
165
166 % detect the round error caused by the non-interger sampling rate
167 % -- added in 20141008
168 if nActiveWindowPoints - length(single_iActiveWindow) == 1
169 single_iActiveWindow(nActiveWindowPoints) = single_iActiveWindow(end) + 1;
170 end
171 %%%%%%%%
172 iActiveWindowTime(i,:) = single_iActiveWindow(1:nActiveWindowPoints);
173 end
174 iMinVarTime = panel_time('GetTimeIndices', Time, MinVarTime);
175
176
177 % ===== LOAD CHANNEL FILE =====
178 % Load channel file
179 ChannelMat = in_bst_channel(sInputs(1).ChannelFile);
180 % Get the channels of interest
181 iChannels = channel_find(ChannelMat.Channel, SensorTypes);
182
183 % ===== LOAD HEAD MODEL =====
184 % Get channel study
185 [sChannelStudy, iChannelStudy] = bst_get('ChannelFile', sInputs(1).ChannelFile);
186 % Load the default head model
187 HeadModelFile = sChannelStudy.HeadModel(sChannelStudy.iHeadModel).FileName;
188 sHeadModel = in_headmodel_bst(HeadModelFile);
189 if (OrientConstr == 3) && (isequal(sHeadModel.HeadModelType,'volume') || isempty(sHeadModel.GridOrient))
190 bst_report('Error', sProcess, sInputs, 'No dipole orientation for cortical constrained beamformer estimation.');
191 return;
192 end
193
194
195 % ===== LOAD ALL INPUT FILES =====
196 % Initialize the covariance matrices
197 MinVarCov = zeros(nChannels, nChannels);
198 nTotalMinVar = zeros(nChannels, nChannels);
199 NAIActiveCov = zeros(nChannels, nChannels, nNAI);
200 nTotalNAIActive = zeros(nChannels, nChannels, nNAI);
201 bst_progress('start', 'Beamformer: LCMV', 'Calulating covariance matrix...', 0, length(sInputs));
202 % Reading all the input files in a big matrix
203 for i = 1:length(sInputs)
204 % Read the file #i
205 DataMat = in_bst_data(sInputs(i).FileName);
206 % Check the dimensions of the recordings matrix in this file
207 if (size(DataMat.F,1) ~= nChannels) || (size(DataMat.F,2) ~= nTime)
208 bst_report('Error', sProcess, sInputs, 'One file has a different number of channels or a different number of time samples.');
209 return;
210 end
211 % Get good channels
212 iGoodChan = find(DataMat.ChannelFlag == 1);
213
214 % Average baseline values
215 FavgActive = mean(DataMat.F(iGoodChan,iMinVarTime), 2);
216 %FavgActive = 0;
217 % Remove average
218 DataActive = bst_bsxfun(@minus, DataMat.F(iGoodChan,iMinVarTime), FavgActive);
219 % Compute covariance for this file
220 fileActiveCov = DataMat.nAvg .* (DataActive * DataActive');
221 % Add file covariance to accumulator
222 MinVarCov(iGoodChan,iGoodChan) = MinVarCov(iGoodChan,iGoodChan) + fileActiveCov;
223 nTotalMinVar(iGoodChan,iGoodChan) = nTotalMinVar(iGoodChan,iGoodChan) + length(iMinVarTime);
224
225 for j = 1:nNAI
226 % Average baseline values
227 FavgNAIActive = mean(DataMat.F(iGoodChan,iActiveWindowTime(j,:)), 2);
228 %FavgNAIActive = 0;
229 % Remove average
230 DataNAIActive = bst_bsxfun(@minus, DataMat.F(iGoodChan,iActiveWindowTime(j,:)), FavgNAIActive);
231 % Compute covariance for this file
232 fileNAIActiveCov = DataMat.nAvg * (DataNAIActive * DataNAIActive');
233 % Add file covariance to accumulator
234 NAIActiveCov(iGoodChan,iGoodChan,j) = NAIActiveCov(iGoodChan,iGoodChan,j) + fileNAIActiveCov;
235 nTotalNAIActive(iGoodChan,iGoodChan,j) = nTotalNAIActive(iGoodChan,iGoodChan,j) + nActiveWindowPoints;
236 end
237 bst_progress('inc',1);
238 end
239
240 % Remove zeros from N matrix
241 nTotalMinVar(nTotalMinVar <= 1) = 2;
242 nTotalNAIActive(nTotalNAIActive <= 1) = 2;
243 % Divide final matrix by number of samples
244 MinVarCov = MinVarCov ./ (nTotalMinVar - 1);
245 NAIActiveCov = NAIActiveCov ./ (nTotalNAIActive - 1);
246
247
248 % ===== PROCESS =====
249 % Extract the covariance matrix of the used channels
250 NoiseCovMat = load(file_fullpath(sChannelStudy.NoiseCov(1).FileName));
251 NoiseCov = NoiseCovMat.NoiseCov(iChannels,iChannels);
252 NAIActiveCov = NAIActiveCov(iChannels,iChannels,:);
253 MinVarCov = MinVarCov(iChannels,iChannels);
254
255 % Normalize the regularization parameter
256 % eigValues = eig(MinVarCov);
257 % Reg_alpha = Reg * max(eigValues);
258 bst_progress('start', 'Beamformer: LCMV', 'Calulating inverse covariance matrix...', 0, nNAI+2);
259 % Calculate the inverse of (Cn+alpha*I)
260 Cn_inv = pinv(NoiseCov);
261 bst_progress('inc',1);
262 Cm_inv = zeros(length(iChannels),length(iChannels),nNAI);
263 for i = 1:nNAI
264 % Calculate the inverse of (Cm+alpha*I)
265 Cm_inv(:,:,i) = pinv(NAIActiveCov(:,:,i));
266 bst_progress('inc',1);
267 end
268
269 MinVarCov_inv = regulatePInv(MinVarCov,Reg);
270 bst_progress('inc',1);
271
272 % Get forward field
273 if (OrientConstr == 3) % constrained LCMV
274 Kernel = bst_gain_orient(sHeadModel.Gain(iChannels,:), sHeadModel.GridOrient);
275 else
276 Kernel = sHeadModel.Gain(iChannels,:);
277 end
278 Kernel(abs(Kernel(:)) < eps) = eps; % Set zero elements to strictly non-zero
279
280
281
282 % ===== CALCULATE NEURAL ACTIVITY INDEX AND FILTERS =====
283 % Different methods
284 switch (OrientConstr)
285 % Unconstrained: maximum eigenvalue as the power
286 case 1
287 nComponents = 3;
288 strConstr = 'Unconstr/svd'; % eig -> svd, corrected in 20141009
289 fcn_lcmv = @(c)lambda1(lcmvPInv(c));
290 % Unconstrained: use trace to calculate the power
291 case 2
292 nComponents = 3;
293 strConstr = 'Unconstr/trace';
294 fcn_lcmv = @(c)trace(lcmvPInv(c));
295 % Constrained
296 case 3
297 nComponents = 1;
298 strConstr = 'Constr';
299 fcn_lcmv = @(c)1./c; % Inverse function
300 end
301 % Number of sources
302 [nChannels, nSources] = size(Kernel);
303 nSources = nSources / nComponents;
304 nBlockProgress = 300;
305 bst_progress('start', 'Beamformer: LCMV', 'Calculating neural activity index...', 0, round(nSources/nBlockProgress));
306 % Initialize variables
307 NAI = zeros(nSources, nNAI);
308 if isSaveFilter
309 spatialFilter = zeros(nSources, nChannels);
310 end
311 % Loop on sources
312 for i = 1:nSources
313 % Indices for unconstrained model
314 iu = nComponents * (i-1) + (1:nComponents);
315 % Obtain gain matrix
316 Gi = Kernel(:,iu);
317 % At each vertex, noise pwr = (G Cn^-1 G')^(-1)
318 NoisePower = fcn_lcmv(Gi' * Cn_inv * Gi);
319 % For each time point
320 for j = 1:nNAI
321 % At each vertex, src pwr = (G Cm^-1 G')^(-1)
322 SourcePower = fcn_lcmv(Gi' * Cm_inv(:,:,j) * Gi);
323 % Calculate the neural activity index:
324 NAI(i,j) = SourcePower / NoisePower;
325 end
326 % Spatial filter
327 if isSaveFilter
328 % Calculate spatial filter W
329 A = Gi' * MinVarCov_inv;
330 W = pinv(A * Gi) * A;
331 % Normalize the spatial filter using the power of control state
332 spatialFilter(iu,:) = W / sqrt(NoisePower);
333 end
334 % Increase by one every block of sources
335 if (mod(i,nBlockProgress) == 0)
336 bst_progress('inc',1);
337 end
338 end
339
340
341 % ===== ASSIGN MAPS AND KERNELS =====
342 bst_progress('start', 'Beamformer: LCMV', 'Saving Results...');
343
344 if (nNAI == 1)
345 ImageGridAmp = [NAI, NAI];
346 TimeIndex = [Time(1), Time(end)];
347 else % Interpolate the NAI to have the same temporal resolution as the data
348 ImageGridAmpOriginal = NAI;
349 ImageGridAmp = zeros(nSources, nTime);
350 for i=1:(nNAI-1)
351 InterpolateTimeWindow = NAIRange(1) + ActiveWindowSize/2 + [(i-1)*NAITResolu, i*NAITResolu];
352 iInterpolateTime = panel_time('GetTimeIndices', Time, InterpolateTimeWindow);
353 nInterpolateTime = length(iInterpolateTime);
354
355 ImageGridAmp(:,iInterpolateTime(1)) = ImageGridAmpOriginal(:,i);
356 ImageGridAmp(:,iInterpolateTime(end)) = ImageGridAmpOriginal(:,i+1);
357
358 for j=2:(nInterpolateTime-1)
359 InterpolatePercentage = (j-1)/(nInterpolateTime-1);
360 ImageGridAmp(:,iInterpolateTime(j)) = ImageGridAmpOriginal(:,i+1)*InterpolatePercentage + ImageGridAmpOriginal(:,i)*(1-InterpolatePercentage);
361 end
362 end
363 TimeIndex = Time;
364 end
365
366
367 % ===== NAI: SAVE THE RESULTS =====
368 % Comment
369 if ~isempty(Comment)
370 Comment = [Comment ': '];
371 end
372 if (nNAI == 1)
373 Comment1 = sprintf('%sLCMV (%s, %d-%dms)', Comment, strConstr, round(NAIRange(1)*1000), round(NAIRange(2)*1000));
374 else
375 Comment1 = sprintf('%sLCMV (%s, %d-%dms, ws:%dms, tr:%dms)', Comment, strConstr, round((NAIRange(1)+ActiveWindowSize/2)*1000), ...
376 round((NAIRange(1) + ActiveWindowSize/2 + (nNAI-1)*NAITResolu)*1000), round(ActiveWindowSize*1000), round(NAITResolu*1000));
377 end
378 % Create a new data file structure
379 ResultsMat = db_template('resultsmat');
380 ResultsMat.ImagingKernel = [];
381 ResultsMat.ImageGridAmp = ImageGridAmp;
382 ResultsMat.Comment = Comment1;
383 ResultsMat.nComponents = 1; % 1 or 3
384 ResultsMat.Function = 'LCMVneuralActivityIndex';
385 ResultsMat.Time = TimeIndex; % Leave it empty if using ImagingKernel
386 ResultsMat.DataFile = [];
387 ResultsMat.HeadModelFile = HeadModelFile;
388 ResultsMat.HeadModelType = sHeadModel.HeadModelType;
389 ResultsMat.ChannelFlag = [];
390 ResultsMat.GoodChannel = iChannels;
391 ResultsMat.SurfaceFile = sHeadModel.SurfaceFile;
392 ResultsMat.GridLoc = []; %sHeadModel.GridLoc
393 % Get the output study (pick the one from the first file)
394 iStudy = sInputs(1).iStudy;
395 % Create a default output filename
396 OutputFiles{1} = bst_process('GetNewFilename', fileparts(sInputs(1).FileName), 'results_LCMV_NAI');
397 % Save on disk
398 save(OutputFiles{1}, '-struct', 'ResultsMat');
399 % Register in database
400 db_add_data(iStudy, OutputFiles{1}, ResultsMat);
401
402
403 % ===== SPATIAL FILTER: SAVE FILE =====
404 if isSaveFilter
405 % Set file comment
406 Comment2 = sprintf('%sLCMV/filter (%s, %d-%dms)', Comment, strConstr, round(MinVarTime(1)*1000), round(MinVarTime(2)*1000));
407 % Create new structure
408 ResultsMat2 = db_template('resultsmat');
409 ResultsMat2.ImagingKernel = spatialFilter;
410 ResultsMat2.ImageGridAmp = [];
411 ResultsMat2.Comment = Comment2;
412 ResultsMat2.Function = 'LCMVspatialfilter';
413 ResultsMat2.Time = []; % Leave it empty if using ImagingKernel
414 ResultsMat2.DataFile = [];
415 ResultsMat2.HeadModelFile = HeadModelFile;
416 ResultsMat2.HeadModelType = sHeadModel.HeadModelType;
417 ResultsMat2.nComponents = nComponents;
418 ResultsMat2.ChannelFlag = [];
419 ResultsMat2.GoodChannel = iChannels;
420 ResultsMat2.SurfaceFile = sHeadModel.SurfaceFile;
421 ResultsMat2.GridLoc = sHeadModel.GridLoc;
422 % Get the output study (pick the one from the first file)
423 iStudy = iChannelStudy;
424 % Create a default output filename (shared kernel)
425 OutputFiles{2} = bst_process('GetNewFilename', fileparts(sInputs(1).ChannelFile), 'results_LCMV_KERNEL');
426 % Save on disk
427 save(OutputFiles{2}, '-struct', 'ResultsMat2');
428 % Register in database
429 db_add_data(iStudy, OutputFiles{2}, ResultsMat2);
430 end
431
432 % profile viewer;
433
434 end
435
436
437
438 %% ===== TRUNCATED PSEUDO-INVERSE =====
439 function X = lcmvPInv(A)
440 % Inverse of 3x3 GCG' in unconstrained beamformers.
441 % Since most head models have rank 2 at each vertex, we cut all the fat and
442 % just take a rank 2 inverse of all the 3x3 matrices
443 [U,S,V] = svd(A,0);
444 S = diag(S);
445 X = V(:,1:2)*diag(1./S(1:2))*U(:,1:2)';
446 end
447
448 function X = regulatePInv(A,Reg)
449 eigValues = eig(A);
450 Reg_alpha = (Reg/100) * max(eigValues);
451 X = inv(A + Reg_alpha*eye(size(A)));
452 % Calculate the inverse of (Cm+alpha*I)
453 % [U,S,V] = svd(A,0);
454 % S = diag(S); % Covariance = Cm = U*S*U'
455 % Si = diag(1 ./ (S + S(1) * Reg / 100)); % 1/(S^2 + lambda I)
456 % X = V*Si*U'; % U * 1/(S^2 + lambda I) * U' = Cm^(-1)
457 end
458
459 function X = lambda1(A)
460 [U,S,V] = svd(A,0);
461 X = S(1,1);
462 end
Attached Files
To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.You are not allowed to attach a file to this page.