-
Notifications
You must be signed in to change notification settings - Fork 0
/
az_detect_events.m
executable file
·188 lines (159 loc) · 6.56 KB
/
az_detect_events.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
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
function [event, hdr] = az_detect_events(varargin)
% AZ_DETECT takes the raw binary recorder files and detects discontinuous
% events due to either trigger events, multiple recordings, or data corruption
%
% event = az_detect_events(FNAME1,FNAME2) returns a 1xM struct holding
% the start/stop sample numbers and times for M triggered events
%
% event = az_detect_events(PREFIX) uses the string PREFIX to specify the
% SRZ file pair
%
% [event, hdr] = az_detect_events(..) also returns the header field
% information for both files
%
% event = az_detect_events(FNAME1,FNAME2,'eventfile') saves event struct
% to the file 'eventfile.mat'
%
% event = az_detect_events(FNAME1,FNAME2,'eventfile','hdrfile') also
% saves the header field information for all samples in 'hdrfile.mat'
%
fprintf('\n\n***********************************************\n')
fprintf('Detecting triggered events in data files\n')
switch nargin
case 1
prefix = varargin{1};
if ~strcmp(prefix(end),'_')
prefix(end+1) = '_';
end
fname1 = [prefix 'side1.srz'];
fname2 = [prefix 'side2.srz'];
case 2
fname1 = varargin{1};
fname2 = varargin{2};
case 3
fname1 = varargin{1};
fname2 = varargin{2};
eventfile = varargin{3};
case 4
fname1 = varargin{1};
fname2 = varargin{2};
eventfile = varargin{3};
hdrfile = varargin{4};
otherwise
error('Incorrect number of input parameters')
end
% get filename prefix from current filename
if ~exist('prefix','var')
prefix = regexp(fname1,'[_\-\ ]');
prefix = fname1(1:prefix(end));
end
% assign eventfile name if not specified
if ~exist('eventfile','var')
eventfile = [prefix 'event.mat'];
end
% assign hdrfile name if not specified
if ~exist('hdrfile','var')
hdrfile = [prefix 'hdr.mat'];
end
% verify files exist
if ~exist(fname1,'file')
error('AZ_DETECT_EVENTS:fnf', 'Could not locate file "%s"', fname1)
end
if ~exist(fname2,'file')
error('AZ_DETECT_EVENTS:fnf', 'Could not locate file "%s"', fname2)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read header content
fprintf('\nReading header information from side 1... \n')
hdr(1) = read_header(fname1);
fprintf('\nReading header information from side 2... \n')
hdr(2) = read_header(fname2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% verify/correct data set alignment, if sync data is available
if isfield(hdr,'sync')
fprintf('\nSynchronizing data sequences... \n')
% match synchronous trigger events to event index in both files
[match,idx(1,:),idx(2,:)] = intersect(hdr(1).sync, hdr(2).sync);
N = numel(match);
M = numel(union(hdr(1).syn, hdr(2).sync));
fprintf(' Found %d of %d synchronous trigger events.\n', N, M);
else
% otherwise, just use first valid block for now and manually align data segments :(
fprintf('\nNo synchronization data was found on auxiliary channels! Assuming perfect alignment exists\n\t(correct this manually by adjusting event index).\n')
% determine number of events in each side
N = [numel(hdr(1).event) numel(hdr(2).event)]-1;
% assign event list and indices
idx = [1:N(1) nan(1,N(2)-N(1)) ; 1:N(2) nan(1,N(1)-N(2))];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% assign event samples and times to struct
event(1:max(N)) = struct('s0',[0 0],'s1',[0 0],'t0',[0 0],'t1',[0 0]);
for n=1:max(N)
event(n).eNum = n; % assign an event number for reference
% repeat for each side
for m=1:2
if (n < numel(hdr(m).event))
event(n).s0(m) = hdr(m).event(idx(m,n));
event(n).s1(m) = hdr(m).event(idx(m,n)+1)-1;
event(n).t0(m) = hdr(m).time(event(n).s0(m));
event(n).t1(m) = hdr(m).time(event(n).s1(m));
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% save results to file
if ~isempty(eventfile)
fprintf('\nSaving event index to "%s"...\n',eventfile)
save(eventfile,'event');
end
if ~isempty(hdrfile)
fprintf('Saving header information to "%s"...\n',hdrfile)
save(hdrfile,'hdr');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Subfunction to read and process header info
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function hdr = read_header(fname)
res = read_SRZ_header_field(fname,[],[2 3 5 7 9]); % can this result be converted directly to (u)int32?
% assign to hdr struct with appropriate names
hdr.count = res(:,1); % Mx1 uint32
hdr.clock = res(:,2); % Mx1 uint32
hdr.gain = res(:,3); % 1x1 uint8
hdr.fs = res(:,4); % 1x1 double
hdr.aux = res(:,5); % Mx1 double
% unwrap clock cycles if overflow occurred
k = find(diff(hdr.clock) < 1); k = [k; length(hdr.clock)];
for kNum = 1:length(k)-1
idx = k(kNum)+1:k(kNum+1);
hdr.clock(idx) = hdr.clock(idx) + kNum*2^32;
end
% convert clock cycles to time relative to start
hdr.time = (hdr.clock - hdr.clock(1)) ./ 106.25e6;
% look for discontinuous data blocks
hdr.block = [1; find(diff(hdr.count)~=1)+1];
fprintf('Found %d data blocks in "%s"\n',numel(hdr.block),fname)
% look for trigger events by difference in clock time - append last sample
hdr.event = [1; find(diff(hdr.clock) > 450)+1; numel(hdr.count)+1];
fprintf('Found %d data segments (trigger events) in "%s"\n', length(hdr.event)-1, fname);
% find index to the first event in each data block
hdr.blockidx = zeros(numel(hdr.block),1);
for i=1:numel(hdr.block)
res = find(hdr.event == hdr.block(i),1);
assert(~isempty(res), 'Start of data block not aligned with any data segment!'); % THIS SHOULD NEVER HAPPEN
hdr.blockidx(i) = res;
end
hdr.blockidx(end+1,1) = numel(hdr.event); % append last segment index
%%% use FFT to look for 10kHz digital data stream on auxiliary channel
% extract digital data from auxiliary channel, if possible
% hdr.aux = hdr.aux - mean(hdr.aux); % subtract DC offset
% % evaluate sync for
% hdr.sync = zeros(size(hdr.event));
%
% % iterate over each event
% for i=1:numel(hdr.event)
% % compare level with 0 to recover DC level square wave
% % synchronize to start pulse - remove prior data
% % sample digital level every X seconds
% % convert the digital serial word to a trigger count
% hdr.sync(i) = res;
% end