[proj,nproj,U] = mne_make_projector(projs,ch_names,bads) proj - The projection operator to apply to the data nproj - How many items in the projector U - The orthogonal basis of the projection vectors (optional) Make an SSP operator projs - A set of projection vectors ch_names - A cell array of channel names bads - Bad channels to exclude
0001 function [proj,nproj,U] = mne_make_projector(projs,ch_names,bads) 0002 % 0003 % [proj,nproj,U] = mne_make_projector(projs,ch_names,bads) 0004 % 0005 % proj - The projection operator to apply to the data 0006 % nproj - How many items in the projector 0007 % U - The orthogonal basis of the projection vectors (optional) 0008 % 0009 % Make an SSP operator 0010 % 0011 % projs - A set of projection vectors 0012 % ch_names - A cell array of channel names 0013 % bads - Bad channels to exclude 0014 % 0015 0016 % 0017 % Author : Matti Hamalainen, MGH Martinos Center 0018 % License : BSD 3-clause 0019 % 0020 % 0021 % Revision 1.11 2008/11/22 17:06:42 msh 0022 % Optionally return U 0023 % 0024 % Revision 1.10 2006/05/05 19:37:47 msh 0025 % Fixed error in mne_make_projector. 0026 % Better detection of small eigenvalues for the projector. 0027 % 0028 % Revision 1.9 2006/05/05 03:50:40 msh 0029 % Added routines to compute L2-norm inverse solutions. 0030 % Added mne_write_inverse_sol_stc to write them in stc files 0031 % Several bug fixes in other files 0032 % 0033 % Revision 1.8 2006/04/26 00:43:22 msh 0034 % Fixed errors in mne_make_projector related to vectors which do not affect the data 0035 % 0036 % Revision 1.7 2006/04/23 15:29:40 msh 0037 % Added MGH to the copyright 0038 % 0039 % Revision 1.6 2006/04/21 14:23:16 msh 0040 % Further improvements in raw data reading 0041 % 0042 % Revision 1.5 2006/04/18 20:44:46 msh 0043 % Added reading of forward solution. 0044 % Use length instead of size when appropriate 0045 % 0046 % Revision 1.4 2006/04/17 15:01:34 msh 0047 % More small improvements. 0048 % 0049 % Revision 1.3 2006/04/15 12:21:00 msh 0050 % Several small improvements 0051 % 0052 % Revision 1.2 2006/04/14 15:49:49 msh 0053 % Improved the channel selection code and added ch_names to measurement info. 0054 % 0055 % Revision 1.1 2006/04/14 03:30:49 msh 0056 % Added mne_make_projector 0057 % 0058 % 0059 0060 me='MNE:mne_make_projector'; 0061 0062 if nargin == 2 0063 bads = []; 0064 elseif nargin ~= 3 0065 error(me,'Incorrect number of arguments'); 0066 end 0067 0068 nchan = length(ch_names); 0069 if nchan == 0 0070 error(me,'No channel names specified'); 0071 end 0072 0073 proj = eye(nchan,nchan); 0074 nproj = 0; 0075 U = []; 0076 % 0077 % Check trivial cases first 0078 % 0079 if isempty(projs) 0080 return; 0081 end 0082 0083 nactive = 0; 0084 nvec = 0; 0085 for k = 1:length(projs) 0086 if projs(k).active 0087 nactive = nactive + 1; 0088 nvec = nvec + projs(k).data.nrow; 0089 end 0090 end 0091 0092 if nactive == 0 0093 return; 0094 end 0095 % 0096 % Pick the appropriate entries 0097 % 0098 vecs = zeros(nchan,nvec); 0099 nvec = 0; 0100 nonzero = 0; 0101 for k = 1:length(projs) 0102 if projs(k).active 0103 one = projs(k); 0104 sel = []; 0105 vecsel = []; 0106 if length(one.data.col_names) ~= length(unique(one.data.col_names)) 0107 error(me,'Channel name list in projection item %d contains duplicate items',k); 0108 end 0109 % 0110 % Get the two selection vectors to pick correct elements from 0111 % the projection vectors omitting bad channels 0112 % 0113 p = 0; 0114 for c = 1:nchan 0115 match = strmatch(ch_names{c},one.data.col_names); 0116 if ~isempty(match) && isempty(strmatch(ch_names{c},bads)) 0117 p = p + 1; 0118 sel(p) = c; 0119 vecsel(p) = match(1); 0120 end 0121 end 0122 % 0123 % If there is something to pick, pickit 0124 % 0125 if ~isempty(sel) 0126 for v = 1:one.data.nrow 0127 vecs(sel,nvec+v) = one.data.data(v,vecsel)'; 0128 end 0129 % 0130 % Rescale for more straightforward detection of small singular values 0131 % 0132 for v = 1:one.data.nrow 0133 onesize = sqrt(vecs(:,nvec+v)'*vecs(:,nvec+v)); 0134 if onesize > 0 0135 vecs(:,nvec+v) = vecs(:,nvec+v)/onesize; 0136 nonzero = nonzero + 1; 0137 end 0138 end 0139 nvec = nvec + one.data.nrow; 0140 end 0141 end 0142 end 0143 % 0144 % Check whether all of the vectors are exactly zero 0145 % 0146 if nonzero == 0 0147 return; 0148 end 0149 % 0150 % Reorthogonalize the vectors 0151 % 0152 [U,S,V] = svd(vecs(:,1:nvec),0); 0153 S = diag(S); 0154 % 0155 % Throw away the linearly dependent guys 0156 % 0157 for k = 1:nvec 0158 if S(k)/S(1) < 1e-2 0159 nvec = k; 0160 break; 0161 end 0162 end 0163 U = U(:,1:nvec); 0164 % 0165 % Here is the celebrated result 0166 % 0167 proj = proj - U*U'; 0168 nproj = nvec; 0169 0170 return; 0171 0172 end