ft_realtime_neuralynxproxy.m

function ft_realtime_neuralynxproxy(cfg)

% FT_REALTIME_NEURALYNXPROXY reads continuous data from a Neuralynx Cheetah
% acquisition system and writes it to a FieldTrip buffer.
%
% The FieldTrip buffer is a network transparent server that allows the acquisition
% client to stream data to it. An analysis client can connect to read the data upon
% request. Multiple clients can connect simultaneously, each analyzing a specific
% aspect of the data concurrently.
%
% Use as
%   ft_realtime_neuralynxproxy(cfg)
%
% The configuration should contain
%   cfg.acquisition          = string, name of computer running the Cheetah software (default = 'fcdc284')
%   cfg.channel              = cell-array, see FT_CHANNELSELECTION (default = 'all')
%
% The target to write the data to is configured as
%   cfg.target.datafile      = string, target destination for the data (default = 'buffer://localhost:1972')
%   cfg.target.dataformat    = string, default is determined automatic
%
% To stop this realtime function, you have to press Ctrl-C
%
% See also FT_REALTIME_SIGNALPROXY, FT_REALTIME_SIGNALVIEWER

% Copyright (C) 2008, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% 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$

% set the defaults
if ~isfield(cfg, 'acquisition'),        cfg.acquisition = 'fcdc284';                      end
if ~isfield(cfg, 'target'),             cfg.target = [];                                  end
if ~isfield(cfg, 'channel'),            cfg.channel = 'all';                              end
if ~isfield(cfg.target, 'datafile'),    cfg.target.datafile = 'buffer://localhost:1972';  end
if ~isfield(cfg.target, 'dataformat'),  cfg.target.dataformat = [];                       end % default is to use autodetection of the output format
if ~isfield(cfg.target, 'eventfile'),   cfg.target.eventfile = 'buffer://localhost:1972'; end
if ~isfield(cfg.target, 'eventformat'), cfg.target.eventformat = [];                      end % default is to use autodetection of the output format

% this requires an external Neuralynx toolbox for the communication with Cheetah
ft_hastoolbox('nlxnetcom', 1);

%% the record pool is used to collect the fragments while they become available
recordPool        = [];
recordStamp       = [];
recordComplete    = [];

%% Initialize
if ~NlxAreWeConnected
  disp(sprintf('Connecting to %s...', cfg.acquisition));
  succeeded = NlxConnectToServer(cfg.acquisition);
  if succeeded ~= 1
    fprintf('FAILED connect to %s. Exiting script.', cfg.acquisition);
    return;
  else
    fprintf('Connected to %s.', cfg.acquisition);
  end
end

%% Identify this program to the server we're connected to.
succeeded = NlxSetApplicationName('My Matlab Script');
if succeeded ~= 1
  disp 'FAILED set the application name'
else
  disp 'PASSED set the application name'
end

%% get a list of all objects in Cheetah, along with their types.
[succeeded, cheetahObjects, cheetahTypes] = NlxGetCheetahObjectsAndTypes;
if succeeded == 0
  disp 'FAILED get cheetah objects and types'
else
  disp 'PASSED get cheetah objects and types'
end

%% open up a stream for all objects
for index = 1:length(cheetahObjects)
  succeeded = NlxOpenStream(cheetahObjects(index));
  if succeeded == 0
    fprintf('FAILED to open stream for %s', char(cheetahObjects(index)));
    break;
  end
end
if succeeded == 1
  disp 'PASSED open stream for all current objects'
end


cscIndex = find(strcmp('CscAcqEnt', cheetahTypes))';
cscCount = length(cscIndex);
evtIndex = find(strcmp('EventAcqEnt', cheetahTypes))';
evtCount = length(evtIndex);

% select a subset of channels
cfg.channel = ft_channelselection(cfg.channel, cheetahObjects(cscIndex));
selIndex    = match_str(cheetahObjects, cfg.channel)';

hdr = []; % FIXME
targethdr = [];

numSamples = 0;
pass = 0;
again = true;

