Non-looped fread of binary file with multiple value types

Hi all,

This is my first post here and I am having quite the trouble with this little program. I am trying to read data out of a proprietary binary generated by a photon counting card. I've got it reading the data correctly but it is super slow. I'm guessing this is because of my implementation of fread, however I can't figure out how to do this without a loop. Any suggestions from you matlab gurus out there? Currently it takes about nine seconds on my computer to read in a 452kB file. I've attached my program with a set of example data. For those of you who don't care to download anything, here is the code excerpt in suspect:

    %  This reads the TTTR event records
    Ofltime = 0;
    Counts = 0;
    outfile = ['t3r.out'];
    fpout = fopen(outfile,'W');
    fprintf(1,'Writing data to %s\n', outfile);
    fprintf(1,'This may take a while...');

    % fprintf(fpout,'      #   timetag    time/s   channel  route\n\n');
    for i=1:NumberOfRecords

      TTTRRecord = fread(fid, 1, 'uint32');
      TimeTag = bitand(TTTRRecord,65535);              %the lowest 16 bits
      Channel = bitand(bitshift(TTTRRecord,-16),4095); %the next 12 bits
      Route   = bitand(bitshift(TTTRRecord,-28),3);    %the next 2 bits
      Valid   = bitand(bitshift(TTTRRecord,-30),1);    %the next bit
      Truetime = (Ofltime + TimeTag) * TTTRGlobclock * 1e-9;

      if Valid
        fprintf(fpout,'%7u  %7u  %11.7f  %5u   %2u', Counts, TimeTag, Truetime, Channel, Route);
        Counts = Counts+1;
      else %this means we have a special record
          if bitand(Channel,2048) 
              Ofltime = Ofltime+65536;
          end;    
          if bitand (Channel,7)  
              fprintf(fpout,'Marker=%u\n',bitand(Channel,7));
          end;    
      end;

      if Valid 
        fprintf(fpout,'\n');
      end;

    end;

Thanks,

Josh

Program zip

https://docs.google.com/open?id=0B84phdetR0DkYktXT2pxWGt2YUE

Answer 1

[Josh Parks]

Hi all,

in case anyone was wondering, this is what I came up with:

