From 831cab7a5fbbccc5ecaf62afd712973712a71918 Mon Sep 17 00:00:00 2001
From: 4gwe <agoue.rau@gmail.com>
Date: Sun, 20 Sep 2026 16:49:46 +0200
Subject: [PATCH 1/8] test: add test for inestimable DCM trials

---
 test/pspm_dcm_test.m | 24 ++++++++++++++++++++++++
 1 file changed, 24 insertions(+)

diff --git a/test/pspm_dcm_test.m b/test/pspm_dcm_test.m
index 0d138e2c..d6eed7da 100644
--- a/test/pspm_dcm_test.m
+++ b/test/pspm_dcm_test.m
@@ -125,6 +125,7 @@ function valid_input(this)
       test_hra1_flex_cs(this);
       test_hra1_flex_cs_missing(this);
       test_hra1_flex_cs_nan(this);
+      test_hra1_flex_cs_inestimable(this);
     end
   end
   methods 
@@ -208,6 +209,29 @@ function test_hra1_flex_cs_nan(this)
 
       end
     end
+    function test_hra1_flex_cs_inestimable(this)
+        % find free filename
+        fn = pspm_find_free_fn(this.modelfile_prfx, '_inestimable.mat');
+
+        [df, ~] = this.get_hra_files(1);
+        [timing, eventnames, trialnames] = this.extract_hra_timings(1);
+
+        model = struct( ...
+            'modelfile', fn, ...
+            'datafile', df, ...
+            'timing', {timing}, ...
+            ... % force the last trial to be marked as inestimable
+            'lasttrialcutoff', 1e6 ...
+            );
+
+        options = struct( ...
+            'dispwin', 0, ...
+            'trlnames', {trialnames}, ...
+            'eventnames', {eventnames} ...
+            );
+
+        this.verifyWarning( @() pspm_dcm(model, options), 'ID:inestimable_trials');
+    end
     function [timing, eventnames, trialnames] = ...
         extract_hra_timings(this, subject)
       [s, c] = this.get_hra_files(subject);

From 0e76205df1af3eb11485d685d94ac3ceb2affb5d Mon Sep 17 00:00:00 2001
From: 4gwe <agoue.rau@gmail.com>
Date: Sun, 20 Sep 2026 17:05:09 +0200
Subject: [PATCH 2/8] Fix HP reward GLM config and runner (#909)

---
 src/pspm_init.m | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/pspm_init.m b/src/pspm_init.m
index 26aff0d8..8bd3ba7d 100644
--- a/src/pspm_init.m
+++ b/src/pspm_init.m
@@ -680,7 +680,7 @@
 defaults.glm(end + 1) = struct(...
     'modality',    'hp', ...
     'modelspec',   'hp_rew', ...
-    'cbf',          struct('fhandle', @bf_hprf_rew, 'args', []), ...
+    'cbf',          struct('fhandle', @pspm_bf_hprf_rew, 'args', []), ...
     'filter',       struct('lpfreq', 0.5, 'lporder', 4,  'hpfreq', 0.015, 'hporder', 4, 'down', 10,   'direction', 'bi'), ...
     'default', 0);
 % GLM for PS (fear-conditioning)

From e5cb7beeefef92fb13ab7a4013397412b4ce30a2 Mon Sep 17 00:00:00 2001
From: 4gwe <agoue.rau@gmail.com>
Date: Sun, 20 Sep 2026 17:08:36 +0200
Subject: [PATCH 3/8] typo

---
 src/pspm_cfg/pspm_cfg_glm_hp_rew.m | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/pspm_cfg/pspm_cfg_glm_hp_rew.m b/src/pspm_cfg/pspm_cfg_glm_hp_rew.m
index f5415e2d..c63a4202 100644
--- a/src/pspm_cfg/pspm_cfg_glm_hp_rew.m
+++ b/src/pspm_cfg/pspm_cfg_glm_hp_rew.m
@@ -1,5 +1,5 @@
 function [glm_hp_rew] = pspm_cfg_glm_hp_rew
-% GLM HP FC
+% GLM HP reward conditioning
 
 % Initialise
 global settings

From 8b7b41509c9a566889327119be86f6a3e746dbd0 Mon Sep 17 00:00:00 2001
From: 4gwe <agoue.rau@gmail.com>
Date: Sun, 20 Sep 2026 17:32:01 +0200
Subject: [PATCH 4/8] test: update

---
 test/pspm_convert_hb2hp_test.m | 15 ++++++++++-----
 1 file changed, 10 insertions(+), 5 deletions(-)

diff --git a/test/pspm_convert_hb2hp_test.m b/test/pspm_convert_hb2hp_test.m
index 5585aceb..1df5ea2d 100644
--- a/test/pspm_convert_hb2hp_test.m
+++ b/test/pspm_convert_hb2hp_test.m
@@ -75,6 +75,8 @@ function basic_conversion(this)
     [sts, infos, data, filestruct] = pspm_load_data(fn);
     this.verifyEqual(sts, 1);
     this.verifyEqual(data{outchannel}.header.chantype, 'hp');
+    this.verifyEqual(data{outchannel}.header.units, 'ms');
+    this.verifyEqual(data{outchannel}.header.sr, sr);
 end
 
 function too_strict_limits(this)
@@ -82,11 +84,6 @@ function too_strict_limits(this)
     sr = 1;
     options = struct('limit_lower', 11, 'limit_upper', 11);
 
-    % options = struct('limit_lower', 10, 'limit_upper', 11);
-    % this.verifyWarning(@() pspm_convert_hb2hp(fn, sr, options), ...
-    % 'ID:too_strict_limits');
-
-    % this.verifyWarningFree(@() pspm_convert_hb2hp(this.input_filename, sr));
     [sts, outchannel] = pspm_convert_hb2hp(fn, sr,options);
     this.verifyEqual(sts, 1);
     [sts, infos, data, filestruct] = pspm_load_data(fn);
@@ -94,7 +91,15 @@ function too_strict_limits(this)
     this.verifyEqual(data{outchannel}.header.chantype, 'hp');
     this.verifyTrue(any(isnan(data{outchannel}.data)))
 end
+function invalid_channel_action(this)
+    fn = this.input_filename;
+    sr = 100;
 
+    options = struct();
+    options.channel_action = 'abc';
+
+    this.verifyWarning( @() pspm_convert_hb2hp(fn, sr, options), 'ID:invalid_input');
+end
 function add_replace_channel_action(this)
     fn = this.input_filename;
     sr = 1;

From ad919d8fac383ef1c88760921ee18def4f99c7ef Mon Sep 17 00:00:00 2001
From: 4gwe <agoue.rau@gmail.com>
Date: Sun, 20 Sep 2026 20:51:09 +0200
Subject: [PATCH 5/8] Fix validation and eye selection in
 pspm_find_valid_fixations

---
 src/pspm_find_valid_fixations.m | 100 +++++++++++++++++++-------------
 1 file changed, 61 insertions(+), 39 deletions(-)

diff --git a/src/pspm_find_valid_fixations.m b/src/pspm_find_valid_fixations.m
index 5e1ef793..a2a60c24 100644
--- a/src/pspm_find_valid_fixations.m
+++ b/src/pspm_find_valid_fixations.m
@@ -20,17 +20,18 @@
 %   the circle around fixation point. Since this usage is currently
 %   considered secondary, it still requires a valid pupil channel as
 %   primary channel, even though unrelated to pupil analysis.
-%   In both usages, valid fixations can be outputted as additional channel.
+%   In both usages, an additional channel marking invalid fixations can be added.
 %   By default, screen centre is assumed as fixation point. If an explicit
 %   fixation point is given, the function assumes that the screen is
 %   perpendicular to the vector from the eye to the fixation point (which
 %   is approximately correct for large enough screen distance).
 % ● Format
-%   [sts, channel_index, fn] = pspm_find_valid_fixations(fn, bitmap, options)
-%   [sts, channel_index, fn] = pspm_find_valid_fixations(fn, circle_degree, distance, unit, options)
+%   [sts, pos_of_channel, fn] = pspm_find_valid_fixations(fn, bitmap, options)
+%   [sts, pos_of_channel, fn] = pspm_find_valid_fixations(fn, circle_degree, distance, unit, options)
 % ● Arguments
-%   *             fn : The actual data file containing the eyelink recording with gaze
-%                      data converted to cm.
+%   *             fn : Data file containing the eye-tracking recording. 
+%                      The corresponding gaze channels must be available in
+%                      distance units for fixation mode, or in distance or pixel units for bitmap mode.
 %   *         bitmap : A nxm matrix of the same size as the display, with 1
 %                      for valid and 0 for invalid gaze points. IMPORTANT: the bitmap has to
 %                      be defined in terms of the eyetracker coordinate system, i.e.
@@ -41,12 +42,13 @@
 %   *       distance : Distance between eye and screen in length units.
 %   *           unit : Unit in which distance is given.
 %   ┌────────options
-%   ├.fixation_point : A nx2 vector containing x and y of the fixation point (with respect
+%   ├.fixation_point : An nx2 matrix containing x and y of the fixation point (with respect
 %   │                  to the given resolution, and in the eyetracker coordinate system).
 %   │                  n should equal either 1 (constant fixation point) or the length
-%   │                  of the actual data. If resolution is not defined the values are
-%   │                  given in percent. Therefore (0.5 0.5) would correspond to the
-%   │                  middle of the screen. Default is (0.5 0.5). Only taken into account
+%   │                  of the actual data. If resolution is not defined, the
+%   │                  coordinates are given as fractions of the screen dimensions. 
+%   │                  Therefore (0.5 0.5) would correspond to the middle of the screen.
+%   │                  Default is (0.5 0.5). Only taken into account
 %   │                  if there is no bitmap.
 %   ├────.resolution : Resolution with which the fixation point is defined (Maximum value
 %   │                  of the x and y coordinates). This can be the screen resolution in
@@ -54,16 +56,18 @@
 %   │                  in cm (e.g. (50 30)). Default is (1 1). Only taken into account
 %   │                  if there is no bitmap.
 %   ├────.screen_dim : Only considered if .plot_gaze_coords is passed; used
-%   │                  plot the gaze data and circle on the actual screen
+%   │                  to plot the gaze data and circle on the actual screen
 %   │                  dimensions rather than using auto scaling. Input
 %   │                  should follow the format: [x_dim, y_dim].
 %   ├.plot_gaze_coords: Define whether to plot the gaze coordinates for visual
 %   │                 inspection of the validation process. Default is false.
-%   ├.channel_action: Define whether to add or replace the data. Default is
-%   │                 'add'. Possible values are 'add' or 'replace'
+%   ├.channel_action: Defines whether the processed channel is added as a new
+%   │                 channel or replaces the selected input channel. Default
+%   │                 is 'add'. Possible values are 'add' or 'replace'.
 %   ├───.add_invalid: [0/1] If this option is enabled, an extra channel will be
-%   │                 written containing information about the valid samples.
-%   │                 Data points equal to 1 correspond to invalid fixation.
+%   │                 added that marks invalid fixation samples.
+%   │                 Data points equal to 1 correspond to invalid fixations,
+%   │                 and data points equal to 0 correspond to valid fixations.
 %   │                 Default is not to add this channel.
 %   └───────.channel: Choose channels in which the data should be set to NaN
 %                     during invalid fixations. This can be a channel
@@ -89,7 +93,7 @@
 % ● Developer
 %   Additional i/o options for recursive calls are not included in the help.
 %   (1) fn can be a data structure as permitted by pspm_load_data,
-%   (2) the output argument pos_of_channels is an index of the channel(s)
+%   (2) the output argument pos_of_channel is an index of the channel(s)
 %   that was/were replaced or added
 %   (3) The third output argument is required for recursive calls
 % ● History
@@ -110,25 +114,27 @@
         ' You have to either pass a bitmap or circle_degree, distance and unit',...
         ' to compute the valid fixations']); return;
 end
