diff --git a/fileio/private/read_neuralynx_ds.m b/fileio/private/read_neuralynx_ds.m
index a5c5a55d56..8bf026c0c9 100644
--- a/fileio/private/read_neuralynx_ds.m
+++ b/fileio/private/read_neuralynx_ds.m
@@ -157,17 +157,42 @@
     end
   end % if ftype_ncs
 
+  % take into account the pause/stop feature + epoch data if discontinuous recording
+  ncs = read_neuralynx_ncs(fname{1});
+  idxPauses = find(ncs.NumValidSamp < recordsize);
+  nValidSamples = ncs.NumValidSamp(idxPauses);
+  nPauses = length(idxPauses);
+  isPause = false(nPauses, 1);
+  for i = 1:nPauses
+      if idxPauses(i) == ncs.NRecords
+          isPause(i) = true;
+      else
+          timeGap = double(diff(ncs.TimeStamp(idxPauses(i):idxPauses(i)+1)));
+          if timeGap > 3 * mode(diff(ncs.TimeStamp)) % threshold to be tuned
+              isPause(i) = true;
+          end
+      end
+  end
+  idxPauses = idxPauses(isPause);
+  nValidSamples = nValidSamples(isPause);
+  nTrials = length(idxPauses);
+  trl = zeros(nTrials, 3);
+  trl(:, 2) = ((idxPauses-1)*recordsize + nValidSamples)';
+  trl(1, 1) = 1;
+  trl(2:end, 1) = (idxPauses(1:end-1)*recordsize+1)';
+
   % construct the header that applies to all channels combined
   hdr.nChans         = length(label);
   hdr.label          = label;
   hdr.filename       = fname;
-  hdr.nTrials        = 1;                           % it is continuous
+  hdr.nTrials        = nTrials;
   hdr.Fs             = SamplingFrequency(1);
-  hdr.nSamplesPre    = 0;                           % it is continuous
-
+  hdr.nSamplesPre    = 0;
+  hdr.trl = trl;
+  
   if ~isempty(ftype_ncs)
     % these elements are only relevant for continuously sampled channels
-    hdr.nSamples           = NRecords(1) * 512;
+    hdr.nSamples           = NRecords(1) * recordsize;
     hdr.FirstTimeStamp     = FirstTimeStamp(1);
     hdr.LastTimeStamp      = LastTimeStamp(1);
     hdr.TimeStampPerSample = TimeStampPerSample(1);
diff --git a/ft_preprocessing.m b/ft_preprocessing.m
index 6f930b091a..f90321c192 100644
--- a/ft_preprocessing.m
+++ b/ft_preprocessing.m
@@ -415,11 +415,15 @@
       trl(1,2) = hdr.nSamples*hdr.nTrials;
       trl(1,3) = -hdr.nSamplesPre;
     else
-      trl = zeros(hdr.nTrials, 3);
-      for i=1:hdr.nTrials
-        trl(i,1) = (i-1)*hdr.nSamples + 1;
-        trl(i,2) = (i  )*hdr.nSamples    ;
-        trl(i,3) = -hdr.nSamplesPre;
+      if isfield(hdr, 'trl')
+        trl = hdr.trl;
+      else
+        trl = zeros(hdr.nTrials, 3);
+	    for i=1:hdr.nTrials
+          trl(i,1) = (i-1)*hdr.nSamples+1;
+          trl(i,2) = (i)*hdr.nSamples;
+          trl(i,3) = -hdr.nSamplesPre;
+	    end
       end
     end
     cfg.trl = trl;
diff --git a/trialfun/ft_trialfun_neuralynx.m b/trialfun/ft_trialfun_neuralynx.m
new file mode 100644
index 0000000000..5b7ae8c315
--- /dev/null
+++ b/trialfun/ft_trialfun_neuralynx.m
@@ -0,0 +1,57 @@
+function trl = ft_trialfun_neuralynx(cfg)
+
+% FT_TRIALFUN_NEURALYNX outputs trials for Neuralynx datasets.
+%
+% Use as
+%   cfg          = [];
+%   cfg.dataset  = 'directory_with_ncs_files';
+%   cfg.trialfun = 'ft_trialfun_neuralynx';
+%   cfg = ft_definetrial(cfg);
+%   data = ft_preprocessing(cfg);
+%
+% A Neuralynx dataset consists of a directory containing as many .ncs files
+% as electrodes. The option cfg.dataset should refer to this directory.
+%
+% Neuralynx data is saved as a series of 512-sample blocks (or "records").
+% When user presses pause or stop buttons, the recording stops at sample i
+% from the current record, and the remaining samples (i+1 to 512) are a
+% copy of samples i+1 to 512 of the previous record. These copied samples
+% are marked as invalid in ncs.NumValidSamples, and have to be discarded as
+% they are not genuine samples and as they introduce redundancies.
+% 
+% Contact ludovic.bellier@berkeley.edu for any questions.
+
+
+gapThreshold = 3; % in number of records (e.g., at fs=8kHz, 3*512 samples represents 192 ms)
+recordsize = 512;
+
+dirname = cfg.dataset;
+if ~any(strcmp(dirname(end), {'/', '\'}))
+    dirname = [dirname filesep];
+end
+
+ls = dir([dirname '*ncs']);
+fname = cellfun(@(x) fullfile(dirname, x), {ls(:).name}, 'un', 0);
+
+ncs = read_neuralynx_ncs(fname{1});
+idxPauses = find(ncs.NumValidSamp < recordsize);
+nValidSamples = ncs.NumValidSamp(idxPauses);
+nPauses = length(idxPauses);
+isPause = false(nPauses, 1);
+for i = 1:nPauses
+    if idxPauses(i) == ncs.NRecords
+        isPause(i) = true;
+    else
+        timeGap = double(diff(ncs.TimeStamp(idxPauses(i):idxPauses(i)+1)));
+        if timeGap > gapThreshold * mode(diff(ncs.TimeStamp))
+            isPause(i) = true;
+        end
+    end
+end
+idxPauses = idxPauses(isPause);
+nValidSamples = nValidSamples(isPause);
+nTrials = length(idxPauses);
+trl = zeros(nTrials, 3);
+trl(:, 2) = ((idxPauses-1)*recordsize + nValidSamples)';
+trl(1, 1) = 1;
+trl(2:end, 1) = (idxPauses(1:end-1)*recordsize+1)';
\ No newline at end of file