function obj = trace(data_src, dt, dy, id, props)
% trace - Load a generic data trace. It can be membrane voltage, current, etc.
%
% Usage:
% obj = trace(data_src, dt, dy, id, props)
%
% Description:
% This object is designed to recognize most data file formats. See the
% data_src parameter below. Traces for specific experimental or simulation
% protocols can extend this class for adding new parameters.
%
% Parameters:
% data_src: Trace data as a column vector OR name of a data file generated by either
% Genesis (.bin, .gbin, .genflac), Neuron, PCDX (.all), Matlab (.mat),
% NeuroShare (.nsn, .nev, .stb, .plx, .nex, .map, .son, .smr,
% .mcd), ASCII (.txt) or Axoclamp (.abf) files. As last
% resort, file will be loaded with Matlab's load() function.
% dt: Time resolution in [s], unless specified in HDF5 or NeuroShare
% file. For Neuron ASCII files, it is used as the file data units.
% dy: y-axis resolution in [ISI (V, A, etc.)], unless specified in HDF5 or NeuroShare file.
% id: Identification string
% props: A structure with any optional properties.
% scale_y: Y-axis scale to be applied to loaded data.
% offset_y: Y-axis offset to be added to loaded and scaled data.
% unit_y: Unit of Y-axis as in 'V' or 'A' (default='V').
% y_label: String to put on Y-axis of plots.
% trace_time_start: Samples in the beginning to discard [dt]
% baseline: Resting potential.
% channel: Channel to read from file Genesis, PCDX, NeuroShare or
% Neuron file, or column in a data vector.
% numTraces: Divide the single column vector of data into this
% many columns by making it a matrix.
% file_type: Specify file type instead of guessing from extension:
% 'genesis': Raw binary files created with Genesis disk_out method.
% 'genesis_flac': Compressed Genesis binary files.
% 'neuron': Binary files created with Neuron's Vector.vwrite method.
% 'neuronascii': Ascii files created from Neuron's Vector objects.
% Uses time step in file to scale given dt (Must be in ms).
% 'pcdx': .ALL data acquisition files from PCDX program.
% 'matlab': Matlab .MAT binary files with matrix data.
% 'neuroshare': One of the vendor formats recognized by
% NeuroShare Windows DLLs. See above and http://neuroshare.org. A
% scale_y value may need to be supplied to get the correct units.
% 'abf': AxoClamp .ABF format read with abf2load from
% Matlab FileExchange.
% 'txt': Generic, space-separated ASCII file.
% file_endian: 'l' for little endian and 'b' for big endian
% (default='n', for native endian). See machineformat option
% in fopen for more info.
% traces: Trace numbers as a numeric array or as a string with
% numeric ranges (e.g., '1 2 5-10 28') for PCDX files.
% spike_finder: Method of finding spikes
% (1 for findFilteredSpikes, 2 for Li Su's
% findspikes, and 3 for Alfonso Delgado Reyes's
% findspikes_old). Methods 2 and 3 require a threshold.
% threshold: Spike finding threshold. For the findspikes method,
% it is either a scalar, or [thres1 thres2] to define
% a range. For findFilteredSpikes it is used on the
% filtered data and the default is 2/3 max amplitude
% of band-passed data, but with a minimum of 15.
% downThreshold: (Only for findFilteredSpikes) Size of the trough
% after the spike peak in filtered data (Default=-2).
% minInit2MaxAmp, minMin2MaxAmp: For spike_shape elimination,
% conditions of minimal allowed values for
% initial point to max point and minimal point to
% max point, respectively (Default=10 for both).
% init_Vm_method: Method of finding spike thresholds during spike
% shape calculation (see spike_shape/spike_shape).
% init_threshold: Spike initiation threshold (deriv or accel).
% (see above methods and implementation in calcInitVm)
% init_lo_thr, init_hi_thr: Low and high thresholds for slope.
% custom_filter: Recommended if sampling rate differs appreciably from 10 kHz.
% If custom_filter == 1, a filter with custom lowpass and highpass
% cutoffs can be specified. This allows for fast and accurate spike
% discrimination. The filter type used is a 2-pole
% butterworth, different than the default high-order
% Cheby2. Creates new prop called 'butterWorth' to
% hold the filter.
% lowPassFreq: If set, it sets a new low pass cutoff for custom filter. Default is 3000Hz
% highPassFreq: If set it sets a new high pass cutoff for custom filter. Default is 50 Hz
% quiet: If 1, reduces the amount of textual description in plots
% and does not add information to id field.
%
% Returns a structure object with the following fields:
% data: The trace column matrix.
% dt, dy, id, props (see above)
%
% General operations on trace objects:
% trace - Construct a new trace object.
% plot - Graph the trace.
% display - Returns and displays the identification string.
% get - Gets attributes of this object and parents.
% subsref - Allows usage of . operator.
% calcAvg - Calculate the average value of the trace period.
% calcMin - Calculate the minimum value of the trace period.
% calcMax - Calculate the minimum value of the trace period.
% periodWhole - Return the whole period of this trace.
% findFilteredSpikes - Use a band-pass filter to smooth the data and
% find spikes afterwards.
% getResults - Calculates a set of tests.
% spike_shape - Build a trace of the average spike shape in here.
% spikes - Build a spikes object with the spikes found here.
%
% Converter methods:
% spikes - Find the spikes and construct a spikes object.
% spike_shape - Construct a spike_shape object from this trace.
%
% Additional methods:
% See methods('trace')
%
% See also: spikes, spike_shape, cip_trace, period
%
% $Id$
%
% Author: Cengiz Gunay <cgunay@emory.edu>, 2004/07/30
% Modified:
% - allow custom filter, Thomas D. Sangrey 2007/12/04
% Copyright (c) 2007 Cengiz Gunay <cengique@users.sf.net>.
% This work is licensed under the Academic Free License ("AFL")
% v. 3.0. To view a copy of this license, please look at the COPYING
% file distributed with this software or visit
% http://opensource.org/licenses/afl-3.0.php.
vs = warning('query', 'verbose');
verbose = strcmp(vs.state, 'on');
if nargin == 0 % Called with no params
obj.data = [];
obj.dt = 1;
obj.dy = 1;
obj.id = '';
obj.props = struct([]);
obj = class(obj,'trace');
elseif isa(data_src,'trace') % copy constructor?
obj = data_src;
else
if ~ exist('props','var')
props = struct;
end
if ~ exist('id','var')
id = '';
end
if isa(data_src, 'char') % filename?
[path, filename, ext] = fileparts(data_src);
ext = lower(ext); % Case insensitive matches for file extension
% Default if native endian of this computer
if ~ isfield(props, 'file_endian')
[c,maxsize,endian] = computer;
props.file_endian = endian;
end
% if file type not specified, use file extension to guess it
if ~ isfield(props, 'file_type')
if strcmpi(ext, '.bin') || strcmpi(ext, '.gbin') % Genesis file
props.file_type = 'genesis';
else
props.file_type = '';
end
end
if strcmpi(props.file_type, 'genesis')
channel = 1; % by default
if isfield(props, 'channel')
channel = props.channel;
end
if ~ isempty(findstr(filename, '_BE_')) || ...
~ isempty(findstr(filename, '_BE.')) || ...
strcmpi(props.file_endian, 'b')
% check if readgenesis.mex* file available
if exist('readgenesis_BE') == 3
% Use big-endian (Mac, Sun) version of readgenesis
data = readgenesis_BE(data_src, channel);
else
% otherwise use native Matlab reader
[data times] = readgenbin(data_src, NaN, NaN, props.file_endian);
data = data(:, channel);
dt = diff(times(1:2)) * 1e-3;
end
else
% check if readgenesis.mex* file available
if exist('readgenesis') == 3
% Use regular (i386 PCs) little-endian version of readgenesis
data = readgenesis(data_src, channel);
else
% otherwise use native Matlab reader
[data times] = readgenbin(data_src, NaN, NaN, props.file_endian);
data = data(:, channel);
dt = diff(times(1:2)) * 1e-3;
end
end
elseif strcmpi(props.file_type, 'genesis_flac') || ...
strcmpi(ext, '.genflac') % Compressed 16-bit genesis file
channel = 1; % by default
if isfield(props, 'channel')
channel = props.channel;
end
data = readgenflac(data_src);
data = data(:, channel);
elseif strcmpi(props.file_type, 'neuron')
[data err] = readNeuronVecBin(data_src, props.file_endian);
if isempty(strfind(err, 'Success'))
error([ 'Failed to load Neuron binary ''' data_src ''' because of ' ...
'error: ' err ]);
end
channel = ':'; % by default
if isfield(props, 'channel')
channel = props.channel;
end
data = data(:, channel);
elseif strcmpi(props.file_type, 'neuronascii')
[data, label] = readNeuronVecAscii(data_src);
dt = dt * diff(data(1:2, 1)); % scale given dt by delta in first column
data = data(:, 2);
id = [ id label];
elseif strcmpi(props.file_type, 'pcdx') || ...
strcmpi(ext, '.all') % PCDX file
%disp('Loading PCDX trace');
data = loadtraces(data_src, props.traces, props.channel, 1);
elseif strcmpi(props.file_type, 'matlab') || ...
strcmpi(ext, '.mat') % MatLab file
s = load(data_src);
fields = fieldnames(s);
data = getfield(s, fields{1}); % Assuming there's only one vector
elseif strcmpi(props.file_type, 'hdf5') || ...
strcmpi(ext, '.hdf5')
% new neurosage file
% reset file type explicitly
props.file_type = 'hdf5';
s1 = ns_read(props.AcquisitionData{props.channel});
% Make sure the data is in a column vector
if size(s1.Y, 1) > size(s1.Y, 2)
data = s1.Y;
else
data = s1.Y';
end
% NeuroShare files with these extensions:
% nsn, .nev, .stb, .plx, .nex, .map, .son, .smr, .mcd
elseif strcmpi(props.file_type, 'neuroshare') || ...
strcmpi(ext, '.nsn') || strcmpi(ext, '.nev') || strcmpi(ext, '.stb') || ...
strcmpi(ext, '.plx') || strcmpi(ext, '.nex') || strcmpi(ext, '.map') || ...
strcmpi(ext, '.son') || strcmpi(ext, '.smr') || strcmpi(ext, ...
'.mcd')
[data, dt, id, props] = trace_loadns(data_src, props);
% assume units in mV? Nothing specified in Matlab loader
if isempty(dy), dy = 1e-3; end
elseif strcmpi(ext, '.abf') || strcmpi(props.file_type, 'abf')
% Load AxoClamp ABF format using MathWorks fileExchange entry
% abf2load
[dt, data, y_units, dy, cell_name] = loadABF(data_src, props);
props.unit_y = y_units;
if ~ isfield(props, 'quiet')
id = [ cell_name id ];
end
elseif strcmpi(ext, '.txt') || strcmpi(props.file_type, 'txt')
data = load(data_src);
data = data(:, getFieldDefault(props, 'channel', ':'));
else
warning(['No matching load function found for file ''' data_src ...
''' or specified type ''' ...
props.file_type '''; defaulting to Matlab load() function.']);
s = load(data_src);
fields = fieldnames(s);
data = getfield(s, fields{1}); % Assuming there's only one vector
data = data(:, getFieldDefault(props, 'channel', ':'));
end
% use the filename as id unless otherwise specified
if ~ exist('id','var') || isempty(id)
id = name;
end
elseif isnumeric(data_src)
data = data_src;
if size(data_src, 2) > size(data_src, 2)
data = data'; % convert to column vector
end
channel = getFieldDefault(props, 'channel', ':');
if ~ strcmp(channel, ':')
data = data(:, channel);
id = [id ', chan ' num2str(channel) ];
end
else
error('Unrecognized data source!');
end
% reshape data if requested
if isfield(props, 'numTraces')
total_len = size(data, 1);
if mod(total_len, props.numTraces) ~= 0
error(['props.numTraces=' num2str(props.numTraces) ...
' does not evenly divide data vector of size=' ...
num2str(size(data, 1))]);
end
data = reshape(data, total_len/props.numTraces, props.numTraces);
end
% Scale the loaded data if desired
if isfield(props, 'scale_y')
% scale HDF5 data only if gain is unspecified in file metadata (by Li Su)
if ~( (isfield(props, 'file_type') && strcmpi(props.file_type, 'hdf5')) ) || ...
(isfield(props, 'AcquisitionData') && ...
strcmpi( ...
props.AcquisitionData{props.channel}.PhysicalDevice, ...
'Unspecified' ...
) ...
)
data = props.scale_y * data;
if verbose
warning(sprintf(['Device unspecified. Trace multiplied by given ' ...
'gain %g'], props.scale_y ));
end
end
end
% Apply offset to data if desired (after scaling?)
if isfield(props, 'offset_y')
data = data + props.offset_y;
end
% Crop the data if desired
if isfield(props, 'trace_time_start')
% data = data(:, props.trace_time_start:end);
data = data(props.trace_time_start:end, :);
end
% Custom filter props
if isfield(props, 'custom_filter') && props.custom_filter == 1
if isfield(props, 'highPassFreq')
hpf = props.highPassFreq;
else
hpf = 50;
end
if isfield(props, 'lowPassFreq')
lpf = props.lowPassFreq;
else
lpf =3000;
end
% make a new filter based on the given dt
% multiply by 2 to find nyquist freq
assert(all(2*[lpf hpf]*dt <= 1), ...
[ 'For custom filter, highPassFreq and lowPassFreq must be <= ' ...
'Nyquist freq (' num2str(1/dt/2) ').' ]);
[b,a] = butter(2,2*hpf*dt,'high');
butterWorth.highPass = struct('b', b, 'a', a);
[b,a] = butter(2,2*lpf*dt,'low');
butterWorth.lowPass = struct('b', b, 'a', a);
props.butterWorth = butterWorth;
end
obj.data = data;
obj.dt = dt;
obj.dy = dy;
obj.id = id;
obj.props = props;
obj = class(obj, 'trace');
end