Pre-processing Data (PRV Signals)

I am trying to recreate a published paper (link) and there are data that needs to be pre-process before calculating the stress rate recognition using physiological data.

My code is as follows:

%% BVP DATA
dataBVP = readmatrix(bvp_directory);
bvp_signal = dataBVP(:, 1);

% Define the original sampling frequency
fs_original = 64; % in Hz

% Define the target sampling frequency
fs_target = 128; % in Hz

% Resample the signal using pchip interpolation
bvptime_original = (0:length(bvp_signal) - 1) / fs_original;
bvptime_target = (0:1/fs_target:bvptime_original(end))';
% bvp_resampled = resample(bvp_signal,128,64);
bvp_resampled = pchip(bvptime_original, bvp_signal, bvptime_target);   

% Find the peaks in the resampled signal
[BVPpeak, BVPpeak_locs] = findpeaks(bvp_resampled, MinPeakDistance = fs_target/2);

% Compute the Pulse-to-Pulse (PP) intervals in seconds
bvp_pp_intervals = diff(BVPpeak_locs)/fs_target;


%% RPPG DATA
rppg_signal = load(rppg_directory);

% Define the original sampling frequency
fs_original = 35; % in Hz

% Define the target sampling frequency
fs_target = 128; % in Hz

% Resample the signal using pchip interpolation
rppgtime_original = (0:length(rppg_signal)-1)/fs_original;
rppgtime_target = (0:1/fs_target:rppgtime_original(end));
% rppg_resampled = resample(rppg_signal,128,35);
rppg_resampled = pchip(rppgtime_original, rppg_signal, rppgtime_target);
% rppg_resampled = interp1(rppgtime_original, rppg_signal, rppgtime_target, 'pchip');

% Find the peaks in the resampled signal
[RPPGpeak, RPPGpeak_locs] = findpeaks(rppg_resampled, MinPeakDistance = fs_target/2);

RPPGpeak_locs128 = RPPGpeak_locs/128;

% Compute the Pulse-to-Pulse (PP) intervals in seconds
rppg_pp_intervals = diff(RPPGpeak_locs)/fs_target;

I am not sure on whether my code is right especially on peak detection and pulse-to-pulse intervals section. Please advised.

Thank you.

Answer 1

[William Rose]

@Jeffery,

It looks to me like your code will do what you want, almost. However, I think your specification of MinPeakDistance=fs_target/2 is misguided, and will not yield the results you want.

I agree with you that Sabour et al. (2022) used the pchip method for interpolation.

Please provide a sample of the blood volume pulse data and remote photoplethysmography data. If you are using data from the UBFC-Phys database, please specifiy which data file (S1 to S56, etc.).

Good luck with your analysis.

[Jeffery]

