Fix TAM z-score normalization

Author: 4gwe

src/pspm_tam.m

@@ -155,7 +155,7 @@
     % Filling up the data and the sampling rates
     y{iFile} = data.data(:);
     sr(iFile) = data.header.sr;
-    
+
     fprintf('.');
 end
 
@@ -174,15 +174,15 @@
   fprintf('\n');
 end
 
-%%  Zscoring the data ->  can be done by pspm_extract_segments ! what is with nans?
-if model.norm
-  fprintf('Zscoring ...\n')
-  n_file = numel(model.datafile);
-  for iFile = 1:n_file
-    % NANZSCORE found in src/ext/VBA/stats&plots
-    [y{iFile},~,~] = nannorm(y{iFile}); %nannorm??? -> nanzscore??
-  end
-end
+%%  Zscoring the data ->  will be done by pspm_extract_segments 
+% if model.norm
+%   fprintf('Zscoring ...\n')
+%   n_file = numel(model.datafile);
+%   for iFile = 1:n_file
+%     % NANZSCORE found in src/ext/VBA/stats&plots
+%     [y{iFile},~,~] = nanzscore(y{iFile}); %nannorm??? -> nanzscore??
+%   end
+% end
 
 %%  Extracting segments
 fprintf('Extracting segments ...\n')
@@ -191,6 +191,7 @@
 extrsgopt.timeunits = model.timeunits;
 extrsgopt.length = model.window;       % segments of 'model.window' time unit long
 extrsgopt.plot = 0;                    % do not plot mean value and std
+extrsgopt.norm = model.norm;
 
 for k=1:n_file
   if strcmpi(model.timeunits, 'markers')
@@ -207,7 +208,6 @@
 
       % In file mode the raw pupil data are loaded again,
       % so normalization must happen inside extract_segments.
-      fileopt.norm = model.norm;
       [lsts, s] = pspm_extract_segments( 'file', model.datafile{k}, model.channel, model.timing{k}, fileopt);
   else
       [lsts, s] = pspm_extract_segments('data', y{k}, sr(k), model.timing{k}, extrsgopt); % wird es richtig gen znormed?

test/pspm_tam_test.m

@@ -675,7 +675,102 @@ function testMarkerTimeunits(testCase, nFiles)
 
 end
 
+function testZscoreNormalization(testCase)
 
+    sr = 10;
+    duration = 50;
+    window = 5;
+
+    t = (0:1/sr:window-1/sr)';
+    response = exp(-((t - 2).^2) / (2 * 0.5^2));
+
+    onsets = [10 20 30];
+
+    y = zeros(duration * sr, 1);
+
+    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;
+    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';
+
+    % This is what we are testing
+    model.norm = 1;
+
+    model.baseline = 0;
+    model.norm_max = 0;
+
+    model.filter = struct( ...
+        'lpfreq', 'none', ...
+        'lporder', 1, ...
+        'hpfreq', 'none', ...
+        'hporder', 1, ...
+        'down', 'none', ...
+        'direction', 'bi');
+
+    options.overwrite = 1;
+
+    [sts, tam] = pspm_tam(model, options);
+
+    testCase.verifyEqual(sts, 1);
+
+    % Expected z-scoring of the COMPLETE continuous time series
+    mu = mean(y, 'omitnan');
+    sigma = std(y, 0, 'omitnan');
+
+    expected = (response - mu) / sigma;
+
+    % TAM subsequently baselines against the first point
+    expected = expected - expected(1);
+
+    testCase.verifyEqual( ...
+        tam.data.Y{1}, ...
+        expected(:), ...
+        'AbsTol', 1e-10);
+
+    testCase.verifyEqual(tam.data.norm, 1);
+
+end