Update Python integration to PsPM UI

src/pspm_cfg/pspm_cfg_pp_heart_data.m

@@ -1,8 +1,7 @@
 
 function pp_heart_data = pspm_cfg_pp_heart_data
 
-% $Id: pspm_cfg_pp_heart_data.m 784 2019-07-08 08:16:46Z esrefo $
-% $Rev: 784 $
+% Updated 27-Mar-2024 by Teddy
 
 % Initialise
 global settings
@@ -192,12 +191,50 @@
     'Heart beat conversion, or directly work on heart beat time stamps, ', ...
     'for example obtained by a pulse oxymeter.']};
 
+%% ppg2hb
+ppg2hb_py_auto          = cfg_const;
+ppg2hb_py_auto.name     = 'Automatically detect Python';
+ppg2hb_py_auto.tag      = 'ppg2hb_py_auto';
+ppg2hb_py_auto.val      = {0};
+ppg2hb_py_auto.help     = {['This only works if a Python environment ',...
+                            'already exists in Matlab, created by ',...
+                            'previous PsPM function calls or manually.']};
+
+ppg2hb_py_path          = cfg_files;
+ppg2hb_py_path.name     = 'Manually define Python';
+ppg2hb_py_path.tag      = 'ppg2hb_py_path';
+ppg2hb_py_path.num      = [1 1];
+ppg2hb_py_path.help     = {'Please specify python executable file on the computer.'};
+
+ppg2hb_py_detect        = cfg_choice;
+ppg2hb_py_detect.name   = 'HeartPy';
+ppg2hb_py_detect.tag    = 'heart_py';
+ppg2hb_py_detect.val    = {ppg2hb_py_auto};
+ppg2hb_py_detect.values = {ppg2hb_py_auto, ppg2hb_py_path};
+ppg2hb_py_detect.help   = {'Mode of detecting python path in the operating system.'};
+
+ppg2hb_classic          = cfg_const;
+ppg2hb_classic.name     = 'Classic';
+ppg2hb_classic.tag      = 'ppg2hb_classic';
+ppg2hb_classic.val      = {0};
+ppg2hb_classic.help     = {'Classic mode.'};
+
+ppg2hb_method           = cfg_choice;
+ppg2hb_method.name      = 'Select the method of converting the data';
+ppg2hb_method.tag       = 'ppg2hb_convert';
+ppg2hb_method.val       = {ppg2hb_classic};
+ppg2hb_method.values    = {ppg2hb_classic, ppg2hb_py_detect};
+ppg2hb_method.help      = {['Convert the PPG data into heart rate by using the ', ...
+                          'selected method.']};
+
 ppg2hb_chan         = pspm_cfg_channel_selector('peripheral pulse oxymetry');
+ppg2hb_chan.help    = {['Number of peripheral pulse oximetry channel ', ...
+     '(default: last peripheral puls oximetry channel)']};
 
 ppg2hb              = cfg_exbranch;
-ppg2hb.name         = 'Convert Peripheral pulse oximetry to Heart Beat';
+ppg2hb.name         = 'Convert peripheral pulse oximetry to Heart Beat';
 ppg2hb.tag          = 'ppg2hb';
-ppg2hb.val          = {ppg2hb_chan};
+ppg2hb.val          = {ppg2hb_chan, ppg2hb_method};
 ppg2hb.help          = {['Convert Peripheral pulse oximetry to ', ...
     'Heart Beat events.']};
 

src/pspm_cfg/pspm_cfg_run_pp_heart_data.m

@@ -1,5 +1,5 @@
 function out = pspm_cfg_run_pp_heart_data(job)
-% Updated on 08-01-2024 by Teddy
+% Updated on 26-03-2024 by Teddy
 fn = job.datafile{1};
 outputs = cell(size(job.pp_type));
 for i = 1:numel(job.pp_type)
@@ -64,6 +64,14 @@
         end
       case 'ppg2hb'
         options = struct();
+        if ~isfield(job.pp_type{i}.ppg2hb.ppg2hb_convert, 'heart_py')
+          options.method = 'classic';
+        else
+          options.method = 'heartpy';
+          if isfield(job.pp_type{i}.ppg2hb.ppg2hb_convert.heart_py, 'ppg2hb_py_path')
+            options.python_path = job.pp_type{i}.ppg2hb.ppg2hb_convert.heart_py.ppg2hb_py_path{1};
+          end
+        end
         options.channel = chan;
         options = pspm_update_struct(options, job, {'channel_action'});
         [sts, winfo] = pspm_convert_ppg2hb(fn, options);