-if numel(varargin{1}) > 1
+if numel(varargin) <= 2 && numel(varargin{1}) > 1
     mode = 'bitmap';
     bitmap = varargin{1};
-    if ~ismatrix(bitmap) || (~isnumeric(bitmap) && ~islogical(bitmap))
-        warning('ID:invalid_input', ['The bitmap must be a matrix and must',...
-            ' contain numeric or logical values.']); return;
+    if ~ismatrix(bitmap) || (~isnumeric(bitmap) && ~islogical(bitmap)) || any(~ismember(bitmap(:), [0 1]))        
+         warning('ID:invalid_input', ['The bitmap must be a numeric or logical matrix containing only 0 and 1.']); return;
     end
     if numel(varargin) < 2
         options = struct();
-        options.mode = 'bitmap';
     else
         options = varargin{2};
-        options.mode = 'bitmap';
+        if ~isstruct(options)
+            warning('ID:invalid_input', 'Options must be a struct.');
+            return;
+        end
     end
+    options.mode = 'bitmap';
 else
     mode = 'fixation';
     if numel(varargin) < 3
         warning('ID:invalid_input', ['Not enough input arguments.', ...
-            ' You have to set circle_degree, distance and unit',...
+            ' You have to provide circle_degree, distance and unit',...
             ' to compute the valid fixations']); return;
     end
     circle_degree = varargin{1};
@@ -146,11 +152,11 @@
             options.mode = 'fixation';
         end
     end
-    if ~isnumeric(circle_degree)
-        warning('ID:invalid_input', 'Circle_degree is not numeric.');
+    if ~(isnumeric(circle_degree) && isscalar(circle_degree) && isreal(circle_degree) && isfinite(circle_degree) && circle_degree >= 0)
+        warning('ID:invalid_input', 'Circle_degree must be a finite, non-negative numeric scalar.');
         return;
-    elseif ~isnumeric(distance)
-        warning('ID:invalid_input', 'Distance is not set or not numeric.');
+    elseif ~(isnumeric(distance) && isscalar(distance) && isreal(distance) && isfinite(distance) && distance > 0)
+        warning('ID:invalid_input', 'Distance must be a finite, positive numeric scalar.');
         return;
     elseif ~ischar(unit)
         warning('ID:invalid_input', 'Unit should be a char.');
@@ -164,6 +170,7 @@
         [nsts,distance] = pspm_convert_unit(distance,unit ,'mm');
         if nsts~=1
             warning('ID:invalid_input', 'Failed to convert distance to mm.');
+            return;
         end
     end
 end
