run_meg_csdm

status = run_meg_csdm(skip_samples,Epoch,MEG_max,EEG_max,EEG_reference,
    add_bad_chan,max_freq,meg_filename);

skip_samples = # of samples from beginning of file to average 
Epoch = epoch length in seconds to use
MEG_max = maximum value (in picotesla) to keep an MEG channel
EEG_max = maximum value (in microvolts) to keep an EEG channel
EEG_reference = reference to use for EEG (defaults to average mastoids)
add_bad_chan = bad_channels to ignore
max_freq = highest frequency to write out 
megfilename = meg_filename (must end in .MEG or .meg) gui o.k. by skipping
 channels are sequenced in the following absolute channel order:
 1 to 148 :: MEG
 149 to 149 + NReference :: reference channels  
 149 +NReference +1 to 149 +NReference + NEeg :: EEG signals 
    Note:  this is not necessarily the same in the csdm file 
        it depends on the EEG reference selected.  if averag
        mastoids is selected this is true but if a bipolar 
        array is selected you mayhave many more channels than the 
        original data file.  The csdm file will however carry 
        an updated NEeg
 149 +NReference +NEeg +1 to  149 + NReference + NEeg +NAnalog :: any
        analog channels

 Downstream analysis and visualization routines require the absolute channel
 numbers rather than the original channel numbers.

