mne_prepare_inverse_operator

 [inv] = mne_prepare_inverse_operator(orig,nave,lambda2,dSPM,sLORETA)

 Prepare for actually computing the inverse

 orig        - The inverse operator structure read from a file
 nave        - Number of averages (scales the noise covariance)
 lambda2     - The regularization factor
 dSPM        - Compute the noise-normalization factors for dSPM?
 sLORETA     - Compute the noise-normalization factors for sLORETA?

0001 function [inv] = mne_prepare_inverse_operator(orig,nave,lambda2,dSPM,sLORETA)
0002 %
0003 % [inv] = mne_prepare_inverse_operator(orig,nave,lambda2,dSPM,sLORETA)
0004 %
0005 % Prepare for actually computing the inverse
0006 %
0007 % orig        - The inverse operator structure read from a file
0008 % nave        - Number of averages (scales the noise covariance)
0009 % lambda2     - The regularization factor
0010 % dSPM        - Compute the noise-normalization factors for dSPM?
0011 % sLORETA     - Compute the noise-normalization factors for sLORETA?
0012 %
0013 
0014 %
0015 %
0016 %   Author : Matti Hamalainen, MGH Martinos Center
0017 %   License : BSD 3-clause
0018 %
0019 %   Revision 1.4  2009/03/09 11:22:23  msh
0020 %   Fixed the noise-normalization factors for the case where eigen leads come
0021 %   already weighted
0022 %
0023 %   Revision 1.3  2008/05/26 10:49:26  msh
0024 %   Update to incorporate the already weighted lead field basis
0025 %
0026 %   Revision 1.2  2006/05/05 19:37:47  msh
0027 %   Fixed error in mne_make_projector.
0028 %   Better detection of small eigenvalues for the projector.
0029 %
0030 %   Revision 1.1  2006/05/05 03:50:40  msh
0031 %   Added routines to compute L2-norm inverse solutions.
0032 %   Added mne_write_inverse_sol_stc to write them in stc files
0033 %   Several bug fixes in other files
0034 %
0035 %
0036 
0037 me='MNE:mne_prepare_inverse_operator';
0038 
0039 global FIFF;
0040 if isempty(FIFF)
0041     FIFF = fiff_define_constants();
0042 end
0043 
0044 if nargin ~= 4 && nargin ~=5
0045     error(me,'Wrong number of arguments');
0046  end
0047 
0048 if nargin == 4
0049    sLORETA = false;
0050 end
0051 
0052 if nave <= 0
0053     error(me,'The number of averages should be positive');
0054 end
0055 fprintf(1,'Preparing the inverse operator for use...\n');
0056 inv = orig;
0057 %
0058 %   Scale some of the stuff
0059 %
0060 scale               = double(inv.nave)/double(nave);
0061 inv.noise_cov.data  = scale*inv.noise_cov.data;
0062 inv.noise_cov.eig   = scale*inv.noise_cov.eig;
0063 inv.source_cov.data = scale*inv.source_cov.data;
0064 %
0065 if inv.eigen_leads_weighted
0066     inv.eigen_leads.data = sqrt(scale)*inv.eigen_leads.data;
0067 end
0068 %
0069 fprintf(1,'\tScaled noise and source covariance from nave = %d to nave = %d\n',inv.nave,nave);
0070 inv.nave = nave;
0071 %
0072 %   Create the diagonal matrix for computing the regularized inverse
0073 %
0074 inv.reginv = inv.sing./(inv.sing.*inv.sing + lambda2);
0075 fprintf(1,'\tCreated the regularized inverter\n');
0076 %
0077 %   Create the projection operator
0078 %
0079 [ inv.proj, ncomp ] = mne_make_projector(inv.projs,inv.noise_cov.names);
0080 if ncomp > 0
0081     
0082     fprintf(1,'\tCreated an SSP operator (subspace dimension = %d)\n',ncomp);
0083 end
0084 %
0085 %   Create the whitener
0086 %
0087 inv.whitener = zeros(inv.noise_cov.dim);
0088 if inv.noise_cov.diag == 0
0089     %
0090     %   Omit the zeroes due to projection
0091     %
0092     nnzero = 0;
0093     for k = ncomp+1:inv.noise_cov.dim
0094         if inv.noise_cov.eig(k) > 0
0095             inv.whitener(k,k) = 1.0/sqrt(inv.noise_cov.eig(k));
0096             nnzero = nnzero + 1;
0097         end
0098     end
0099     %
0100     %   Rows of eigvec are the eigenvectors
0101     %
0102     inv.whitener = inv.whitener*inv.noise_cov.eigvec;
0103     fprintf(1,'\tCreated the whitener using a full noise covariance matrix (%d small eigenvalues omitted)\n', ...
0104         inv.noise_cov.dim - nnzero);
0105 else
0106     %
0107     %   No need to omit the zeroes due to projection
0108     %
0109     for k = 1:inv.noise_cov.dim
0110         inv.whitener(k,k) = 1.0/sqrt(inv.noise_cov.data(k));
0111     end
0112     fprintf(1,'\tCreated the whitener using a diagonal noise covariance matrix (%d small eigenvalues discarded)\n',ncomp);
0113 end
0114 %
0115 %   Finally, compute the noise-normalization factors
0116 %
0117 if dSPM || sLORETA
0118     noise_norm = zeros(inv.eigen_leads.nrow,1);
0119     if dSPM
0120        fprintf(1,'\tComputing noise-normalization factors (dSPM)...');
0121        noise_weight = inv.reginv;
0122     else
0123        fprintf(1,'\tComputing noise-normalization factors (sLORETA)...');
0124        noise_weight = inv.reginv.*sqrt((1 + inv.sing.*inv.sing/lambda2));
0125     end
0126     if inv.eigen_leads_weighted
0127        for k = 1:inv.eigen_leads.nrow
0128           one = inv.eigen_leads.data(k,:).*noise_weight';
0129           noise_norm(k) = sqrt(one*one');
0130        end
0131     else
0132        for k = 1:inv.eigen_leads.nrow
0133           one = sqrt(inv.source_cov.data(k))*(inv.eigen_leads.data(k,:).*noise_weight');
0134           noise_norm(k) = sqrt(one*one');
0135        end
0136     end
0137     %
0138     %   Compute the final result
0139     %
0140     if inv.source_ori == FIFF.FIFFV_MNE_FREE_ORI
0141         %
0142         %   The three-component case is a little bit more involved
0143         %   The variances at three consequtive entries must be squeared and
0144         %   added together
0145         %
0146         %   Even in this case return only one noise-normalization factor
0147         %   per source location
0148         %
0149         noise_norm = sqrt(mne_combine_xyz(noise_norm));
0150         %
0151         %   This would replicate the same value on three consequtive
0152         %   entries
0153         %
0154         %   noise_norm = kron(sqrt(mne_combine_xyz(noise_norm)),ones(3,1));
0155     end
0156     inv.noisenorm = diag(sparse(1./abs(noise_norm)));
0157     fprintf(1,'[done]\n');
0158 else
0159     inv.noisenorm = [];
0160 end
0161 
0162 return;
0163 
0164 end