@@ -76,4 +84,4 @@
     end
   end
 end
-out = {fn};
+out = {fn};

src/pspm_convert_ppg2hb.m

@@ -10,17 +10,19 @@
 %   [ sts, outinfo ] = pspm_convert_ppg2hb( fn, options )
 % ● Arguments
 %                 fn: file name with path
-%   ┌─────── options
-%   ├───────.channel: [optional, numeric/string, default: 'ppg', i.e. last 
+%            channel: ppg channel number, default: last ppg channel
+%   ┌────────options: struct with following possible fields
+%   ├────────.method: 'classic' (default) or 'heartpy'.
+%   ├───────.channel: [optional, numeric/string, default: 'ppg', i.e. last
 %   │                 PPG channel in the file]
 %   │                 Channel type or channel ID to be preprocessed.
-%   │                 Channel can be specified by its index (numeric) in the 
+%   │                 Channel can be specified by its index (numeric) in the
 %   │                 file, or by channel type (string).
 %   │                 If there are multiple channels with this type, only
 %   │                 the last one will be processed. If you want to
 %   │                 process several PPG channels in a PsPM file,
 %   │                 call this function multiple times with the index of
-%   │                 each channel.  In this case, set the option 
+%   │                 each channel.  In this case, set the option
 %   │                 'channel_action' to 'add',  to store each
 %   │                 resulting 'hb' channel separately.
 %   ├───.diagnostics: [true/FALSE]
@@ -29,15 +31,23 @@
 %   │                 Defines whether the interpolated
 %   │                 data should be added or the corresponding channel
 %   │                 should be replaced.
-%   └───────────.lsm: [integer]
-%                     large spikes mode compensates for large spikes
-%                     while generating template by removing the [integer]
-%                     largest percentile of spikes from consideration.
+%   ├───────.missing: allows to specify missing (e. g. artefact) epochs in the
+%   │                 data file. See pspm_get_timing for epoch definition. This
+%   │                 must always be specified in SECONDS. These epochs will be
+%   │                 set to 0 for the detection.
+%   │                 Default: no missing values
+%   ├───────────.lsm: [integer] for method 'classic'
+%   │                 large spikes mode compensates for large spikes
+%   │                 while generating template by removing the [integer]
+%   │                 largest percentile of spikes from consideration.
+%   └───.python_path: [char] for method 'heartpy'
+%                     The path where python can be found. Mandatory if
+%                     python environment is not yet set up
 % ● History
 %   Introduced in PsPM 3.1
 %   Written in 2016 by Samuel Gerster (University of Zurich)
 %                      Tobias Moser (University of Zurich)
-%   Maintained in 2022 by Teddy Chao (UCL)
+%   Updated in 2024 by Dominik Bach/Uzay Gokay (Uni Bonn)
 
 %% Initialise
 global settings
@@ -51,9 +61,13 @@
 % -------------------------------------------------------------------------
 if nargin < 1
   warning('ID:invalid_input', 'No input. Don''t know what to do.'); return;
-elseif nargin < 2 
-    options = struct();
+elseif ~ischar(fn)
+  warning('ID:invalid_input', 'Need file name string as first input.'); return;
+elseif nargin < 2
+  options = struct();
+  options.channel = 'ppg';
 end
+
 options = pspm_options(options, 'convert_ppg2hb');
 if options.invalid
   return
@@ -71,96 +85,152 @@
 %% Large spikes mode
 %--------------------------------------------------------------------------
 ppg = data.data;
-% large spike mode
-if options.lsm
-  fprintf('Entering large spikes mode. This might take some time.');
-  % Look for all peaks lower than 200 bpm (multiple of two in heart rate
-  %  to compensate for absolute value and therefore twice as mani maxima)
-  [pks,pis] = findpeaks(abs(ppg),...
-    'MinPeakDistance',30/200*data.header.sr);
-  % 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);
-  % define large spikes index as last lsm percentile (or as adapted above)
-  lsi = pks>quantile(pks,q);
-  %define a minimum peak prominence 2/3 of non large spikes range (more
-  %or less)
-  minProm = max(pks(~lsi))*2/3;
-  % save indexes of large spikes for later removal while generating
-  % template
-  lsi = pis(lsi);
-  fprintf('   done.\n');
-else
-  minProm = range(ppg)/3;
-end
+sr = data.header.sr;
 