@@ -194,15 +201,25 @@
     elseif strcmpi(mode, 'bitmap')
         channels_correct_units = find(~contains(channelunits_list, 'degree'));
     end
+    
     gazedata = struct('infos', alldata.infos, 'data', {alldata.data(channels_correct_units)});
 
     [sts_gaze, gaze_x, gaze_y, eye] = pspm_load_gaze(gazedata, data.header.chantype);
 
     if sts_gaze < 1
         warning('ID:invalid_input', ['Unable to perform gaze ', ...
-          'validation. Cannot find gaze channels with distance ',...
-          'unit values. Maybe you need to convert them with ', ...
-          'pspm_convert_gaze()']);
+            'validation. Cannot find corresponding gaze channels ', ...
+            'with compatible units. Maybe you need to convert them ', ...
+            'with pspm_convert_gaze().']);
+        return;
+    end
+
+    if numel(gaze_x.data) ~= numel(gaze_y.data) ||  numel(gaze_x.data) ~= numel(data.data)
+        warning('ID:invalid_input', ['Pupil and corresponding gaze channels must have ', 'the same number of data points.']);
+        return;
+    end
+    if gaze_x.header.sr ~= gaze_y.header.sr || gaze_x.header.sr ~= data.header.sr
+        warning('ID:invalid_input', ['Pupil and corresponding gaze channels must have ', 'the same sampling rate.']);
         return;
     end
 
@@ -213,13 +230,17 @@
         case 'fixation'
             % expand fixation point to size of data
             fix_point = options.fixation_point;
-            if size(fix_point, 1) == 1
+            if size(fix_point, 2) ~= 2
+                warning('ID:invalid_input', ...
+                    'Fixation point must have exactly two columns for x and y coordinates.');
+                return;
+            elseif size(fix_point, 1) == 1
                 fix_point = repmat(fix_point(:)', numel(gaze_x.data), 1);
-            elseif size(fix_point, 1) ~= numel(gaze_x)
+            elseif size(fix_point, 1) ~= numel(gaze_x.data)    
                 warning('ID:invalid_input', ['Fixation point has wrong ', ...
                     'dimensions - it should be 1x2 or nx2 where n is the ', ...
                     'number of gaze data points.']);
-                return
+                return;
             end
 
             % normalise fixation point to fraction of full screen
@@ -398,16 +419,17 @@
     if (rsts(1) < 1 && rsts(2) < 1)
         return;
     elseif (rsts(1) < 1 || rsts(2) < 1)
+        successful_idx = find(rsts > 0);
         pos_of_channel(rsts < 1) = [];
+
+        [~, best_eye] = pspm_find_eye(channels{successful_idx});
+        alldata.infos.source.best_eye = best_eye;
     else
-        % update best eye
-        eye_stat = Inf(1,numel(alldata.infos.source.eyesObserved));
-        for i = 1:numel(alldata.infos.source.eyesObserved)
-            e_stat = alldata.infos.source.chan_stats(pos_of_channel);
-            eye_stat(i) = max(cellfun(@(x) x.nan_ratio, e_stat));
-        end
+        eye_stat = cellfun(@(x) x.nan_ratio, alldata.infos.source.chan_stats(pos_of_channel));
+        
         [~, min_idx] = min(eye_stat);
-        alldata.infos.source.best_eye = lower(alldata.infos.source.eyesObserved(min_idx));
+        [~, best_eye] = pspm_find_eye(channels{min_idx});
+        alldata.infos.source.best_eye = best_eye;
     end
 end
 

From eac1b789e767813adf927eb325e7c8012341e4cd Mon Sep 17 00:00:00 2001
From: 4gwe <agoue.rau@gmail.com>
Date: Tue, 22 Sep 2026 12:38:48 +0200
Subject: [PATCH 6/8] nanmean -> mean([tmp_data.sem],2, 'omitnan')

---
 src/pspm_tam.m | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/pspm_tam.m b/src/pspm_tam.m
index f5008a21..98749c74 100644
--- a/src/pspm_tam.m
+++ b/src/pspm_tam.m
@@ -162,7 +162,7 @@
 
 % Checking if the sampling rate is the same for all samples.
 if n_file > 1 && any(diff(sr) > 0)
-  if model.filter.down > min(sr)) ||...                                    % if filter.down is less than the minimal sr
+  if (model.filter.down > min(sr)) ||...                                    % if filter.down is less than the minimal sr
       strcmpi(model.filter.down,'none')                                    % if filter.down is none
     model.filter.down = min(sr);
     fprintf('\nSampling rate differs between sessions. Data will be downsampled.\n')
@@ -239,19 +239,19 @@
 if exist('std_exp_cond','var')
   tmp_data = [segm{:,std_exp_cond.ind}];
 
-  std_exp_cond.data = nanmean([tmp_data.mean],2);
-  std_exp_cond.std = nanmean([tmp_data.std],2);
-  std_exp_cond.sem = nanmean([tmp_data.sem],2);
+  std_exp_cond.data = mean([tmp_data.mean],2, 'omitnan'); %nanmean([tmp_data.mean],2);
+  std_exp_cond.std =  mean([tmp_data.std],2, 'omitnan');
+  std_exp_cond.sem =  mean([tmp_data.sem],2, 'omitnan');
 end
 
 for i=1:n_exp_cond
 
   tmp_data = [segm{:,i}];
 
-  tmp_data_new.data = nanmean([tmp_data.mean],2);
-  tmp_data_new.std = nanmean([tmp_data.std],2);
-  tmp_data_new.sem = nanmean([tmp_data.sem],2);
-  tmp_data_new.t = nanmean([tmp_data.t],2);
+  tmp_data_new.data = mean([tmp_data.mean],2, 'omitnan');
+  tmp_data_new.std = mean([tmp_data.std],2, 'omitnan');
+  tmp_data_new.sem = mean([tmp_data.sem],2, 'omitnan');
+  tmp_data_new.t = mean([tmp_data.t],2, 'omitnan');
 
   % Subtracting the standard experimental condition
   if exist('std_exp_cond','var') && i~=std_exp_cond.ind

From ad8973f64f9a70c003fa308be63601c94934999c Mon Sep 17 00:00:00 2001
From: 4gwe <agoue.rau@gmail.com>
Date: Wed, 23 Sep 2026 17:36:35 +0200
Subject: [PATCH 7/8] Fix TAM filtering defaults, segment extraction, and
 outputs

---
 src/pspm_check_model.m |  11 +++-
 src/pspm_glm.m         |   2 +-
 src/pspm_init.m        |   8 +--
 src/pspm_tam.m         |  48 ++++++++++-------
 test/pspm_tam_test.m   | 115 +++++++++++++++++++++++++++++++++++++++++
 5 files changed, 158 insertions(+), 26 deletions(-)
 create mode 100644 test/pspm_tam_test.m

diff --git a/src/pspm_check_model.m b/src/pspm_check_model.m
index 320d4c2b..768d28b5 100644
--- a/src/pspm_check_model.m
+++ b/src/pspm_check_model.m
@@ -519,8 +519,15 @@
     end
 end
 
