Fix pspm_convert_ppg2hb HeartPy handling and extend test coverage

.github/workflows/main.yml

@@ -20,8 +20,12 @@ jobs:
         with:
           python-version: '3.10' 
 
-      - name: Install bioread
-        run: python -m pip install bioread==3.0.1
+      - name: Install Python dependencies
+        run: python -m pip install bioread==3.0.1 heartpy==1.2.7
+
+      - name: Verify Python dependencies
+        run: python -c "import bioread, heartpy; print('Python dependencies available')"
+
 
       - name: Detect python
         run: python 'test/pspm_py_find_location.py'

src/pspm_convert_ppg2hb.m

@@ -53,6 +53,7 @@
 %   Written in 2016 by Samuel Gerster (University of Zurich)
 %                      Tobias Moser (University of Zurich)
 %   Updated in 2024 by Dominik Bach/Uzay Gokay (Uni Bonn)
+%   Updated in 2026 by Bernhard von Raußendorf (Uni Bonn)
 
 %% Initialise
 global settings
@@ -118,7 +119,7 @@
   psts = pspm_check_python_modules('heartpy');
   if psts < 1, return; end
 
-  filtered_ppg = py.heartpy.filter_signal(ppg, ...
+  filtered_ppg = py.heartpy.filter_signal(ppg.', ...
     pyargs('cutoff', [1,20], ...
     'filtertype',  'bandpass', ...
     'sample_rate', sr, ...
@@ -127,18 +128,20 @@
   try
     tup = py.heartpy.process(filtered_ppg, pyargs('sample_rate', sr));
     wd = tup{1};
-    m = tup{2};
+
     py_peak_list =  py.array.array('d',(wd{'peaklist'}));
-    py_removed =  py.array.array('d',(wd{'removed_beats'}));
-    peak_list = double(py_peak_list) ;
-    rejected_peaks = double(py_removed);
+    binary_peak_list = double(py.array.array('d', wd{'binary_peaklist'}));
+
+    peak_list = double(py_peak_list);
+    peak_list = peak_list(logical(binary_peak_list));
+
     msg = sprintf(['Heart beat detection from PPG with HeartPy ',...
       'HB-timeseries added to data on %s'],...
       datetime("today"));
     hb = peak_list(:) / sr;
-  catch
-    msg = sprintf('HeartPy did not find any heart beats on %s', datetime("today"));
-    hb = [];
+catch ME
+    warning('ID:heartpy_error', 'HeartPy processing failed: %s', ME.message);
+    return;
   end
 else
   %% large spike mode
@@ -148,6 +151,7 @@
     %  to compensate for absolute value and therefore twice as mani maxima)
     [pks,pis] = findpeaks(abs(ppg),...
       'MinPeakDistance',30/200*sr);
+    if isempty(pks),warning('ID:NoPulse', 'No pulse found, nothing done.');return;end % to protect from 0/0 -> NaN
     % Ensure at least one spike is removed by adapting quantil to realistic
     % values, given number of detected spikes
     q = floor(length(pks)*(1-options.lsm/100))/length(pks);
@@ -230,14 +234,15 @@
   [~,hb] = findpeaks(ppg_corr/max(ppg_corr),...
     sr,...
     'MinPeakdistance',min_pulse_period/sr);
+  msg = sprintf('Heart beat detection from ppg with cross correlation HB-timeseries added to data on %s', datetime("today"));
+
 end
 
 
 %% Prepare output and save
 %--------------------------------------------------------------------------
 % save data
 fprintf('Saving data.');
-msg = sprintf('Heart beat detection from ppg with cross correlation HB-timeseries added to data on %s', datetime("today"));
 
 newdata.data = hb(:);
 newdata.header.sr = 1;

src/pspm_options.m

@@ -108,7 +108,8 @@
     options = autofill_channel_action(options);
     options = autofill(options, 'channel',                'ppg',     '*Int*Char'        );
     options = autofill(options, 'diagnostics',            0,          1                 );
-    options = autofill(options, 'lsm',                    0,          [0,100-10^-10]    );
+    options = autofill(options, 'lsm',                    0,         'integer>=', 0     );
+    options = autofill(options, 'lsm',                    0,          '<', 100          );
     options = autofill(options, 'method',                 'classic',  {'classic', 'heartpy'});
     options = autofill(options, 'python_path',            '',         '*Char'           );
     options = autofill(options, 'missing',                {},         '*Char*Num'       );

test/pspm_convert_hb2hp_test.m

@@ -49,19 +49,19 @@ function invalid_input(this)
   this.verifyWarning(@() pspm_convert_hb2hp('abc'), 'ID:invalid_input');
   % Verify not a string sample rate
   this.verifyWarning(@() pspm_convert_hb2hp('abc','abc'), 'ID:invalid_input');
-  % Verify not a numeric channel
+  % Verify options is not a struct
   this.verifyWarning(@() pspm_convert_hb2hp('abc',2,'abc'), 'ID:invalid_input');
 
-  % 
-  % % Verify that call of pspm_load_data fails
-  % this.verifyWarning(@() pspm_convert_hb2hp(files{1},100), 'ID:nonexistent_file');
-  % % Verify that interpolation does not have enough points
-  % % this.verifyWarning(@() pspm_convert_hb2hp(files{2}, 100), 'ID:too_strict_limits');
-  % % Verify that call of pspm_write_channel fails
-  % options.channel_action = 'abc';
-  % this.verifyWarning(@() pspm_convert_hb2hp(files{3},100,[],options),'ID:invalid_input');
-  % %options.channel_action = 'add';
-  % %this.verifyWarningFree(@()pspm_convert_hb2hp(files{1},100,[],options));
+  fn = this.input_filename;
+
+  options.channel_action = 'abc';
+  this.verifyWarning( @() pspm_convert_hb2hp(fn, 100, options), 'ID:invalid_input');
+
+  options = struct('limit_lower', 2, 'limit_upper', 1);
+  this.verifyWarning( @() pspm_convert_hb2hp(fn, 100, options), 'ID:invalid_input');
+
+  options = struct('limit_lower', -1);
+  this.verifyWarning( @() pspm_convert_hb2hp(fn, 100, options), 'ID:invalid_input');
 
 end
 
@@ -77,62 +77,196 @@ function basic_conversion(this)
     this.verifyEqual(data{outchannel}.header.chantype, 'hp');
 end
 
+function selects_hb_channel_by_index(this)
+    fn = [tempname, '.mat'];
+
+    channels = { ...
+        struct('chantype', 'hb'), ...
+        struct('chantype', 'hb') ...
+        };
+
+    testdata = pspm_testdata_gen(channels, 5);
+
+    % Channel 1 -> HP = 1000 ms
+    testdata.data{1}.data = ...
+        [0.5; 1.5; 2.5; 3.5; 4.5];
+
+    % Channel 2 -> HP = 500 ms
+    testdata.data{2}.data = ...
+        [0.5; 1.0; 1.5; 2.0; 2.5; 3.0; 3.5; 4.0; 4.5];
+
+    save(fn, '-struct', 'testdata');
+    cleanup = onCleanup(@() delete(fn));
+
+    options.channel = 1;
+    options.channel_action = 'add';
+
+    [sts, outchannel] = pspm_convert_hb2hp(fn, 2, options);
+    this.verifyEqual(sts, 1);
+
+    [sts, ~, data] = pspm_load_data(fn);
+    this.verifyEqual(sts, 1);
+
+    hp = data{outchannel};
+
+    % Explicitly selecting channel 1 must produce 1000 ms.
+    this.verifyEqual(hp.data, 1000 * ones(size(hp.data)), 'AbsTol', 1e-10);
+end
+
 function too_strict_limits(this)
     fn = this.input_filename;
     sr = 1;
-    options = struct('limit_lower', 11, 'limit_upper', 11);
+    options = struct('limit_lower', 10, '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);
     this.verifyEqual(sts, 1);
     this.verifyEqual(data{outchannel}.header.chantype, 'hp');
-    this.verifyTrue(any(isnan(data{outchannel}.data)))
+    this.verifyTrue(all(isnan(data{outchannel}.data)));
+end
+
+function too_strict_limits_warning(this)
+    fn = this.input_filename;
+    sr = 1;
+
+    options = struct('limit_lower', 10, 'limit_upper', 11);
+    this.verifyWarning(@() pspm_convert_hb2hp(fn, sr, options), 'ID:too_strict_limits');
 end
 
 function add_replace_channel_action(this)
     fn = this.input_filename;
     sr = 1;
     options.channel_action = 'add';
 
-    [sts, infos, ~, filestruct] = pspm_load_data(fn);
+    [sts, ~, ~, filestruct] = pspm_load_data(fn);
     this.verifyEqual(sts, 1);
-    % display(infos.history)
 
-    % add 1st hp 
-    [sts, outchannel] = pspm_convert_hb2hp(fn, sr, options);
+    original_num_channels = filestruct.numofchan;
+
+    % add 1st hp
+    [sts, hp1] = pspm_convert_hb2hp(fn, sr, options);
     this.verifyEqual(sts, 1);
-    [sts, ~, data, ~ ] = pspm_load_data(fn);
+
+    [sts, ~, data, ~] = pspm_load_data(fn);
     this.verifyEqual(sts, 1);
-    this.verifyEqual(data{outchannel}.header.chantype, 'hp');
-    this.verifyEqual(filestruct.numofchan + 1 , numel(data));
-    % display(infos.history)
+    this.verifyEqual(data{hp1}.header.chantype, 'hp');
+    this.verifyEqual(original_num_channels + 1, numel(data));
 
     % add 2nd hp
-    [sts, outchannel] = pspm_convert_hb2hp(fn, sr, options);
+    [sts, hp2] = pspm_convert_hb2hp(fn, sr, options);
     this.verifyEqual(sts, 1);
-    [sts, ~, data, ~ ] = pspm_load_data(fn);
 
+    [sts, ~, data, ~] = pspm_load_data(fn);
     this.verifyEqual(sts, 1);
-    this.verifyEqual(data{outchannel}.header.chantype, 'hp');
-    this.verifyEqual(filestruct.numofchan + 2 , numel(data));
+    this.verifyEqual(data{hp2}.header.chantype, 'hp');
+    this.verifyEqual(original_num_channels + 2, numel(data));
 
+    % Make sure the two added channels are different channels
+    this.verifyNotEqual(hp1, hp2);
 
-    % replace hp (last)
-    options.channel_action = 'replace';    
-    [sts, outchannel] = pspm_convert_hb2hp(fn, sr, options);
+    % replace last hp
+    options.channel_action = 'replace';
+
+    [sts, hp_replaced] = pspm_convert_hb2hp(fn, sr, options);
+    this.verifyEqual(sts, 1);
+
+    [sts, ~, data, ~] = pspm_load_data(fn);
     this.verifyEqual(sts, 1);
-    [sts, infos, data, ~] = pspm_load_data(fn);
+
+    % Number of channels must stay the same
+    this.verifyEqual(original_num_channels + 2, numel(data));
+
+    % The last HP channel should have been replaced
+    this.verifyEqual(hp_replaced, hp2);
+
+    % Replaced channel is still HP
+    this.verifyEqual(data{hp_replaced}.header.chantype, 'hp');
+end
+
+function known_heart_period_values(this)
+    fn = [tempname, '.mat'];
+
+    duration = 5;
+    sr = 2;
+
+    channels = {struct('chantype', 'hb')};
+    testdata = pspm_testdata_gen(channels, duration);
+    testdata.data{1}.data = [0.5; 1.5; 3.0; 4.0];
+
+    save(fn, '-struct', 'testdata');
+    cleanup = onCleanup(@() delete(fn));
 
-    % infos.history
+    % Convert HB -> HP.
+    [sts, outchannel] = pspm_convert_hb2hp(fn, sr);
     this.verifyEqual(sts, 1);
-    this.verifyEqual(data{outchannel}.header.chantype, 'hp');
-    this.verifyEqual(filestruct.numofchan + 2 , numel(data));
+
+    [sts, ~, data, ~] = pspm_load_data(fn);
+    this.verifyEqual(sts, 1);
+
+    hp = data{outchannel};
+
+    % Expected interpolated HP signal.
+    expected_hp = [ ...
+        1000;
+        1000;
+        1000;
+        1166.666666666667;
+        1333.333333333333;
+        1500;
+        1250;
+        1000;
+        1000;
+        1000
+        ];
+
+    % Verify actual numerical conversion.
+    this.verifyEqual(hp.data, expected_hp, 'AbsTol', 1e-10);
+
+    % Verify metadata.
+    this.verifyEqual(hp.header.chantype, 'hp');
+    this.verifyEqual(hp.header.units, 'ms');
+    this.verifyEqual(hp.header.sr, sr);
+end
+
+function selects_last_hb_channel_by_default(this)
+    fn = [tempname, '.mat'];
+
+    channels = { ...
+        struct('chantype', 'hb'), ...
+        struct('chantype', 'hb') ...
+        };
+
+    testdata = pspm_testdata_gen(channels, 5);
+
+    % Channel 1: IBI = 1 s -> HP = 1000 ms
+    testdata.data{1}.data = [0.5; 1.5; 2.5; 3.5; 4.5];
+
+    % Channel 2: IBI = 0.5 s -> HP = 500 ms
+    testdata.data{2}.data = ...
+        [0.5; 1.0; 1.5; 2.0; 2.5; 3.0; 3.5; 4.0; 4.5];
+
+    save(fn, '-struct', 'testdata');
+    cleanup = onCleanup(@() delete(fn));
+
+    sr = 2;
+    options.channel_action = 'add';
+
+    [sts, outchannel] = pspm_convert_hb2hp(fn, sr, options);
+
+    this.verifyEqual(sts, 1);
+
+    [sts, ~, data] = pspm_load_data(fn);
+    this.verifyEqual(sts, 1);
+
+    hp = data{outchannel};
+
+    this.verifyEqual(hp.header.chantype, 'hp');
+    this.verifyEqual(hp.header.units, 'ms');
+    this.verifyEqual(hp.header.sr, sr);
+
+    % Default 'hb' must select the LAST HB channel.
+    this.verifyEqual(hp.data, 500 * ones(size(hp.data)), 'AbsTol', 1e-10);
 
 end
 end