while again
  pass = pass + 1;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % get the events
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  for objectIndex = evtIndex(:)'
    objectToRetrieve = cheetahObjects{objectIndex};
    [succeeded,  timeStampArray, eventIDArray, ttlValueArray, eventStringArray, numRecordsReturned, numRecordsDropped ] = NlxGetNewEventData(objectToRetrieve);

    if ~succeeded
      fprintf('FAILED to get new new events from stream %s on pass %d', objectToRetrieve, pass);
      break;

    elseif numRecordsReturned==0
      disp(sprintf('FAILED to get new events from stream %s on pass %d', objectToRetrieve, pass));
      break;

    else
      disp(sprintf('Retrieved %d new events from %s with %d dropped.', numRecordsReturned, objectToRetrieve, numRecordsDropped));

      if ~isempty(targethdr) && isfield(targethdr, 'FirstSampleTimeStamp') && isfield(targethdr, 'TimeStampPerSample')
        % convert to event array
        event = [];
        for i=1:numRecordsReturned
          event(i).type     = 'ttl';
          event(i).value    = ttlValueArray(i);
          event(i).sample   = double(uint64(timeStampArray(i))-targethdr.FirstSampleTimeStamp)/targethdr.TimeStampPerSample;
          event(i).offset   = 0;
          event(i).duration = 0;
        end
        % write the events
        ft_write_event(cfg.target.eventfile, event, 'eventformat', cfg.target.eventformat);
      end
    end % if succeeded
  end % for all event streams

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % get the continuous data
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  objectData = [];

  for objectIndex = cscIndex(selIndex)
    objectToRetrieve = cheetahObjects{objectIndex};
    % determine the type of acquisition entity we are currently indexed to and call the appropriate function for that type
    [succeeded, dataArray, timeStampArray, channelNumberArray, samplingFreqArray, numValidSamplesArray, numRecordsReturned, numRecordsDropped ] = NlxGetNewCSCData(objectToRetrieve);

    if ~succeeded
      fprintf('FAILED to get new data for CSC stream %s on pass %d', objectToRetrieve, pass);
      break;

    elseif numRecordsReturned==0
      disp(sprintf('FAILED to get new data for CSC stream %s on pass %d', objectToRetrieve, pass));
      break;

    else
      disp(sprintf('Retrieved %d CSC records for %s with %d dropped.', numRecordsReturned, objectToRetrieve, numRecordsDropped));

      objectData(objectIndex).dataArray            = reshape(dataArray, 512, numRecordsReturned);
      objectData(objectIndex).timeStampArray       = timeStampArray;
      objectData(objectIndex).channelNumberArray   = channelNumberArray;
      objectData(objectIndex).samplingFreqArray    = samplingFreqArray;
      objectData(objectIndex).numValidSamplesArray = numValidSamplesArray;
      objectData(objectIndex).numRecordsReturned   = numRecordsReturned;
      objectData(objectIndex).numRecordsDropped    = numRecordsDropped;

      selectValid = (numValidSamplesArray==512);
      if ~all(selectValid)
        ft_warning('some samples were not valid')
      end

      % remove the records that are somehow corrupted
      objectData(objectIndex).dataArray            = objectData(objectIndex).dataArray(:,selectValid);
      objectData(objectIndex).timeStampArray       = objectData(objectIndex).timeStampArray(selectValid);
      objectData(objectIndex).channelNumberArray   = objectData(objectIndex).channelNumberArray(selectValid);
      objectData(objectIndex).samplingFreqArray    = objectData(objectIndex).samplingFreqArray(selectValid);
      objectData(objectIndex).numValidSamplesArray = objectData(objectIndex).numValidSamplesArray(selectValid);
      objectData(objectIndex).numRecordsValid      = sum(selectValid);
      objectData(objectIndex).numRecordsReturned   = numRecordsReturned;
      objectData(objectIndex).numRecordsDropped    = numRecordsDropped;

    end % if succeeded

  end % for objectIndex

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % create a FieldTrip compatible header structure
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if isempty(hdr) && ~isempty(objectData) && ~isempty(objectData(1))
    hdr.Fs                 = objectData(1).samplingFreqArray(1);  % sampling frequency
    hdr.nChans             = length(selIndex);                    % number of channels
    hdr.nSamples           = 0;                                   % number of samples per trial
    hdr.nSamplesPre        = 0;                                   % number of pre-trigger samples in each trial
    hdr.nTrials            = 1;                                   % number of trials, 1 means continuous
    hdr.label              = cheetahObjects(cscIndex(selIndex));  % cell-array with labels of each channel
    hdr.FirstTimeStamp     = nan;
    hdr.TimeStampPerSample = nan;
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % reorganize the objects into records that each contain a Nchan X Nsamples
  % data matrix
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  for i=1:numel(objectData)
    nrec = numel(objectData(i).channelNumberArray);
    for j=1:nrec

      if ~isempty(recordPool)
        % The temporary variable is to work around a MATLAB bug in version 7.5 on windows
        % see http://www.mathworks.com/support/bugreports/details.html?rp=425978
        tmp = objectData(i).timeStampArray(j);
        sel = find(recordStamp==tmp);
      else
        sel = [];
      end

      if isempty(sel)
        % create a new record in the pool
        recordComplete(end+1)       = false;
        recordStamp(end+1)          = objectData(i).timeStampArray(j);
        recordPool(end+1).timeStamp = objectData(i).timeStampArray(j);
        recordPool(end  ).chanCount = 0;
        recordPool(end  ).skipCount = 0;
        recordPool(end  ).present   = false(cscCount, 1);
        recordPool(end  ).dat       = zeros(cscCount, 512);
        sel = numel(recordPool);
      end

      % add this data to the selected record
      chan = objectData(i).channelNumberArray(j) + 1; % should be offset 1 instead of 0
      recordPool(sel).present(chan) = true;
      recordPool(sel).dat(chan,:)   = objectData(i).dataArray(:,j);
      recordPool(sel).chanCount     = recordPool(sel).chanCount + 1;
      recordComplete(sel)           = all(recordPool(sel).present(selIndex));
    end
  end % for all objects that were acuired

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % write the records that are complete
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if any(recordComplete)
    last = find(recordComplete, 1, 'last');

    % give some feedback and fix the recordComplete vector
    fprintf('writing %d records\n', last);