function [bvp_pp_intervals, bvp_resampled] = F_resampleBVP(directory,show)
    % Load the signal from the CSV file
    dataBVP = readmatrix(directory);
    bvp_signal = dataBVP(:, 1);

    % Define the original sampling frequency
    fs_original = 64; % in Hz

    % Define the target sampling frequency
    fs_target = 128; % in Hz

    % Resample the signal using pchip interpolation
    bvptime_original = (0:length(bvp_signal) - 1) / fs_original;
    bvptime_target = (0:1/fs_target:bvptime_original(end))';
    % bvp_resampled = resample(bvp_signal,128,64);
    bvp_resampled = pchip(bvptime_original, bvp_signal, bvptime_target);   

    % Find the peaks in the resampled signal
    [BVPpeak, BVPpeak_locs] = findpeaks(bvp_resampled, MinPeakDistance = fs_target/2);

    % Compute the Pulse-to-Pulse (PP) intervals in seconds
    bvp_pp_intervals = diff(BVPpeak_locs)/fs_target;

    % x = diff(BVPpeak_locs128);

    % Input PRV signal
    bvpPRV = bvp_pp_intervals; % Replace with your actual PRV signal
    
    % Parameters
    N = length(bvpPRV); % PRV signal length
    L = 51; % Median filter length
    T = 0.15; % Threshold value
    
    % Median filtering
    % medPP = medfilt1(PRV, L); % Median filter with 51 samples
    
    % Mirror flip at boundaries
    boundary = ((L-1)/2) + 1;
    bvpPP_flip_left = flip(bvpPRV(2:boundary));
    bvpPP_flip_right = flip(bvpPRV(N-((L-1)/2):N-1));
    
    % Tranpose
    % bvpPRV = bvpPRV';
    % bvpPP_flip_right = bvpPP_flip_right';
    % bvpPP_flip_left = bvpPP_flip_left';
    
    bvpPPnew = [bvpPP_flip_left', bvpPRV', bvpPP_flip_right'];
    bvpmedPP = medfilt1(bvpPPnew, L);
    
    assignin('base', 'bvpPPnew', bvpPPnew);

    % Threshold-based filtering
    M = N + L - 1; % Length of PPnew
    threshold = T * max(abs(bvpmedPP)); % Threshold value
    
    % Filtering
    % for i = boundary:M-((L-1)/2)
    %         disp(i);
    %         disp(bvpPPnew(i));
    %         disp(bvpmedPP(i));
    % 
    %     if abs(bvpPPnew(i) - bvpmedPP(i)) > threshold
    %         disp(bvpPPnew(i));
    %         disp(bvpmedPP(i));
    %         bvpPPnew(i) = bvpmedPP(i);
    %     end
    % end
    
    for i = boundary:M-((L-1)/2)
        if abs(bvpPPnew(i) - bvpmedPP(i)) > threshold
            bvpPPnew(i) = bvpmedPP(i);
        end
    end

    % Extract filtered PRV samples
    filtered_bvpPRV = bvpPPnew(boundary:M-(L-1)/2);
    
    % Display the filtered PRV signal
    disp(filtered_bvpPRV);
    
    if show
        figure;
        hold on;
        plot(bvp_pp_intervals,'b');
        plot(filtered_bvpPRV,'r');
        legend("Original BVP","Filtered BVP");
        xlim([0 length(bvp_pp_intervals)])
        ylim([0.5, 1.4]);
        title('Filtered Contact Pulse Rate Variability (cPRV) signal');
    end
end

This is the full function to resample, extract PRV signal and do PRV filtering for my BVP signal.

The result is as so:

UBFC-Phys (Subject 18 - T2)

[William Rose]

@Jeffery,

[Edit: adjusted the discussion of findpeaks()]

you have posted a longer code fragment and more of the paper by Sabour et al., 2022. As I said in my earlier answer, I would be wary of using

findpeaks(bvp_resampled, MinPeakDistance = fs_target/2);

when detecting peaks. Since you have not specified a sampling rate when calling findpeaks(), the MinPeakDistance is in units of samples. Therefore your command will result in finding peaks that are separated by at least half a second, and it will reject peaks corresponding to a pulse rate of greater than 120 bpm. Do you really want to do that?

I do not have the "Kubios HRV (ver. 3.3) users's guide" which Sabour cites, so I cannot comment on whether your method correctly implements the median filter. I am not sure if Sabour et al. apply the median filter on the 128 Hz signal before peak detection, or they use the median filter on the HRV sequence after peak detection.

The reflection of the signal at the ends, before median filtering, is a widely used technique, but there are different opinions about how to do it. Should the extension be flipped, or upside down and flipped? You use flipped. Maybe Sabour et al. do the same, but I am not sure. Upside down and flipped has the advantage of producing a signal whose derivative is approximately continuous at the boundary. Consider the example below. The red trace is flipped, the blue trace is upside down and flipped. Which is better?

t=0:1/40:1; tExt=1:1/40:1.25;
x=sin(2*pi*t)+1;
x1=fliplr(x(end-10:end)); x2=-fliplr(x(end-10:end))+2*x(end);
plot(t,x,'-k.',tExt,x1,'-r.',tExt,x2,'-b.')

If you wish to discuss further, please conatct me using the envelope icon which appears when you click the WR icon at the top of my posts.

Answer 2

[Star Strider]

The Signal Processing Toolbox resample function is specifically intended for signal processing and incorporates anti-aliasing filters. So it would be best to use it instead. I usually use:

[sr,tr] = resample(s,t,Fs);

where ‘s’ and ‘t’ are the original signal and time vectors, and ‘sr’ and ‘tr’ are the resampled versions, with ‘Fs’ the desired resampled sampling frequency. The advantage is that this will produce a uniformly-sampled signal regardless of whether the original signal was uniformly sampled. Uniform sampling is required for most signal processing (the nufft and nufftn funcitons being the only exceptions).

I do not know what the characteristics of your signsls are, however it is relatively straightforward to filter baseline variation and high frequency noise in MATLAB if you need to do that. I usually use the bandpass funciton with the 'ImpulseResponse','iir' name-value pair to produce the most efficient filter.

I would calculate the peak-to-peak intervals (and I assume these are what are usually considered to be pulse plethysmograph (PPG) signals) by first converting the locations to time:

bvp_pp_intervals = diff(tr(BVPpeak_locs));

and similarly for the others. There are several ways of being certain that findpeaks returns the information you want. I usually include the 'MinPeakProminence' name-value pair to keep from returning minor peaks and noise as identified peaks.

(I usually use the linspace function to calculate the time vectors, since it’s usually easier.)

The one thing I would caution is using PPG to quantify heart rate variability (HRV). It is not good for that, because HRV analysis is only reliable using EKG data. This is because the QRS complexes used to analyse HRV have to be preceeded by a normal P-wave in order to be useful. The PPG data do not provide that information. .

.

[William Rose]

@Jeffery,

I agree with @Star Strider that resample() has advantages over interp1(). I also agree with @Star Strider that heart rate variability studies should use the ECG. If I received a manuscript to review which called it "heart rate variability", I would remind the authors that if they did not analyze the ECG, it's not heart rate variability.

An increasing number of psychophysics studies estimate pulse rate from video and/or photoplethysmography (PPG) recordings, because video and PPG is less expensive and is considered less invasive (although doing an ECG hardly seems invasive to me). Therefore the term "pulse rate variability" (PRV) is used by Sabour et al. (2022), whose work you are seeking to understand and replicate.

Sabour et al. impose mental stress on their subjects. They wish to use PRV to detect mental stress. Sabour et al. claim that "the same features can be extracted from the two signals" (PRV and HRV). They cite Schäfer & Vagedes (2013) to support this claim. But Schafer & Vagedes 2013, in their review of many studies of HRV and PRV, concluded "Physical activity and some mental stressors seem to impair the agreement of PRV and HRV, often to an inacceptable extent." This raises methodological questions about studies which use PRV in lieu of HRV.

I commend to Sabour et al for collecting the data and making it freely available. This allows you and others to test methods for detecting stress from PPG data and video recordings.

Since you want to start by replicating Sabour et al., 2022, it makes sense to do what they did. They say they used "shape-preserving piecewise cubic interpolation" to resample. So I think it is reasonable for you to use

bvpResamp=interp1(tOrig,bvpOrig,tResamp,"pchip")

or

bvpResamp=pchip(tOrig,bvpOrig,tResamp);

which yield the same results.

[Jeffery]

Thank you very much! I did use it to resample my bvp signals. However, I am not really sure on what I did for peak detection, pulse-to-pulse intervals and PRV filtering is correct. Once again, thank you very much, I appreciate your help @Star Strider and @William Rose.

[Star Strider]

The SPT resample function resamples the entire waveform, and offers 'pchip' as one interpolation method. I have several posts here that analyse PPG signals, including filtering them and creating ensembles and ensemble averages. As to your PRV method being ‘correct’, it depends on what you want to get from it. If you want to detect HRV substituting it for EKG analysis, than I would simply not use it. If on the other hand you are looking for arrythmias, that might be worthwhile, understanding that the vascular tree is a ‘transmission line’ with its own impedance dynamics. I would use PPG only to monitor haemoglobin oxygenation and the peripheral pulse waveform (including the presence and relative timing of the dicrotic notch), since that is all that it actually measures. Everything else about it is inference.