Skip to content

Commit

Permalink
Implement documented and input checked version of pspm_ecg2hb_fmri
Browse files Browse the repository at this point in the history 
git-svn-id: https://svn.code.sf.net/p/pspm/svn/trunk@732 23ec9670-3a1e-4a12-9187-a0ba8bce340d
  • Loading branch information
@eozd
eozd committed Jun 20, 2019
1 parent 87473db commit 75e2cb8
Show file tree
Hide file tree
Showing 2 changed files with 137 additions and 8 deletions.
29 changes: 22 additions & 7 deletions amri_eegfmri/amri_eeg_rpeak.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,13 @@
% to
% [~,imax]=max(weights' .* ccorr(anarrowrange));
% so that a regular expectation is calculated.
% - Remove the second
% kTEO.k = 1;
% line right before constructing TEO.
% - Add ecg_lowcutoff and ecg_highcutoff.
% - Make all algorithm parameters modifiable using
% keyword arguments.
% - Remove number of input argument checks

%%
function r_peaks = amri_eeg_rpeak(ecg,varargin)
Expand All @@ -95,6 +102,9 @@
kTEO.lowcutoff = 8; % low cutoff
kTEO.highcutoff= 40; % high cutoff

ecg_lowcutoff = 0.5;
ecg_highcutoff = 40;

thres.mincc = 0.5; % cross corr to the template
thres.maxrpa = 3; % max relative peak amplitude
thres.minrpa = 0.4; % min relative peak amplitude
Expand All @@ -108,11 +118,6 @@
% Collect keyword-value pairs
% *************************************************************************

if (nargin> 2 && rem(nargin,2) == 1)
printf('amri_eeg_cbc(): Even number of input arguments???')
return
end

for i = 1:2:size(varargin,2) % for each Keyword
Keyword = varargin{i};
Value = varargin{i+1};
Expand All @@ -126,6 +131,17 @@
end
elseif strcmpi(Keyword,'whatisy')
whatisy=Value;
elseif strcmpi(Keyword, 'teoparams')
kTEO.k = Value(1);
kTEO.lowcutoff = Value(2);
kTEO.highcutoff = Value(3);
elseif strcmpi(Keyword, 'ecgcutoff')
ecg_lowcutoff = Value(1);
ecg_highcutoff = Value(2);
elseif strcmpi(Keyword, 'mincc')
thres.mincc = Value(1);
elseif strcmpi(Keyword, 'minrpa')
thres.minrpa = Value(1);
end
end

Expand All @@ -138,11 +154,10 @@

% bandpass filtering from 10 to 40 Hz, in order to suppress T-waves
if strcmpi(filter_method,'ifft')
ecg.data = amri_sig_filtfft(ecg.data,ecg.srate,0.5,40);
ecg.data = amri_sig_filtfft(ecg.data,ecg.srate, ecg_lowcutoff, ecg_highcutoff);
E = amri_sig_filtfft(ecg.data,ecg.srate,kTEO.lowcutoff,kTEO.highcutoff);
end

kTEO.k=1;
% construct a complex lead by k-Teager energy operator
if flag_verbose>0
fprintf(['amri_eeg_rpeak(): construct k-TEO complex (k=' ...
Expand Down
116 changes: 115 additions & 1 deletion pspm_ecg2hb_fmri.m
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,13 @@
function [sts, out_channel] = pspm_ecg2hb_fmri(fn, options)
% pspm_ecg2hb_fmri performs R-peak detection from an ECG signal using the steps
% decribed in R-peak detection section of [1]. This function uses a modified
% version of the amri_eeg_rpeak.m code that can be obtained from [2]. Modified
% version with a list of changes made is shipped with PsPM under amri_eegfmri
% directory.
%
% Once the R-peaks are computed, according to the option 'channel_action',
% it will either replace an existing heartbeat channel or add it as a new
% channel to the provided file.
%
% FORMAT: [sts, out_channel] = pspm_ecg2hb_fmri(fn)
% [sts, out_channel] = pspm_ecg2hb_fmri(fn, options)
Expand All @@ -21,6 +30,51 @@
% each channel. Further, use 'add' mode to store each
% resulting 'heartbeat' channel separately.
%
% signal_to_use: ['ecg'/'teo'/'auto'] Choose which signal will be used
% as the input to the core R-peak detection steps. When
% 'ecg', filtered ECG signal will be used. When 'teo',
% Teager Enery Operator will be applied to the filtered
% ECG signal before feeding it to R-peak finding part.
% When 'auto', the option that results in the higher
% maximal autocorrelation will be used.
% (Default: 'auto')
%
% heartbeat_rate: [numeric] Minimum and maximum heartbeat rates (BPM)
% to use in the algorithm. Must be a numeric array of
% length 2, i.e. [min_bpm max_bpm].
% (Default: [20 200])
% (Unit: beats per minute)
%
% ecg_bandpass: [numeric] Minimum and maximum frequencies to use
% during bandpass filter of the raw ECG signal to
% construct filtered ECG signal.
% (Default: [0.5 40])
% (Unit: Hz)
%
% teo_bandpass: [numeric] Minimum and maximum frequencies to use
% during bandpass filter of filtered ECG signal to
% construct TEO input signal.
% (Default: [8 40])
% (Unit: Hz)
%
% teo_order: [numeric] Order of the TEO operator. Must be integer.
% For a discrete time signal x(t) and order k,
% TEO[x(t); k] is defined as
%
% TEO[x(t); k] = x(t)x(t) - x(t-k)x(t+k)
%
% (Default: 1)
%
% min_cross_corr: [numeric] Minimum cross correlation between a candidate
% R-peak and the found template such that the candidate is
% classified as an R-peak.
% (Default: 0.5)
%
% min_relative_amplitude:
% [numeric] Minimum relative peak amplitude of a candidate
% R-peak such that it is classified as an R-peak.
% (Default: 0.4)
%
% channel_action: ['add'/'replace'] Defines whether corrected data
% should be added or the corresponding preprocessed
% channel should be replaced. Note that 'replace' mode
Expand Down Expand Up @@ -58,9 +112,62 @@
if ~isfield(options, 'channel_action')
options.channel_action = 'replace';
end
if ~isfield(options, 'signal_to_use')
options.signal_to_use = 'auto';
end
if ~isfield(options, 'heartbeat_rate')
options.heartbeat_rate = [20 200];
end
if ~isfield(options, 'ecg_bandpass')
options.ecg_bandpass = [0.5 40];
end
if ~isfield(options, 'teo_bandpass')
options.teo_bandpass = [8 40];
end
if ~isfield(options, 'teo_order')
options.teo_order = 1;
end
if ~isfield(options, 'min_cross_corr')
options.min_cross_corr = 0.5;
end
if ~isfield(options, 'min_relative_amplitude')
options.min_relative_amplitude = 0.4;
end

% input checks
% -------------------------------------------------------------------------
if ~ismember(options.channel_action, {'add', 'replace'})
warning('ID:invalid_input', 'Option channel_action must be either ''add'' or ''replace''');
return;
end
if ~ismember(options.signal_to_use, {'ecg', 'teo', 'auto'})
warning('ID:invalid_input', 'Option signal_to_use must be one of ''ecg'',''teo'' or ''auto''');
return;
end
if ~isnumeric(options.heartbeat_rate) || any(options.heartbeat_rate <= 0)
warning('ID:invalid_input', 'Option heartbeat_rate must contain positive numbers');
return;
end
if ~isnumeric(options.ecg_bandpass) || any(options.ecg_bandpass <= 0)
warning('ID:invalid_input', 'Option ecg_bandpass must contain positive numbers');
return;
end
if ~isnumeric(options.teo_bandpass) || any(options.teo_bandpass <= 0)
warning('ID:invalid_input', 'Option teo_bandpass must contain positive numbers');
return;
end
if ~isnumeric(options.teo_order) || options.teo_order <= 0 || mod(options.teo_order, 1) ~= 0
warning('ID:invalid_input', 'Option teo_order must be a positive integer');
return;
end
if ~isnumeric(options.min_cross_corr)
warning('ID:invalid_input', 'Option min_cross_corr must be numeric');
return;
end
if ~isnumeric(options.min_relative_amplitude)
warning('ID:invalid_input', 'Option min_relative_amplitude must be numeric');
return;
end

% load
% -------------------------------------------------------------------------
Expand All @@ -74,7 +181,14 @@
addpath(pspm_path('amri_eegfmri'));
ecg.data = data{1}.data;
ecg.srate = data{1}.header.sr;
heartbeats{1}.data = amri_eeg_rpeak(ecg);
heartbeats{1}.data = amri_eeg_rpeak(ecg, ...
'WhatIsY', options.signal_to_use, ...
'PulseRate', options.heartbeat_rate, ...
'TEOParams', [options.teo_order options.teo_bandpass], ...
'ECGcutoff', options.ecg_bandpass, ...
'mincc', options.min_cross_corr, ...
'minrpa', options.min_relative_amplitude ...
);
rmpath(pspm_path('amri_eegfmri'));

% save
Expand Down

0 comments on commit 75e2cb8

Please sign in to comment.