function sds = sync_datasets(dss, node_names, PLOTRANGE, ONLYDATA)
  % load datasets
  len = length(dss);
  sds = {};
  first_ts = [];
  last_ts = [];
  for i=1:len
    ds = read_syncdata(dss{i});
    first_ts(i) = ds(1:1);
    last_ts(i) = ds(end,1);   
    sds{i} = ds;
  end
  
  % get first and last common timestamp
  ts0 = max(first_ts);
  tsend = min(last_ts);
  
  for i=1:len
    ds = sds{i}; 
    ds = ds(ds(:,1) > ts0 + (ds(2,1) - ds(1,1)) / 2 & ds(:,1) < tsend, :); % drop non-common sections
    ds(:,1) = ds(:,1) - ts0; % substract ts0 from each dataset's timestamp array
    sds{i} = ds;
  end
  
  % reference dataset and timescale
  ds_ref = sds{1};
  ts_ref = ds_ref(:,1);
  
  % create plot
  if (~ONLYDATA)
    figure('Position', [100 100 1000 1000])
  else
    figure('Position', [100 100 1000 400])
  end
  clf
  
  % plot samples
  if (~ONLYDATA)
    subplot(3,1,1);
  end
  
  marks = ["-x", "-o"];
  
  %PLOTRANGE = 1017:1020;
  
  VOLT_PER_BIN = 42.80E-06;
  
  % get beginning indices
  for i=1:len
    ds = sds{i};
    ts = ds(:,1);
    
    s = size(ds);
    
    for k = 2:s(2)
    %for k = 2:2
      if (ONLYDATA)
        mark = "-";
      else
        mark = marks(k-1);
      end
    
      plot(ts(PLOTRANGE), ds(PLOTRANGE,k) * VOLT_PER_BIN, mark);
      hold on
    end
  end
  
  grid on
  
  if (ONLYDATA)
    %xlabel("Idő [s]");
    xlabel("Time [s]")
  end
  
  %ylabel("Feszültség [V]");
  ylabel("Voltage [V]");
  xlim([ts_ref(PLOTRANGE(1)) ts_ref(PLOTRANGE(end))]);
    
  legend_lines = {};
  
  for i=1:length(node_names)
    legend_lines{i,1} = node_names{i,1};
  end
  
  legend(legend_lines);
  
  if (ONLYDATA)
    return;
  end
  
  % plot signal difference
  subplot(3,1,2);
  
  % normalize signals
  ds = sds{1};
  ds_norm_range{1} = ds(PLOTRANGE,2:end) ./ max(abs(ds(PLOTRANGE,2:end)));
  legend_lines = {};
  
  for i=2:length(node_names)
    ds = sds{i};
    ds_norm_range{i} = ds(PLOTRANGE,2:end) ./ max(abs(ds(PLOTRANGE,2:end)));
    
    for ch=1:2 
      legend_lines{2 * (i - 2) + ch, 1} = strcat(node_names{i,1}, "-", node_names{1,1}, " CH", num2str(ch));
    end
    
    ts = ds(:,1);
    diff_data = ((ds_norm_range{i} - ds_norm_range{1})); % ./ ds_norm_range{1});
    
%     for l=1:length(diff_data)
%       diff_data(l,:) = diff_data(l,:) ./ ds_norm_range{1}(l,:);
%     end
    
    mark = marks(k-1);
    plot(ts(PLOTRANGE), diff_data, mark);
  end
  
  legend(legend_lines);
  %ylabel("Erősítéshibával kompenzált különbség")
  ylabel("Gain error compensated difference")
  xlim([ts_ref(PLOTRANGE(1)) ts_ref(PLOTRANGE(end))]);
  
  grid on
  
  % plot timestamp errors
  subplot(3,1,3);
  
  legend_lines = {};
  
  for i = 2:len
    ts_err = ts(PLOTRANGE) - ts_ref(PLOTRANGE);
    plot(ts_ref(PLOTRANGE), ts_err * 1E+09, "-o");
    hold on
    
    legend_lines{i - 1, 1} = strcat(node_names{i,1}, "-", node_names{1,1});
  end
  
  grid on
  xlabel("Time [s]")
  %xlabel("Idő [s]");
  ylabel("Time error [ns]")
  %ylabel("Időhiba [ns]");
  legend(legend_lines);
  xlim([ts_ref(PLOTRANGE(1)) ts_ref(PLOTRANGE(end))]);
  
end