forked from open-ephys/analysis-tools
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathload_open_ephys_data_faster.m
167 lines (156 loc) · 6.85 KB
/
load_open_ephys_data_faster.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
function [data, timestamps, info] = load_open_ephys_data_faster(filename)
%
% [data, timestamps, info] = load_open_ephys_data(filename)
%
% Loads continuous, event, or spike data files into Matlab.
%
% Inputs:
%
% filename: path to file
%
%
% Outputs:
%
% data: either an array continuous samples (in microvolts),
% a matrix of spike waveforms (in microvolts),
% or an array of event channels (integers)
%
% timestamps: in seconds
%
% info: structure with header and other information
%
%
%
% DISCLAIMER:
%
% Both the Open Ephys data format and this m-file are works in progress.
% There's no guarantee that they will preserve the integrity of your
% data. They will both be updated rather frequently, so try to use the
% most recent version of this file, if possible.
%
%
%
% ------------------------------------------------------------------
%
% Copyright (C) 2014 Open Ephys
%
% ------------------------------------------------------------------
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% <http://www.gnu.org/licenses/>.
%
[~,~,filetype] = fileparts(filename);
if ~any(strcmp(filetype,{'.events','.continuous','.spikes'}))
error('File extension not recognized. Please use a ''.continuous'', ''.spikes'', or ''.events'' file.');
end
fid = fopen(filename);
fseek(fid,0,'eof');
filesize = ftell(fid);
NUM_HEADER_BYTES = 1024;
fseek(fid,0,'bof');
hdr = fread(fid, NUM_HEADER_BYTES, 'char*1');
info = getHeader(hdr);
if isfield(info.header, 'version')
version = info.header.version;
else
version = 0.0;
end
switch filetype
case '.events'
bStr = {'timestamps' 'sampleNum' 'eventType' 'nodeId' 'eventId' 'data' 'recNum'};
bTypes = {'int64' 'uint16' 'uint8' 'uint8' 'uint8' 'uint8' 'uint16'};
bRepeat = {1 1 1 1 1 1 1};
dblock = struct('Repeat',bRepeat,'Types', bTypes,'Str',bStr);
if version < 0.2, dblock(7) = []; end
if version < 0.1, dblock(1).Types = 'uint64'; end
case '.continuous'
SAMPLES_PER_RECORD = 1024;
bStr = {'ts' 'nsamples' 'recNum' 'data' 'recordMarker'};
bTypes = {'int64' 'uint16' 'uint16' 'int16' 'uint8'};
bRepeat = {1 1 1 SAMPLES_PER_RECORD 10};
dblock = struct('Repeat',bRepeat,'Types', bTypes,'Str',bStr);
if version < 0.2, dblock(3) = []; end
if version < 0.1, dblock(1).Types = 'uint64'; dblock(2).Types = 'int16'; end
case '.spikes'
num_channels = info.header.num_channels;
num_samples = 40;
bStr = {'eventType' 'timestamps' 'timestamps_software' 'source' 'nChannels' 'nSamples' 'sortedId' 'electrodeID' 'channel' 'color' 'pcProj' 'samplingFrequencyHz' 'data' 'gain' 'threshold' 'recordingNumber'};
bTypes = {'uint8' 'int64' 'int64' 'uint16' 'uint16' 'uint16' 'uint16' 'uint16' 'uint16' 'uint8' 'float32' 'uint16' 'uint16' 'float32' 'uint16' 'uint16'};
bRepeat = {1 1 1 1 1 1 1 1 1 3 2 1 num_channels*num_samples num_channels num_channels 1};
dblock = struct('Repeat',bRepeat,'Types', bTypes,'Str',bStr);
if version < 0.4, dblock(7:12) = []; dblock(8).Types = 'uint16'; end
if version == 0.3, dblock = [dblock(1), struct('Repeat',1,'Types','uint32','Str','ts'), dblock(2:end)]; end
if version < 0.3, dblock(2) = []; end
if version < 0.2, dblock(9) = []; end
if version < 0.1, dblock(2).Types = 'uint64'; end
end
blockBytes = str2double(regexp({dblock.Types},'\d{1,2}$','match', 'once')) ./8 .* cell2mat({dblock.Repeat});
numIdx = floor((filesize - NUM_HEADER_BYTES)/sum(blockBytes));
switch filetype
case '.events'
timestamps = segRead('timestamps')./info.header.sampleRate;
info.sampleNum = segRead('sampleNum');
info.eventType = segRead('eventType');
info.nodeId = segRead('nodeId');
info.eventId = segRead('eventId');
data = segRead('data');
if version >= 0.2, info.recNum = segRead('recNum'); end
case '.continuous'
if nargout>1
info.ts = segRead('ts');
end
info.nsamples = segRead('nsamples');
if ~all(info.nsamples == SAMPLES_PER_RECORD)&& version >= 0.1, error('Found corrupted record'); end
if version >= 0.2, info.recNum = segRead('recNum'); end
data = segRead_data('data', 'b').*info.header.bitVolts; % read in continuousdata
if nargout>1 % do not create timestamp arrays unless they are requested
timestamps = nan(size(data));
current_sample = 0;
for record = 1:length(info.ts)
timestamps(current_sample+1:current_sample+info.nsamples(record)) = info.ts(record):info.ts(record)+info.nsamples(record)-1;
current_sample = current_sample + info.nsamples(record);
end
end
case '.spikes'
timestamps = segRead('timestamps')./info.header.sampleRate;
info.source = segRead('source');
info.samplenum = segRead('nSamples');
info.gain = permute(reshape(segRead('gain'), num_channels, numIdx), [2 1]);
info.thresh = permute(reshape(segRead('threshold'), num_channels, numIdx), [2 1]);
if version >= 0.4, info.sortedId = segRead('sortedId'); end
if version >= 0.2, info.recNum = segRead('recordingNumber'); end
data = permute(reshape(segRead('data'), num_samples, num_channels, numIdx), [3 1 2]);
data = (data-32768)./ permute(repmat(info.gain/1000,[1 1 num_samples]), [1 3 2]);
end
fclose(fid);
function seg = segRead_data(segName, mf)
%% This function is specifically for reading continuous data.
% It keeps the data in int16 precision, which can drastically decrease
% memory consumption
if nargin == 1, mf = 'l'; end
segNum = find(strcmp({dblock.Str},segName));
fseek(fid, sum(blockBytes(1:segNum-1))+NUM_HEADER_BYTES, 'bof');
seg = fread(fid, numIdx*dblock(segNum).Repeat, [sprintf('%d*%s', ...
dblock(segNum).Repeat,dblock(segNum).Types) '=>int16'], sum(blockBytes) - blockBytes(segNum), mf);
end
function seg = segRead(segName, mf)
if nargin == 1, mf = 'l'; end
segNum = find(strcmp({dblock.Str},segName));
fseek(fid, sum(blockBytes(1:segNum-1))+NUM_HEADER_BYTES, 'bof');
seg = fread(fid, numIdx*dblock(segNum).Repeat, sprintf('%d*%s', ...
dblock(segNum).Repeat,dblock(segNum).Types), sum(blockBytes) - blockBytes(segNum), mf);
end
end
function info = getHeader(hdr)
eval(char(hdr'));
info.header = header;
end