mne_ex_evoked_grad_amp

   function [res] = mne_ex_evoked_grad_amp(inname,bmin,bmax,outname)

   Compute the magnitude of the tangential gradient at each
   sensor location using the planar gradiometer data and
   optionally output the result to a fif file.

   inname      The input file name. All average data sets are
               read and processed
   bmin,bmax   Baseline limits in seconds
   outname     Optional output file name


   Function returns the data which was or would have been written
   to the file

0001 function [res] = mne_ex_evoked_grad_amp(inname,bmin,bmax,outname)
0002 %
0003 %   function [res] = mne_ex_evoked_grad_amp(inname,bmin,bmax,outname)
0004 %
0005 %   Compute the magnitude of the tangential gradient at each
0006 %   sensor location using the planar gradiometer data and
0007 %   optionally output the result to a fif file.
0008 %
0009 %   inname      The input file name. All average data sets are
0010 %               read and processed
0011 %   bmin,bmax   Baseline limits in seconds
0012 %   outname     Optional output file name
0013 %
0014 %
0015 %   Function returns the data which was or would have been written
0016 %   to the file
0017 %
0018 
0019 %
0020 %   Author : Matti Hamalainen, MGH Martinos Center
0021 %   License : BSD 3-clause
0022 %
0023 
0024 me='MNE:mne_ex_evoked_grad_amp';
0025 
0026 global FIFF;
0027 if isempty(FIFF)
0028     FIFF = fiff_define_constants();
0029 end
0030 
0031 if nargin == 1
0032     do_baseline = false;
0033 elseif (nargin == 3 || nargin == 4)
0034     do_baseline = true;
0035 else
0036     error(me,'Wrong number of arguments');
0037 end
0038 %
0039 %   Read the data
0040 %
0041 data = fiff_read_evoked_all(inname);
0042 %
0043 %   Figure out the planar gradiometer pairs
0044 %
0045 pairs = zeros(data.info.nchan,2);
0046 npair = 0;
0047 k = 1;
0048 while k < data.info.nchan
0049     %
0050     %   First check the coil types
0051     %
0052     coil1 = data.info.chs(k).coil_type;
0053     coil2 = data.info.chs(k+1).coil_type;
0054     if (coil1 == coil2 && ...
0055             (coil1 == 2 || coil1 == 3012 || coil1 == 3013))
0056         one = data.info.ch_names{k};
0057         two = data.info.ch_names{k+1};
0058         lastone = one(length(one));
0059         lasttwo = two(length(two));
0060         %
0061         %   Then the channel names
0062         %
0063         if (strcmp(one(1:3),'MEG') && strcmp(one(1:3),'MEG'))
0064             if (strcmp(one(1:(length(one)-1)),two(1:(length(two)-1))) && ...
0065                     ((lastone == '2' && lasttwo == '3') || ...
0066                     (lastone == '3' && lasttwo == '2')))
0067                 npair = npair + 1;
0068                 pairs(npair,1) = k;
0069                 pairs(npair,2) = k+1;
0070                 k = k + 1;
0071             end
0072         end
0073     end
0074     k = k + 1;
0075 end
0076 
0077 if npair == 0
0078     error(me,'No planar gradiometers in these data');
0079 end
0080 %
0081 %   Compute the amplitudes
0082 %
0083 fprintf(1,'Computing the amplitudes');
0084 if do_baseline
0085     fprintf(1,' (Baseline = %7.1f ... %7.1f ms)',1000*bmin,1000*bmax);
0086 end
0087 fprintf(1,'...');
0088 for k = 1:length(data.evoked)
0089     epochs = data.evoked(k).epochs;
0090     %
0091     %  Setup baseline limits
0092     %
0093     if do_baseline
0094         b1 = double(data.info.sfreq*bmin - data.evoked(k).first);
0095         b2 = double(data.info.sfreq*bmax - data.evoked(k).first);
0096         if b1 < 1
0097             b1 = 1;
0098         end
0099         if b2 > size(epochs,2)
0100             b2 = size(epochs,2)
0101         end
0102     else
0103         b1 = 1;
0104         b2 = 1;
0105     end
0106     %
0107     %   Go through all pairs
0108     %
0109     for p = 1:npair
0110         one = pairs(p,1);
0111         two = pairs(p,2);
0112         if b2 > b1
0113             base1 = sum(epochs(one,b1:b2))/(b2-b1);
0114             base2 = sum(epochs(two,b1:b2))/(b2-b1);
0115             epochs(one,:) = sqrt((epochs(one,:)-base1).*(epochs(one, ...
0116                 :)-base1)+(epochs(two,:)-base2).*(epochs(two,:)-base2));
0117         else
0118             epochs(one,:) = sqrt(epochs(one,:).*epochs(one, ...
0119                 :)+epochs(two,:).*epochs(two,:));
0120         end
0121     end
0122     data.evoked(k).epochs = epochs;
0123     fprintf(1,'.');
0124 end
0125 fprintf(1,'[done]\n');
0126 %
0127 %   Compose the selection name list
0128 %
0129 pairs = pairs(1:npair,:);
0130 for k = 1:npair
0131     ch_sel_names{k} = data.info.ch_names{pairs(k,1)};
0132 end
0133 %
0134 %   Omit MEG channels but include others
0135 %
0136 for p = 1:data.info.nchan
0137     if (data.info.chs(p).kind ~= FIFF.FIFFV_MEG_CH)
0138         k = k + 1;
0139         ch_sel_names{k} = data.info.ch_names{p};
0140     end
0141 end
0142 %
0143 %   Modify the bad channel list
0144 %
0145 if ~isempty(data.info.bads)
0146     nbad = length(data.info.bads);
0147     for k = 1:npair
0148         one = data.info.ch_names{pairs(k,1)};
0149         two = data.info.ch_names{pairs(k,2)};
0150         %
0151         %   If one channel of the planar gradiometer is marked bad,
0152         %   add the other to the bad channel list
0153         %
0154         if (~isempty(strmatch(one,data.info.bads)) && ...
0155                 isempty(strmatch(two,data.info.bads)))
0156             nbad = nbad + 1;
0157             data.info.bads{nbad} = two;
0158         elseif (isempty(strmatch(one,data.info.bads)) && ...
0159                 ~isempty(strmatch(two,data.info.bads)))
0160             nbad = nbad + 1;
0161             data.info.bads{nbad} = one;
0162         end
0163     end
0164 end
0165 %
0166 %   Do the picking
0167 %
0168 res = fiff_pick_channels_evoked(data,ch_sel_names);
0169 %
0170 %   Optionally write an output file
0171 %
0172 if nargin == 4
0173     fiff_write_evoked(outname,res);
0174     fprintf(1,'Wrote %s\n',outname);
0175 end