run_riv_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_riv_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 csdmfilename1 = [megfilename(1:size(megfilename,2) - 4) '.1.csdm'];
0099 csdmfilename2 = [megfilename(1:size(megfilename,2) - 4) '.2.csdm'];
0100 csdmfilenamer1 = [megfilename(1:size(megfilename,2) - 4) '.r1.csdm'];
0101 csdmfilenamer2 = [megfilename(1:size(megfilename,2) - 4) '.r2.csdm'];
0102 csdm_fid1 = fopen(csdmfilename1,'wb');
0103 csdm_fid2 = fopen(csdmfilename2,'wb');
0104 csdm_fidr1 = fopen(csdmfilenamer1,'wb');
0105 csdm_fidr2 = fopen(csdmfilenamer2,'wb');
0106 add_bad_chan = [add_bad_chan bad_sensors];
0107 csdmversion = -1;
0108 csdm_NMeg = 148;
0109 %
0110 % note: this has to be generalize dwitht he appropriate routine.
0111 %    to allow for bipolar pairs
0112 csdm_NEeg = NEeg;
0113 csdm_NReference = NReference;
0114 csdm_NAnalog = NAnalog;
0115 
0116 csdm_NChan = csdm_NMeg + csdm_NEeg + csdm_NReference + csdm_NAnalog;
0117 
0118 if csdm_NChan ~= max(analog_channels);
0119     error('channel identification mismatch');
0120 end;
0121 stack1 = zeros(1,csdm_NChan);
0122 stack2 = zeros(1,csdm_NChan);
0123 stackr1 = zeros(1,csdm_NChan);
0124 stackr2 = zeros(1,csdm_NChan);
0125 csdm_NChan_file = (csdm_NChan*(csdm_NChan+1))/2;
0126 %write out csdm header version
0127 fwrite(csdm_fid1,csdmversion,'integer*2');
0128 fwrite(csdm_fid2,csdmversion,'integer*2');
0129 fwrite(csdm_fidr1,csdmversion,'integer*2');
0130 fwrite(csdm_fidr2,csdmversion,'integer*2');
0131 avgcsdm1 = zeros(numsamp/4,csdm_NChan_file);
0132 avgcsdm2 = zeros(numsamp/4,csdm_NChan_file);
0133 avgcsdmr1 = zeros(numsamp/4,csdm_NChan_file);
0134 avgcsdmr2 = zeros(numsamp/4,csdm_NChan_file);
0135 rand_wt1 = rand(1,NWindow);
0136 rand_wt2 = ones(1,NWindow) - rand_wt1;
0137 for i = 1:NWindow
0138     if i == 1
0139         trialdata = rd_meg_allch(meg_fid,header_length,NChan,numsamp,skip_samples+1);
0140     else
0141         trialdata = rd_meg_allch(meg_fid,header_length,NChan,numsamp);
0142     end;
0143     trialdata2 = trialdata(1:2:numsamp-1,:);
0144     trialdata_ord = fix_trial_order(trialdata2,meg_channels,meg_indices,eeg_channels,eeg_indices,analog_channels,analog_indices,reference_channels,reference_indices);
0145     mask = artifact_edit_meg(trialdata_ord,meg_channels,eeg_channels,MEG_max,EEG_max,add_bad_chan);
0146     if sum(mask) > 100
0147         if ~isempty(eeg_channels)
0148         trialdata_ord = reference_eeg_for_meg(trialdata_ord,eeg_channels,mask,EEG_reference);
0149         end
0150         wt1 = size(find(response(skip_samples+(i-1)*numsamp+1:skip_samples+i*numsamp)== 1),2)/numsamp;
0151         wt2 = size(find(response(skip_samples+(i-1)*numsamp+1:skip_samples+i*numsamp)== 2),2)/numsamp;
0152         csdm_trial = csdm(trialdata_ord);
0153         avgcsdm1 = avgcsdm1 + wt1*csdm_trial;
0154         avgcsdm2 = avgcsdm2 + wt2*csdm_trial;
0155         avgcsdmr1 = avgcsdmr1 + rand_wt1*csdm_trial;
0156         avgcsdmr2 = avgcsdmr2 + rand_wt2*csdm_trial;
0157         stack1 = stack1 + wt1*mask;
0158         stack2 = stack2 + wt2*mask;
0159         stackr1 = stackr1 + rand_wt1*mask;
0160         stackr2 = stackr2 + rand_wt2*mask;
0161 end;
0162 end; 
0163 icount = 1;
0164 good_cross1 = zeros(1,size(avgcsdm1,2));
0165 good_cross2 = zeros(1,size(avgcsdm2,2));
0166 good_crossr1 = zeros(1,size(avgcsdm1,2));
0167 good_crossr2 = zeros(1,size(avgcsdm1,2));
0168 for ichan = 1:csdm_NChan
0169     for jchan = ichan:csdm_NChan
0170         good_cross1(icount) = min([stack1(ichan) stack1(jchan)]);
0171         if good_cross1(icount) == 0
0172             good_cross1(icount) = 1;
0173         end;
0174         good_cross2(icount) = min([stack2(ichan) stack2(jchan)]);
0175         if good_cross2(icount) == 0
0176             good_cross2(icount) = 1;
0177         end;
0178         good_crossr1(icount) = min([stackr1(ichan) stackr1(jchan)]);
0179         if good_crossr1(icount) == 0
0180             good_crossr1(icount) = 1;
0181         end;
0182         good_crossr2(icount) = min([stackr2(ichan) stackr2(jchan)]);
0183         if good_crossr2(icount) == 0
0184             good_crossr2(icount) = 1;
0185         end;
0186         icount = icount+1;
0187     end;
0188 end;
0189 avgcsdm1 = avgcsdm1./(ones(size(avgcsdm1,1),1)*good_cross1);
0190 avgcsdm2 = avgcsdm2./(ones(size(avgcsdm1,1),1)*good_cross2);
0191 avgcsdmr1 = avgcsdmr1./(ones(size(avgcsdm1,1),1)*good_crossr1);
0192 avgcsdmr2 = avgcsdmr2./(ones(size(avgcsdm1,1),1)*good_crossr2);
0193 
0194 NEpoch_used1 = fix(max(good_cross1)); 
0195 NEpoch_used2 = fix(max(good_cross2)); 
0196 NEpoch_usedr1 = fix(max(good_crossr1)); 
0197 NEpoch_usedr2 = fix(max(good_crossr2)); 
0198 
0199 
0200 more_bad_chan1 = find(stack1 < 0.75*NEpoch_used1);
0201 more_bad_chan2 = find(stack2 < 0.75*NEpoch_used2);
0202 more_bad_chanr1 = find(stackr1 < 0.75*NEpoch_usedr1);
0203 more_bad_chanr2 = find(stackr2 < 0.75*NEpoch_usedr2);
0204 
0205 add_bad_chan1 = [add_bad_chan more_bad_chan1];
0206 add_bad_chan2 = [add_bad_chan more_bad_chan2];
0207 add_bad_chanr1 = [add_bad_chan more_bad_chanr1];
0208 add_bad_chanr2 = [add_bad_chan more_bad_chanr2];
0209 
0210 fwrite(csdm_fid1,csdm_NChan_file,'integer*2');
0211 fwrite(csdm_fid1,csdm_NMeg,'integer*2');
0212 fwrite(csdm_fid1,csdm_NReference,'integer*2');
0213 fwrite(csdm_fid1,csdm_NEeg,'integer*2');
0214 fwrite(csdm_fid1,csdm_NAnalog,'integer*2');
0215 NBad_chan = size(add_bad_chan1,2);
0216 fwrite(csdm_fid1,NBad_chan,'integer*2');
0217 fwrite(csdm_fid1,add_bad_chan1,'integer*2');
0218 disp('the total set of bad channels were:')
0219 add_bad_chan1
0220 fwrite(csdm_fid1,NEpoch_used1,'integer*2');
0221 fwrite(csdm_fid1,Epoch,'real*4');
0222 max_freq = fix(max_freq*Epoch)+1;
0223 fwrite(csdm_fid1,max_freq,'integer*2');
0224 fwrite(csdm_fid1,real(avgcsdm1(1:max_freq,:))','real*4');
0225 fwrite(csdm_fid1,imag(avgcsdm1(1:max_freq,:))','real*4');
0226 
0227 fwrite(csdm_fid2,csdm_NChan_file,'integer*2');
0228 fwrite(csdm_fid2,csdm_NMeg,'integer*2');
0229 fwrite(csdm_fid2,csdm_NReference,'integer*2');
0230 fwrite(csdm_fid2,csdm_NEeg,'integer*2');
0231 fwrite(csdm_fid2,csdm_NAnalog,'integer*2');
0232 NBad_chan = size(add_bad_chan2,2);
0233 fwrite(csdm_fid2,NBad_chan,'integer*2');
0234 fwrite(csdm_fid2,add_bad_chan2,'integer*2');
0235 disp('the total set of bad channels were:')
0236 add_bad_chan2
0237 fwrite(csdm_fid2,NEpoch_used2,'integer*2');
0238 fwrite(csdm_fid2,Epoch,'real*4');
0239 max_freq = fix(max_freq*Epoch)+1;
0240 fwrite(csdm_fid2,max_freq,'integer*2');
0241 fwrite(csdm_fid2,real(avgcsdm2(1:max_freq,:))','real*4');
0242 fwrite(csdm_fid2,imag(avgcsdm2(1:max_freq,:))','real*4');
0243 
0244 fwrite(csdm_fidr1,csdm_NChan_file,'integer*2');
0245 fwrite(csdm_fidr1,csdm_NMeg,'integer*2');
0246 fwrite(csdm_fidr1,csdm_NReference,'integer*2');
0247 fwrite(csdm_fidr1,csdm_NEeg,'integer*2');
0248 fwrite(csdm_fidr1,csdm_NAnalog,'integer*2');
0249 NBad_chan = size(add_bad_chanr1,2);
0250 fwrite(csdm_fidr1,NBad_chan,'integer*2');
0251 fwrite(csdm_fidr1,add_bad_chanr1,'integer*2');
0252 disp('the total set of bad channels were:')
0253 add_bad_chanr1
0254 fwrite(csdm_fidr1,NEpoch_usedr1,'integer*2');
0255 fwrite(csdm_fidr1,Epoch,'real*4');
0256 max_freq = fix(max_freq*Epoch)+1;
0257 fwrite(csdm_fidr1,max_freq,'integer*2');
0258 fwrite(csdm_fidr1,real(avgcsdmr1(1:max_freq,:))','real*4');
0259 fwrite(csdm_fidr1,imag(avgcsdmr1(1:max_freq,:))','real*4');
0260 
0261 fwrite(csdm_fidr2,csdm_NChan_file,'integer*2');
0262 fwrite(csdm_fidr2,csdm_NMeg,'integer*2');
0263 fwrite(csdm_fidr2,csdm_NReference,'integer*2');
0264 fwrite(csdm_fidr2,csdm_NEeg,'integer*2');
0265 fwrite(csdm_fidr2,csdm_NAnalog,'integer*2');
0266 NBad_chan = size(add_bad_chanr2,2);
0267 fwrite(csdm_fidr2,NBad_chan,'integer*2');
0268 fwrite(csdm_fidr2,add_bad_chanr2,'integer*2');
0269 disp('the total set of bad channels were:')
0270 add_bad_chanr2
0271 fwrite(csdm_fidr2,NEpoch_usedr2,'integer*2');
0272 fwrite(csdm_fidr2,Epoch,'real*4');
0273 max_freq = fix(max_freq*Epoch)+1;
0274 fwrite(csdm_fidr2,max_freq,'integer*2');
0275 fwrite(csdm_fidr2,real(avgcsdmr2(1:max_freq,:))','real*4');
0276 fwrite(csdm_fidr2,imag(avgcsdmr2(1:max_freq,:))','real*4');
0277 fclose('all');
0278 status = 1;
0279 
0280 
0281 
0282 
0283 
0284 
0285 
0286 
0287 
0288