-if ~isfield(model.filter, 'down') || ~isnumeric(model.filter.down)
-    warning('ID:invalid_input', 'Filter structure needs a numeric ''down'' field.'); return;
+if ~isfield(model.filter, 'down')
+    warning('ID:invalid_input', 'Filter structure needs a ''down'' field.'); return;
+end
+
+valid_down = (isnumeric(model.filter.down) && isscalar(model.filter.down) && isfinite(model.filter.down) && model.filter.down > 0) ...
+          || (ischar(model.filter.down) &&   strcmpi(model.filter.down, 'none'));
+
+if ~valid_down    
+    warning('ID:invalid_input',  'Filter structure needs a positive numeric ''down'' field or ''none''.'); return;
 end
 
 
diff --git a/src/pspm_glm.m b/src/pspm_glm.m
index a9426544..4f5eb401 100644
--- a/src/pspm_glm.m
+++ b/src/pspm_glm.m
@@ -1,4 +1,4 @@
-    function [sts, glm] = pspm_glm(model, options)
+function [sts, glm] = pspm_glm(model, options)
 % ● Description
 %   pspm_glm specifies a within-subject general linear convolution model 
 %   (GLM) of predicted signals and calculates amplitude estimates for these 
diff --git a/src/pspm_init.m b/src/pspm_init.m
index 8bd3ba7d..aa50b2a2 100644
--- a/src/pspm_init.m
+++ b/src/pspm_init.m
@@ -756,13 +756,15 @@
   'cif', struct('fhandle', @pspm_bf_ldrf_gm, 'args', [0, 2.76 , 0.09 , 0.31],...     % input function & default parameters
   'lb', [0,0,0,0], 'ub', [0,Inf,Inf,Inf]),...                          % & the lower/upper bounds
   'filter', struct('lpfreq', 'none', 'lporder', 0,  ...                               % default filter
-  'hpfreq', 'none', 'hporder', 0, 'down', 0, 'direction', 'bi'));
+  'hpfreq', 'none', 'hporder', 0, 'down', 'none', 'direction', 'bi'));
 defaults.tam(2) = struct(...
   'modality',  'pupil', ...
   'modelspec', 'constriction',...
   'cbf', struct('fhandle', @pspm_bf_lcrf_gm, 'args', [0.2, 3.24 , 0.18 , 0.43]),...
-  'cif', struct('fhandle', @pspm_bf_lcrf_gm, 'args', [0, 2.76 , 0.09 , 0.31], 'lb', [0,0,0,0], 'ub', [0,Inf,Inf,Inf]),...
-  'filter', struct('lpfreq', 'none', 'lporder', 0, 'hpfreq', 'none', 'hporder', 0, 'down', 0, 'direction', 'bi'));
+  'cif', struct('fhandle', @pspm_bf_lcrf_gm, 'args', [0, 2.76 , 0.09 , 0.31], ...
+  'lb', [0,0,0,0], 'ub', [0,Inf,Inf,Inf]),...
+  'filter', struct('lpfreq', 'none', 'lporder', 0, ...
+  'hpfreq', 'none', 'hporder', 0, 'down', 'none', 'direction', 'bi'));
 %% 9 FIRST LEVEL settings
 % 9.1 allowed first level model types
 defaults.first = {'glm', 'sf', 'dcm', 'tam'};
diff --git a/src/pspm_tam.m b/src/pspm_tam.m
index 98749c74..74e177d6 100644
--- a/src/pspm_tam.m
+++ b/src/pspm_tam.m
@@ -1,4 +1,4 @@
-function tam = pspm_tam(model, options)
+function [sts, tam] = pspm_tam(model, options)
 % ● Description
 %   TAM stands for Trial Average Model and allows to fit models on trial-averaged data.
 %   pspm_tam starts by extracting and averaging signal segments of length `model.window`
@@ -161,10 +161,11 @@
 oldsr = sr;
 
 % Checking if the sampling rate is the same for all samples.
-if n_file > 1 && any(diff(sr) > 0)
-  if (model.filter.down > min(sr)) ||...                                    % if filter.down is less than the minimal sr
-      strcmpi(model.filter.down,'none')                                    % if filter.down is none
-    model.filter.down = min(sr);
+if n_file > 1 && any(diff(sr) ~= 0)    % not any(diff(sr) > 0) 
+  if (ischar(model.filter.down) && strcmpi(model.filter.down, 'none'))  || ...    % if filter.down is none         
+     (isnumeric(model.filter.down) && model.filter.down > min(sr))                            % if filter.down is less than the minimal sr
+    
+      model.filter.down = min(sr);
     fprintf('\nSampling rate differs between sessions. Data will be downsampled.\n')
   end
 else
@@ -194,14 +195,14 @@
     extrsgopt.marker_chan = markers(k);
   end
 
-  [lsts, s] = pspm_extract_segments('manual', y(k), sr(k), model.timing(k), extrsgopt);
+  [lsts, s] = pspm_extract_segments('data', y{k}, sr(k), model.timing{k}, extrsgopt); % 'markers'???
   if lsts<1, warning('ID:error_extract_segments','An error occured in pspm_extract_segments.'); return; end
 
   for i=1:n_exp_cond
-    tmp_data.mean = s.segments{i,1}.mean;
-    tmp_data.std = s.segments{i,1}.std;
-    tmp_data.sem = s.segments{i,1}.sem;
-    tmp_data.t = s.segments{i,1}.t;
+    tmp_data.mean = s.segments{i,1}.mean; % why not eg 50x1 
+    tmp_data.std = s.segments{i,1}.std;% why not eg 50x1
+    tmp_data.sem = s.segments{i,1}.sem;% why not eg 50x1
+    tmp_data.t = s.segments{i,1}.t; 
     % a cell array of struct and of size (n_file x n_exp_cond) where each
     % line correspond to a given file and each column to an
     % experimental condition
@@ -228,8 +229,9 @@
 end
 
 % changing the sampling rate
-sr = model.filter.down*ones(size(sr));
-
+if isnumeric(model.filter.down)
+    sr = model.filter.down * ones(size(sr)); % needs to be tested
+end
 
 %% Determining mean values
 fprintf('Preparing for fitting ...\n')
@@ -271,7 +273,7 @@
     tmp_data_new.sem = tmp_data_new.sem + tmp_data_new.sem(tmp_max_ind); % the error adds up
   end
 
-  mean{1,i} = tmp_data_new;
+  mean_data{1,i} = tmp_data_new;
 
   clear tmp_data tmp_data_new tmp_max tmp_max_ind
 end
@@ -280,7 +282,7 @@
 fprintf('Fitting ...\n')
 
 for i=1:n_exp_cond
-  raw_y = mean{1,i}.data;
+  raw_y = mean_data{1,i}.data;
 
   n = model.window;
   td = n / length(raw_y);
@@ -299,7 +301,7 @@
   % Minimization of RSS
   warning off all
   [~, fitted{1,i}.optargs, fitted{1,i}.fval, sts, fmincon_output] = ...