%     for i=1:last
%       if recordComplete(i)
%         fprintf('writing complete record with timestamp %d\n', recordPool(i).timeStamp);
%       else
%         fprintf('writing incomplete record with timestamp %d\n', recordPool(i).timeStamp);
%       end
%     end

    % all records up to the last complete one will be written
    recordComplete(1:last) = true;

    if ~isempty(last)
      % concatenate all data in a single Nchans X Nsamples array
      % only use the selected channels, the other channels will all be zero
      alldata = zeros(length(selIndex), last*512);
      for i=1:last
        begsample = (i-1)*512+1;
        endsample = (i  )*512;
        alldata(:,begsample:endsample) = recordPool(i).dat(selIndex,:);
      end

      if numSamples==0
        % write the header and write data of this record to the buffer
        ft_write_data(cfg.target.datafile, alldata, 'dataformat', cfg.target.dataformat, 'header', hdr, 'append', false);
        targethdr.FirstSampleTimeStamp = uint64(recordPool(1).timeStamp); % this is only defined once
      else
        % write the data of this record to the buffer
        ft_write_data(cfg.target.datafile, alldata, 'dataformat', cfg.target.dataformat, 'append', true);
      end
      targethdr.TimeStampPerSample   = double(uint64(recordPool(i).timeStamp) - targethdr.FirstSampleTimeStamp)/numSamples;
      numSamples                     = numSamples + size(recordPool(i).dat,2);

      clear alldata
    end
  end

  % remove the records that have been written, i.e. keep all others
  recordPool      = recordPool(~recordComplete);
  recordStamp     = recordStamp(~recordComplete);
  recordComplete  = recordComplete(~recordComplete);

  fprintf('pass %d, there are %d records remaining in the pool\n', pass, numel(recordPool));

end % while again