Remove trend and detect peaks in a photoplethysmogram(PPG) signal

Question

Kahuthaman Silvadorai 2018-2-5

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/380879-remove-trend-and-detect-peaks-in-a-photoplethysmogram-ppg-signal

评论： Patricia 2023-1-20

There seems to be a sinusoidal like trend in my PPG signal. I believe this is due to respiration. I would like to detect the PPG peaks(maxima) in this signal after removing the trend or even better detect the peaks without removing the trends. Please help to provide a sample code as I am kind of new to signal processing for bio-signals.

clear all;
clc;
readData = csvread('p3_sit.csv',0,0);
ppg_head = readData(:,1).';
fs = 960; 
rec_time_minutes = (length(ppg_head)-1)/60000 %recording time calculation
wn=10/(fs/2); %lowpass 10Hz for ppg
[b,a] = butter(8,wn); 
ppg_head_data = filter(b,a,ppg_head);
N1 = length(ppg_head_data);
t1 = [0:N1-1]/fs;
%{comment here %}
figure(1)
plot(t1,ppg_head_data);
xlabel('second');ylabel('Volts');title('PPG Head Signal')

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2018-2-5

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/380879-remove-trend-and-detect-peaks-in-a-photoplethysmogram-ppg-signal#answer_303556

在 MATLAB Online 中打开

If you want to eliminate the baseline wander and offset, change to a highpass filter with a different cutoff frequency:

wn=5/(fs/2); %lowpass 10Hz for ppg
[b,a] = butter(8,wn,'high'); 
ppg_head_data = filter(b,a,ppg_head);

There are much better ways to design and implement filters. Here is one:

t = linspace(0, numel(ppg_head), numel(ppg_head))/fs;
Fs = 960;                                                   % Sampling Frequency (Hz)
Fn = Fs/2;                                                  % Nyquist Frequency (Hz)
Wp = 24.5/Fn;                                               % Stopband Frequency (Normalised)
Ws = 25.0/Fn;                                               % Passband Frequency (Normalised)
Rp =   1;                                                   % Passband Ripple (dB)
Rs =  50;                                                   % Stopband Ripple (dB)
[n,Ws] = cheb2ord(Wp,Ws,Rp,Rs);                             % Filter Order
[z,p,k] = cheby2(n,Rs,Ws,'high');                           % Filter Design, Sepcify Bandstop
[sos,g] = zp2sos(z,p,k);                                    % Convert To Second-Order-Section For Stability
figure(2)
freqz(sos, 2^16, Fs)                                        % Filter Bode Plot
ppg_head_data = filtfilt(sos, g, ppg_head);                 % Filter Signal

Use the Signal Processing Toolbox findpeaks (link) function on your filtered data to locate the peaks:

[pks,locs] = findpeaks(ppg_head_data, 'MinPeakHeight',2E+5);

You can use the ‘locs’ index vector to refer to the peaks and times in the original unfiltered signal.

21 个评论
显示 19更早的评论隐藏 19更早的评论

Star Strider 2018-2-6

在 MATLAB Online 中打开

This works:

Fs = 960;                                                   % Sampling Frequency (Hz)
t = linspace(0, numel(ppg_head), numel(ppg_head))/Fs;
Ts = mean(diff(t));                                         % Sampling Interval
Fs = 1/Ts;                                                  % Sampling Frequency
Fn = Fs/2;                                                  % Nyquist Frequency
L = numel(t);                                               % Signal Length
ppg_head = ppg_head - mean(ppg_head);
FTsig = fft(ppg_head)/L;                                    % Fourier Transform
Fv = linspace(0, 1, fix(L/2)+1)*Fn;                         % Frequency Vector
Iv = 1:numel(Fv);                                           % Index Vector
figure(1)
plot(Fv, abs(FTsig(Iv))*2)
grid
% axis([0  3    ylim])
Fs = 960;                                                   % Sampling Frequency (Hz)
Fn = Fs/2;                                                  % Nyquist Frequency (Hz)
Wp = [1.5 10.00]/Fn;                                        % Passband Frequency (Normalised)
Ws = [1.0 10.50]/Fn;                                        % Stopband Frequency (Normalised)
Rp =   1;                                                   % Passband Ripple (dB)
Rs =  50;                                                   % Stopband Ripple (dB)
[n,Ws] = cheb2ord(Wp,Ws,Rp,Rs);                             % Filter Order
[z,p,k] = cheby2(n,Rs,Ws);                                  % Filter Design, Sepcify Bandpass
[sos,g] = zp2sos(z,p,k);                                    % Convert To Second-Order-Section For Stability
figure(2)
freqz(sos, 2^16, Fs)                                        % Filter Bode Plot
ppg_head_data = filtfilt(sos, g, ppg_head);                 % Filter Signal
[pks,locs] = findpeaks(ppg_head_data, 'MinPeakHeight',1E+5, 'MinPeakDist',350);
figure(3)
plot(t, ppg_head_data, '-b')
hold on
plot(t(locs), pks, '^r')
hold off
axis([0  20.6    ylim])