-    evalc('fmincon(RSS,model.if.args,[],[],[],[],model.if.lb,model.if.ub)');
+    evalc('fmincon(RSS,model.if.args,[],[],[],[],model.if.lb,model.if.ub)'); % lb ub for opt.
   warning on all
   if sts == 0
     warning('ID:fmincon',['During the fitting process, ''fmincon'' exceeded', ...
@@ -320,7 +322,7 @@
   % Calculating the predicted signal that will be included in the output structure
   fitted{1,i}.data = predicted_y(fitted{1,i}.optargs);
   % Cutting away tail
-  tmp_y = mean{1,i}.data;
+  tmp_y = mean_data{1,i}.data;
   fitted{1,i}.data(size(tmp_y,1)+1:end) = [];
 
 end
@@ -336,9 +338,15 @@
 tam.bf            = model.bf;
 tam.if            = model.if;
 
+% Check if filtered
+no_lp = ischar(model.filter.lpfreq) && strcmpi(model.filter.lpfreq, 'none');
+no_hp = ischar(model.filter.hpfreq) && strcmpi(model.filter.hpfreq, 'none');
+no_down = ischar(model.filter.down) && strcmpi(model.filter.down, 'none');
+filtered = ~(no_lp && no_hp && no_down);
+
 % Collecting fitting data
-tmp_mean = [mean{1,:}];
-tam.data.Y        = {tmp_mean.data};
+tmp_mean = [mean_data{1,:}];
+tam.data.Y        = {tmp_mean.data}; % measured data
 tam.data.X        = {tmp_mean.t};
 tam.data.std      = {tmp_mean.std};
 tam.data.sem      = {tmp_mean.sem};
@@ -358,12 +366,12 @@
 tmp_fitted = [fitted{1,:}];
 tam.fit.Y         = {tmp_fitted.data};
 tam.fit.X         = {tmp_mean.t};
-tam.fit.rss       = {tmp_fitted.fval};  % RSS (residual sum square)
+tam.fit.rss       = {tmp_fitted.fval};  
 tam.fit.args      = {tmp_fitted.optargs};
 tam.fit.sr        = num2cell(sr(:).');
 
 tam.infos.duration     = model.window;
-tam.infos.durationinfo = 'duration in seconds';
+tam.infos.durationinfo = 'duration in seconds'; % not allways true! -> timeunits='samples'
 
 tam.timing        = model.timing;
 
diff --git a/test/pspm_tam_test.m b/test/pspm_tam_test.m
new file mode 100644
index 00000000..5a93f0a5
--- /dev/null
+++ b/test/pspm_tam_test.m
@@ -0,0 +1,115 @@
+classdef pspm_tam_test < matlab.unittest.TestCase
+
+methods (Test)
+
+function testSingleConditionTrialAverage(testCase)
+
+    %% Parameters
+    sr = 10;
+    duration = 50;
+    window = 5;
+
+    %% Known trial response
+    t = (0:1/sr:window-1/sr)';
+
+    response = exp( ...
+        -((t - 2).^2) / (2 * 0.5^2));
+
+    %% Continuous pupil signal
+    y = zeros(duration * sr, 1);
+
+    onsets = [10 20 30];
+
+    for iTrial = 1:numel(onsets)
+
+        startSample = round(onsets(iTrial) * sr) + 1;
+        stopSample = startSample + numel(response) - 1;
+
+        y(startSample:stopSample) = ...
+            y(startSample:stopSample) + response;
+
+    end
+
+    %% Timing
+    timing.names = {'condition_1'};
+    timing.onsets = {onsets};
+    timing.durations = {zeros(size(onsets))};
+
+    %% Create PsPM pupil file
+    data = cell(1,1);
+
+    data{1}.header = struct( ...
+        'chantype', 'pupil', ...
+        'sr', sr, ...
+        'units', 'a.u.');
+
+    data{1}.data = y(:);
+
+    infos.duration = duration;
+    infos.durationinfo = 'Duration in seconds';
+    infos.source = struct();
+
+    datafile = [tempname '.mat'];
+    save(datafile, 'data', 'infos');
+    
+
+
+
+    testCase.addTeardown(@() deleteIfExists(datafile));
+
+    %% Verify test file itself
+    [sts, loaded] = pspm_load_channel( datafile, 'pupil', 'pupil');
+
+    testCase.verifyEqual(sts, 1);
+    testCase.verifyEqual(loaded.header.sr, sr);
+    testCase.verifyEqual(loaded.data(:), y(:));
+
+
+    %% Configure TAM model
+    modelfile = [tempname '.mat'];
+    testCase.addTeardown(@() deleteIfExists(modelfile));
+
+    model = struct();
+
+    model.modelfile = modelfile;
+    model.datafile = {datafile};
+    model.timing = {timing};
+
+    model.timeunits = 'seconds';
+    model.window = window;
+
+    model.modelspec = 'dilation';
+    model.modality = 'pupil';
+    model.channel = 'pupil';
+
+    % Keep preprocessing simple for this test
+    model.norm = 0;
+    model.baseline = 0;
+    model.norm_max = 0;
+    % model.std_exp_cond = 'none';
+
+    options = struct();
+    options.overwrite = 1;
+
+    %% Run TAM
+    [sts, tam] = pspm_tam(model, options);
+
+    %% Basic output checks
+    testCase.verifyTrue(isstruct(tam));
+    testCase.verifyTrue(isfield(tam, 'data'));
+    testCase.verifyTrue(isfield(tam.data, 'Y'));
+    testCase.verifyEqual(numel(tam.data.Y), 1);
+end
+
+end
+
+end
+
+
+function deleteIfExists(filename)
+
+if exist(filename, 'file')
+delete(filename);
+end
+
+end
\ No newline at end of file

From 8cf179d1c68dece3e9482dd4e212fe19cc55eae1 Mon Sep 17 00:00:00 2001
From: 4gwe <agoue.rau@gmail.com>
Date: Thu, 24 Sep 2026 15:19:28 +0200
Subject: [PATCH 8/8] Fix TAM preprocessing and standard condition handling

---
 src/pspm_check_model.m |  26 ++-
 src/pspm_tam.m         |  69 +++++---
 test/pspm_tam_test.m   | 357 +++++++++++++++++++++++++++++++++++++----
 3 files changed, 381 insertions(+), 71 deletions(-)

diff --git a/src/pspm_check_model.m b/src/pspm_check_model.m
index 768d28b5..caf7c52a 100644
--- a/src/pspm_check_model.m
+++ b/src/pspm_check_model.m
@@ -441,22 +441,20 @@
     ' index corresponding to it.'];
   if ~isfield(model,'std_exp_cond')
     model.std_exp_cond = 'none';
-  elseif ~ischar(model.std_exp_cond) && ~isnumeric(model.std_exp_cond)
-    warning('ID:invalid_input',std_cond_war_msg); return;
   elseif ischar(model.std_exp_cond)
-    tmp_ind = cellfun(@(x) strcmpi(model.std_exp_cond,x), model.timing{1}.names);
-    if ~any(tmp_ind)
-      warning('ID:invalid_input',std_cond_war_msg); return;
-    end
-    std_exp_cond.name = model.std_exp_cond;
-    std_exp_cond.ind = find(tmp_ind);
+      % 'none' is explicitly allowed
+      if ~strcmpi(model.std_exp_cond, 'none')
+          tmp_ind = strcmpi(model.std_exp_cond, model.timing{1}.names);
+          if ~any(tmp_ind); warning('ID:invalid_input', std_cond_war_msg);   return;    end
+      end
   elseif isnumeric(model.std_exp_cond)
-    if model.std_exp_cond < 1 || ...
-        model.std_exp_cond > numel(model.timing{1}.names)
-      warning('ID:invalid_input',std_cond_war_msg); return;
-    end
-    std_exp_cond.name = model.timing{1}.names(model.std_exp_cond);
-    std_exp_cond.ind = model.std_exp_cond;
+      % Must be one valid scalar integer index
+      if ~isscalar(model.std_exp_cond) || ~isfinite(model.std_exp_cond) || ...
+          model.std_exp_cond < 1 || model.std_exp_cond > numel(model.timing{1}.names)
+          warning('ID:invalid_input', std_cond_war_msg); return;
+      end
+  else
+    warning('ID:invalid_input',std_cond_war_msg); return;  
   end
   clear std_cond_war_msg tmp_ind
 
diff --git a/src/pspm_tam.m b/src/pspm_tam.m
index 74e177d6..e2d51b77 100644
--- a/src/pspm_tam.m
+++ b/src/pspm_tam.m
@@ -12,7 +12,7 @@
 %   ├─────.timing:  a multiple condition file name (single session) OR
 %   │               a cell array of multiple condition file names OR
 %   │               a struct (single session) with fields .names, .onsets,
-%   │               and (optional) .durations OR
+%   │               and (optional) duration OR
 %   │               a cell array of struct OR
 %   │               a struct with fields 'markerinfos', 'markervalues',
 %   │               'names' OR
@@ -83,6 +83,7 @@
 %   Introduced In PsPM 4.2
 %   Written in 2020 by Ivan Rojkov (University of Zurich)
 %   Maintained in 2022 by Teddy
+%   Maintained in 2026 by Bernhard A. von Raußendorf
 % ● Developer
 %   The fitting process is a residual least square minimisation where the
 %   predicted value is calculated as following:
@@ -126,6 +127,18 @@
     return
 end
 
+% 2.3 check ....
+std_exp_cond = [];
+
+if ~(ischar(model.std_exp_cond) && strcmpi(model.std_exp_cond, 'none'))
+    if ischar(model.std_exp_cond)
+        std_exp_cond.ind = find(strcmpi( model.std_exp_cond, model.timing{1}.names), 1);
+    else
+        std_exp_cond.ind = model.std_exp_cond;
+    end
+    std_exp_cond.name = model.timing{1}.names{std_exp_cond.ind};
+end
+
 %% Loading files
 
 fprintf('Computing Trial Average Model: %s \n', model.modelfile);
@@ -163,7 +176,7 @@
 % Checking if the sampling rate is the same for all samples.
 if n_file > 1 && any(diff(sr) ~= 0)    % not any(diff(sr) > 0) 
   if (ischar(model.filter.down) && strcmpi(model.filter.down, 'none'))  || ...    % if filter.down is none         
-     (isnumeric(model.filter.down) && model.filter.down > min(sr))                            % if filter.down is less than the minimal sr
+     (isnumeric(model.filter.down) && model.filter.down > min(sr))                % if filter.down is less than the minimal sr
     
       model.filter.down = min(sr);
     fprintf('\nSampling rate differs between sessions. Data will be downsampled.\n')
@@ -202,7 +215,7 @@
     tmp_data.mean = s.segments{i,1}.mean; % why not eg 50x1 
     tmp_data.std = s.segments{i,1}.std;% why not eg 50x1
     tmp_data.sem = s.segments{i,1}.sem;% why not eg 50x1
-    tmp_data.t = s.segments{i,1}.t; 
+    tmp_data.t = s.segments{i,1}.t; % 
     % a cell array of struct and of size (n_file x n_exp_cond) where each
     % line correspond to a given file and each column to an
     % experimental condition
@@ -215,33 +228,44 @@
 
 %% Downsample the data
 % if a filter was specified or if the data differ in sr
-fprintf('Filtering ...\n')
-for i = 1:n_exp_cond
-    for k = 1:n_file
 
-        model.filter.sr = sr(k);
+fprintf('Filtering ...\n') % maybe only if there is a filtering?
 
-        [lsts, segm{k,i}, ~] = structfun(@(x) pspm_prepdata(x, model.filter),segm{k,i},'UniformOutput',false);
-        if any(structfun(@(x) x<1,lsts)), warning('ID:error_prepdata','An error occured in pspm_prepdata.'); return; end
 
-        clear new_sr lsts
+fields = {'mean', 'std', 'sem'};
+for i = 1:n_exp_cond
+    for k = 1:n_file
+        model.filter.sr = sr(k); % adds the corresponding sr to the filer
+        
+        for f = 1:numel(fields)
+            field = fields{f};
+            [lsts, segm{k,i}.(field), new_sr(k)] = pspm_prepdata(segm{k,i}.(field), model.filter);
+            if lsts < 1; warning('ID:error_prepdata', 'An error occured in pspm_prepdata.');  return; end
+        end
+        [llsts, segm{k,i}.t , new_sr(k) ] = pspm_downsample(segm{k,i}.t,  sr(k), new_sr(k) ); % does not need to be filtered!
+        if llsts < 1; warning('ID:error_downsample', 'An error occured in pspm_downsample.');  return; end
+
     end
 end
 
-% changing the sampling rate
-if isnumeric(model.filter.down)
-    sr = model.filter.down * ones(size(sr)); % needs to be tested
-end
+sr = new_sr; % maybe test upstream if all new_sr are the same!
+% % changing the sampling rate
+% if isnumeric(model.filter.down)
+%     sr = model.filter.down * ones(size(new_sr)); % needs to be tested why not the values of prepdata?
+% end
+
+clear new_sr
 
 %% Determining mean values
 fprintf('Preparing for fitting ...\n')
 
 baseline_index = floor(sr(1)*model.baseline)+1;
 
-if exist('std_exp_cond','var')
+
+if ~isempty(std_exp_cond)
   tmp_data = [segm{:,std_exp_cond.ind}];
 
-  std_exp_cond.data = mean([tmp_data.mean],2, 'omitnan'); %nanmean([tmp_data.mean],2);
+  std_exp_cond.data = mean([tmp_data.mean],2, 'omitnan');
   std_exp_cond.std =  mean([tmp_data.std],2, 'omitnan');
   std_exp_cond.sem =  mean([tmp_data.sem],2, 'omitnan');
 end
@@ -256,7 +280,7 @@
   tmp_data_new.t = mean([tmp_data.t],2, 'omitnan');
 
   % Subtracting the standard experimental condition
-  if exist('std_exp_cond','var') && i~=std_exp_cond.ind
+  if ~isempty(std_exp_cond) && i~=std_exp_cond.ind
     tmp_data_new.data = tmp_data_new.data - std_exp_cond.data;
     tmp_data_new.std = tmp_data_new.std + std_exp_cond.std;      % the error adds up
     tmp_data_new.sem = tmp_data_new.sem + std_exp_cond.sem;      % the error adds up
@@ -266,11 +290,12 @@
   tmp_data_new.data = tmp_data_new.data - tmp_data_new.data(baseline_index);
 
   % Dividing by the max value
-  if model.norm_max
-    [tmp_max,tmp_max_ind] = max(tmp_data_new.data);
+  if model.norm_max 
+    [tmp_max,tmp_max_ind] = max(tmp_data_new.data); % this does not give us the first peak!!!
+    % maybe add     if tmp_max == 0 
     tmp_data_new.data = tmp_data_new.data/tmp_max;
-    tmp_data_new.std = tmp_data_new.std + tmp_data_new.std(tmp_max_ind); % the error adds up
-    tmp_data_new.sem = tmp_data_new.sem + tmp_data_new.sem(tmp_max_ind); % the error adds up
+    tmp_data_new.std = tmp_data_new.std + tmp_data_new.std(tmp_max_ind); % the error adds up . why not /abs(tmp_max)
+    tmp_data_new.sem = tmp_data_new.sem + tmp_data_new.sem(tmp_max_ind); % the error adds up . why not /abs(tmp_max)
   end
 
   mean_data{1,i} = tmp_data_new;
@@ -355,7 +380,7 @@
 tam.data.normd  = model.norm;
 tam.data.norm     = model.norm;
 
-if exist('std_exp_cond','var')
+if ~isempty(std_exp_cond)
   tam.data.std_exp_cond.name  = std_exp_cond.name;
   tam.data.std_exp_cond.ind   = std_exp_cond.ind;
 else
diff --git a/test/pspm_tam_test.m b/test/pspm_tam_test.m
index 5a93f0a5..171526ab 100644
--- a/test/pspm_tam_test.m
+++ b/test/pspm_tam_test.m
@@ -1,48 +1,207 @@
 classdef pspm_tam_test < matlab.unittest.TestCase
 
+
+properties (TestParameter)
+    nFiles = struct( ...
+        'oneFile', 1, ...
+        'twoFiles', 2);
+end
+
 methods (Test)
 
-function testSingleConditionTrialAverage(testCase)
+function testSingleConditionTrialAverage(testCase, nFiles)
 
-    %% Parameters
     sr = 10;
     duration = 50;
     window = 5;
 
-    %% Known trial response
     t = (0:1/sr:window-1/sr)';
+    response = exp(-((t - 2).^2) / (2 * 0.5^2));
 
-    response = exp( ...
-        -((t - 2).^2) / (2 * 0.5^2));
-
-    %% Continuous pupil signal
     y = zeros(duration * sr, 1);
-
     onsets = [10 20 30];
 
     for iTrial = 1:numel(onsets)
-
         startSample = round(onsets(iTrial) * sr) + 1;
         stopSample = startSample + numel(response) - 1;
 
         y(startSample:stopSample) = ...
             y(startSample:stopSample) + response;
+    end
+
+    timing.names = {'condition_1'};
+    timing.onsets = {onsets};
+    timing.durations = {zeros(size(onsets))};
+
+    datafiles = cell(1, nFiles);
+
+    for k = 1:nFiles
 
+        data = cell(1,1);
+        data{1}.header = struct( ...
+            'chantype', 'pupil', ...
+            'sr', sr, ...
+            'units', 'a.u.');
+
+        data{1}.data = y(:);
+
+        infos.duration = duration;
+        infos.durationinfo = 'Duration in seconds';
+        infos.source = struct();
+
+        datafiles{k} = [tempname '.mat'];
+
+        save(datafiles{k}, 'data', 'infos');
+
+        testCase.addTeardown( ...
+            @() deleteIfExists(datafiles{k}));
+    end
+
+    modelfile = [tempname '.mat'];
+    testCase.addTeardown(@() deleteIfExists(modelfile));
+
+    model = struct();
+    model.modelfile = modelfile;
+
+    model.datafile = datafiles;
+
+    % Same timing definition for each file/session
+    model.timing = repmat({timing}, 1, nFiles);
+
+    model.timeunits = 'seconds';
+    model.window = window;
+    model.modelspec = 'dilation';
+    model.modality = 'pupil';
+    model.channel = 'pupil';
+    model.norm = 0;
+    model.baseline = 0;
+    model.norm_max = 0;
+
+    options.overwrite = 1;
+
+    [sts, tam] = pspm_tam(model, options);
+
+    testCase.verifyEqual(sts, 1);
+
+    expected = response - response(1);
+
+    testCase.verifyEqual( ...
+        tam.data.Y{1}, ...
+        expected(:), ...
+        'AbsTol', 1e-10);
+
+    % Sampling rate should stay 10 Hz for every session
+    testCase.verifyEqual( ...
+        tam.data.sr, ...
+        repmat({sr}, 1, nFiles));
+
+end
+
+function testDownsamplingTo5Hz(testCase)
+    sr = 10;
+    duration = 50;
+    window = 5;
+
+    % 5 seconds at 10 Hz -> 50 samples
+    t = (0:1/sr:window-1/sr)';
+    response = exp(-((t - 2).^2) / (2 * 0.5^2));
+
+    y = zeros(duration * sr, 1);
+    onsets = [10 20 30];
+
+    for iTrial = 1:numel(onsets)
+        startSample = round(onsets(iTrial) * sr) + 1;
+        stopSample = startSample + numel(response) - 1;
+
+        y(startSample:stopSample) = y(startSample:stopSample) + response;
     end
 
-    %% Timing
     timing.names = {'condition_1'};
     timing.onsets = {onsets};
     timing.durations = {zeros(size(onsets))};
 
-    %% Create PsPM pupil file
     data = cell(1,1);
+    data{1}.header = struct( 'chantype', 'pupil', 'sr', sr, 'units', 'a.u.');
+    data{1}.data = y(:);
+
+    infos.duration = duration;
+    infos.durationinfo = 'Duration in seconds';
+    infos.source = struct();
+
+    datafile = [tempname '.mat'];
+    save(datafile, 'data', 'infos');
+    testCase.addTeardown(@() deleteIfExists(datafile));
+
+    modelfile = [tempname '.mat'];
+    testCase.addTeardown(@() deleteIfExists(modelfile));
+
+    model = struct();
+    model.modelfile = modelfile;
+    model.datafile = {datafile};
+    model.timing = {timing};
+    model.timeunits = 'seconds';
+    model.window = window;
+    model.modelspec = 'dilation';
+    model.modality = 'pupil';
+    model.channel = 'pupil';
+    model.norm = 0;
+    model.baseline = 0;
+    model.norm_max = 0;
+
+    % Explicitly request downsampling from 10 Hz to 5 Hz no filtering
+    model.filter = struct( 'lpfreq', 'none', 'lporder', 0, 'hpfreq', 'none', 'hporder', 0, 'down', 5, 'direction', 'bi');
+
+    options = struct();
+    options.overwrite = 1;
+
+    [sts, tam] = pspm_tam(model, options);
+
+    testCase.assertEqual(sts, 1);
+
+    % 5 seconds at 5 Hz -> 25 samples
+    expectedSize = [25 1];
+
+    testCase.verifySize(tam.data.Y{1},   expectedSize); % mean
+    testCase.verifySize(tam.data.std{1}, expectedSize);
+    testCase.verifySize(tam.data.sem{1}, expectedSize);
+    testCase.verifySize(tam.data.X{1},   expectedSize); % t
+
+    testCase.verifyEqual(tam.data.sr{1}, 5);
+
+end
+
+function testLowpassFilterWithoutDownsampling(testCase)
+
+    sr = 10;
+    duration = 50;
+    window = 5;
 
+    t = (0:1/sr:window-1/sr)';
+
+    % Signal with a slower response plus high-frequency component
+    response = exp(-((t - 2).^2) / (2 * 0.5^2)) ...
+        + 0.1 * sin(2*pi*4*t);
+
+    y = zeros(duration * sr, 1);
+    onsets = [10 20 30];
+
+    for iTrial = 1:numel(onsets)
+        startSample = round(onsets(iTrial) * sr) + 1;
+        stopSample = startSample + numel(response) - 1;
+
+        y(startSample:stopSample) = ...
+            y(startSample:stopSample) + response;
+    end
+
+    timing.names = {'condition_1'};
+    timing.onsets = {onsets};
+    timing.durations = {zeros(size(onsets))};
+
+    data = cell(1,1);
     data{1}.header = struct( ...
         'chantype', 'pupil', ...
         'sr', sr, ...
         'units', 'a.u.');
-
     data{1}.data = y(:);
 
     infos.duration = duration;
@@ -51,58 +210,186 @@ function testSingleConditionTrialAverage(testCase)
 
     datafile = [tempname '.mat'];
     save(datafile, 'data', 'infos');
-    
+    testCase.addTeardown(@() deleteIfExists(datafile));
+
+    modelfile = [tempname '.mat'];
+    testCase.addTeardown(@() deleteIfExists(modelfile));
+
+    model = struct();
+    model.modelfile = modelfile;
+    model.datafile = {datafile};
+    model.timing = {timing};
+    model.timeunits = 'seconds';
+    model.window = window;
+    model.modelspec = 'dilation';
+    model.modality = 'pupil';
+    model.channel = 'pupil';
+    model.norm = 0;
+    model.baseline = 0;
+    model.norm_max = 0;
 
+    model.filter = struct( ...
+        'lpfreq', 1, ...
+        'lporder', 1, ...
+        'hpfreq', 'none', ...
+        'hporder', 0, ...
+        'down', 'none', ...
+        'direction', 'bi');
 
+    options.overwrite = 1;
 
-    testCase.addTeardown(@() deleteIfExists(datafile));
+    [sts, tam] = pspm_tam(model, options);
 
-    %% Verify test file itself
-    [sts, loaded] = pspm_load_channel( datafile, 'pupil', 'pupil');
+    testCase.assertEqual(sts, 1);
 
-    testCase.verifyEqual(sts, 1);
-    testCase.verifyEqual(loaded.header.sr, sr);
-    testCase.verifyEqual(loaded.data(:), y(:));
+    % No downsampling -> still 50 samples
+    testCase.verifySize(tam.data.Y{1},   [50 1]);
+    testCase.verifySize(tam.data.std{1}, [50 1]);
+    testCase.verifySize(tam.data.sem{1}, [50 1]);
+    testCase.verifySize(tam.data.X{1},   [50 1]);
+
+    % Sampling rate must stay unchanged
+    testCase.verifyEqual(tam.data.sr{1}, 10);
+
+    % Filter flag
+    testCase.verifyTrue(tam.data.filtered);
+
+    % Time vector must NOT have been filtered
+    expected_t = (1:50)' / sr;
+
+    testCase.verifyEqual( tam.data.X{1}, expected_t, 'AbsTol', 1e-12);
+
+    % Filtering should actually change the signal
+    unfiltered = response - response(1);
 
+    testCase.verifyGreaterThan( norm(tam.data.Y{1} - unfiltered), 1e-6);
 
-    %% Configure TAM model
+end
+
+
+function testTwoSessions(this)
+
+    sr = 10;
+    duration = 50;
+    window = 5;
+
+    t = (0:1/sr:window-1/sr)';
+
+    response1 = exp(-((t - 2).^2) / (2 * 0.5^2));
+    response2 = 2 * response1;
+
+    onsets = [10 20 30];
+
+    % ----- Session 1 -----
+    y1 = zeros(duration * sr, 1);
+
+    for iTrial = 1:numel(onsets)
+        startSample = round(onsets(iTrial) * sr) + 1;
+        stopSample = startSample + numel(response1) - 1;
+
+        y1(startSample:stopSample) = ...
+            y1(startSample:stopSample) + response1;
+    end
+
+    % ----- Session 2 -----
+    y2 = zeros(duration * sr, 1);
+
+    for iTrial = 1:numel(onsets)
+        startSample = round(onsets(iTrial) * sr) + 1;
+        stopSample = startSample + numel(response2) - 1;
+
+        y2(startSample:stopSample) = ...
+            y2(startSample:stopSample) + response2;
+    end
+
+    % Timing for both sessions
+    timing1.names = {'condition_1'};
+    timing1.onsets = {onsets};
+    timing1.durations = {zeros(size(onsets))};
+
+    timing2 = timing1;
+
+    % Temporary files
+    datafile1 = [tempname '.mat'];
+    datafile2 = [tempname '.mat'];
     modelfile = [tempname '.mat'];
-    testCase.addTeardown(@() deleteIfExists(modelfile));
 
-    model = struct();
+    this.addTeardown(@() deleteIfExists(datafile1));
+    this.addTeardown(@() deleteIfExists(datafile2));
+    this.addTeardown(@() deleteIfExists(modelfile));
 
+    % Save session 1
+    data = cell(1,1);
+    data{1}.header = struct( ...
+        'chantype', 'pupil', ...
+        'sr', sr, ...
+        'units', 'a.u.');
+    data{1}.data = y1;
+
+    infos.duration = duration;
+    infos.durationinfo = 'Duration in seconds';
+    infos.source = struct();
+
+    save(datafile1, 'data', 'infos');
+
+    % Save session 2
+    data{1}.data = y2;
+    save(datafile2, 'data', 'infos');
+
+    % Model
+    model = struct();
     model.modelfile = modelfile;
-    model.datafile = {datafile};
-    model.timing = {timing};
+
+    model.datafile{1} = datafile1;
+    model.datafile{2} = datafile2;
+
+    model.timing{1} = timing1;
+    model.timing{2} = timing2;
 
     model.timeunits = 'seconds';
     model.window = window;
-
     model.modelspec = 'dilation';
     model.modality = 'pupil';
     model.channel = 'pupil';
 
-    % Keep preprocessing simple for this test
     model.norm = 0;
     model.baseline = 0;
     model.norm_max = 0;
-    % model.std_exp_cond = 'none';
 
-    options = struct();
+    model.filter = struct( ...
+        'lpfreq', 'none', ...
+        'lporder', 1, ...
+        'hpfreq', 'none', ...
+        'hporder', 1, ...
+        'down', 'none', ...
+        'direction', 'bi');
+
     options.overwrite = 1;
 
-    %% Run TAM
+    % Run TAM
     [sts, tam] = pspm_tam(model, options);
 
-    %% Basic output checks
-    testCase.verifyTrue(isstruct(tam));
-    testCase.verifyTrue(isfield(tam, 'data'));
-    testCase.verifyTrue(isfield(tam.data, 'Y'));
-    testCase.verifyEqual(numel(tam.data.Y), 1);
-end
+    this.verifyEqual(sts, 1);
+
+    % Session average:
+    % (response1 + response2) / 2 = 1.5 * response1
+    expected = (response1 + response2) / 2;
+
+    % TAM baseline correction at first datapoint
+    expected = expected - expected(1);
+
+    this.verifyEqual( tam.data.Y{1}, expected(:), 'AbsTol', 1e-10);
+
+    % Both sessions keep 10 Hz
+    this.verifyEqual(tam.data.sr, {10, 10});
+
+    % Still one experimental condition
+    this.verifyEqual(numel(tam.data.Y), 1);
 
 end
 
+
+end
 end