phase_scatter

status =phase_scatter(freqs,obs,array,lines,phase_type,phasefname);

freqs = frequencies to plot, e.g. [3 30] limits the plot at 3 and 30 Hz  
obs = observations in the power file to plot (defaults to all)
array = either an array name,e.g., 'all','hemisphere','quadrant','frontback',
    or 'oned' or a list of channels
 Restrictions nobs <= 3. If nobs == 1 and multiple frequencies 
 are listed then freqs < = 3. If nobs > 1 obs go on same plot one for
 each freq.  if nobs == 1 all the freqs go on the same plot
lines = 'black' or 'color'
phase_type = type of phase angle toplot either 'cos' or 'angles'
phasefname = phasefilename (can be gui by skipping this param)

note also draws in solid black phase due to uncorrelated sources
in the four spheres model.  If 'laplacian' is a reference the curve
for laplacians is drawn as well.

0001 function status =phase_scatter(freqs,obs,array,lines,phase_type,phasefname);
0002 %status =phase_scatter(freqs,obs,array,lines,phase_type,phasefname);
0003 %
0004 %freqs = frequencies to plot, e.g. [3 30] limits the plot at 3 and 30 Hz
0005 %obs = observations in the power file to plot (defaults to all)
0006 %array = either an array name,e.g., 'all','hemisphere','quadrant','frontback',
0007 %    or 'oned' or a list of channels
0008 % Restrictions nobs <= 3. If nobs == 1 and multiple frequencies
0009 % are listed then freqs < = 3. If nobs > 1 obs go on same plot one for
0010 % each freq.  if nobs == 1 all the freqs go on the same plot
0011 %lines = 'black' or 'color'
0012 %phase_type = type of phase angle toplot either 'cos' or 'angles'
0013 %phasefname = phasefilename (can be gui by skipping this param)
0014 %
0015 %note also draws in solid black phase due to uncorrelated sources
0016 %in the four spheres model.  If 'laplacian' is a reference the curve
0017 %for laplacians is drawn as well.
0018 if nargin < 6
0019     [fid] = get_fid('rb');
0020 else
0021     fid = fopen(phasefname,'rb');
0022 end;
0023 version = fread(fid,1,'int16');
0024 if version ~= -2
0025     error('this is not a phase file');
0026 end;
0027  [nfiles,obs_labels,Epoch,Window_Length,NEpoch,Nbad_chan,bad_chan,ref_flag,reference,NChan,NFreq] = rd_anal_hdr(fid);
0028 if nargin == 0
0029     freqs = [0:40];
0030     obs = [1:nfiles];
0031     array = 'all';
0032     lines = 'black';
0033     phase_type = 'angles';
0034 end;
0035 if nargin == 1
0036     obs = [1:nfiles];
0037     array = 'all';
0038     lines = 'black';
0039     phase_type = 'angles'
0040 end;        
0041 if nargin == 2
0042     array = 'all';
0043     lines = 'black';
0044     phase_type = 'angles'
0045 end;    
0046 if nargin == 3
0047     lines = 'black';
0048     phase_type = 'angles';
0049 end;
0050 if nargin == 4
0051     phase_type = 'angles';
0052 end;
0053 if isempty(lines)
0054     lines = 'black';
0055 end
0056 if isempty(array)
0057     array = 'all';
0058 end
0059 if isempty(obs)
0060     obs = [1:nfiles];
0061 end;
0062 if isempty(freqs)
0063     freqs = [0:40];
0064 end;
0065 if isempty(phase_type)
0066     phase_type = 'angles';
0067 end;
0068 if size(obs,2) > 3
0069     fclose('all')
0070     error('too many observations > 3 to be plotted in a single plot')
0071 end;
0072 if size(obs,2) == 1
0073     if size(freqs,2) > 3 
0074         error('too many frequencies for a single plot')
0075     end;
0076 end;
0077 ch_pair_indices;
0078 if strcmp(lines,'black')
0079 line_type = ['w- ';'w--';'w-.'];
0080 scat_type = ['wx';'wo';'w.'];
0081 else    
0082 line_type = ['r-';'b-';'w-'];
0083 scat_type = ['rx';'bo';'w.'];
0084 end;
0085 iline = 1;
0086 [xelec,yelec,zelec] = electrodes(129);
0087 twod_dist = twod_pos(xelec,yelec,zelec,9.2);
0088 phase_array = real(get_eeg_coherence(fid,freqs,obs,nfiles,NChan,NFreq,Epoch));
0089 if strcmp(phase_type,'cos')
0090     phase_array = cos(phase_array);
0091     phasemax = 1;
0092 else
0093     phasemax = pi;
0094 end
0095 if isstr(array)
0096     [arrayg,gnames] = arrays(array,[]);
0097 else
0098     arrayg = array;
0099 end
0100 if size(obs,2) > 1
0101     for g = 1:size(arrayg,1)
0102             for i= 1:size(freqs,2)
0103                 figure
0104                 hold on
0105                 iline = 1;
0106                 
0107                 for j = 1:size(obs,2)
0108                 
0109                     if isstr(array)
0110                         [arrayh,hnames] = arrays(array,bad_chan(obs(j),find(bad_chan(obs(j),:))));
0111                     else
0112                         arrayh = array;
0113                     end
0114                 pairels = arrayh(g,find(arrayh(g,:)));
0115                 npairs = size(pairels,2)*(size(pairels,2)+1)/2;
0116                 pairs = zeros(1,npairs);
0117                 ipairs = 1;
0118                 for ip = 1:size(pairels,2)
0119                     for jp=ip:size(pairels,2)
0120                         pairs(ipairs) = chpair(pairels(ip),pairels(jp));
0121                         ipairs = ipairs + 1;
0122                     end;
0123                 end;
0124                 pairnames = hnames(g,:);
0125                 hold on, plot(real(twod_dist(pairs)),phase_array((j-1)*size(freqs,2)+i,pairs),scat_type(iline,:))
0126                 p_coh = polyfit(twod_dist(pairs),phase_array((j-1)*size(freqs,2)+i,pairs),4);
0127                 coh_fit = real(polyval(p_coh,[0:1:fix(max(twod_dist(pairs)))]));
0128                 hold on, plot([0:1:max(twod_dist(pairs))],coh_fit,line_type(iline,:)), axis([0 max(twod_dist(pairs))+1 -phasemax*1.1 phasemax*1.1])
0129                 title(['Phase versus Distance:  @' int2str(freqs(i)) ' Hz ' hnames(g,:)]);
0130                 text([max(twod_dist(pairs))-2.5],[1-0.1*j],[scat_type(iline,2) ' ' deblank(setstr(obs_labels(obs(j),:)))])
0131                 xlabel('Interelectrode Distance (cm)')
0132                 ylabel('Phase')
0133                  iline = iline+ 1;
0134             end
0135         end
0136     end
0137 else
0138     if isstr(array)
0139         [arrayh,hnames] = arrays(array,bad_chan(obs,find(bad_chan(obs,:))));
0140     else
0141         arrayh = array;
0142     end
0143     for g = 1:size(arrayh,1)
0144         iline = 1;
0145         figure
0146         hold on
0147         
0148         for j = 1:size(freqs,2)
0149             pairels = arrayh(g,find(arrayh(g,:)));
0150             
0151             npairs = size(pairels,2)*(size(pairels,2)+1)/2;
0152             pairs = zeros(1,npairs);
0153             ipairs = 1;
0154             for ip = 1:size(pairels,2)
0155                 for jp=ip:size(pairels,2)
0156                     pairs(ipairs) = chpair(pairels(ip),pairels(jp));
0157                     ipairs = ipairs + 1;
0158                 end;
0159             end;
0160             chpair(pairels(1),pairels(2))
0161             pairnames = hnames(g,:);
0162             
0163             hold on, plot(real(twod_dist(pairs)),phase_array(j,pairs),scat_type(iline,:))
0164             p_coh = polyfit(real(twod_dist(pairs)),phase_array(j,pairs),4);
0165             coh_fit = real(polyval(p_coh,[0:1:fix(max(twod_dist(pairs)))]));
0166             hold on, plot([0:1:max(twod_dist(pairs))],coh_fit,line_type(iline,:)), axis([0 max(twod_dist(pairs))+1 -phasemax*1.1 phasemax*1.1])
0167             title(['Phase versus Distance: ' deblank(setstr(obs_labels(obs(1),:))) ' ' hnames(g,:)]);
0168             text([max(twod_dist(pairs))-2.5],[1-0.1*j],[scat_type(iline,2) ' ' int2str(freqs(j)) ' Hz'])
0169             xlabel('Interelectrode Distance (cm)')
0170             ylabel('Phase')
0171             iline = iline+ 1;
0172         end
0173     end
0174 end
0175 status = 1;
0176 
0177 
0178 
0179 
0180 
0181 
0182 
0183 
0184 
0185 
0186 
0187 
0188 
0189 
0190 
0191 
0192 
0193 
0194 
0195 
0196 
0197 
0198 
0199 
0200 
0201 
0202 
0203 
0204