%% Data Processing

    %  This reads the TTTR event records
    outfile = [filename(1:length(filename)-4),'.txt'];
    fprintf(1,'Analyzing data from %s\n', outfile);
    fprintf(1,'This may take a few seconds...');

    Ofltime = 0;
    Counts = 0;
    TTTRData = fread(fid, [1 NumberOfRecords], 'uint32');
    TTTRProcessed = zeros(NumberOfRecords,5);

    %read values out of bitrecord for each photon and assimilate
    %accordinlgy
    position = ftell(fid)
    ValidA = bitand(bitshift(TTTRData(:),-30),ones(length(TTTRData),1)) == 1;
    CountsA = cumsum(ValidA);
    ChannelA = bitand(bitshift(TTTRData(:),-16),ones(length(TTTRData),1).*4095);
    %special record
    OfltimeA = cumsum(bitand(ChannelA(:),ones(length(TTTRData),1).*2048)./2048.*65536.*(~ValidA));
    TimeTagA = bitand(TTTRData(:),ones(length(TTTRData),1).*65535);
    TrueTimeA = (OfltimeA + TimeTagA).*TTTRGlobclock.*1e-9;
    RouteA = bitand(bitshift(TTTRData(:),-28),ones(length(TTTRData),1).*3);

    %zero unimportant data tags reduction
    ChannelA = ChannelA.*ValidA;
    TimeTagA = TimeTagA.*ValidA;
    RouteA = RouteA.*ValidA;
    TrueTimeA = TrueTimeA.*ValidA;
    CountsA = CountsA.*ValidA;

    %concatenate values for single matrix printing/saving
    TTTRProcessed = horzcat(CountsA,TimeTagA,TrueTimeA,ChannelA,RouteA);

    %removes non-data rows
    TTTRProcessed = TTTRProcessed(any(TTTRProcessed,2),:);

    if printTXT == 1
        fprintf(1,'Writing data to %s\n', outfile);
        fpout = fopen(outfile,'W');
        fprintf(fpout,'%7u\t%7u\t%11.7f\t%5u\t%2u\n',TTTRProcessed.');
        fclose(fid);
        fclose(fpout);
        fprintf(1,'\n%u photon records written to %s\n\n',max(CountsA),outfile);
    end

Answer 2

[Walter Roberson]

• you can read the entire file into memory and process it once it is in memory
• if Valid is not true, you can skip calculation of Route and so on.
• You can extract Valid with a single bitand() without doing a shift, as you only care if it is zero or non-zero
• assign the result of bitand(Channel,7) to a variable so you do not need to recalculate it for printing purposes

Answer 3

[Josh Parks]

Thanks walter, but is there a way to completely vectorize the data read in? This is what I've come up with so far, but it still takes 2 minutes to run on large files because of the loop

    %  This reads the TTTR event records
    Ofltime = 0;
    Counts = 0;
    outfile = ['t3r.out'];
    fpout = fopen(outfile,'W');
    fprintf(1,'Writing data to %s\n', outfile);
    fprintf(1,'This may take a while...');
    TTTRData = fread(fid, [1 NumberOfRecords], 'uint32');
    TTTRProcessed = zeros(NumberOfRecords,5);

    % fprintf(fpout,'      #   timetag    time/s   channel  route\n\n');
    for i=1:NumberOfRecords
      TTTRRecord = TTTRData(i);
      Valid   = bitand(bitshift(TTTRRecord,-30),1);    %the next bit
      Channel = bitand(bitshift(TTTRRecord,-16),4095); %the next 12 bits
      if Valid
        TimeTag = bitand(TTTRRecord,65535);              %the lowest 16 bits
        Route   = bitand(bitshift(TTTRRecord,-28),3);  %the next 2 bits 

        Truetime = (Ofltime + TimeTag) * TTTRGlobclock * 1e-9;
        TTTRProcessed(i,1) = Counts;
        TTTRProcessed(i,2) = TimeTag;
        TTTRProcessed(i,3) = Truetime;
        TTTRProcessed(i,4) = Channel;
        TTTRProcessed(i,5) = Route;
%         fprintf(fpout,'%7u  %7u  %11.7f  %5u   %2u', Counts, TimeTag, Truetime, Channel, Route);
          Counts = Counts+1;
      else %this means we have a special record
          if bitand(Channel,2048) 
              Ofltime = Ofltime+65536;
          end;    
          if bitand(Channel,7)  
              TTTRProcessed(i,1) = -1;
              TTTRProcessed(i,2) = -1;
              TTTRProcessed(i,3) = -1;
              TTTRProcessed(i,4) = -1;
              TTTRProcessed(i,5) = -1;
          end;    
      end;
    end;

    %removes non-data rows
    TTTRProcessed = TTTRProcessed(any(TTTRProcessed,2),:);

    fprintf(fpout,'%7u\t%7u\t%11.7f\t%5u\t%2u\n',TTTRProcessed.');

    fclose(fid);
    fclose(fpout);
    fprintf(1,'\n%u photon records written to %s\n',Counts,outfile);

[Walter Roberson]

bitand() and bitshift() both work on arrays, so you can do all those extractions at once (possibly including TimeTag and Route). After that you could use logical indexing to select the Valid entries for printing.

[Jan]

Try if a multiplication is much faster than bitshift(). This has been the case at least in older Matlab versions, as far as I remember in R2009a.

Answer 4

[Jan]

Small but easy improvement:

TTTRRecord = TTTRData(i);
Valid   = bitand(bitshift(TTTRRecord,-30),1);    %the next bit
Channel = bitand(bitshift(TTTRRecord,-16),4095); %the next 12 bits
if Valid
  ...

To:

TTTRRecord = TTTRData(i);
Valid = bitand(TTTRRecord), 1073741824);    % 2^30
if Valid
  Channel = bitand(bitshift(TTTRRecord,-16),4095); %the next 12 bits
  ...

Perhaps:

    Channel = bitand(TTTRRecord/65536, 4095); %the next 12 bits

Another standard procedure ist to move all repeated calculations out of the loop, e.g.:

c1 = TTTRGlobclock * 1e-9;
...
     Truetime = (Ofltime + TimeTag) * c1;

Another improvement: Instead of 5 single calls use:

TTTRProcessed(i,1:5) = -1;
c1 = TTTRGlobclock * 1e-9;
TTTRProcessed = zeros(NumberOfRecords,5);
for i=1:NumberOfRecords
  TTTRRecord = TTTRData(i);
  Valid = bitand(TTTRRecord), 1073741824);        % 2^30
  if Valid
    Channel = bitand(TTTRRecord / 65536, 4095);   % next 12 bits
    TimeTag = bitand(TTTRRecord, 65535);          % lowest 16 bits
    Route   = bitand(TTTRRecord / 268435456), 3); % next 2 bits 
    Truetime = (Ofltime + TimeTag) * c1;

      TTTRProcessed(i,1) = Counts;
      TTTRProcessed(i,2) = TimeTag;
      TTTRProcessed(i,3) = Truetime;
      TTTRProcessed(i,4) = Channel;
      TTTRProcessed(i,5) = Route;
      Counts = Counts + 1;
    else %this means we have a special record
      if bitand(Channel,2048) 
        Ofltime = Ofltime + 65536;
      end
      if bitand(Channel, 7)
        TTTRProcessed(i, 1:5) = -1;
      end    
    end
  end