You may need to experiment with the findpeaks arguments. I only looked at a section of the signal in detail.

632541 2021-4-21

Hi Star Strider ,

In the code ,

FTsig = fft(ppg_head)/L;

Can you please tell why divided by L?

What does dividing number of samples signify ?

And how is it related to signal amplitude?

Star Strider 2021-4-21

@Kaveri A — It normalises it by the signal length. (It’s also discussed in the fft documentation.)

请先登录，再进行评论。

Answer 2

Peter H Charlton 2022-11-25

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/380879-remove-trend-and-detect-peaks-in-a-photoplethysmogram-ppg-signal#answer_1111198

编辑：Peter H Charlton 2022-12-28

在 MATLAB Online 中打开

Since the previous answer, a new toolbox has become available to detect peaks in a photoplethysomogram signal:

https://uk.mathworks.com/matlabcentral/fileexchange/107644-ppg-beats

After downloading the toolbox and adding it to the Matlab path (as described here), the toolbox can be used to detect peaks in this signal as follows (this is a slight modification of one of the tutorials in the toolbox documentation, here):

% setup using the code provided in the question
readData = csvread('p3_sit.csv',0,0);
ppg_head = readData(:,1).';
fs = 960; 
wn=10/(fs/2); %lowpass 10Hz for ppg
[b,a] = butter(8,wn); 
ppg_head_data = filter(b,a,ppg_head);
% format data as required by the toolbox
S.v = ppg_head_data(:);   % PPG samples in a column vector 
S.fs = fs;                % sampling frequency in Hz
% detect PPG beats
beat_detector = 'MSPTD';  % specify a beat detection algorithm to use
[peaks, onsets, mid_amps] = detect_ppg_beats(S, beat_detector);  % detect beats
% plot the result
figure('Position', [20,20,1000,350])     % Setup figure
subplot('Position', [0.05,0.17,0.92,0.82])
t = [0:length(S.v)-1]/S.fs;             % Make time vector
plot(t, S.v, 'b'), hold on,             % Plot PPG signal
plot(t(peaks), S.v(peaks), 'or'),       % Plot detected beats
ftsize = 20;                            % Tidy up plot
set(gca, 'FontSize', ftsize, 'YTick', [], 'Box', 'off');
ylabel('PPG', 'FontSize', ftsize),
xlabel('Time (s)', 'FontSize', ftsize)

This produces the following plot, which shows beats detected in the signal:

Zooming in to a 10-second segment:

xlim([10, 20]); h = findobj(gca,'Type','line'); set(h, 'LineWidth', 2)

This plot shows the detected peaks (red circles).

Disclaimer: I'm one of the authors of this PPG beat detection toolbox.

10 个评论
显示 8更早的评论隐藏 8更早的评论

Patricia 2023-1-10

在 MATLAB Online 中打开

Thank you very much for the answer!

Excuse me, it is one of my first times using MATLAB and I am a bit lost. May I ask you something else? I have two signals, one coming from a camera (it is remote PPG, 32FPS) and another one coming from a conventional sensor (so it is typical PPG signal, fs=60) and I have to process them and find the peaks and IBIs on each signal. As I told you it is one of my first times working with these signals and I studied the toolbox to use it but I was wondering if the previous code would also work for my case. I have resampled (downsample) the PPG signal to the remote PPG signal frequency (60 to 32) to have them aligned and be able to compare them but once I try to use the lowpass filter you propose on the previous answers I see no diference on the output signal, moreover, I guess I need more pre-processing than what I am already doing since I think the algorithm is not detecting the peaks correctly (at least on the conventional ppg signal) because of my pre-processing. (And last question, can I use this code for both of the signals? I think it would work with both but I still have that question in mind)