0001 function status = run_meg_csdm(skip_samples,Epoch,MEG_max,EEG_max,EEG_reference,add_bad_chan,max_freq,megfilename);
0002 %status = run_meg_csdm(skip_samples,Epoch,MEG_max,EEG_max,EEG_reference,
0003 %    add_bad_chan,max_freq,meg_filename);
0004 %
0005 %skip_samples = # of samples from beginning of file to average
0006 %Epoch = epoch length in seconds to use
0007 %MEG_max = maximum value (in picotesla) to keep an MEG channel
0008 %EEG_max = maximum value (in microvolts) to keep an EEG channel
0009 %EEG_reference = reference to use for EEG (defaults to average mastoids)
0010 %add_bad_chan = bad_channels to ignore
0011 %max_freq = highest frequency to write out
0012 %megfilename = meg_filename (must end in .MEG or .meg) gui o.k. by skipping
0013 % channels are sequenced in the following absolute channel order:
0014 % 1 to 148 :: MEG
0015 % 149 to 149 + NReference :: reference channels
0016 % 149 +NReference +1 to 149 +NReference + NEeg :: EEG signals
0017 %    Note:  this is not necessarily the same in the csdm file
0018 %        it depends on the EEG reference selected.  if averag
0019 %        mastoids is selected this is true but if a bipolar
0020 %        array is selected you mayhave many more channels than the
0021 %        original data file.  The csdm file will however carry
0022 %        an updated NEeg
0023 % 149 +NReference +NEeg +1 to  149 + NReference + NEeg +NAnalog :: any
0024 %        analog channels
0025 %
0026 % Downstream analysis and visualization routines require the absolute channel
0027 % numbers rather than the original channel numbers.
0028 
0029 
0030 if nargin < 2
0031     error('insufficient arguments');
0032 end;
0033 
0034 if nargin == 2
0035     MEG_max =  1;
0036     EEG_max = 1;
0037     EEG_reference = 'mastoids';
0038     add_bad_chan = [];
0039     max_freq = 50;
0040 end;
0041 
0042 if nargin == 3
0043     EEG_max = 1;
0044     EEG_reference = 'mastoids';
0045     add_bad_chan = [];
0046     max_freq = 50;
0047 end;
0048 
0049 if nargin == 4
0050     EEG_reference = 'mastoids'
0051     add_bad_chan = [];
0052     max_freq = 50;
0053 end;
0054 
0055 if nargin == 5
0056     add_bad_chan = [];
0057     max_freq = 50;    
0058 end;
0059 
0060 if nargin == 6
0061     max_freq = 50;
0062 end;
0063 
0064 if isempty(MEG_max)
0065     MEG_max = 1;
0066 end
0067 
0068 if isempty(EEG_max)
0069     EEG_max = 1;
0070 end;
0071 
0072 if isempty(EEG_reference)
0073     EEG_reference = 'mastoids';
0074 end;
0075 
0076 if isempty(add_bad_chan)
0077     add_bad_chan = [];
0078 end;
0079 
0080 if isempty(max_freq)
0081     max_freq = 50;
0082 end;
0083 
0084 if nargin < 8
0085     [fid,megfilename] = get_fid('rb');
0086     fclose(fid);
0087 end;
0088 
0089 meg_fid = open_file_w_byte_order(megfilename,1);
0090 
0091 [version, NChan, NMeg, NEeg, NReference, NBad_Sensors, NTrigger, NResponse, NUtility, NAnalog, Samp_Rate, NData_Epoch, Names, NSamp, trigger, response, header_length] = rd_meg_hdr(meg_fid);
0092 
0093 [meg_indices,meg_channels, bad_sensors,eeg_indices,eeg_channels,analog_indices,analog_channels,reference_indices,reference_channels] = parse_names(Names,NMeg,NEeg,NBad_Sensors,NAnalog,NReference);
0094 numsamp = 4*fix(Epoch*Samp_Rate/4);
0095 Epoch = numsamp/Samp_Rate;
0096 NWindow = fix((NSamp-skip_samples)/numsamp);
0097 
0098 csdmfilename = [megfilename(1:size(megfilename,2) - 4) '.csdm'];
0099 
0100 csdm_fid = fopen(csdmfilename,'wb');
0101 add_bad_chan = [add_bad_chan bad_sensors];
0102 csdmversion = -1;
0103 csdm_NMeg = 148;
0104 %
0105 % note: this has to be generalize dwitht he appropriate routine.
0106 %    to allow for bipolar pairs
0107 csdm_NEeg = NEeg;
0108 csdm_NReference = NReference;
0109 csdm_NAnalog = NAnalog;
0110 
0111 csdm_NChan = csdm_NMeg + csdm_NEeg + csdm_NReference + csdm_NAnalog;
0112 
0113 if csdm_NChan ~= max(analog_channels);
0114     error('channel identification mismatch');
0115 end;
0116 stack = zeros(1,csdm_NChan);
0117 csdm_NChan_file = (csdm_NChan*(csdm_NChan+1))/2;
0118 %write out csdm header version
0119 
0120 fwrite(csdm_fid,csdmversion,'integer*2');
0121 avgcsdm = zeros(numsamp/4,csdm_NChan_file);
0122 for i = 1:NWindow
0123     if i == 1
0124         trialdata = rd_meg_allch(meg_fid,header_length,NChan,numsamp,skip_samples+1);
0125     else
0126         trialdata = rd_meg_allch(meg_fid,header_length,NChan,numsamp);
0127     end;
0128     trialdata2 = trialdata(1:2:numsamp-1,:);
0129     trialdata_ord = fix_trial_order(trialdata2,meg_channels,meg_indices,eeg_channels,eeg_indices,analog_channels,analog_indices,reference_channels,reference_indices);
0130     mask = artifact_edit_meg(trialdata_ord,meg_channels,eeg_channels,MEG_max,EEG_max,add_bad_chan);
0131     if sum(mask) > 100
0132         if ~isempty(eeg_channels)
0133         trialdata_ord = reference_eeg_for_meg(trialdata_ord,eeg_channels,mask,EEG_reference);
0134         end
0135         avgcsdm = avgcsdm + csdm(trialdata_ord);
0136         stack = stack + mask;
0137         
0138     end;
0139 end; 
0140 icount = 1;
0141 good_cross = zeros(1,size(avgcsdm,2));
0142 for ichan = 1:csdm_NChan
0143     for jchan = ichan:csdm_NChan
0144         good_cross(icount) = min([stack(ichan) stack(jchan)]);
0145         if good_cross(icount) == 0
0146             good_cross(icount) = 1;
0147         end;
0148         icount = icount+1;
0149     end;
0150 end;
0151 avgcsdm = avgcsdm./(ones(size(avgcsdm,1),1)*good_cross);
0152 
0153 NEpoch_used = max(good_cross); 
0154 
0155 more_bad_chan = find(stack < 0.75*NEpoch_used);
0156 
0157 add_bad_chan = [add_bad_chan more_bad_chan];
0158 
0159 fwrite(csdm_fid,csdm_NChan_file,'integer*2');
0160 fwrite(csdm_fid,csdm_NMeg,'integer*2');
0161 fwrite(csdm_fid,csdm_NReference,'integer*2');
0162 fwrite(csdm_fid,csdm_NEeg,'integer*2');
0163 fwrite(csdm_fid,csdm_NAnalog,'integer*2');
0164 NBad_chan = size(add_bad_chan,2);
0165 fwrite(csdm_fid,NBad_chan,'integer*2');
0166 fwrite(csdm_fid,add_bad_chan,'integer*2');
0167 disp('the total set of bad channels were:')
0168 add_bad_chan
0169 fwrite(csdm_fid,NEpoch_used,'integer*2');
0170 fwrite(csdm_fid,Epoch,'real*4');
0171 max_freq = fix(max_freq*Epoch)+1;
0172 fwrite(csdm_fid,max_freq,'integer*2');
0173 fwrite(csdm_fid,real(avgcsdm(1:max_freq,:))','real*4');
0174 fwrite(csdm_fid,imag(avgcsdm(1:max_freq,:))','real*4');
0175 fclose('all');
0176 status = 1;
0177 
0178 
0179 
0180 
0181 
0182 
0183 
0184 
0185