make_coherence

[coh_mat] = make_coherence(csdmfiles,outfname,win, cf_list, cf_flag, bad_chans);

  this function makes coherence files derived files from csdm files

  csdmfiles = the names of the csdm files must be stored in an array where each row 
              is a filename.   
  outfname = root name of output file.  The number of the channel or frequency with
  respect to which the coherence matrix is computed is appended to the
  root name.  _f3 indicates frequency bin 3 of the original csdm.
    _c69 indicates channel 69.
            
 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

    cf_list = list of channels or frequency bins
    cf_flag = 'chan' if coherence is computed relative to one channel.  Resulting matrix
                 is Npoints X NChan
             = 'freq' if coherence is computed relative to all channels, but only for a 
                single frequency.  Resulting matrix is NChan X NChan

    bad_chans = list of bad channels -- one row for each row of csdmfiles

0001 function [coh_mat] = make_coherence(csdmfiles,outfname,win, cf_list, cf_flag, bad_chans);
0002 %[coh_mat] = make_coherence(csdmfiles,outfname,win, cf_list, cf_flag, bad_chans);
0003 %
0004 %  this function makes coherence files derived files from csdm files
0005 %
0006 %  csdmfiles = the names of the csdm files must be stored in an array where each row
0007 %              is a filename.
0008 %  outfname = root name of output file.  The number of the channel or frequency with
0009 %  respect to which the coherence matrix is computed is appended to the
0010 %  root name.  _f3 indicates frequency bin 3 of the original csdm.
0011 %    _c69 indicates channel 69.
0012 %
0013 % win =     which window (observation number) to use in each file
0014 %        if win == [], all observations are used
0015 %        if length(win) > 0 and there is more than one csdmfile specified, a
0016 %           separate file is created for each window, with NObs in each file
0017 %            corresponding to the number of input csdmfiles
0018 %
0019 %    cf_list = list of channels or frequency bins
0020 %    cf_flag = 'chan' if coherence is computed relative to one channel.  Resulting matrix
0021 %                 is Npoints X NChan
0022 %             = 'freq' if coherence is computed relative to all channels, but only for a
0023 %                single frequency.  Resulting matrix is NChan X NChan
0024 %
0025 %    bad_chans = list of bad channels -- one row for each row of csdmfiles
0026 
0027 %  01/20/96 PJ Started modification tracking and modified for propagation of bad_chans
0028 %  01/27/96 PJ Combined make_coh_by_freq and make_coh_by_chan
0029 
0030 status = -1;
0031 format compact
0032 
0033 %
0034 %  Check input arguments
0035 %
0036 
0037 if (nargin<5)
0038   error('Too few arguments to make_coherence()')
0039 elseif nargin<6
0040   bad_chans = [];
0041 else
0042   if size(bad_chans,1) ~= size(csdmfiles,1)
0043     error('csdmfile list and bad_chans lists don''t match');
0044   end
0045 end;
0046 
0047 if ~(strcmp(cf_flag, 'freq') | strcmp(cf_flag, 'chan'))
0048   error(['Invalid cf_flag: ' cf_flag]);
0049 end
0050 
0051 ave_hdr_offsets_v;
0052 
0053 %
0054 % Read in copy of the csdm header
0055 %
0056 
0057 [csdm_fhdr,csdm_chdr]=get_csdm_hdr(csdmfiles(1,:));
0058 
0059 if win == []
0060   win = 1:csdm_chdr(1,NObs);
0061 end
0062 
0063 
0064 %
0065 % Copy csdm headers to average header
0066 %
0067 
0068 [ave_fhdr,ave_chdr,ave_ename,ave_cnames,ave_fcom,ave_ftext] = csdmhdr2avehdr(csdmfiles, win);
0069 
0070 % Correct the NPoints field for each cell header
0071 % This must be done independently in each m-file since in some cases the
0072 % resulting matrices are NChan X NChan instead of NPoints X NChan
0073 
0074 for c = 1:csdm_fhdr(NCells)
0075   if strcmp(cf_flag, 'chan')
0076     ave_chdr(c,NPoints) = csdm_chdr(c,NPoints)/2;
0077   else
0078      ave_chdr(c,NPoints) = ave_fhdr(NChan);
0079   end 
0080 end
0081 
0082 %
0083 % Determine the number of output files corresponding to windows
0084 %
0085 
0086 if size(csdmfiles,1) > 1 & length(win) > 1
0087   nwin_files = length(win);
0088 else
0089   nwin_files = 1;
0090 end
0091 
0092 nitems = size(cf_list,2);
0093 nwin = length(win);
0094 
0095 %
0096 % Get fileids for all of the resulting coherence files
0097 % The file ids are stored in a nitems X nwin matrix
0098 %
0099 
0100 new_root = outfname;
0101 if strcmp(cf_flag, 'freq')
0102   new_root = [new_root '_cf'];  % cf = coherence by frequency
0103 else
0104   new_root = [new_root '_cc'];    % cc = coherence by channel
0105 end
0106 
0107 fid_mat = zeros(nitems,nwin_files);
0108 
0109 for iref = 1:nitems
0110     iref_str = int2str(cf_list(iref));
0111     outfname1 = [new_root iref_str];
0112     if nwin_files == 1
0113       fid_mat(iref, 1) = fopen(outfname1,'wb');
0114       wt_ave_hdr_v(fid_mat(iref,1),ave_fhdr,ave_chdr,ave_ename,ave_cnames,ave_fcom,ave_ftext);
0115       fseek(fid_mat(iref, 1), 0, 'eof');
0116     else
0117       for fref = 1:nwin_files
0118         text_ref = int2str(fref);
0119         outfname1 = [new_root iref_str '_w' text_ref];
0120         fid_mat(iref, fref) = fopen(outfname1,'wb');
0121         wt_ave_hdr_v(fid_mat(iref,fref),ave_fhdr,ave_chdr,ave_ename,ave_cnames,ave_fcom,ave_ftext);
0122 %        fseek(fid_mat(iref, fref), 0, 'eof');
0123         fclose(fid_mat(iref,fref));
0124       end
0125     end
0126 end;
0127 ch_pair_indices;
0128 
0129 ch_str = [];
0130 for c = 1:length(cf_list)
0131   ch_str = [ch_str int2str(cf_list(c)) ' '];
0132 end
0133 if strcmp(cf_flag, 'freq')
0134   disp(['Computing all channel coherences for frequency bins ' ch_str]);
0135 else
0136   disp(['Computing coherences for all frequency bins with respect to chans: ' ch_str]);
0137 end
0138 
0139 [cell_data_offsets]=get_csdm_cell_off(csdm_fhdr, csdm_chdr);
0140 
0141 for icell = 1:csdm_fhdr(NCells)
0142   disp(['Working on cell ' int2str(icell)]);
0143   npts = ave_chdr(icell,NPoints);
0144   nchan = ave_fhdr(NChan);
0145     
0146   for ifile = 1:size(csdmfiles,1)
0147     disp(['  Input file: ' csdmfiles(ifile,:)]);
0148     for iobs = 1:nwin
0149       disp(['  Observation ' int2str(iobs) '  (' int2str(win(iobs)) ')']);
0150 
0151         % read in the csdm matrix
0152         [csdm_data, csdm_fhdr, csdm_chdr] = read_csdm(csdmfiles(ifile,:), icell, win(iobs));
0153         
0154         % compute coherence matrix successively for each item in cf_list
0155         for iref = 1:nitems
0156           if strcmp(cf_flag, 'freq')
0157             disp(['    Frequency band ' int2str(cf_list(iref))]);
0158           else
0159             disp(['    Channel ' int2str(cf_list(iref))]);
0160           end
0161           
0162           % Size matrices
0163           if strcmp(cf_flag, 'freq')
0164             power_data = zeros(1,nchan);
0165           else
0166             power_data = zeros(npts, nchan);
0167           end
0168           
0169           coh_data = zeros(npts,nchan);
0170           numerator = zeros(npts, nchan);
0171     
0172           % Generate the matrix of autospectra
0173           if strcmp(cf_flag, 'freq')                    
0174               for ichan = 1:nchan
0175                 power_data(ichan) = csdm_data(cf_list(iref),chpair(ichan,ichan));
0176             end;          
0177 
0178              % Put in NaNs where there are bad_chans -- will be zeroed later
0179             if bad_chans ~= []
0180                  power_data(bad_chans(ifile,find(bad_chans(ifile,:)))) = ones(size(find(bad_chans(ifile,:)))) * NaN;
0181               end
0182               power_data = power_data' * power_data;
0183           else
0184             for ichan = 1:nchan 
0185               power_data(:,ichan) = csdm_data(1:npts,chpair(ichan,ichan));
0186             end
0187             if bad_chans ~= []
0188                 power_data(:,bad_chans(ifile,find(bad_chans(ifile,:)))) = ones(size(power_data(:,bad_chans(ifile,find(bad_chans(ifile,:)))))) * NaN;
0189             end;
0190           end
0191           
0192           % Make sure the autospectra matrix is valid
0193           testreal = isreal(sum(power_data));
0194           if testreal ~= 1
0195             power_data(1,[6 54 96])
0196             error('complex numbers in power data matrix')
0197           end;
0198             
0199           %
0200           % PERFORM THE COHERENCE CALCULATION
0201           %
0202           
0203           % COHERENCE BY FREQUENCY CASE
0204           if strcmp(cf_flag, 'freq')
0205              % Generate the numerator for the coherence calc
0206               if strcmp(computer, 'MAC2')
0207                 for rchan = 1:nchan
0208               numerator(rchan,:) = ppc_cmplx_mult(csdm_data(cf_list(iref),chpair(rchan, 1:nchan)),conj(csdm_data(cf_list(iref),chpair(rchan, 1:nchan))));
0209              end
0210             else
0211               for rchan = 1:nchan
0212             numerator(rchan,:) = csdm_data(cf_list(iref),chpair(rchan, 1:nchan)).*conj(csdm_data(cf_list(iref),chpair(rchan, 1:nchan)));
0213                 end
0214               end
0215           
0216               % Compute the coherence
0217              coh_data = numerator./power_data;
0218         
0219           else    % COHERENCE BY CHANNEL CASE
0220               % Generate the numerator for the coherence calc
0221             if strcmp(computer, 'MAC2')
0222               for ichan = 1:nchan
0223                 numerator(:,ichan) = ppc_cmplx_mult(csdm_data(1:npts,chpair(cf_list(iref), ichan)),conj(csdm_data(1:npts,chpair(cf_list(iref), ichan))));
0224               end
0225             else
0226               for ichan = 1:nchan
0227                 numerator(:,ichan) = csdm_data(1:npts,chpair(cf_list(iref), ichan)).*conj(csdm_data(1:npts,chpair(cf_list(iref), ichan)));
0228               end
0229             end
0230             
0231             % Compute the coherence
0232             coh_data = numerator./((power_data(1:npts,cf_list(iref))*ones(1,nchan)) .* power_data);
0233           end
0234           
0235           %
0236           % END OF COHERENCE CALC
0237           %
0238           
0239           % Replace NaNs with zeros
0240           coh_data(find(isnan(coh_data))) = zeros(size(find(isnan(coh_data))));
0241           
0242           %
0243           % Error check various aspects of the coherence matrix
0244           %
0245           if strcmp(cf_flag, 'freq')
0246               testcoh = sum(diag(coh_data));
0247             if testcoh ~= (nchan - length(find(bad_chans(ifile,:))))
0248                 if ~isempty(findstr(csdmfiles(ifile,:), '2D'))
0249                 disp(['You are a moron who failed to detect bad chans' int2str(testcoh)])
0250                   else
0251                 disp(['Undetected Bad Chan ' int2str(testcoh)]);
0252                 end
0253             end
0254           else
0255             testcoh = sum(coh_data(1:npts,cf_list(iref)));
0256             if testcoh ~= npts
0257               disp(['Your character is suspect ' int2str(testcoh)]);
0258             end
0259           end
0260           
0261           if ~isreal(coh_data)
0262             error('Coherence matrix is not real')
0263           end
0264                     
0265           coh_mat = coh_data;
0266           
0267           %
0268           % Scale coherence matrix to have representation as int16
0269           %
0270           
0271           coh_data = coh_data * 500;
0272 
0273           iref_str = int2str(cf_list(iref));
0274           if nwin_files == 1
0275             outfname1 = [new_root iref_str];
0276                   fid_mat = fopen(outfname1,'ab');
0277                  fseek(fid_mat, 0, 'eof');
0278                   num_written = fwrite(fid_mat,coh_data','int16');
0279             fclose(fid_mat);
0280           else
0281               for fref = 1:nwin_files
0282                     text_ref = int2str(fref);
0283                         outfname1 = [new_root iref_str '_w' text_ref];
0284                         fid_mat = fopen(outfname1,'ab');
0285                         fseek(fid_mat, 0, 'eof');
0286                 num_written = fwrite(fid_mat,coh_data','int16');
0287                   if num_written~= size(coh_data,1)*size(coh_data,2)
0288                     error('write failed')
0289                   end;
0290 
0291                         fclose(fid_mat);
0292                 end
0293           end
0294 %          if nwin_files > 1
0295 %            num_written = fwrite(fid_mat(iref, iobs),coh_data','int16');
0296 %          else
0297 %            num_written = fwrite(fid_mat(iref, 1),coh_data','int16');
0298 %          end
0299             
0300           if num_written~= size(coh_data,1)*size(coh_data,2)
0301             error('write failed')
0302           end;
0303         end;  %for iref =
0304       end  % for iobs
0305     end; %for ifile =
0306 end; % for icell
0307 fclose('all');
0308 
0309 format loose
0310 status = 0;
0311 return
0312