mne_ex_compute_inverse

 [res] = mne_ex_compute_inverse(fname_data,setno,fname_inv,nave,lambda2,dSPM,sLORETA)

 An example on how to compute a L2-norm inverse solution
 Actual code using these principles might be different because 
 the inverse operator is often reused across data sets.


 fname_data  - Name of the data file
 setno       - Data set number
 fname_inv   - Inverse operator file name
 nave        - Number of averages (scales the noise covariance)
               If negative, the number of averages in the data will be
               used
 lambda2     - The regularization factor
 dSPM        - do dSPM?
 sLORETA     - do sLORETA?

0001 function [res] = mne_ex_compute_inverse(fname_data,setno,fname_inv,nave,lambda2,dSPM,sLORETA)
0002 %
0003 % [res] = mne_ex_compute_inverse(fname_data,setno,fname_inv,nave,lambda2,dSPM,sLORETA)
0004 %
0005 % An example on how to compute a L2-norm inverse solution
0006 % Actual code using these principles might be different because
0007 % the inverse operator is often reused across data sets.
0008 %
0009 %
0010 % fname_data  - Name of the data file
0011 % setno       - Data set number
0012 % fname_inv   - Inverse operator file name
0013 % nave        - Number of averages (scales the noise covariance)
0014 %               If negative, the number of averages in the data will be
0015 %               used
0016 % lambda2     - The regularization factor
0017 % dSPM        - do dSPM?
0018 % sLORETA     - do sLORETA?
0019 %
0020 
0021 %
0022 %
0023 %   Author : Matti Hamalainen, MGH Martinos Center
0024 %   License : BSD 3-clause
0025 %
0026 %   Revision 1.2  2008/05/26 10:49:26  msh
0027 %   Update to incorporate the already weighted lead field basis
0028 %
0029 %   Revision 1.1  2006/05/05 03:50:40  msh
0030 %   Added routines to compute L2-norm inverse solutions.
0031 %   Added mne_write_inverse_sol_stc to write them in stc files
0032 %   Several bug fixes in other files
0033 %
0034 %
0035 
0036 me='MNE:mne_ex_compute_inverse';
0037 
0038 global FIFF;
0039 if isempty(FIFF)
0040    FIFF = fiff_define_constants();
0041 end
0042 
0043 if nargin ~= 6 && nargin ~= 7
0044    error(me,'Incorrect number of arguments'); 
0045 end
0046 
0047 if nargin == 6
0048    sLORETA = false;
0049 end
0050 %
0051 %   Read the data first
0052 %
0053 data = fiff_read_evoked(fname_data,setno);
0054 %
0055 %   Then the inverse operator
0056 %
0057 inv = mne_read_inverse_operator(fname_inv);
0058 %
0059 %   Set up the inverse according to the parameters
0060 %
0061 if nave < 0
0062     nave = data.evoked.nave;
0063 end
0064 inv = mne_prepare_inverse_operator(inv,nave,lambda2,dSPM,sLORETA);
0065 %
0066 %   Pick the correct channels from the data
0067 %
0068 data = fiff_pick_channels_evoked(data,inv.noise_cov.names);
0069 fprintf(1,'Picked %d channels from the data\n',data.info.nchan);
0070 fprintf(1,'Computing inverse...');
0071 %
0072 %   Simple matrix multiplication followed by combination of the
0073 %   three current components
0074 %
0075 %   This does all the data transformations to compute the weights for the
0076 %   eigenleads
0077 %
0078 trans = diag(sparse(inv.reginv))*inv.eigen_fields.data*inv.whitener*inv.proj*double(data.evoked(1).epochs);
0079 %
0080 %   Transformation into current distributions by weighting the eigenleads
0081 %   with the weights computed above
0082 %
0083 if inv.eigen_leads_weighted
0084    %
0085    %     R^0.5 has been already factored in
0086    %
0087    fprintf(1,'(eigenleads already weighted)...');
0088    sol   = inv.eigen_leads.data*trans;
0089 else
0090    %
0091    %     R^0.5 has to factored in
0092    %
0093    fprintf(1,'(eigenleads need to be weighted)...');
0094    sol   = diag(sparse(sqrt(inv.source_cov.data)))*inv.eigen_leads.data*trans;
0095 end
0096    
0097 if inv.source_ori == FIFF.FIFFV_MNE_FREE_ORI
0098     fprintf(1,'combining the current components...');
0099     sol1 = zeros(size(sol,1)/3,size(sol,2));
0100     for k = 1:size(sol,2)
0101         sol1(:,k) = sqrt(mne_combine_xyz(sol(:,k)));
0102     end
0103     sol = sol1;
0104 end
0105 if dSPM
0106     fprintf(1,'(dSPM)...');
0107     sol = inv.noisenorm*sol;
0108 elseif sLORETA
0109     fprintf(1,'(sLORETA)...');
0110     sol = inv.noisenorm*sol;
0111 end
0112 res.inv   = inv;
0113 res.sol   = sol;
0114 res.tmin  = double(data.evoked(1).first)/data.info.sfreq;
0115 res.tstep = 1/data.info.sfreq;
0116 fprintf(1,'[done]\n');
0117 
0118 return;
0119 end
0120