Main Data History
Show Index Toggle 0 comments
  •  Quick Edit
  • Homework 3

    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

    \label{fig_datac1} Subset of channels from data1.