To construct the L2 MNE operator, we used the formula found in the lecture slides: \[P_{MNE}=RL^T(LRL^T + \lambda C_n)^{-1}\]
The following Matlab script was used in the simple case without depth weighting:
load ex3_data.mat
%% MNE
L2=L*L';
snr = sqrt(n_trial);
lambda = trace(L2)/(trace(Cn)*snr^2);
ev = eigs(L2+lambda*Cn, rank(data1));
tol = ev(end)
p_mne=L'*pinv(L2+lambda*Cn, tol);
src1 = p_mne*data1;
src2 = p_mne*data2;
Here an identity matrix was used for the source covariance matrix.
Next, to include depth weighting, we modified the source covariance matrix:
%% Depth-Weighted MNE
W_i = diag(sqrt(sum(L,1).^2./306));
lambda = trace(L*W_i*L')/(trace(Cn)*snr^2);
ev = eigs(L*W_i*L'+lambda*Cn, rank(data1));
tol = ev(end)
pw_mne=W_i*L'*pinv(L*W_i*L'+lambda*Cn,tol);
srcw1 = pw_mne*data1;
srcw2 = pw_mne*data2;
Finally, we constructed a beamformer operator for the data:
%% Beamformer
N = 306;
M = 5124;
ev = eigs(Cd, rank(data1));
tol = ev(end)
p_bf = (pinv(Cd,tol)*L)';
denominator = diag(p_bf*L);
p_bf = p_bf./repmat(denominator,1, N);
srcb1 = p_bf*data1;
srcb2 = p_bf*data2;
using the following formula from the lecture slides \[P_{BF,\theta} = \dfrac{(C_\mathrm{d}^{-1} L_\theta)^T}{L_\theta^T C_\mathrm{d}^{-1} L_\theta},\] where \(L_\theta\) is the gain vector for source point \(\theta\).
We wanted to visualize the source estimates as a function of time. We achieved this with the following script:
%% Plots
close all
vis_surface_data(srcb1(:,115), 0.1, max(srcb1(:)), anat_decim)
%%
close all
figure(3)
hold on
vv = var(data1,[], 2);
plot(timeaxis,data1(vv > 0.3*max(vv),:))
plot([1,1]*timeaxis(115),ylim(gca()), 'k--')
plot([1,1]*timeaxis(140),ylim(gca()), 'k--')
plot([1,1]*timeaxis(190),ylim(gca()), 'k--')
xlabel('time [s]')
axis tight
%%
close all
figure(4)
hold on
vv = var(data2,[], 2);
plot(timeaxis,data2(vv > 0.3*max(vv),:))
plot([1,1]*timeaxis(125),ylim(gca()), 'k--')
plot([1,1]*timeaxis(190),ylim(gca()), 'k--')
xlabel('time [s]')
axis tight