make_factor

status = make_factor(csdmfiles,outfname,win,binmin,binmax,nfactor,bad_chan,output);

    this function makes factor 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 = root name of output file factor # is appended
    win = which windows (observation number) to use in each file 
    COMING SOON    ifwin = -1 generate one file for each observation
            if win = 0 generate one file fo reach subject
    freqlist = frequency list to analyze
    nfactor = #of factor files to create
    bad_chan = badchannel list 
    output = 'cos' - cosine of phase 'ang' phase angle defaults to 'cos'

0001 function status = make_factor(csdmfiles,outfname,win,freqlist,nfactor,bad_chan,output);
0002 %status = make_factor(csdmfiles,outfname,win,binmin,binmax,nfactor,bad_chan,output);
0003 %
0004 %    this function makes factor 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 = root name of output file factor # is appended
0008 %    win = which windows (observation number) to use in each file
0009 %    COMING SOON    ifwin = -1 generate one file for each observation
0010 %            if win = 0 generate one file fo reach subject
0011 %    freqlist = frequency list to analyze
0012 %    nfactor = #of factor files to create
0013 %    bad_chan = badchannel list
0014 %    output = 'cos' - cosine of phase 'ang' phase angle defaults to 'cos'
0015 %
0016 
0017 if ~(nargin==6|nargin==7|nargin == 5)
0018     error('improper arguments')
0019 end;
0020 
0021 if nargin == 5
0022     bad_chan = [];
0023     output = 'cos';
0024 elseif nargin == 6
0025     output = 'cos';
0026 end;
0027 
0028 
0029 %
0030 % read in a csdm header in order to make an output file header
0031 %
0032 ave_hdr_offsets_v;
0033 %csdm_fid = fopen(csdmfiles(1,:),'r');
0034 %[csdm_fhdr,csdm_chdr,ename,czeros,cgains,cnames,fcom,ftext, coff]=rd_csdm_hdr_v(csdm_fid);
0035 %fclose(csdm_fid);
0036 %[cell_data_offsets]=get_csdm_cell_off(csdm_fhdr, csdm_chdr);%
0037 %
0038 %if win > csdm_chdr(1,NObs)
0039 %    error('window out of range')
0040 %end;
0041 %
0042 %ave_fhdr = csdm_fhdr;
0043 % Copy average file subject specifics array
0044 %ave_chdr = concat_csdm_chdr(csdmfiles, win);
0045 
0046 %if csdm_fhdr(NChan) == 8385
0047 %    ave_fhdr(NChan) = 129;
0048 %elseif csdm_fhdr(NChan) == 2145
0049 %    ave_fhdr(NChan) = 65;
0050 %elseif csdm_fhdr(NChan) == 33153;
0051 %    ave_fhdr(NChan) = 257;
0052 %else
0053 %    error('Nchan is unacceptable in this csdm file')
0054 %end;
0055 %
0056 % Read in copy of the csdm header
0057 %
0058 
0059 [csdm_fhdr,csdm_chdr]=get_csdm_hdr(csdmfiles(1,:));
0060 
0061 if win == []
0062   win = 1:csdm_chdr(1,NObs);
0063 end
0064 
0065 %
0066 % Copy csdm headers to average header
0067 %
0068 
0069 [ave_fhdr,ave_chdr,ave_ename,ave_cnames,ave_fcom,ave_ftext] = csdmhdr2avehdr(csdmfiles, win);
0070 
0071 % Correct the NPoints field for each cell header
0072 % This must be done independently in each m-file since in some cases the
0073 % resulting matrices are NChan X NChan instead of NPoints X NChan
0074 
0075 for c = 1:csdm_fhdr(NCells)
0076         ave_chdr(c,NPoints) = size(freqlist,2);
0077 end
0078 
0079 %
0080 % Determine the number of output files
0081 %
0082 
0083 if size(csdmfiles,1) > 1 & length(win) > 1
0084   nwin_files = length(win);
0085 else
0086   nwin_files = 1;
0087 end
0088 
0089 nwin = length(win);
0090 
0091 %
0092 % Write out the average file header -- same header if multiple files
0093 %
0094 
0095 out_root = outfname;
0096 for i = 1:nfactor
0097     for f = 1:nwin_files
0098       if nwin_files > 1 & win ~= []
0099             outfname1 = [out_root '_mag_' int2str(win(f)) '_ft_' int2str(i)];
0100             outfname2 = [outfname '_' output '_' int2str(win(f)) '_ft_' int2str(i)];
0101         else
0102         outfname1 = [out_root '_mag_ft_' int2str(i)];
0103             outfname2 = [out_root '_' output '_ft_' int2str(i)];
0104              end
0105     
0106         ave_fid = fopen(outfname1,'wb');
0107         ave_fid2 = fopen(outfname2,'wb');
0108         wt_ave_hdr_v(ave_fid,ave_fhdr,ave_chdr,ave_ename,ave_cnames,ave_fcom,ave_ftext);
0109         wt_ave_hdr_v(ave_fid2,ave_fhdr,ave_chdr,ave_ename,ave_cnames,ave_fcom,ave_ftext);
0110 
0111       fclose(ave_fid);
0112       fclose(ave_fid2);
0113     end
0114 end;
0115 
0116 ch_pair_indices;
0117 
0118 [cell_data_offsets]=get_csdm_cell_off(csdm_fhdr, csdm_chdr);
0119 
0120 file_mask = ones(size(csdmfiles,1),ave_fhdr(NChan));
0121 if (nargin == 7|nargin == 8)
0122     for i=1:size(bad_chan,1)
0123         bchan = find(bad_chan(i,:));
0124         if bchan ~= []
0125             file_mask(i,bad_chan(bchan)) = zeros(1,size(bchan,2));
0126         end;
0127     end;
0128 end;
0129 %for i=1:nfactor
0130 %    outfname1 = [outfname '_mag_' int2str(i)]
0131 %    outfname2 = [outfname '_' output '_' int2str(i)]
0132 %    avefid(i) = fopen(outfname1,'wb');
0133 %    avefid(i+nfactor) = fopen(outfname2, 'wb');
0134 %    % Write out the average file header
0135 %    wt_ave_hdr_v(avefid(i),ave_fhdr,ave_chdr,ave_ename,ave_cnames,ave_fcom,ave_ftext);
0136 %    fseek(avefid(i), 0, 'eof');
0137 %        wt_ave_hdr_v(avefid(i+nfactor),ave_fhdr,ave_chdr,ave_ename,ave_cnames,ave_fcom,ave_ftext);
0138 %        fseek(avefid(i+nfactor), 0, 'eof');
0139 %    fclose(avefid(i));
0140 %    fclose(avefid(i+nfactor));
0141 %end;
0142 
0143 ch_pair_indices;
0144 outfname3 = [outfname '_eig'];
0145 eig_fid = fopen(outfname3,'wb');
0146 fwrite(eig_fid,[csdm_fhdr(NCells) size(csdmfiles,1) ave_chdr(1,NPoints) nwin],'int16');
0147 for icell = 1:csdm_fhdr(NCells)
0148     for ifile = 1:size(csdmfiles,1)
0149        for f = 1:nwin
0150         [csdm_data, csdm_fhdr, csdm_chdr] = read_csdm(csdmfiles(ifile,:), icell, win(f));
0151             
0152         csdm_mat = zeros(ave_fhdr(NChan));
0153         for jfreq = 1:size(freqlist,2)
0154             ifreq= freqlist(jfreq);
0155             for ichan =1:ave_fhdr(NChan)
0156                 for jchan = ichan:ave_fhdr(NChan);
0157                     csdm_mat(ichan,jchan) = csdm_data(ifreq,chpair(ichan,jchan));
0158                     csdm_mat(jchan,ichan) = conj(csdm_mat(ichan,jchan));
0159                 end;
0160             end;
0161             good_chan = find(file_mask(ifile,:));
0162             csd = zeros(size(good_chan,2));
0163             csd = csdm_mat(good_chan,good_chan);
0164         
0165             [Vec_csd, val_csd] = eig(csd,'nobalance');
0166 %            val_csd([1:3 127:129],[1:3 127:129])
0167             if ~isreal(val_csd)
0168 %                disp('eigenvalues are not real, fixing')
0169 %                save val_csdm
0170 %                val_csd(size(good_chan,2),size(good_chan,2))
0171                 val_csd = real(val_csd);
0172             end;
0173             vec_csd = vec_csd';
0174             abs_vec_csd = 500*abs(vec_csd);
0175             fwrite(eig_fid,[icell ifreq size(val_csd,1) f],'int16');
0176 %            fwrite(eig_fid,ifile,'int16');
0177 %            fwrite(eig_fid,ifreq,'int16');
0178 %            fwrite(eig_fid,size(val_csd,1),'i
0179                 fwrite(eig_fid,diag(val_csd),'float');
0180         
0181             phase_data = vec_csd;
0182                         if output == 'cos'
0183                                 phase_data = 500*real(phase_data)./abs(phase_data);
0184                         else
0185                                 real_sign = sign(real(phase_data));
0186                                 phase_data = angle(phase_data);
0187                                 ang_sign = sign(phase_data);
0188                                 for iph = 1:size(phase_data,1)
0189                                         for jph = size(phase_data,2)
0190                                                 if real_sign(iph,jph) < 0
0191                                                         phase_data(iph,jph) = phase_data(iph,jph) + pi*real_sign(iph,jph)*ang_sign(iph,jph);
0192                                                 end;
0193                                         end;
0194                                 end;
0195 
0196 
0197                                 phase_data = 500*phase_data;
0198                         end
0199             phase_data(find(isnan(phase_data))) = zeros(size(find(isnan(phase_data))));
0200 
0201             if ~isreal(phase_data)
0202                           error('Phase matrix is not real')
0203                     end
0204 %            if nwin_files > 1 & win ~= []
0205 %                outfname1 = [out_root '_mag_' int2str(win(f))];
0206 %                    outfname2 = [outfname '_' output '_' int2str(win(f))];
0207 %              else
0208 %                   outfname1 = [out_root '_mag'];
0209 %                outfname2 = [out_root '_' output];
0210 %              end
0211 
0212 
0213             for ifac = 1:nfactor
0214                 mag = zeros(1, ave_fhdr(NChan));
0215                 ph = zeros(1,ave_fhdr(NChan));
0216                 mag(1,good_chan) = abs_vec_csd(ifac,:);
0217                 ph(1,good_chan) = phase_data(ifac,:);
0218 %                    outfname1 = [outfname '_mag_' int2str(ifac)];
0219 %                    outfname2 = [outfname '_' output '_' int2str(ifac)];
0220                             if nwin_files > 1 & win ~= []
0221                                 outfname1 = [out_root '_mag_' int2str(win(f)) '_ft_' int2str(ifac)];
0222                                 outfname2 = [outfname '_' output '_' int2str(win(f)) '_ft_' int2str(ifac)];
0223                             else
0224                                    outfname1 = [out_root '_mag_ft_' int2str(ifac)];
0225                                    outfname2 = [out_root '_' output '_ft_' int2str(ifac)];
0226                             end
0227 
0228                     avefid1 = fopen(outfname1,'ab');
0229                     avefid2= fopen(outfname2, 'ab');
0230 
0231                 num_written = fwrite(avefid1, mag','int16');
0232                 if num_written~= size(mag,2)
0233                     error('write out croaked')
0234                 end;
0235                 num_written = fwrite(avefid2, ph','int16');
0236                             if num_written~= size(ph,2)
0237                                     error('write out croaked')
0238                             end;
0239                 fclose(avefid1);
0240                 fclose(avefid2);
0241             end    
0242         end;
0243       end;
0244     end; %for ifile =
0245 end; % for icell
0246 fclose('all');
0247 
0248 
0249