Skip to content

Commit

Permalink
Added SimMEEG 21a to master branch
Browse files Browse the repository at this point in the history
  • Loading branch information
@branelab
branelab committed Nov 9, 2021
1 parent a9141cf commit ac37a47
Show file tree
Hide file tree
Showing 100 changed files with 19,088 additions and 0 deletions.
86 changes: 86 additions & 0 deletions BRANELab_calc_cov.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,86 @@
function [R,N,Rbar,Nbar,Rinv,Ninv]=BRANELab_calc_cov(data,act_samps, ctrl_samps)
%function [R,N,Rbar,Nbar,Rinv,Ninv]=BRANELab_calc_cov(data,act_samps, ctrl_samps)
% INPUT:
% data = post-stimulus singla + noise matrix (samples x channels x trials)
% act_samps = samples in active interval used for calculating signal-to-noise ratio.
% ctrl_samps = samples in control interval used for calculating signal-to-noise ratio.
%
% OUTPUT:
% R = covariance matrix for active interval to be used in "BRANElab_LCMV_Scalar_beamformer.m"
% N = covariance matrix for control interval to be used in "BRANElab_LCMV_Scalar_beamformer.m"
% Rinv = inverted covariance matrix for active interval to be used in "BRANElab_LCMV_Scalar_beamformer.m"
% Ninv = inverted covariance matrix for control interval to be used in "BRANElab_LCMV_Scalar_beamformer.m"
% Rbar = covariance matrix for active interval for multi-source event-related beamformer covariance Rbar =((b-mean(b))*(b-mean(b))'). see Moiseev et al. 2011 NeuroImage.
% Nbar = covariance matrix for control interval for multi-source event-related beamformer covariance Rbar =((b-mean(b))*(b-mean(b))'). see Moiseev et al. 2011 NeuroImage.
%
% written by A. Herdman, UBC June 13, 2015

%% NOTES
% tikhonov = Tikhonov regularization as a percent (default = 10%). NOTE: This is only used if Rinv or Ninv are badly scaled (not invertible). ;
% Using pinv(R) instead of svd(R) and regularization

%%
%if nargin<4; tikhonov=10; end
RN_type = {'Active Interval' 'Control Interval'};
RN_name={'Rinv' 'Ninv'};
act_samps(act_samps==0)=1;
for c=1:2; % loop trhough for active and control intervals
% first baselining data to mean for each trial;
if c==1;
xdata = data(act_samps,:,:) - repmat( nanmean(data(act_samps,:,:),1), [length(act_samps) 1 1]);
elseif c==2;
xdata = data(ctrl_samps,:,:) - repmat( nanmean(data(ctrl_samps,:,:),1), [length(ctrl_samps) 1 1]);
end
[num_samps,num_chans,num_trials]=size(xdata);

% fprintf('calculate covariance then average across trials: %s\n',RN_type{c});
cov_dat=nan(num_chans,num_chans,num_trials);
for t=1:num_trials;
F=squeeze(xdata(:,:,t))';
RN_data(:,:,t)=cov(F');
%RN_data(:,:,t)=(F*F')./ ((ones(size(F,1),size(F,1))*num_samps)-1); %This is equivalent to cov(F') in line above.
end
cov_dat=squeeze(nanmean(cov_dat,3)); % averaging across trials
RN_data=squeeze(nanmean(RN_data,3)); % averaging across trials

% calculating Rbar for multi-source event-related beamformer covariance Rbar =((b-mean(b))*(b-mean(b))'). see Moiseev et al. 2011 NeuroImage.
% fprintf('average across trials then calculate covariance of average: %s\n',RN_type{c});
F=squeeze(nanmean(xdata,3)); % averaging across trials to get evoked data then calculating covariance
Rbar_data=cov(F,1);
% Code below gives almost exactly same result as above Rbar_data=cov(F,1). difference is rounding errors close to the floating-poinjt precision of the system.
% for s=1:size(F,1); Rbar_data2(s,:,:) = ( F(s,:) -mean(F))'*(F(s,:)-mean(F)); end;
% Rbar_data2=squeeze(nanmean(Rbar_data2));

% if c==1; Rinv=cov_dat; R=RN_data; elseif c==2; Ninv=cov_dat; N=RN_data; end
if c==1; % for act_int
R=RN_data;
Rinv=inv(RN_data); % covariance matrix should never be near singular for BRANELAb's LCMV beamformers - if it is then there is a problem.
Rbar = Rbar_data;
if rcond(Rinv)<=eps;
fprintf('\nWARNING: Inverted Matrix "%s" is badly scaled and is near singular to working precision!\nPlease view the covariance matrices using surf(%s)\n',RN_name{c},RN_name{c});
Rinv=pinv(R); % pseudo-inverse of covariance (similar to what is done in NutMEG
% fprintf('Inverting R matrix using Tikhonov regularization: %.f%%\n',tikhonov);
% % code from Brainstorm3 for inverting a matrix with regularization. NOT IMPLEMENTED because matrix should be invertible.
% [U,S] = svd(R); % Covariance = 1/(m-1) * F * F'
% S = diag(S); % Covariance = Cm = U*S*U'
% Si = diag(1 ./ (S + S(1) * tikhonov / 100)); % 1/(S^2 + lambda I)
% Rinv = U*Si*U'; % U * 1/(S^2 + lambda I) * U' = Cm^(-1)
end
elseif c==2; % for ctrl_int
N=RN_data;
Ninv=inv(RN_data); % covariance matrix should never be near singular for BRANELAb's LCMV beamformers - if it is then there is a problem.
Nbar = Rbar_data;
if rcond(Ninv)<=eps;
fprintf('\nWARNING: Inverted Matrix "%s" is badly scaled and is near singular to working precision!\nPlease view the covariance matrices using surf(%s)\n',RN_name{c},RN_name{c});
Ninv=pinv(N); % pseudo-inverse of covariance (similar to what is done in NutMEG
% fprintf('Inverting N matrix using Tikhonov regularization: %.f%%\n',tikhonov);
% % code from Brainstorm3 for inverting a matrix with regularization. NOT IMPLEMENTED because matrix should be invertible.
% [U,S] = svd(N); % Covariance = 1/(m-1) * F * F'
% S = diag(S); % Covariance = Cm = U*S*U'
% Si = diag(1 ./ (S + S(1) * tikhonov / 100)); % 1/(S^2 + lambda I)
% Ninv = U*Si*U'; % U * 1/(S^2 + lambda I) * U' = Cm^(-1)
end
end
end


57 changes: 57 additions & 0 deletions BRANELab_find_peak_voxel_thresh.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,57 @@
function [peak_voxel,p_idx]=BRANELab_find_peak_voxel_thresh(voxel_vals,thresh_val,thresh_limit);
%function [peak_voxel,p_idx]=BRANELab_find_peak_voxel_thresh(voxel_vals,thresh_val,thresh_limit);
%
% This program will step through the voxels and find the peak absolute amplitude
% within a region separated by the thresh_limit value
%
% INPUT:
% voxel_vals = xyz coordinates of each voxel (M x N = voxels x [x y z amp ...] ) must be at least an M x 4 matrix
% thresh_val = value of voxel amplitude to threshold data before peak searching
% thresh_limit = minimum distance between peak voxels (distance thus depends on voxel step size). This creates a bounding box with thresh limit distance.
%
% OUTPUT:
% [peak_voxel] = peak voxels with all the columns in original voxel_vals
%
% written by A. Herdman June 30,2005
%
% modified by A. Herdman June 6, 2015 - updated description to make more sense.
%

voxel_abs=voxel_vals;
voxel_abs(:,4)=abs(voxel_vals(:,4));

if length(thresh_val)==1;
[lv_idx]=find(voxel_abs(:,4)>thresh_val);
elseif length(thresh_val)==2;
[lv_idx1]=find(voxel_abs(:,4)>thresh_val(2));
[lv_idx2]=find(voxel_abs(:,4)<thresh_val(1));
lv_idx=[lv_idx1;lv_idx2];
end

p_idx=[];
if isempty(lv_idx)~=1

voxel_lv=voxel_abs(lv_idx,:);
voxel_org=voxel_vals(lv_idx,:);

p=0;

for n=1:size(voxel_lv,1)
vox_roi=voxel_lv(n,:);
vox_roi_org=voxel_org(n,:);

[roi_idx,z,zz]=find((voxel_lv(:,1)< vox_roi(1,1)+thresh_limit & voxel_lv(:,1)>vox_roi(1,1)-thresh_limit)...
& (voxel_lv(:,2)< vox_roi(1,2)+thresh_limit & voxel_lv(:,2)>vox_roi(1,2)-thresh_limit)...
& (voxel_lv(:,3)< vox_roi(1,3)+thresh_limit & voxel_lv(:,3)>vox_roi(1,3)-thresh_limit));

if vox_roi(1,4)==max(voxel_lv(roi_idx,4))
p=p+1;
%fprintf(1,'Number of peaks = %.f\n',p)
peak_voxel(p,:)=vox_roi_org;
p_idx=[p_idx n];
end
end
else
peak_voxel=[];
end
p_idx=lv_idx(p_idx);
102 changes: 102 additions & 0 deletions BRANElab_LCMV_beamformer_SPA.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,102 @@
function [SPA]=BRANElab_LCMV_beamformer_SPA(H,data,act_samps,ctrl_samps,loc_flag,noise_alpha)
%function [wts,Hscalar,ori,P,RNcov,null_thresh]=BRANElab_LCMV_beamformer_SPA(H,data,act_samps,ctrl_samps,loc_flag)
% This program will find voxels with maximal amplitudes from scalar beamformer
%
% INPUT:
% H = lead field [channels x 3 orientations x voxels]
% data = channel data used to find voxels with maximum signal-to-noise. [samples x channels x trials] .
% act_samps = samples in active interval used for calculating signal-to-noise ratio.
% ctrl_samps = samples in control interval used for calculating signal-to-noise ratio.
% loc_flag = (0) returns 'pseudoZ' for single-source or 'MPZ' for multisource.
% (1) returns 'pseudoZ_er' for single-source or 'MER' for multisource.
% (2) returns 'pseudoZ_rer' for single-source or 'rMER' for multisource.
% (3) returns 'activity_index' for single-source or 'MAI' for multisource.
% Note: single-source localizer returned if Href=[]; or single-source localizer returned if Href=[channels x voxels].
%
% OUTPUT:
% wts = source weights [channels x voxels];
% Hscalar = scalar lead-field matrix based on single- or mult-source
% ori = orientation of voxels [voxels x (X,Y,Z)];
% P.type = possible localizers - formulae from Moiseev et al., 2011 .
% pseudoZ = (H'*Rinv*H)/(H'*Rinv*N*Rinv*H);
% pseudoZ_er = (H'*Rinv*Rbar*Rinv*H)/(H'*Rinv*N*Rinv*H)
% pseudoZ_rer = (H'*Rinv*Rbar*Rinv*H)/(H'*Rinv*H)
% activity_index = (H'*Ninv*H)/(H'*Rinv*H)
% MPZ = trace((H'*Rinv*H)/(H'*Rinv*N*Rinv*H))-n, where n = number of multiple sources.
% MER = trace((H'*Rinv*Rbar*Rinv*H)/(H'*Rinv*N*Rinv*H))-n, where n = number of multiple sources.
% rMER = trace((H'*Rinv*Rbar*Rinv*H)/(HF'*Rinv*HF))-n, where n = number of multiple sources.
% MAI = trace((H'*Ninv*H)/(HF'*Rinv*HF))-n, where n = number of multiple sources.
% P.img = image for localizer values for the active interval [voxels x 1]
% P_ctrl.img = image for localizer values for the control interval [voxels x 1]
% P_log.img = normalized power image relative to ctrl_int: P_log.img=20*log(P.img./P_ctrl.img);
%
% RNcov = R, Rinv, N, Ninv, Rbar, Nbar covariance matrices.
% null_thresh = localizer threshold based on LCMV of ctrl interval by splitting ctrl samps into half active and half control.
%
% written by A. Herdman Jun 13, 2015
%

%% beamformer based on following equation
% W = R^-1 * H * ( H' * R^-1 * H)^-1 , where R = covariance matrix [sensor x sensor] and H = Ledfield [sensor x voxel].


fprintf('Performing Single-Source beamforming\n');

%% LCMV
% calculating covarince matrices for signal+noise (Rinv) and noise (Ninv
xs=round(size(ctrl_samps,2)/2);

%% checking rank of data. If it doesn't match the num_chans then white noise will be added until rank sufficient.
R=[]; t=0;
% r2=rank(squeeze(data(:,:,1)));
[R,~,~,~,~,~]=BRANELab_calc_cov(data,act_samps,ctrl_samps(xs+1:end));
r2=rank(R);
xstd=min(min(std(data)))/size(data,2);
while r2<size(H,1)
t=t+1;
data=data+(xstd*randn(size(data)));
[R,~,~,~,~,~]=BRANELab_calc_cov(data,act_samps,ctrl_samps(xs+1:end));
r2=rank(R);
fprintf('Rank = %.f\tChannel Count=%.f\n',r2,size(data,2));
if 100*(t*xstd)/max(max(std(data)))>10; fprintf('WARNING! Insufficient Rank in Data or Data is empty\n'); SPA=[]; return; end % stopping infinite loop if added noise exceeds 10%
% fprintf('Rank = %.f\n',r2);
fprintf('Percent white noise added for sufficient rank = %.3f %%\n\n',100*(t*xstd)/max(max(std(data))));
end

%[R,N,Rbar,Nbar,Rinv,Ninv]=BRANELab_calc_cov(data,act_samps,ctrl_samps);
[R,N,Rbar,Nbar,Rinv,Ninv]=BRANELab_calc_cov(data,act_samps,ctrl_samps(xs+1:end));
RNcov.R=R; RNcov.Rinv=Rinv; RNcov.Rbar=Rbar; RNcov.N=N; RNcov.Ninv=Ninv; RNcov.Nbar=Nbar;

H_idx=1:size(H,3); [wts,Hscalar,ori,P]=bl_lcmv_scalar_v4(H,H_idx,[],RNcov,[],loc_flag);


%% Noise estimate from control samples
% splitting ctrl interval in half to calculate null distribution of noise
%xs=round(size(ctrl_samps,2)/2);
[nR,nN,nRbar,nNbar,nRinv,nNinv]=BRANELab_calc_cov(data,ctrl_samps(1:xs),ctrl_samps(xs+1:end));
nRNcov.R=nR; nRNcov.Rinv=nRinv; nRNcov.Rbar=nRbar; nRNcov.N=nN; nRNcov.Ninv=nNinv; nRNcov.Nbar=nNbar;

% localizer image for ctrl_samps.
if loc_flag==0; Rx=nR; Nx=nN; bm_type2='pseudoZ'; % pseudoZ
elseif loc_flag==1; Rx=nRbar; Nx=nN; bm_type2='pseudoZ_er'; % pseudoZ_er
elseif loc_flag==2; Rx=nRbar; Nx=nR; bm_type2='pseudoZ_rer'; % pseudoZ_er
elseif loc_flag==3; Rx=nN; Nx=nR; bm_type2='acitvity_index'; % pseudoZ_er2
end

% First calculation of noise - like SPA
H_idx=1:size(H,3); [~,~,~,nP]=bl_lcmv_scalar_v4(H,H_idx,[],nRNcov,[],loc_flag);
nd=sort(nP.img); null_thresh=nd(ceil(length(nd)*noise_alpha));

% OUTPUTS
SPA.wts=wts;
SPA.Hscalar=Hscalar;
SPA.ori=ori;
SPA.P=P;
SPA.P_ctrl=nP;
SPA.RNcov=RNcov;
SPA.null_thresh=null_thresh;
% SPA.P_log.type='Log normalized Map';
% SPA.P_log.img=20*log(P.img./nP.img);
% SPA.P_SNR.img=SPA.P.img./SPA.P_ctrl.img;


Loading

0 comments on commit ac37a47

Please sign in to comment.