mne_make_compensator

 [comp] = mne_make_compensator(info,from,to,exclude_comp_chs)

 info              - measurement info as returned by the fif reading routines
 from              - compensation in the input data
 to                - desired compensation in the output
 exclude_comp_chs  - exclude compensation channels from the output (optional)

0001 function [comp] = mne_make_compensator(info,from,to,exclude_comp_chs)
0002 %
0003 % [comp] = mne_make_compensator(info,from,to,exclude_comp_chs)
0004 %
0005 % info              - measurement info as returned by the fif reading routines
0006 % from              - compensation in the input data
0007 % to                - desired compensation in the output
0008 % exclude_comp_chs  - exclude compensation channels from the output (optional)
0009 %
0010 
0011 %
0012 % Create a compensation matrix to bring the data from one compensation
0013 % state to another
0014 %
0015 
0016 
0017 %
0018 %   Author : Matti Hamalainen, MGH Martinos Center
0019 %   License : BSD 3-clause
0020 %
0021 %
0022 %   Revision 1.8  2006/05/03 19:09:03  msh
0023 %   Fixed even more compatibility issues.
0024 %
0025 %   Revision 1.7  2006/04/23 15:29:40  msh
0026 %   Added MGH to the copyright
0027 %
0028 %   Revision 1.6  2006/04/18 21:11:31  msh
0029 %   Added option to omit compensation channels from the compensator
0030 %
0031 %   Revision 1.5  2006/04/18 20:44:46  msh
0032 %   Added reading of forward solution.
0033 %   Use length instead of size when appropriate
0034 %
0035 %   Revision 1.4  2006/04/17 11:52:15  msh
0036 %   Added coil definition stuff
0037 %
0038 %   Revision 1.3  2006/04/15 12:21:00  msh
0039 %   Several small improvements
0040 %
0041 %   Revision 1.2  2006/04/14 15:49:49  msh
0042 %   Improved the channel selection code and added ch_names to measurement info.
0043 %
0044 %   Revision 1.1  2006/04/12 10:51:19  msh
0045 %   Added projection writing and compensation routines
0046 %
0047 %
0048 
0049 me='MNE:mne_make_compensator';
0050 
0051 global FIFF;
0052 if isempty(FIFF)
0053     FIFF = fiff_define_constants();
0054 end
0055 
0056 if nargin == 3
0057     exclude_comp_chs = false;
0058 elseif nargin ~= 4
0059     error(me,'Incorrect number of arguments');
0060 end
0061 
0062 if from == to
0063     comp = zeros(info.nchan,info.nchan);
0064     return;
0065 end
0066 
0067 if from == 0
0068     C1 = zeros(info.nchan,info.nchan);
0069 else
0070     try
0071         C1 = make_compensator(info,from);
0072     catch
0073         error(me,'Cannot create compensator C1 (%s)',mne_omit_first_line(lasterr));
0074     end
0075     if isempty(C1)
0076         error(me,'Desired compensation matrix (kind = %d) not found',from);
0077     end
0078 end
0079 if to == 0
0080     C2 = zeros(info.nchan,info.nchan);
0081 else
0082     try
0083         C2 = make_compensator(info,to);
0084     catch
0085         error(me,'Cannot create compensator C2 (%s)',mne_omit_first_line(lasterr));
0086     end
0087     if isempty(C2)
0088         error(me,'Desired compensation matrix (kind = %d) not found',to);
0089     end
0090 end
0091 %
0092 %   s_orig = s_from + C1*s_from = (I + C1)*s_from
0093 %   s_to   = s_orig - C2*s_orig = (I - C2)*s_orig
0094 %   s_to   = (I - C2)*(I + C1)*s_from = (I + C1 - C2 - C2*C1)*s_from
0095 %
0096 comp = eye(info.nchan,info.nchan) + C1 - C2 - C2*C1;
0097 
0098 if exclude_comp_chs
0099     pick  = zeros(info.nchan, 1);
0100     npick = 0;
0101     for k = 1:info.nchan
0102         if info.chs(k).kind ~= FIFF.FIFFV_REF_MEG_CH
0103             npick = npick + 1;
0104             pick(npick) = k;
0105         end
0106     end
0107     if npick == 0
0108         error(me,'Nothing remains after excluding the compensation channels');
0109     end
0110     comp = comp(pick(1:npick),:);
0111 end
0112 
0113 return;
0114 
0115     function this_comp = make_compensator(info,kind)
0116 
0117         for k = 1:length(info.comps)
0118             if info.comps(k).kind == kind
0119                 this_data = info.comps(k).data;
0120                 %
0121                 %   Create the preselector
0122                 %
0123                 presel  = zeros(this_data.ncol,info.nchan);
0124                 for col = 1:this_data.ncol
0125                     c = strmatch(this_data.col_names{col},info.ch_names,'exact');
0126                     if isempty(c)
0127                         error(me,'Channel %s is not available in data',this_data.col_names{col});
0128                     elseif length(c) > 1
0129                         error(me,'Ambiguous channel %s',mat.col_names{col});
0130                     end
0131                     presel(col,c) = 1.0;
0132                 end
0133                 %
0134                 %   Create the postselector
0135                 %
0136                 postsel = zeros(info.nchan,this_data.nrow);
0137                 for c = 1:info.nchan
0138                     row = strmatch(info.ch_names{c},this_data.row_names,'exact');
0139                     if length(row) > 1
0140                         error(me,'Ambiguous channel %s', info.ch_names{c});
0141                     elseif length(row) == 1
0142                         postsel(c,row) = 1.0;
0143                     end
0144                 end
0145                 this_comp = postsel*this_data.data*presel;
0146                 return;
0147             end
0148         end
0149         this_comp = [];
0150         return;
0151     end
0152 
0153 end