Non-looped fread of binary file with multiple value types

2 views (last 30 days)
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

Answers (4)

Josh Parks
Josh Parks on 27 Sep 2012
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

Walter Roberson
Walter Roberson on 18 Sep 2012
  • 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

Josh Parks
Josh Parks on 20 Sep 2012
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);
  2 Comments
Walter Roberson
Walter Roberson on 20 Sep 2012
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
Jan on 20 Sep 2012
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.

Sign in to comment.


Jan
Jan on 20 Sep 2012
Edited: Jan on 20 Sep 2012
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

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!