make_power_sc

power_data = make_power_sc(csdmfiles,outfname,win,output);

    this function makes power files derived files from csdm files
    the names of the csdm files must be stored in an array where each row 
        is a filename - csdmfiles.   
    outfname = name of output file
    win =     which window (observation number) to use in each file 
            if win == [], all observations are used
            if length(win) > 0 and there is more than one csdmfile specified, a
           separate file is created for each window, with NObs in each file
            corresponding to the number of input csdmfiles
    output = 'amp' for amplitude, 'pow' for power, and 'log' for log
        'norm' for by space normalization,  
        
    defaults to 'amp' and bad_chan = []

0001 function scale_values = make_power_sc(csdmfiles,outfname,win,output,bad_chan);
0002 %power_data = make_power_sc(csdmfiles,outfname,win,output);
0003 %
0004 %    this function makes power files derived files from csdm files
0005 %    the names of the csdm files must be stored in an array where each row
0006 %        is a filename - csdmfiles.
0007 %    outfname = name of output file
0008 %    win =     which window (observation number) to use in each file
0009 %            if win == [], all observations are used
0010 %            if length(win) > 0 and there is more than one csdmfile specified, a
0011 %           separate file is created for each window, with NObs in each file
0012 %            corresponding to the number of input csdmfiles
0013 %    output = 'amp' for amplitude, 'pow' for power, and 'log' for log
0014 %        'norm' for by space normalization,
0015 %
0016 %    defaults to 'amp' and bad_chan = []
0017 %
0018 
0019 % Original version by RS
0020 % Modification history:
0021 %  Introduced get_csdm_hdr and csdmhdr2avehdr -- PJ
0022 %
0023 %  02/05/96 PJ added 'norm' option which returns proportion variance
0024 %              explained scaled by 500.  This eliminates the risk of int16 overflow
0025 %              as when using 'pow'.
0026 
0027 if ~(nargin==3|nargin ==5)
0028     error('improper arguments')
0029 end;
0030 
0031 if nargin == 3
0032     output ='amp';
0033     bad_chan = [];
0034 elseif nargin == 4
0035     bad_chan = [];
0036 end
0037 
0038 ave_hdr_offsets_v;
0039 
0040 %
0041 % Read in copy of the csdm header
0042 %
0043 
0044 [csdm_fhdr,csdm_chdr]=get_csdm_hdr(csdmfiles(1,:));
0045 
0046 if win == []
0047   win = 1:csdm_chdr(1,NObs);
0048 end
0049 
0050 %
0051 % Copy csdm headers to average header
0052 %
0053 
0054 [ave_fhdr,ave_chdr,ave_ename,ave_cnames,ave_fcom,ave_ftext] = csdmhdr2avehdr(csdmfiles, win);
0055 
0056 % Correct the NPoints field for each cell header
0057 % This must be done independently in each m-file since in some cases the
0058 % resulting matrices are NChan X NChan instead of NPoints X NChan
0059 
0060 for c = 1:csdm_fhdr(NCells)
0061     ave_chdr(c,NPoints) = csdm_chdr(c,NPoints)/2;
0062 end
0063 
0064 %
0065 % Determine the number of output files
0066 %
0067 
0068 if size(csdmfiles,1) > 1 & length(win) > 1
0069   nwin_files = length(win);
0070 else
0071   nwin_files = 1;
0072 end
0073 
0074 nwin = length(win);
0075 
0076 %
0077 % Write out the average file header -- same header if multiple files
0078 %
0079 
0080 out_root = outfname;
0081 for f = 1:nwin_files
0082   if nwin_files > 1 & win ~= []
0083     outfname1 = [out_root '_' int2str(win(f))];
0084   else
0085     outfname1 = out_root;
0086   end
0087 end
0088   
0089   ave_fid = fopen(outfname1,'wb');
0090   wt_ave_hdr_v(ave_fid,ave_fhdr,ave_chdr,ave_ename,ave_cnames,ave_fcom,ave_ftext);
0091   fclose(ave_fid);
0092 end
0093 if ave_fhdr(NChan) == 129
0094 ch_pair_indices;
0095 else
0096 ch_pair_indices_65;
0097 end;
0098 
0099 [cell_data_offsets]=get_csdm_cell_off(csdm_fhdr, csdm_chdr);
0100 if ~strcmp(output,'norm')
0101     scale_values = 0;
0102 end;
0103 is = 1;
0104 for icell = 1:csdm_fhdr(NCells)
0105   disp(['Working on cell ' int2str(icell)]);
0106     for ifile = 1:size(csdmfiles,1)
0107       disp(['  Reading from file: ' csdmfiles(ifile,:)]);
0108       for iwin = 1:nwin
0109         disp(['    Window ' int2str(iwin) '  (' int2str(win(iwin)) ')']);
0110         [csdm_data, csdm_fhdr, csdm_chdr] = read_csdm(csdmfiles(ifile,:), icell, win(iwin));
0111         npts = ave_chdr(icell,NPoints);
0112         power_data = zeros(ave_chdr(icell,NPoints),ave_fhdr(NChan));
0113         for ichan = 1:ave_fhdr(NChan)
0114             power_data(:,ichan) = csdm_data(:,chpair(ichan,ichan));
0115         end;
0116         
0117         % Reality check
0118         testreal = isreal(sum(power_data));
0119         if testreal ~= 1
0120             power_data(11,73)
0121             error('complex numbers')
0122         end;
0123         
0124         if strcmp(output,'amp')
0125              power_data = 500*sqrt(power_data);
0126         elseif strcmp(output,'pow')    
0127             power_data = 500*power_data;
0128           while any(max(max(power_data)) > 2^15)
0129             disp('Warning: Int16 overflow -- scaling data down by factor of 500');
0130             power_data = power_data / 500;
0131           end
0132         elseif strcmp(output,'norm')
0133           row_av = sum(power_data')'./((ave_fhdr(NChan)-size(find(bad_chan(ifile,:)),2))*ones(size(power_data,1),1));
0134           power_data(1:npts,:) = 500*power_data(1:npts,:) ./ (row_av*ones(1,size(power_data,2)));
0135           scale_values(is,:) = row_av';
0136           is = is + 1;
0137         else  % default to log10
0138             power_data = 500*log10(power_data);
0139         end
0140 
0141         if nwin_files > 1
0142           outfname1 = [out_root '_' int2str(win(iwin))];
0143                 else
0144           outfname1 = out_root;
0145         end
0146 
0147         ave_fid = fopen(outfname1, 'ab');
0148         fseek(ave_fid, 0, 'eof');
0149 
0150         num_written = fwrite(ave_fid,power_data','int16');
0151         if num_written~= size(power_data,1)*size(power_data,2)
0152             error('write failed')
0153         end;
0154 
0155         fclose(ave_fid);
0156         
0157       end; % for iwin =
0158     end; %for ifile =
0159 end; % for icell
0160 fclose('all');
0161 
0162