-%% Create template
-%--------------------------------------------------------------------------
-fprintf('Creating template. This might take some time.');
-% Find prominent peaks for a max heart rate of 200 bpm
-[~,pis] = findpeaks(data.data,...
-  'MinPeakDistance',60/200*data.header.sr,...
-  'MinPeakProminence',minProm);
-
-if options.lsm
-  % Remove large spikes from
-  [~,lsi_in_pis,~] = intersect(pis,lsi);
-  pis(lsi_in_pis) = [];
+
+% process missing data
+nan_index = isnan(ppg);
+if ~isempty(nan_index)
+  ppg(nan_index) = 0;
 end
 
-% handle possible errors
-if isempty(pis),warning('ID:NoPulse', 'No pulse found, nothing done.');return;end
-if length(pis)==1,warning('ID:OnePulse', 'Only one pulse found, unable to calculate min_pulse_period.');return;end
-
-% get pulse period lower limit (assumed onset) as 30% of smalest period
-% before detected peaks
-min_pulse_period = min(diff(pis));
-period_index_lower_bound = floor(pis(2:end-1)-.3*min_pulse_period);
-fprintf('...');
-
-% Create template from mean of peak time-locked ppg pulse periods
-pulses = cell2mat(arrayfun(@(x) data.data(x:x+min_pulse_period),period_index_lower_bound','un',0));
-template = mean(pulses,2);
-fprintf('done.\n');
-
-% handle diagnostic plots relevant to template building
-if options.diagnostics
-  t_template = (0:length(template)-1)'/data.header.sr;
-  t_pulses = repmat(t_template,1,length(pis)-2);
-  figure
-  plot(t_pulses,pulses,'--')
-  set(gca,'NextPlot','add')
-  plot(t_template,template,'k','lineWidth',3)
-  xlabel('time [s]')
-  ylabel('Amplitude')
-  title('Generated ppg template (bk) and pulses used (colored)')
+if ~isempty(options.missing)
+  [sts, missing] = pspm_get_timing('epochs', options.missing, 'seconds');
+  if sts < 1, return; end
+  index   = pspm_epochs2logical(missing, numel(ppg), sr);
+  if (sum(index) > 0)
+    ppg(find(index)) = 0; % sometimes the logical indexing does not work even though index contains only 0 and 1 and is of correct size
+  end
 end
 
-%% Cross correlate the signal with the template and find peaks
-%--------------------------------------------------------------------------
-fprintf('Applying template.');
-ppg_corr = xcorr(data.data,template)/sum(template);
-% Truncate ppg_xcorr and realigne it so the max correlation corresponds to
-% templates peak and not beginning of template.
-ppg_corr = ppg_corr(length(data.data)-floor(.3*min_pulse_period):end-floor(.3*min_pulse_period));
-if options.diagnostics
-  t_ppg = (0:length(data.data)-1)'/data.header.sr;
-  figure
-  if length(t_ppg) ~= length(ppg_corr)
-    length(t_ppg)
+% -------------------------------------------------------------------------
+%% Heartpy
+if strcmpi(options.method, 'heartpy')
+  % initialise python
+  if isempty(options.python_path)
+    psts = pspm_check_python;
+  else
+    psts = pspm_check_python(options.python_path);
+  end
+  if psts < 1, return; end
+  psts = pspm_check_python_modules('heartpy');
+  if psts < 1, return; end
+
+  filtered_ppg = py.heartpy.filter_signal(ppg, ...
+    pyargs('cutoff', [1,20], ...
+    'filtertype',  'bandpass', ...
+    'sample_rate', sr, ...
+    'order', 3));
+  filtered_ppg = double(py.array.array('d',(filtered_ppg)));
+  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);
+    msg = sprintf(['Heart beat detection from PPG with cross correlation ',...
+      'HB-timeseries added to data on %s'],...
+      date);
+    hb = peak_list(:) / sr;
+  catch
+    msg = sprintf('HeartPy did not find any heart beats on %s', date);
+    hb = [];
+  end
+else
+  %% large spike mode
+  if options.lsm
+    fprintf('Entering large spikes mode. This might take some time.');
+    % Look for all peaks lower than 200 bpm (multiple of two in heart rate
+    %  to compensate for absolute value and therefore twice as mani maxima)
+    [pks,pis] = findpeaks(abs(ppg),...
+      'MinPeakDistance',30/200*sr);
+    % 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);
+    % define large spikes index as last lsm percentile (or as adapted above)
+    lsi = pks>quantile(pks,q);
+    %define a minimum peak prominence 2/3 of non large spikes range (more
+    %or less)
+    minProm = max(pks(~lsi))*2/3;
+    % save indexes of large spikes for later removal while generating
+    % template
+    lsi = pis(lsi);
+    fprintf('   done.\n');
+  else
+    minProm = range(ppg)/3;
   end
