plot_eeg_power

status = plot_eeg_power(freqs,obs,chans,lines,power_type,powerfname);

freqs = frequencies to plot, e.g. [3 30] limits the plot to 3 to 30 Hz  
obs = observations in the power file to plot (defaults to all)
chans = channels to plot (defaults to 10_20);
 nobs <= 5. 
lines = 'black' or 'color'
power_type = can be either 'power', 'amplitude', or 'logpower'
powerfname = powerfilename (can be gui by skipping or passing a blank)

0001 function status = plot_eeg_power(freqs,obs,chans,lines,power_type,powerfname);
0002 %status = plot_eeg_power(freqs,obs,chans,lines,power_type,powerfname);
0003 %
0004 %freqs = frequencies to plot, e.g. [3 30] limits the plot to 3 to 30 Hz
0005 %obs = observations in the power file to plot (defaults to all)
0006 %chans = channels to plot (defaults to 10_20);
0007 % nobs <= 5.
0008 %lines = 'black' or 'color'
0009 %power_type = can be either 'power', 'amplitude', or 'logpower'
0010 %powerfname = powerfilename (can be gui by skipping or passing a blank)
0011 %
0012 if nargin < 6
0013     [fid] = get_fid('rb');
0014 else
0015     fid = fopen(powerfname,'rb');
0016 end;
0017 version = fread(fid,1,'int16');
0018 if version ~= -3
0019     error('this is not a power file');
0020 end;
0021  [nfiles,obs_labels,Epoch,Window_Length,NEpoch,Nbad_chan,bad_chan,ref_flag,reference,NChan,NFreq] = rd_anal_hdr(fid);
0022 if nargin == 0
0023     freqs = [0 40];
0024     obs = [1:nfiles];
0025     chans = [9 11 23 25 34 37 46 58 60 62 71 84 86 97 105 109 122 124 129];
0026     power_type = 'amplitude';
0027     lines = 'black';
0028 end;
0029 if nargin == 1
0030     obs = [1:nfiles];
0031     chans = [9 11 23 25 34 37 46 58 60 62 71 84 86 97 105 109 122 124 129];
0032     power_type = 'amplitude';
0033     lines = 'black';
0034 end;        
0035 if nargin == 2
0036     chans = [9 11 23 25 34 37 46 58 60 62 71 84 86 97 105 109 122 124 129];
0037     power_type = 'amplitude';
0038     lines = 'black';
0039 end;    
0040 if nargin == 3
0041     power_type = 'amplitude';
0042     lines = 'black';
0043 end;
0044 if nargin == 4
0045     power_type = 'amplitude';
0046 end;
0047 if isempty(lines)
0048     lines = 'black';
0049 end
0050 if isempty(power_type)
0051     power_type = 'amplitude';
0052 end;
0053 if isempty(chans)
0054     
0055     chans = [9 11 23 25 34 37 46 58 60 62 71 84 86 97 105 109 122 124 129];
0056 end
0057 if isempty(obs)
0058     obs = [1:nfiles];
0059 end;
0060 if isempty(freqs)
0061     freqs = [0 40];
0062 end;
0063 if size(obs,2) > 5
0064     fclose('all')
0065     error('too many observations > 5 to be plotted in a single plot')
0066 end;
0067 if size(chans,2) > 1
0068 layout = ceil(sqrt(size(chans,2)+1));
0069 else
0070 layout =1;
0071 end
0072 if strcmp(lines,'black')
0073 line_type = ['k- ';'k--';'k-.';'kx ';'ko '];
0074 else    
0075 line_type = ['k-';'r-';'b-';'g-';'m-';'c-'];
0076 end;
0077 figure
0078 powermax = 0;
0079 obs_mask = zeros(1,nfiles);
0080 obs_mask(obs) = ones(1,size(obs,2));
0081 maxpower = 0;
0082 iline = 1;
0083 for i = 1:nfiles
0084     power = fread(fid,[NChan(i),NFreq(i)],'float');
0085     if obs_mask(i) > 0
0086     power = power';
0087     if strcmp(power_type,'amplitude')
0088         power = sqrt(power);
0089     end;
0090     binmin = freqs(1)*Epoch(i) + 1;
0091     binmax = freqs(2)*Epoch(i) + 1;
0092     freqstep = 1/Epoch(i);
0093     freq = [freqs(1):freqstep:freqs(2)];
0094     for ic = 1:size(chans,2) 
0095         if strcmp(power_type,'logpower')
0096             hold on,subplot(layout,layout,ic),  semilogy(freq,power(binmin:binmax,chans(ic)),line_type(iline,:))
0097             if layout < 3
0098             title(['Channel ' int2str(chans(ic))])
0099             xlabel('Frequency(Hz)')
0100             ylabel('Power')
0101             end;
0102         else
0103             hold on, subplot(layout,layout,ic), plot(freq,power(binmin:binmax,chans(ic)),line_type(iline,:))
0104             if layout < 3
0105                 title(['Channel ' int2str(chans(ic))])
0106                 xlabel('Frequency(Hz)')
0107                 if strcmp(power_type,'amplitude')
0108                     ylabel('Amplitude');
0109                 else
0110                     ylabel('Power')
0111                 end                
0112             end;
0113         end    
0114     end
0115     mpower = max(max(power(binmin:binmax,chans)));
0116     maxpower = max([maxpower mpower]);
0117     iline = iline + 1;
0118     end;
0119 end    
0120 for ic = 1:size(chans,2)
0121     subplot(layout,layout,ic), axis([0 freqs(2) 0 1.05*maxpower]);
0122 end;
0123 if layout < 2
0124     for i = 1:size(obs,2)
0125         hold on, plot([0.6*freqs(2) 0.7*freqs(2)],[(1-0.1*i)*maxpower  (1-0.1*i)*maxpower],line_type(i,:))
0126         text(0.75*freqs(2),(1-0.1*i)*maxpower,deblank(setstr(obs_labels(obs(i),:))));
0127     end;
0128 end
0129 if layout > 2
0130     for i = 1:size(chans,2)
0131         hold on, subplot(layout,layout,i), text(0.8*freqs(2),0.8*maxpower,int2str(chans(i)));
0132     end;
0133 end;
0134 if layout > 1
0135     for i = 1:size(obs,2)
0136         hold on, subplot(layout,layout,layout*layout),plot([0.2 1],[1 1]*(1.1-0.2*i),line_type(i,:)), axis([0 1.5 0 1])
0137         hold on, subplot(layout,layout,layout*layout),text(1.2,(1.1-0.2*i),deblank(setstr(obs_labels(obs(i),:))))
0138         hold on, subplot(layout,layout,layout*layout),axis('off')
0139     end;
0140     hold on, subplot(layout,layout,layout*layout),axis('off');
0141 end;
0142 status = 1;
0143 
0144 
0145         
0146 
0147 
0148