/MatlabCode/branches/Greg's Branch/Analysis/nex2stro/nex2stro.m
MATLAB | 481 lines | 316 code | 71 blank | 94 comment | 41 complexity | 7981a80205a1c6b8faf587e93a7cad7d MD5 | raw file
Possible License(s): GPL-2.0, AGPL-1.0
- function stro = nex2stro(fileName)
- %
- % EXAMPLE: stro = nex2stro(fileName);
- %
- %This engine will take nex files and return a (S)ummary
- %(T)rials (R)aster (O)ther (stro) structure. 'filename' should be the
- %complete path to the nex file.
- %
- % cah 1/21/08 started development...
- % gdlh/cah 2/14/08 replaced global variable architecture with nested
- % functions.
-
- %start the wait dialog and, initialize all the appropriate structures,
- %then open a .nex file if no input argument was given.
- waitdlg = waitbar(0, 'Opening a .nex file');
- [stro, ecodes, spikes, anlg, waves, events] = setupStructures();
- if (nargin == 0 || isempty(fileName))
- [fname, pathname] = uigetfile('*.nex', 'Select a NeuroExplorer file');
- stro.sum.fileName = strcat(pathname, fname);
- else
- stro.sum.fileName = fileName;
- end
- % Unpack the .nex file. No arguments are needed b/c all the local
- %variables are accessible to the nested functions.
- unpackNexFile();
- cleanupAnalogSignals(); %determine if the start times are identical for each channel
- %determine which paradigm file to use, then eval the appropriate header
- %file. Be sure to return all the appropriate variables (for scoping
- %purposes). the feval statement will crash if one of the outputs is
- %undefined in paradigmFile
- paradigmID = dat2num(ecodes, 8999, 4000, 'int', 0);
- stro.sum.paradigmID = paradigmID{1}(1);
- paradigmFile = paradigmLibrary(stro.sum.paradigmID);
- [C, trialParamsHeader, exptParamsHeader, additionalRasterFields, badLogic, otherInstructions] = feval(paradigmFile);
- %Fill up some of the stro.sum contents
- getExperimentalParams; %iterates over the exptParamsHeader reconstructing those elements
- stro.sum.rasterCells = makeRasterHeader(additionalRasterFields); %defined in the 'paradigmFile'.
- stro.sum.trialFields = trialParamsHeader; %defined in the 'paradigmFile'
- stro.sum.otherFields = otherInstructions;
- stro.sum.analog.sigid = anlg.names;
- stro.sum.analog.storeRates = anlg.ADFrequency;
- stro.sum.analog.ADtoMV = anlg.ADtoMV;
- stro.sum.waves.wf_id = waves.names;
- stro.sum.waves.storeRates = waves.WFrequency;
- %iterate by trials filling up the .raster and .trial fields of stro.
- nTrials = length(events.start);
- stro.ras = cell(nTrials, length(stro.sum.rasterCells));
- stro.trial = nan(nTrials, size(trialParamsHeader, 2));
- stro.sum.absTrialNum = nan(nTrials,1);
- goodTrial = 0; %the counter of good trials
- for b = 1:nTrials;
- trialCodes = checkTrial(events.start(b), events.stop(b), badLogic);
- if trialCodes;
- goodTrial = goodTrial + 1;
- rasterByTrial(events.start(b), events.stop(b), goodTrial, b, additionalRasterFields, trialCodes);
- indexByTrial(trialCodes, trialParamsHeader, goodTrial, C.VALOFFSET);
- otherByTrial(trialCodes, otherInstructions, goodTrial, C);
- stro.sum.absTrialNum(goodTrial) = b; %a mapping b/w actual and good trial nums
- end
- if ~(rem(b, 3)) %update the waitdlg every three trials
- waitbar(b./nTrials, waitdlg, ' ...Processing trials');
- end
- end
- %now cleanup the bad trials (empty rows) of stro.trial, stro.ras, and
- %stro.sum.absTrialNum
- waitbar(1, waitdlg, ' ... cleaning up');
- stro.ras((goodTrial+1):end, :) = [];
- stro.trial((goodTrial+1):end, :) = [];
- stro.sum.absTrialNum(goodTrial+1:end) = [];
- %close up the wait dialog
- close(waitdlg);
- %****************************************%
- % END OF MAIN LOOP. SUBFUNCTIONS BELOW %
- %****************************************%
-
- %%
- function [stro, ecodes, spikes, anlg, waves, events] = setupStructures();
- %this function is necessary because later on we'll test the length of
- %spikes.names and anlg.names to iterate over neuronal and anlg data. If
- %we decide not to record from an anlg channel, setting the field
- %anlg.name to an empty vector (i.e., []) will allow this program to
- %function properly. ditto for spike channels. In addition, these
- %variables are global, so previous definitions of the variables could
- %intrude here.
-
- ecodes = [];
- spikes = struct('names', {{}}, ...
- 'timestamps', {{}});
- anlg = struct('names', {{}}, ...
- 'ADFrequency', {{}}, ...
- 'ADtoMV', {{}}, ...
- 'fragStarts', {{}}, ...
- 'nPtsPerFrag', {{}}, ...
- 'data', {{}});
- waves = struct('names', {{}}, ...
- 'WFrequency', {{}}, ...
- 'timestamps', {{}}, ...
- 'waveforms', {{}});
- events = struct();
-
- stro = struct('sum', struct(...
- 'fileName', [], ...
- 'date', date, ...
- 'paradigmID', [], ...
- 'waves', struct(...
- 'wf_id', [], ...
- 'storeRates', []), ...
- 'analog', struct( ...
- 'sigid', [], ...
- 'storeRates', [], ...
- 'ADtoMV', []),...
- 'trialFields', {{}}, ...
- 'rasterCells', {{}}, ...
- 'otherFields', {{}}, ...
- 'exptParams', struct(), ...
- 'absTrialNum', [], ...
- 'concat', []), ...
- 'trial', {{}}, ...
- 'ras', {{}}, ...
- 'other', {{}});
- end
- %%
- function unpackNexFile()
- %open differently for mac and pc:
- if isunix
- fid = fopen(stro.sum.fileName, 'r', 'l');
- elseif ispc
- fid = fopen(stro.sum.fileName, 'r');
- end
- if(fid == -1)
- error('Unable to open file')
- return
- end
- %open the file and make sure that it's a neuroexplorer file
- magic = fread(fid, 1, 'int32');
- if magic ~= 827868494
- error 'The file is not a valid .nex file'
- end
- %read in the header info
- version = fread(fid, 1, 'int32');
- comment = deblank(char(fread(fid, 256, 'char')'));
- freq = fread(fid, 1, 'double'); %acquisition rate for things other than continuous variables
- tbeg = fread(fid, 1, 'int32') ./ freq;
- tend = fread(fid, 1, 'int32') ./ freq;
- nvar = fread(fid, 1, 'int32');
- fseek(fid, 260, 'cof'); % skip location of next header and padding
- %now read the file putting the appropriate 'type' of Nex info in either
- %the ecodes, spikes, or analog vectors. For now ignore the other 'types'
- neuronCount = 0;
- contSigCount = 0;
- waveCount = 0;
- for i=1:nvar
- type = fread(fid, 1, 'int32');
- varVersion = fread(fid, 1, 'int32');
- name = deblank(char(fread(fid, [1 64], 'char')));
- offset = fread(fid, 1, 'int32');
- n = fread(fid, 1, 'int32');
- dummy = fread(fid, 32, 'char'); %sikps what's commented out below
- %WireNumber = fread(fid, 1, 'int32');
- %UnitNumber = fread(fid, 1, 'int32');
- %Gain = fread(fid, 1, 'int32');
- %Filter = fread(fid, 1, 'int32');
- %XPos = fread(fid, 1, 'double');
- %YPos = fread(fid, 1, 'double');
- WFrequency = fread(fid, 1, 'double'); % wf sampling fr.
- ADtoMV = fread(fid, 1, 'double'); % coeff to convert from AD values to Millivolts.
- NPointsWave = fread(fid, 1, 'int32'); % number of points in each wave
- NMarkers = fread(fid, 1, 'int32'); % how many values are associated with each marker
- MarkerLength = fread(fid, 1, 'int32'); % how many characters are in each marker value
- MVOfffset = fread(fid, 1, 'double'); % coeff to shift AD values in Millivolts: mv = raw*ADtoMV+MVOfffset
- filePosition = ftell(fid);
- switch type
- case 0 % neuron
- waitbar(0, waitdlg,' ...adding spike data');
- neuronCount = neuronCount+1;
- spikes.names(neuronCount) = {name};
- fseek(fid, offset, 'bof');
- spikes.timestamps{neuronCount} = fread(fid, [n 1], 'int32')./freq;
- fseek(fid, filePosition, 'bof');
- case 1 % event (start/stop)
- waitbar(0, waitdlg, sprintf(' ...adding %s times', name));
- fseek(fid, offset, 'bof');
- ts = fread(fid, [n 1], 'int32')./freq;
- events = setfield(events, lower(name), ts);
- fseek(fid, filePosition, 'bof');
- case 2 % interval
- %ignore this type for now
- fprintf('Ignoring the ''interval'' type \n');
- case 3 % waveform
- waveCount = waveCount+1;
- waves.names(waveCount) = {name};
- waves.WFrequency{waveCount} = WFrequency;
- fseek(fid, offset, 'bof');
- waves.timestamps{waveCount} = fread(fid, [n 1], 'int32')./freq;
- waves.waveforms{waveCount} = (fread(fid, [NPointsWave n], 'int16').*ADtoMV + MVOfffset)'; %transpose to make trials go down columns and time accross rows
- fseek(fid, filePosition, 'bof');
- case 4 % population vector
- %ignore this type for now
- fprintf('Ignoring the ''population vector'' type \n');
- case 5 % continuous variable (i.e. eye signals)
- waitbar(0, waitdlg, sprintf(' ...adding analog data: %s', name));
- contSigCount = contSigCount+1;
- anlg.names(contSigCount) = {name};
- anlg.ADtoMV(contSigCount) = {ADtoMV};
- anlg.ADFrequency{contSigCount} = WFrequency;
- fseek(fid, offset, 'bof');
- % get the start times of each analog signal snippet
- anlg.fragStarts{contSigCount} = fread(fid, [n 1], 'int32')./freq;
- %get the number of points sammpled during each of the snippets
- nFrag = fread(fid, [n 1], 'int32');
- nFrag(n+1) = NPointsWave;
- anlg.nPtsPerFrag{contSigCount} = nFrag;
- %now bring in the AtoD data for the entire recording
- anlg.data{contSigCount} = fread(fid, [NPointsWave 1], 'int16').*ADtoMV + MVOfffset;
- fseek(fid, filePosition, 'bof');
- case 6 % marker (i.e. ecodes)
- waitbar(0, waitdlg, ' ...adding ecodes');
- fseek(fid, offset, 'bof');
- timeStamps = fread(fid, [n 1], 'int32') ./ freq;
- %check to make sure that there are actually markers to
- %retreive before trying
- if (NMarkers ~= 1)
- error('bad number of markers in file');
- end
- %check to make sure that you're about to read in strobed
- %codes
- mname = char(fread(fid, [1 64], 'char'));
- if ~strncmp('DIO', mname, 3)
- error('unknown marker name')
- end
- %now collect the marker (ecode) data. convert them to
- %numeric representations and return them along with timestamps
- markers = deblank(fscanf(fid, '%c', [MarkerLength n])'); %transpose, then take off the trailing blank
- markers = markers-48; %convert to numeric
- powersOfTen = 10.^[(MarkerLength-2) : -1 : 0]';
- ecodes = [ecodes; [timeStamps, markers*powersOfTen]];
- %return to where you started in the file!
- fseek(fid, filePosition, 'bof');
- otherwise
- fprintf('unknown variable type <%s> \n', type);
- end
- dummy = fread(fid, 60, 'char'); %skip some junk to get to the next block of data
- end
- fclose(fid);
- end
- %%
- function trialCodes = checkTrial(trialStart, trialStop, badLogic)
- %pull out the relavent ecodes
- ind = (ecodes(:,1) >= trialStart) & (ecodes(:,1) <= trialStop);
- trialCodes = ecodes(ind,:);
-
- %test the ecodes to determine if the trial is good or bad. 'badLogic'
- %should be defined in the paradigm header file.
- bad = eval(badLogic);
- %return an empty vector if it's a bad trial
- if bad
- trialCodes = [];
- end
- end
- %%
- function rasterByTrial(start, stop, gtCounter, totTrialCounter, additionalRasterFields, trialCodes)
- %fill up the rows of stro.ras trial by trial. treat sike data, anlg
- %sigs, and other independently. Format is:
- % <neuron 1> <neuron 2> ... <waveForm1>... <anlg 1> <anlg 2> ... <atime> <other>
-
- %each neuron's spikes will go into it's own column:
- if ~isempty(spikes.names)
- for ind = 1:length(spikes.names);
- stro.ras{gtCounter, ind} = spikes.timestamps{ind}((spikes.timestamps{ind} >= start) & (spikes.timestamps{ind} <= stop));
- end
- else
- ind = 0;
- end
- %iterate over the waveform data (if available):
- for a = 1:length(waves.names);
- ind = ind + 1;
- whichSpikes = ((waves.timestamps{a} >= start) & (waves.timestamps{a} <= stop));
- stro.ras{gtCounter, ind} = waves.waveforms{a}(whichSpikes, :);
- end
- %iterate over anlg channels
- for a = 1:length(anlg.names);
- ind = ind + 1;
- startInd = anlg.nPtsPerFrag{a}(totTrialCounter) + 1;
- stopInd = anlg.nPtsPerFrag{a}(totTrialCounter+1);
- stro.ras{gtCounter, ind} = anlg.data{a}(startInd:stopInd);
- end
-
- %explicitly add the anlg start time. All the analog channels should
- %have the same number of start times, so indexing this way is o.k.
- if ~isempty(anlg.fragStarts)
- ind = ind + 1;
- stro.ras{gtCounter, ind} = anlg.fragStarts{1}(totTrialCounter);
- end
- %now look for additional elements as specified by
- %additionalRasterFields (which is defined in the 'paradigmFile'
- for a = 1:size(additionalRasterFields, 2);
- ind = ind + 1;
- nums = dat2num(trialCodes(:,2), [additionalRasterFields{4, a}], C.VALOFFSET, additionalRasterFields{2, a}, additionalRasterFields{3, a});
- stro.ras{gtCounter, ind} = nums{1};
- end
- end
- %%
- function indexByTrial(trialCodes, trialParamsHeader, goodTrial, offset)
- %deal with the ints, longs, and doubles
- types = {'int', 'long', 'float', 'double'};
- for a = 1:length(types);
- typeInd = strmatch(types{a}, strvcat(trialParamsHeader{2,:}));
- if typeInd
- nums = dat2num(trialCodes(:,2), [trialParamsHeader{4, typeInd}], offset, types{a}, 0);
- stro.trial(goodTrial, typeInd) = deal([nums{:}]);
- end
- end
- %now deal with the 'time' variables
- timeInd = strmatch('time', strvcat(trialParamsHeader{2,:}));
- timeInd_code = [trialParamsHeader{4, timeInd}];
- for a = 1:length(timeInd)
- time = trialCodes((trialCodes(:,2) == timeInd_code(a)), 1); %there will be an error if numel(time) > 1
- if ~isempty(time)
- stro.trial(goodTrial, timeInd(a)) = time;
- end
- end
- end
- %%
- function otherByTrial(trialCodes, otherInstructions, gtCounter, codeStruct)
- ind = 0;
- for a = 1:length(otherInstructions);
- ind = ind+1;
- stro.other{gtCounter, a} = feval(otherInstructions{2,a}, trialCodes, codeStruct);
- end
-
- end
- %%
- function rasterCells = makeRasterHeader(additionalRasterFields)
- %start with the spike data then add the anlg data. this means that
- %stro.raster will be composed of columns:
- % <neuron 1> <neuron 2> ... <neuron 3> <anlg 1> <anlg 2> ... <anlg 3> <atime> <other>
-
- rasterCells = {};
- if ~isempty(spikes.names)
- for ind = 1:length(spikes.names)
- rasterCells{ind} = spikes.names{ind};
- end
- else
- ind = 0; %incase there isn't any spike data
- end
- %iterate over waveforms (if they exist)
- for a = 1:length(waves.names);
- ind = ind + 1;
- rasterCells{ind} = waves.names{a};
- end
- %now add the analog channels
- for a = 1:length(anlg.names)
- ind = ind + 1;
- rasterCells{ind} = anlg.names{a};
- end
- %explictly add the 'anlgStart' field
- if ~isempty(anlg.fragStarts)
- ind = ind + 1;
- rasterCells{ind} = 'anlgStartTime';
- end
- %now add the additional raster headers (from the paradigm header file)
- for a = 1:size(additionalRasterFields, 2)
- ind = ind + 1;
- rasterCells{ind} = additionalRasterFields{1,a};
- end
- end
- %%
- function getExperimentalParams()
- %deal with the entries into the stro.sum.exptParams one by one.
- for a = 1:size(exptParamsHeader, 2);
- nums = dat2num(ecodes, [exptParamsHeader{4, a}], C.VALOFFSET, exptParamsHeader{2, a}, exptParamsHeader{3, a});
- stro.sum.exptParams = setfield(stro.sum.exptParams, exptParamsHeader{1, a}, nums{1});
- end
- end
- %%
- function cleanupAnalogSignals()
- %check to make sure that the number of analog fragments is the same for
- %each channel. ditto for the number of AtoD points per fragment
-
- %Ensure an equal number of starts and stops (i.e. check to
- %make sure that data acquisition was not arrested part way through a
- %trial)
- excessTrial = (length(events.start) - length(events.stop));
- if (excessTrial == 1)
- events.start(end) = [];
- fprintf(' Unequal number of STARTs and STOPs. Deleted the terminal START.\n');
- elseif (excessTrial == -1 && (events.stop(1) < events.start(1)))
- events.stop(1) = [];
- fprintf(' Unequal number of STARTs and STOPs. Deleted the initial STOP.\n');
- elseif (excessTrial ~= 0)
- error(' Unequal number of STARTS and STOPS.');
- end
-
- %now determine if there are the same number of trial starts/stops
- %as there are anlg starts/stops. delete any anlg data that occurs
- %before the first trial.
- if isempty(anlg.fragStarts)
- return
- end
- err = abs(events.start(1) - anlg.fragStarts{1});
- startIdx = find(err == min(err));
- for a = 1:length(anlg.names)
- anlg.fragStarts{a}(1:(startIdx-1)) = [];
- anlg.nPtsPerFrag{a}(1:(startIdx-1)) = [];
- end
-
- if numel(anlg.names) > 1
- if ~(isequal(anlg.fragStarts{1:end}))
- error('frag starts are unequal')
- end
- if ~(isequal(anlg.nPtsPerFrag{1:end}))
- error('points per frag are unequal')
- end
- end
- end
- end