-  plot(t_ppg,ppg_corr,t_ppg,data.data)
-  xlabel('time [s]')
-  ylabel('Amplitude')
-  title('ppg cross-corelated with template and ppg')
-  legend('ppg (X) template','ppg')
+
+  %% Create template
+  %--------------------------------------------------------------------------
+  fprintf('Creating template. This might take some time.');
+  % Find prominent peaks for a max heart rate of 200 bpm
+  [~,pis] = findpeaks(ppg,...
+    'MinPeakDistance',60/200*sr,...
+    'MinPeakProminence',minProm);
+
+  if options.lsm
+    % Remove large spikes from
+    [~,lsi_in_pis,~] = intersect(pis,lsi);
+    pis(lsi_in_pis) = [];
+  end
+
+  % handle possible errors
+  if isempty(pis),warning('ID:NoPulse', 'No pulse found, nothing done.');return;end
+  if length(pis)==1,warning('ID:OnePulse', 'Only one pulse found, unable to calculate min_pulse_period.');return;end
+
+  % get pulse period lower limit (assumed onset) as 30% of smalest period
+  % before detected peaks
+  min_pulse_period = min(diff(pis));
+  period_index_lower_bound = floor(pis(2:end-1)-.3*min_pulse_period);
+  fprintf('...');
+
+  % Create template from mean of peak time-locked ppg pulse periods
+  pulses = cell2mat(arrayfun(@(x) ppg(x:x+min_pulse_period),period_index_lower_bound','un',0));
+  template = mean(pulses,2);
+  fprintf('done.\n');
+
+  % handle diagnostic plots relevant to template building
+  if options.diagnostics
+    t_template = (0:length(template)-1)'/sr;
+    t_pulses = repmat(t_template,1,length(pis)-2);
+    figure
+    plot(t_pulses,pulses,'--')
+    set(gca,'NextPlot','add')
+    plot(t_template,template,'k','lineWidth',3)
+    xlabel('time [s]')
+    ylabel('Amplitude')
+    title('Generated ppg template (bk) and pulses used (colored)')
+  end
+
+  %% Cross correlate the signal with the template and find peaks
+  %--------------------------------------------------------------------------
+  fprintf('Applying template.');
+  ppg_corr = xcorr(ppg,template)/sum(template);
+  % Truncate ppg_xcorr and realigne it so the max correlation corresponds to
+  % templates peak and not beginning of template.
+  ppg_corr = ppg_corr(length(ppg)-floor(.3*min_pulse_period):end-floor(.3*min_pulse_period));
+  if options.diagnostics
+    t_ppg = (0:length(ppg)-1)'/sr;
+    figure
+    if length(t_ppg) ~= length(ppg_corr)
+      length(t_ppg)
+    end
+    plot(t_ppg,ppg_corr,t_ppg,ppg)
+    xlabel('time [s]')
+    ylabel('Amplitude')
+    title('ppg cross-corelated with template and ppg')
+    legend('ppg (X) template','ppg')
+  end
+  % Get peaks that are at least one template width appart. These are the best
+  % correlation points.
+  [~,hb] = findpeaks(ppg_corr/max(ppg_corr),...
+    sr,...
+    'MinPeakdistance',min_pulse_period/sr);
+  fprintf('   done.\n');
 end
-% Get peaks that are at least one template width appart. These are the best
-% correlation points.
-[~,hb] = findpeaks(ppg_corr/max(ppg_corr),...
-  data.header.sr,...
-  'MinPeakdistance',min_pulse_period/data.header.sr);
-fprintf('   done.\n');
+
 
 %% Prepare output and save
 %--------------------------------------------------------------------------