plot_meg_data

plot_meg_data(samples, chans, freqs,Lo_Pass,Hi_Pass,
              megfname);
plots meg data and its spectra for a given range of samples and channels. 
artifact editing on the basis of amplitude criterion in pico (10^-12) Tesla
can be used to blank out bad channels or tune editing parameters.

samples = [min_samp max_samp]
chans = [chan1 chan2 chan3]
freqs = [freq_min freq_max]
Lo_Pass = low_pass filter frequency (defaults to 50)
Hi_Pass = hi_pass filter frequency (defaults to 0)
megfname = meg_filename (if not provided a gui comes up)

0001 function status = plot_meg_data(samples, chans, freqs, Lo_Pass,Hi_Pass,megfname);
0002 
0003 %plot_meg_data(samples, chans, freqs,Lo_Pass,Hi_Pass,
0004 %              megfname);
0005 %plots meg data and its spectra for a given range of samples and channels.
0006 %artifact editing on the basis of amplitude criterion in pico (10^-12) Tesla
0007 %can be used to blank out bad channels or tune editing parameters.
0008 %
0009 %samples = [min_samp max_samp]
0010 %chans = [chan1 chan2 chan3]
0011 %freqs = [freq_min freq_max]
0012 %Lo_Pass = low_pass filter frequency (defaults to 50)
0013 %Hi_Pass = hi_pass filter frequency (defaults to 0)
0014 %megfname = meg_filename (if not provided a gui comes up)
0015 %
0016 if nargin < 1
0017     error('samples not specified');
0018 end;
0019 
0020 if nargin == 1
0021     chans = 1;
0022     freqs = [0 50];
0023     Lo_Pass = 50;
0024     Hi_Pass = 1;
0025 end;
0026 if nargin == 2
0027     freqs = [0 50];
0028     Lo_Pass = 50;
0029     Hi_Pass = 1;
0030 end;
0031 if nargin == 3
0032     Lo_Pass = 50;
0033     Hi_Pass = 1;
0034 end;
0035 if nargin == 4
0036     Hi_Pass = 1;
0037 end;
0038 
0039 if nargin < 6
0040     [meg_fid, megfname]  = get_fid('rb');
0041     fclose(meg_fid);
0042 end;
0043 
0044 meg_fid = open_file_w_byte_order(megfname,1);
0045 
0046 if isempty(Hi_Pass);
0047     Hi_Pass= 1;
0048 end;
0049 if isempty(Lo_Pass)
0050     Lo_Pass = 50;
0051 end;
0052 if isempty(freqs);
0053     freqs = [0 40];
0054 end
0055 
0056 if isempty(chans)
0057     chans = 1;
0058 end;
0059 [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);
0060 [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);
0061 trialdata = rd_meg_allch(meg_fid,header_length,NChan,samples(2) - samples(1)+1,samples(1));
0062 trialdata_ord = fix_trial_order(trialdata,meg_channels,meg_indices,eeg_channels,eeg_indices,analog_channels,analog_indices,reference_channels,reference_indices);
0063 trialdata_ord = zeromean(trialdata_ord);
0064 max_trialdata = max(max(abs(trialdata_ord(:,chans))));
0065 max_trialdata = max_trialdata*1.1;
0066 step = 1/Samp_Rate;
0067 time = [step*samples(1):step:step*samples(2)];
0068 [blow,alow] = butter(10,Lo_Pass/(fix(Samp_Rate)/2));
0069 [bhigh,ahigh] = butter(2,[Hi_Pass]/(fix(Samp_Rate)/2),'high');
0070 for i = 1:size(chans,2)
0071     trialdata_ord(:,chans(i)) = filtfilt(blow,alow,trialdata_ord(:,chans(i)));
0072     trialdata_ord(:,chans(i)) = filtfilt(blow,alow,trialdata_ord(:,chans(i)));
0073 end;
0074 figure
0075 for i = 1:size(chans,2)
0076     hold on,plot(time,trialdata_ord(:,chans(i))+(i-1)*max_trialdata*ones(size(trialdata_ord,1),1),'k-')
0077     hold on, plot([min(time)-0.25 max(time)+0.25],[(i-1)*max_trialdata (i-1)*max_trialdata],'k--')
0078     hold on, text(max(time)+0.4,(i-1)*max_trialdata,int2str(chans(i)));
0079 end
0080 total_samp = samples(2)-samples(1)+1;
0081 xlabel('Time (seconds)')
0082 ylabel('PicoTesla')
0083     if ~isempty(freqs)
0084         figure
0085         if freqs(1) == 0
0086             freqs(1) = 0.1;
0087         end;
0088         fft_trial = abs(fft(trialdata_ord(:,chans)))/(total_samp+1);
0089         binmin = ceil(freqs(1)*(total_samp/(fix(Samp_Rate))));
0090         binmax = fix(freqs(2)*(total_samp/(fix(Samp_Rate))));
0091 
0092         max_fft = max(max(fft_trial(binmin:binmax,:)));
0093         frequency = [0:fix(size(fft_trial,1)/2)-1]/(total_samp/fix(Samp_Rate));
0094         for i = 1:size(chans,2)
0095             hold on,plot(frequency,fft_trial(1:size(frequency,2),i)+(i-1)*max_fft*ones(size(frequency,2),1),'k-')
0096             hold on,axis([freqs(1) freqs(2) 0 size(chans,2)*max_fft ])
0097             hold on, plot([min(freqs)-1 max(freqs)*1.1],[(i-1)*max_fft (i-1)*max_fft],'k--')
0098             hold on, text(max(freqs)*0.75,(i-0.5)*max_fft,int2str(chans(i)));
0099                         
0100         end
0101         xlabel('Frequency(Hz)')
0102         ylabel('Amplitude');
0103     end    
0104     
0105     
0106 
0107 
0108 
0109 
0110