make_power

power_data = make_power(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
        defaults to 'amp'

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