Excuse me for the long question but I am willing to learn more and be able to finish this project! :)

I am attaching the code so far here and the data on the question, thank you in advance!

code:

%remote signal
Signals = load('RPPG.mat')
RPPGSignal = Signals.ProcessedData.RPPG;
fs1= 32;
%sensor signal
MonHearthRate=readtable('PPG_sensor.csv');
heartRate= table2array(MonHearthRate);
PPGSignal= heartRate.';
fs2=60;
%% resampling (downsample the signal coming from the sensor to rPPG (60Hz to 32HZ)
PPGSignal_resampled = (resample(PPGSignal,8,15)).'; %8/15 * 60Hz = 32Hz
PPGSignal_resampled = PPGSignal_resampled.';
%% Preparing the PPG signal
ppg_noOff = detrend(PPGSignal_resampled);
%% Detect peaks on noisy PPG Signals - Only on sensor PPG signal
% setup 
fs = 32; 
wn=10/(fs/2); %lowpass 10Hz for ppg
[b,a] = butter(8,wn); 
ppg_head_data = filter(b,a,y);
% format data as required by the toolbox
S.v = y(:);   % PPG samples in a column vector 
S.fs = fs;                % sampling frequency in Hz
% detect PPG beats
beat_detector = 'MSPTD';  % specify a beat detection algorithm to use
[peaks, onsets, mid_amps] = detect_ppg_beats(S, beat_detector);  % detect beats
% plot the result
figure('Position', [20,20,1000,350])     % Setup figure
subplot('Position', [0.05,0.17,0.92,0.82])
t = [0:length(S.v)-1]/S.fs;             % Make time vector
plot(t, S.v, 'b'), hold on,             % Plot PPG signal
plot(t(peaks), S.v(peaks), 'or'),       % Plot detected beats
ftsize = 20;                            % Tidy up plot
set(gca, 'FontSize', ftsize, 'YTick', [], 'Box', 'off');
ylabel('PPG', 'FontSize', ftsize),
xlabel('Time (s)', 'FontSize', ftsize)

Peter H Charlton 2023-1-10

编辑：Peter H Charlton 2023-1-12

在 MATLAB Online 中打开

This code seems to work. Change 'do_remote' to zero to process the PPG sensor signal.

%% Load data
%remote signal
Signals = load('RPPG.mat');
RPPGSignal = Signals.ProcessedData.RPPG;
fs1= 32;
%sensor signal
MonHearthRate=readtable('PPG_sensor.csv');
heartRate= table2array(MonHearthRate);
PPGSignal= heartRate.';
fs2 = 60;
%% format data as required by the toolbox
do_remote = 1;
if do_remote
    S.v = RPPGSignal;
    S.fs = fs1;                % sampling frequency in Hz
else
    S.v = PPGSignal;
    S.fs = fs2;
end
%% detect PPG beats
beat_detector = 'MSPTD';  % specify a beat detection algorithm to use
[peaks, onsets, mid_amps] = detect_ppg_beats(S, beat_detector);  % detect beats
%% plot the result
figure('Position', [20,20,1000,350])     % Setup figure
subplot('Position', [0.05,0.17,0.92,0.82])
t = [0:length(S.v)-1]/S.fs;             % Make time vector
plot(t, S.v, 'b'), hold on,             % Plot PPG signal
plot(t(peaks), S.v(peaks), 'or'),       % Plot detected beats
ftsize = 20;                            % Tidy up plot
set(gca, 'FontSize', ftsize, 'YTick', [], 'Box', 'off');
ylabel('PPG', 'FontSize', ftsize),
xlabel('Time (s)', 'FontSize', ftsize)

Patricia 2023-1-12

Thank you! :) And would you recommed to use the preprocess_ppg_signal function as a method to prepare the signal? Or just the bandpass filtering between 0.67 and 8 Hz?

Patricia 2023-1-20

Hi again and excuse me for reopening the question. I was trying the performance assessment of the detector using the 'capnobase' option and once I get it working I do not get any result as a return, stops at this point: