-
Notifications
You must be signed in to change notification settings - Fork 1
/
find_pulse_response_endpoints.m
executable file
·377 lines (315 loc) · 15.3 KB
/
find_pulse_response_endpoints.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
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
function [indResponseStart, indResponseEnd, hasJump, ...
indBeforePulseStart, indBeforePulseEnd] = ...
find_pulse_response_endpoints (vectors, siMs, varargin)
%% Returns the start and end indices of the first pulse response (from pulse start to 20 ms after pulse ends by default) from vector(s)
% Usage: [indResponseStart, indResponseEnd, hasJump, ...
% indBeforePulseStart, indBeforePulseEnd] = ...
% find_pulse_response_endpoints (vectors, siMs, varargin)
%
% Explanation:
% TODO
%
% Examples:
% load_examples;
% [idxStart, idxEnd] = find_pulse_response_endpoints(myPulseResponse1a, siMs)
% [idxStart, idxEnd] = find_pulse_response_endpoints(myPulseResponse1a, siMs, 'ResponseLength', 0)
%
% Outputs:
% indResponseStart - indices of pulse response start
% specified as a positive integer (or NaN) vector
% indResponseEnd - indices of pulse response end
% specified as a positive integer (or NaN) vector
% hasJump - whether there was a significant "jump" detected
% specified as a logical vector
% indBeforePulseStart - indices of pulse start
% specified as a positive integer (or NaN) vector
% indBeforePulseEnd - indices of pulse end
% specified as a positive integer (or NaN) vector
%
% Arguments:
% vectors - vector(s) containing a pulse response
% Note: If a cell array, each element must be a vector
% If a non-vector array, each column is a vector
% must be a numeric array or a cell array of numeric vectors
% siMs - sampling interval(s) in ms
% must be a positive vector
% varargin - 'PulseVectors': vectors that contains the pulse itself
% Note: If a cell array, each element must be a vector
% If a non-vector array, each column is a vector
% must be a numeric array or a cell array of numeric vectors
% default == [] (not used)
% - 'SameAsPulse': whether always the same as
% the current pulse endpoints
% must be numeric/logical 1 (true) or 0 (false)
% default == true
% - 'ResponseLengthMs': length of the pulse response
% after pulse endpoint in ms
% must be a nonnegative scalar
% default = 20 ms
% - 'BaselineLengthMs': length of the pulse response
% before pulse start in ms
% must be a nonnegative scalar
% default = 0 ms
%
% Requires:
% cd/create_time_vectors.m
% cd/find_first_jump.m
% cd/find_pulse_endpoints.m
% cd/fit_first_order_response.m TODO
% cd/iscellnumeric.m
% cd/match_format_vector_sets.m
%
% Used by:
% cd/compute_initial_slopes.m
% cd/filter_and_extract_pulse_response.m
% cd/parse_pulse_response.m
% File History:
% 2018-08-13 Adapted from compute_average_initial_slopes.m
% 2018-09-17 Changed required arguement tVecCpr to siMs
% 2018-09-17 Added optional parameters SameAsPulse and ResponseLengthMs
% 2018-09-23 Added the optional parameter baselineLengthSamples
% 2018-10-09 Renamed isUnbalanced -> hasJump and improved documentation
% 2018-12-15 Now returns NaN if there is no pulse
% 2018-12-15 Now prevents the endpoints from exceeding bounds
% 2018-12-17 Now allows siMs to be a vector
% 2019-09-09 Now allows PulseVectors to be omitted
%% Hard-coded parameters
nSamplesPerJump = 2; % number of samples apart for calculating jump
signal2Noise = 10; % signal to noise ratio for a jump
noiseWindowSize = 5; % number of samples to measure baseline noise
%% Default values for optional arguments
pulseVectorsDefault = []; % don't use pulse vectors by default
sameAsPulseDefault = true; % use pulse endpoints by default
responseLengthMsDefault = 20; % a response of 20 ms by default
baselineLengthMsDefault = 0; % don't include a baseline by default
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Deal with arguments
% Check number of required arguments
if nargin < 2
error(create_error_for_nargin(mfilename));
end
% Set up Input Parser Scheme
iP = inputParser;
iP.FunctionName = mfilename;
% Add required inputs to the Input Parser
addRequired(iP, 'vectors', ... % vectors
@(x) assert(isnumeric(x) || iscellnumeric(x), ...
['vectors must be either a numeric array', ...
'or a cell array of numeric arrays!']));
addRequired(iP, 'siMs', ...
@(x) validateattributes(x, {'numeric'}, {'positive', 'vector'}));
% Add parameter-value pairs to the Input Parser
addParameter(iP, 'PulseVectors', pulseVectorsDefault, ...
@(x) isnumeric(x) || iscell(x) && all(cellfun(@isnumeric, x)) );
addParameter(iP, 'SameAsPulse', sameAsPulseDefault, ...
@(x) validateattributes(x, {'logical', 'numeric'}, {'binary'}));
addParameter(iP, 'ResponseLengthMs', responseLengthMsDefault, ...
@(x) validateattributes(x, {'numeric'}, {'nonnegative', 'scalar'}));
addParameter(iP, 'BaselineLengthMs', baselineLengthMsDefault, ...
@(x) validateattributes(x, {'numeric'}, {'nonnegative', 'scalar'}));
% Read from the Input Parser
parse(iP, vectors, siMs, varargin{:});
pulseVectors = iP.Results.PulseVectors;
sameAsPulse = iP.Results.SameAsPulse;
responseLengthMs = iP.Results.ResponseLengthMs;
baselineLengthMs = iP.Results.BaselineLengthMs;
%% Preparation
% Match up pulseVectors with vectors and make sure they are both cell arrays
[pulseVectors, vectors] = ...
match_format_vector_sets(pulseVectors, vectors, 'ForceCellOutputs', true);
% Make sure numel(siMs) is not greater than numel(vectors)
if numel(siMs) > numel(vectors)
error('Number of sampling intervals cannot exceeed the number of vectors!');
end
% Convert siMs to a cell array and match up size with vectors
[siMs, vectors] = ...
match_format_vector_sets(num2cell(siMs), vectors);
%% Do the job
[indResponseStart, indResponseEnd, hasJump, indBeforePulseStart, indBeforePulseEnd] = ...
cellfun(@(x, y, z) fpre_helper(x, y, z, sameAsPulse, responseLengthMs, ...
baselineLengthMs, nSamplesPerJump, signal2Noise, noiseWindowSize), ...
vectors, siMs, pulseVectors);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [idxResponseStart, idxResponseEnd, hasJump, ...
idxPulseStart, idxPulseEnd] = ...
fpre_helper (vector, siMs, pulseVector, sameAsPulse, ...
responseLengthMs, baselineLengthMs, ...
nSamplesPerJump, signal2Noise, noiseWindowSize)
% Compute the total number of samples
nSamples = length(vector);
% Compute the length of the pulse response in samples
cprLengthSamples = floor(responseLengthMs / siMs);
% Compute the baseline length in samples
baselineLengthSamples = floor(baselineLengthMs / siMs);
% Find the start and end points of the pulse
if ~isempty(pulseVector)
% If provided, use the trace containg the pulse
[idxPulseStart, idxPulseEnd] = find_pulse_endpoints(pulseVector);
% If an index is not detected, return warning message
if isnan(idxPulseStart) || isnan(idxPulseEnd)
fprintf('The pulse could not be detected!\n');
idxResponseStart = NaN;
idxResponseEnd = NaN;
hasJump = false;
idxPulseStart = NaN;
idxPulseEnd = NaN;
return
end
else
% If not, estimate by fitting to a moving pulse_response,
% then look for inflection points
% Construct a time vector
timeVec = create_time_vectors(nSamples, 'SamplingIntervalMs', siMs, ...
'TimeUnits', 'ms');
% Shift the vector to start at 0
vectorShifted = vector - vector(1);
% TODO: Pull out as fit_first_order_response.m
%%
% Compute the total duration
totalDuration = timeVec(end) - timeVec(1) + siMs;
% Estimate amplitude
[absAmplitudeEstimate, idxMaxAbs] = max(abs(vectorShifted));
amplitudeEstimate = vectorShifted(idxMaxAbs);
% Find the first point that reached at least 1/4 of amplitude
idxFirstDip = find(abs(vectorShifted) > absAmplitudeEstimate / 4, ...
1, 'first');
timeFirstDip = idxFirstDip * siMs;
% Find the last point that reached at least 3/4 of amplitude
idxFirstReturn = find(abs(vectorShifted) > absAmplitudeEstimate * 3 / 4, ...
1, 'last');
% Estimate the duration
durationEstimate = (idxFirstReturn - idxFirstDip) * siMs;
tauEstimate = durationEstimate;
% Set up fitting parameters for a first order response
[eqForm, aFittype, coeffInit, coeffLower, coeffUpper] = ...
fit_setup_first_order_response('TotalDuration', totalDuration, ...
'AmplitudeEstimate', amplitudeEstimate, ...
'TauEstimate', tauEstimate, ...
'DurationEstimate', durationEstimate, ...
'DelayEstimate', timeFirstDip);
% Fit the first order response
[fitObject, goodnessOfFit, algorithmInfo] = ...
fit(timeVec, vectorShifted, aFittype, 'StartPoint', coeffInit, ...
'Lower', coeffLower, 'Upper', coeffUpper);
% Parse the results
fitResults = parse_fitobject(fitObject);
%%
% Extract fitted parameters
coeffNames = fitResults.coeffNames;
coeffValues = fitResults.coeffValues;
duration = coeffValues(strcmp(coeffNames, 'c'));
delay = coeffValues(strcmp(coeffNames, 'd'));
% Get the pulse times
timePulseStart = delay;
timePulseEnd = delay + duration;
% Get the closest pulse indices
idxPulseStart = floor(timePulseStart / siMs);
idxPulseEnd = floor(timePulseEnd / siMs);
end
% Use windows straddling the the start/end points of the current pulse
% as regions for finding the start/end points of the current pulse response
% Note: this will always cause idxPulseStart/idxPulseEnd to be the first
% index for checking, and will check noiseWindowSize more points
idxRegion1Start = max([idxPulseStart - noiseWindowSize, 1]);
idxRegion1End = min([idxPulseStart + noiseWindowSize, nSamples]);
responseRegion1 = vector(idxRegion1Start:idxRegion1End);
idxRegion2Start = max([idxPulseEnd - noiseWindowSize, 1]);
idxRegion2End = min([idxPulseEnd + noiseWindowSize, nSamples]);
responseRegion2 = vector(idxRegion2Start:idxRegion2End);
% Initialize hasJump as false
hasJump = false;
% Find the start point of the current pulse response
% by detecting the 'first jump' in the region of interest #1
% If it doesn't exist, use the start point of the current pulse
[~, idxTemp1] = ...
find_first_jump(responseRegion1, 'NSamplesPerJump', nSamplesPerJump, ...
'Signal2Noise', signal2Noise, ...
'NoiseWindowSize', noiseWindowSize);
if ~isempty(idxTemp1) && ~sameAsPulse
idxResponseStart = (idxRegion1Start - 1) + idxTemp1 - baselineLengthSamples;
hasJump = true;
else
idxResponseStart = idxPulseStart - baselineLengthSamples;
end
% Find the end point of the current pulse response
% by detecting the 'first jump' in the region of interest #2
% If it doesn't exist, use the end point of the current pulse
% plus cprLengthSamples
[~, idxTemp2] = ...
find_first_jump(responseRegion2, 'NSamplesPerJump', nSamplesPerJump, ...
'Signal2Noise', signal2Noise, ...
'NoiseWindowSize', noiseWindowSize);
if ~isempty(idxTemp2) && ~sameAsPulse
idxResponseEnd = (idxRegion2Start - 1) + idxTemp2 + cprLengthSamples;
hasJump = true;
else
idxResponseEnd = idxPulseEnd + cprLengthSamples;
end
% Make sure the indices are within bounds
idxResponseStart = max([idxResponseStart, 1]);
idxResponseEnd = min([idxResponseEnd, nSamples]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%{
OLD CODE:
% Compute slope right after current pulse start
idxFirst1 = idxPulseStart;
idxLast1 = idxPulseStart + nSamples - 1;
startSlope1 = compute_slope(tvecCpr, vectors, idxFirst1, idxLast1);
% Compute slope right after current pulse end
idxFirst2 = idxPulseEnd;
idxLast2 = idxPulseEnd + nSamples - 1;
endSlope2 = compute_slope(tvecCpr, vectors, idxFirst2, idxLast2);
% Compute slope right after current pulse start
idxFirst3 = idxCpStart2;
idxLast3 = idxCpStart2 + nSamples - 1;
startSlope3 = compute_slope(tvecCpr, vectors, idxFirst3, idxLast3);
% Compute slope right after current pulse end
idxFirst4 = idxCpEnd2;
idxLast4 = idxCpEnd2 + nSamples - 1;
endSlope4 = compute_slope(tvecCpr, vectors, idxFirst4, idxLast4);
% Choose the more negative of the start slopes
% and the more positive of the end slopes
startSlope = min([startSlope1, startSlope3]);
endSlope = max([endSlope2, endSlope4]);
addRequired(iP, 'nSamples', ...
@(x) validateattributes(x, {'numeric'}, {'scalar'}));
parse(iP, tvecCpr, vectors, pulseVector, nSamples);
function [avgSlope, startSlope, endSlope, indsUsedForPlot] = find_pulse_response_endpoints (tvecCpr, vectors, pulseVector, nSamples)
% Note: function calls are slower
% /home/Matlab/Adams_Functions/compute_slope.m
startSlope = compute_slope(tvecCpr, vectors, idxFirst1, idxLast1);
endSlope = compute_slope(tvecCpr, vectors, idxFirst2, idxLast2);
% Crop the voltage trace
vvecCropped = vectors((idxResponseStart + 1):end);
% tvecCpr - time vector of the current pulse response in ms
% must be a numeric vector
addRequired(iP, 'tvecCpr', ...
@(x) validateattributes(x, {'numeric'}, {'vector'}));
parse(iP, tvecCpr, vectors, varargin{:});
@(x) validateattributes(x, {'numeric'}, {'vector'}));
% varargin - 'PulseVector': vector that contains the pulse itself
% must be a numeric vector
% default == [] (not used)
pulseVectorDefault = []; % don't use pulse vector by default
addParameter(iP, 'PulseVector', pulseVectorDefault, ...
@(x) validateattributes(x, {'numeric'}, {'vector'}));
pulseVector = iP.Results.PulseVector;
idxResponseStart = zeros(nVectors, 1);
idxResponseEnd = zeros(nVectors, 1);
hasJump = zeros(nVectors, 1);
idxPulseStart = zeros(nVectors, 1);
idxPulseEnd = zeros(nVectors, 1);
% cd/match_array_counts.m
[pulseVectors, vectors] = ...
match_array_counts(pulseVectors, vectors, 'ForceCellOutputs', true);
if isempty(idxPulseStart) || isempty(idxPulseEnd)
idxResponseStart = [];
idxResponseEnd = [];
hasJump = [];
idxPulseStart = [];
idxPulseEnd = [];
end
if ~isempty(idxTemp1) && ~sameAsPulse
if ~isempty(idxTemp2) && ~sameAsPulse
%}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%