/external/fieldtrip/trunk/fileio/private/read_neuralynx_dma.m
MATLAB | 280 lines | 175 code | 16 blank | 89 comment | 40 complexity | 74089d797fb0187d136ab5e8bbe77247 MD5 | raw file
Possible License(s): GPL-2.0, GPL-3.0
- function [dat] = read_neuralynx_dma(filename, begsample, endsample, channel);
- % READ_NEURALYNX_DMA reads specified samples and channels data from a Neuralynx DMA log file
- %
- % Use as
- % [hdr] = read_neuralynx_dma(filename)
- % [dat] = read_neuralynx_dma(filename, begsample, endsample)
- % [dat] = read_neuralynx_dma(filename, begsample, endsample, chanindx)
- %
- % The channel specification can be a vector with indices, or a single string with the value
- % 'all', 'stx', 'pid', 'siz', 'tsh', 'tsl',
- % 'cpu', 'ttl', 'x01', ..., 'x10'
- %
- % This function returns the electrophysiological data in AD units
- % and not in uV. You should look up the details of the headstage and
- % the Neuralynx amplifier and scale the values accordingly.
- % Copyright (C) 2005-2008, Robert Oostenveld
- %
- % This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
- % for the documentation and details.
- %
- % FieldTrip 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.
- %
- % FieldTrip 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.
- %
- % You should have received a copy of the GNU General Public License
- % along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
- %
- % $Id: read_neuralynx_dma.m 2885 2011-02-16 09:41:58Z roboos $
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % The data is simply a stream of constant size records, the size of the
- % record is dependant on the number of Input Boards in the system. 8 in your
- % case, there are 32 - 32 bit integers from each Input Board.
- %
- % // Digital Lynx Packets look like this for a 1 board system:
- % // stx (or SOP) 0x00000800
- % // Packet ID 0x00000001
- % // siz 0x0000002A (# A/D data ints + #extra wds = 32+10)
- % // TimeStamp Hi 1 32 bit word
- % // TimeStamp Low 1 32 bit word
- % // cpu 1 32 bit word
- % // ttl 1 32 bit word
- % // 10 extras 10 32 bit words
- % // A/D data 32 32 bit words per board
- % // CRC 1 32 bit XOR of the entire packet including stx
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- needhdr = (nargin==1);
- needdat = (nargin>=2);
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % since Cheetah 4.80 the DMS log files have a 16384 byte header
- % determine whether that is the case in this file
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- fid = fopen(filename, 'rb', 'ieee-le');
- buf = fread(fid, 10, 'uchar=>uchar');
- if all(buf(1:4)==35)
- % the file starts with "####"
- hdroffset = 16384;
- elseif all(buf>31 & buf<127)
- % the file does not start with "####" but the header seems to be ascii anyway
- hdroffset = 16384;
- else
- hdroffset = 0;
- end
- if hdroffset
- % read the ascii header from the file, it does not contain much usefull information
- orig = neuralynx_getheader(filename);
- end
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % since Cheetah 5.0 the DMA log file does not always start with a complete
- % packet but can start with an incomplete "junk" packet after the header.
- % The junk due to that incomplete packet should then also be skipped.
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- if hdroffset
- % skip the standard ascii header
- fseek(fid, hdroffset, 'bof');
- % read a single block
- buf = fread(fid, 274, 'uint32=>uint32');
- % search for the first known values of a proper packet
- junk = find((buf(1:(end-1))==2048) & (buf(2:(end-0))==1), 1) - 1;
- if ~isempty(junk)
- % add the additional junk samples to the header offset
- hdroffset = hdroffset + junk*4;
- end
- end
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % since modifying our hardware from 256 to 32 channel (19 March 2007), the
- % number of channels in a DMA log file can vary. Start with the assumption
- % that the file is from the 256 channel system with 8 boards that we have
- % in Nijmegen and verify the CRC
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- nboards = 8;
- nchanboard = 32;
- blocksize = nboards*nchanboard+18; % in 4 byte words
- fseek(fid, hdroffset, 'bof');
- % read a single block
- buf = fread(fid, blocksize, 'uint32=>uint32');
- % determine the number of boards by systematically checking the CRC value
- while (nboards>0)
- if neuralynx_crc(buf(1:blocksize-1))==buf(blocksize)
- % the CRC is correct for this number of boards
- break
- else
- nboards = nboards - 1; % decrease by one
- blocksize = nboards*nchanboard+18; % in 4 byte words
- end
- end
- fclose(fid);
- if ~nboards
- error('could not determine the number of boards and channels');
- end
- % deal with ascii and numeric input for the channel selection
- if nargin>3 && ~ischar(channel)
- chanindx = channel;
- channel = 'selection';
- elseif nargin>3 && ischar(channel)
- switch channel
- case 'stx'
- chanindx = 1;
- case 'pid'
- chanindx = 2;
- case 'siz'
- chanindx = 3;
- case 'tsh'
- chanindx = 4;
- case 'tsl'
- chanindx = 5;
- case 'cpu'
- chanindx = 6;
- case 'ttl'
- chanindx = 7;
- case 'x01'
- chanindx = 8;
- case 'x02'
- chanindx = 9;
- case 'x03'
- chanindx = 10;
- case 'x04'
- chanindx = 11;
- case 'x05'
- chanindx = 12;
- case 'x06'
- chanindx = 13;
- case 'x07'
- chanindx = 14;
- case 'x08'
- chanindx = 15;
- case 'x09'
- chanindx = 16;
- case 'x10'
- chanindx = 17;
- case 'crc'
- chanindx = blocksize;
- case 'all'
- chanindx = 1:blocksize;
- otherwise
- error('unknown value in channel');
- end
- elseif nargin<4
- channel = 'all';
- chanindx = 1:blocksize;
- end
- if needhdr
- % these have to be hardcoded, since the DMA logfile does not contain header information
- hdr = [];
- if hdroffset
- % remember the ascii header from the file, it does not contain much usefull information
- hdr.orig = orig;
- end
- hdr.Fs = 32556; % sampling frequency
- hdr.nChans = nboards*nchanboard; % number of analog channels
- hdr.nSamples = inf; % number of samples per trial
- hdr.nSamplesPre = 0; % number of pre-trigger samples in each trial
- hdr.nTrials = 1; % number of trials
- hdr.label = {}; % cell-array with labels of each channel
- for i=1:hdr.nChans
- chanlabel{i} = sprintf('csc%03d', i);
- end
- statuslabel = {
- 'stx'
- 'pid'
- 'siz'
- 'tsh'
- 'tsl'
- 'cpu'
- 'ttl'
- 'x01'
- 'x02'
- 'x03'
- 'x04'
- 'x05'
- 'x06'
- 'x07'
- 'x08'
- 'x09'
- 'x10'
- };
- hdr.label = cat(1, statuslabel(:), chanlabel(:), {'crc'}); % concatenate all channel labels
- hdr.nChans = length(hdr.label); % this includes all channels
- % determine the length of the file, expressed in samples
- fid = fopen(filename, 'rb', 'ieee-le');
- fseek(fid, 0, 'eof');
- hdr.nSamples = floor((ftell(fid)-hdroffset)/(blocksize*4));
- % determine the timestamp of the first and of the last sample
- datoffset = 0*blocksize;
- fseek(fid, hdroffset + datoffset, 'bof');
- buf = fread(fid, blocksize, 'uint32=>uint32');
- beg_stx = buf(1);
- beg_tsh = buf(4);
- beg_tsl = buf(5);
- datoffset = (hdr.nSamples-1)*blocksize*4;
- fseek(fid, hdroffset + datoffset, 'bof');
- buf = fread(fid, blocksize, 'uint32=>uint32');
- end_stx = buf(1);
- end_tsh = buf(4);
- end_tsl = buf(5);
- fclose(fid);
- % check that the junk at the beginning was correctly detected
- if (beg_stx~=2048)
- error('problem with STX at the begin of the file');
- end
- % check that the file is truely continuous, i.e. no gaps with junk in between
- if (end_stx~=2048)
- error('problem with STX at the end of the file');
- end
- % combine the two uint32 words into a uint64
- hdr.FirstTimeStamp = timestamp_neuralynx(beg_tsl, beg_tsh);
- hdr.LastTimeStamp = timestamp_neuralynx(end_tsl, end_tsh);
- hdr.TimeStampPerSample = double(hdr.LastTimeStamp-hdr.FirstTimeStamp)./(hdr.nSamples-1); % this should be double, since it can be fractional
- % only return the header information
- dat = hdr;
- dat.orig.Offset = hdroffset;
- elseif length(chanindx)>1
- % read the desired samples for all channels
- fid = fopen(filename, 'rb', 'ieee-le');
- datoffset = (begsample-1)*blocksize*4;
- fseek(fid, hdroffset + datoffset, 'bof'); % skip the header and samples that do not have to be read
- % read the data
- dat = fread(fid, [blocksize (endsample-begsample+1)], 'int32=>int32');
- fclose(fid);
- if size(dat,2)<(endsample-begsample+1)
- error('could not read all samples');
- end
- if ~strcmp(channel, 'all')
- % select the subset of desired channels
- dat = dat(chanindx,:);
- end
- elseif length(chanindx)==1
- % read the desired samples for only one channel
- fid = fopen(filename, 'rb', 'ieee-le');
- datoffset = (begsample-1)*blocksize*4;
- fseek(fid, hdroffset + datoffset, 'bof'); % skip the header and samples that do not have to be read
- fseek(fid, (chanindx-1)*4, 'cof'); % skip to the channel of interest
- % Note that the last block with 274*4 bytes can sometimes be incomplete, which is relevant when endsample=inf
- dat = fread(fid, [1 (endsample-begsample+1)], 'int32=>int32', (nboards*32+18-1)*4);
- if size(dat,2)<(endsample-begsample+1)
- error('could not read all samples');
- end
- fclose(fid);
- end